|
| 1 | +import numpy as np |
| 2 | +import scipy as sp |
| 3 | +import scipy.linalg as la |
| 4 | +import quantecon as qe |
| 5 | +import matplotlib.pyplot as plt |
| 6 | +from scipy.stats import norm, lognorm |
| 7 | + |
| 8 | + |
| 9 | + |
| 10 | +class AMF_LSS_VAR: |
| 11 | + """ |
| 12 | + 此类将加性(乘性)泛函转换为QuantEcon线性状态空间系统。 |
| 13 | + """ |
| 14 | + |
| 15 | + def __init__(self, A, B, D, F=None, nu=None, x_0=None): |
| 16 | + # 解包必需的元素 |
| 17 | + self.nx, self.nk = B.shape |
| 18 | + self.A, self.B = A, B |
| 19 | + |
| 20 | + # 检查D的维度(从标量情况扩展) |
| 21 | + if len(D.shape) > 1 and D.shape[0] != 1: |
| 22 | + self.nm = D.shape[0] |
| 23 | + self.D = D |
| 24 | + elif len(D.shape) > 1 and D.shape[0] == 1: |
| 25 | + self.nm = 1 |
| 26 | + self.D = D |
| 27 | + else: |
| 28 | + self.nm = 1 |
| 29 | + self.D = np.expand_dims(D, 0) |
| 30 | + |
| 31 | + # 为加性分解创建空间 |
| 32 | + self.add_decomp = None |
| 33 | + self.mult_decomp = None |
| 34 | + |
| 35 | + # 设置F |
| 36 | + if not np.any(F): |
| 37 | + self.F = np.zeros((self.nk, 1)) |
| 38 | + else: |
| 39 | + self.F = F |
| 40 | + |
| 41 | + # 设置nu |
| 42 | + if not np.any(nu): |
| 43 | + self.nu = np.zeros((self.nm, 1)) |
| 44 | + elif type(nu) == float: |
| 45 | + self.nu = np.asarray([[nu]]) |
| 46 | + elif len(nu.shape) == 1: |
| 47 | + self.nu = np.expand_dims(nu, 1) |
| 48 | + else: |
| 49 | + self.nu = nu |
| 50 | + |
| 51 | + if self.nu.shape[0] != self.D.shape[0]: |
| 52 | + raise ValueError("nu的维度与D不一致!") |
| 53 | + |
| 54 | + # 初始化模拟器 |
| 55 | + self.x_0 = x_0 |
| 56 | + |
| 57 | + # 构建大型状态空间表示 |
| 58 | + self.lss = self.construct_ss() |
| 59 | + |
| 60 | + def construct_ss(self): |
| 61 | + """ |
| 62 | + 这将创建可以传递到quantecon LSS类的状态空间表示。 |
| 63 | + """ |
| 64 | + |
| 65 | + # 提取有用信息 |
| 66 | + nx, nk, nm = self.nx, self.nk, self.nm |
| 67 | + A, B, D, F, nu = self.A, self.B, self.D, self.F, self.nu |
| 68 | + |
| 69 | + if self.add_decomp: |
| 70 | + nu, H, g = self.add_decomp |
| 71 | + else: |
| 72 | + nu, H, g = self.additive_decomp() |
| 73 | + |
| 74 | + # 用0和1填充lss矩阵的辅助块 |
| 75 | + nx0c = np.zeros((nx, 1)) |
| 76 | + nx0r = np.zeros(nx) |
| 77 | + nx1 = np.ones(nx) |
| 78 | + nk0 = np.zeros(nk) |
| 79 | + ny0c = np.zeros((nm, 1)) |
| 80 | + ny0r = np.zeros(nm) |
| 81 | + ny1m = np.eye(nm) |
| 82 | + ny0m = np.zeros((nm, nm)) |
| 83 | + nyx0m = np.zeros_like(D) |
| 84 | + |
| 85 | + # 为LSS构建A矩阵 |
| 86 | + # 状态顺序为:[1, t, xt, yt, mt] |
| 87 | + A1 = np.hstack([1, 0, nx0r, ny0r, ny0r]) # 1的转移 |
| 88 | + A2 = np.hstack([1, 1, nx0r, ny0r, ny0r]) # t的转移 |
| 89 | + A3 = np.hstack([nx0c, nx0c, A, nyx0m.T, nyx0m.T]) # x_{t+1}的转移 |
| 90 | + A4 = np.hstack([nu, ny0c, D, ny1m, ny0m]) # y_{t+1}的转移 |
| 91 | + A5 = np.hstack([ny0c, ny0c, nyx0m, ny0m, ny1m]) # m_{t+1}的转移 |
| 92 | + Abar = np.vstack([A1, A2, A3, A4, A5]) |
| 93 | + |
| 94 | + # 为LSS构建B矩阵 |
| 95 | + Bbar = np.vstack([nk0, nk0, B, F, H]) |
| 96 | + |
| 97 | + # 为LSS构建G矩阵 |
| 98 | + # 观测顺序为:[xt, yt, mt, st, tt] |
| 99 | + G1 = np.hstack([nx0c, nx0c, np.eye(nx), nyx0m.T, nyx0m.T]) # x_{t}的选择器 |
| 100 | + G2 = np.hstack([ny0c, ny0c, nyx0m, ny1m, ny0m]) # y_{t}的选择器 |
| 101 | + G3 = np.hstack([ny0c, ny0c, nyx0m, ny0m, ny1m]) # 鞅的选择器 |
| 102 | + G4 = np.hstack([ny0c, ny0c, -g, ny0m, ny0m]) # 平稳部分的选择器 |
| 103 | + G5 = np.hstack([ny0c, nu, nyx0m, ny0m, ny0m]) # 趋势的选择器 |
| 104 | + Gbar = np.vstack([G1, G2, G3, G4, G5]) |
| 105 | + |
| 106 | + # 为LSS构建H矩阵 |
| 107 | + Hbar = np.zeros((Gbar.shape[0], nk)) |
| 108 | + |
| 109 | + # 构建LSS类型 |
| 110 | + if not np.any(self.x_0): |
| 111 | + x0 = np.hstack([1, 0, nx0r, ny0r, ny0r]) |
| 112 | + else: |
| 113 | + x0 = self.x_0 |
| 114 | + |
| 115 | + S0 = np.zeros((len(x0), len(x0))) |
| 116 | + lss = qe.lss.LinearStateSpace(Abar, Bbar, Gbar, Hbar, mu_0=x0, Sigma_0=S0) |
| 117 | + |
| 118 | + return lss |
| 119 | + |
| 120 | + def additive_decomp(self): |
| 121 | + """ |
| 122 | + 返回鞅分解的值 |
| 123 | + - nu : Y的无条件均值差异 |
| 124 | + - H : (线性)鞅分量的系数 (kappa_a) |
| 125 | + - g : 平稳分量g(x)的系数 |
| 126 | + - Y_0 : 应该是X_0的函数(目前设置为0.0) |
| 127 | + """ |
| 128 | + I = np.eye(self.nx) |
| 129 | + A_res = la.solve(I - self.A, I) |
| 130 | + g = self.D @ A_res |
| 131 | + H = self.F + self.D @ A_res @ self.B |
| 132 | + |
| 133 | + return self.nu, H, g |
| 134 | + |
| 135 | + def multiplicative_decomp(self): |
| 136 | + """ |
| 137 | + 返回乘性分解的值(示例5.4.4) |
| 138 | + - nu_tilde : 特征值 |
| 139 | + - H : Jensen项的向量 |
| 140 | + """ |
| 141 | + nu, H, g = self.additive_decomp() |
| 142 | + nu_tilde = nu + (.5)*np.expand_dims(np.diag(H @ H.T), 1) |
| 143 | + |
| 144 | + return nu_tilde, H, g |
| 145 | + |
| 146 | + |
| 147 | + |
| 148 | + |
| 149 | + |
| 150 | +def future_moments(amf_future, N=25): |
| 151 | + """ |
| 152 | + 计算未来时刻的矩 |
| 153 | + """ |
| 154 | + nx, nk, nm = amf_future.nx, amf_future.nk, amf_future.nm |
| 155 | + nu_tilde, H, g = amf_future.multiplicative_decomp() |
| 156 | + |
| 157 | + # 分配空间(nm是加性泛函的数量) |
| 158 | + mbounds = np.zeros((3, N)) |
| 159 | + sbounds = np.zeros((3, N)) |
| 160 | + ybounds = np.zeros((3, N)) |
| 161 | + #mbounds_mult = np.zeros((3, N)) |
| 162 | + #sbounds_mult = np.zeros((3, N)) |
| 163 | + |
| 164 | + # 模拟所需的时长 |
| 165 | + moment_generator = amf_future.lss.moment_sequence() |
| 166 | + tmoms = next(moment_generator) |
| 167 | + |
| 168 | + # 提取总体矩 |
| 169 | + for t in range (N-1): |
| 170 | + tmoms = next(moment_generator) |
| 171 | + ymeans = tmoms[1] |
| 172 | + yvar = tmoms[3] |
| 173 | + |
| 174 | + # 每个加性泛函的上下界 |
| 175 | + yadd_dist = norm(ymeans[nx], np.sqrt(yvar[nx, nx])) |
| 176 | + ybounds[:2, t+1] = yadd_dist.ppf([0.1, .9]) |
| 177 | + ybounds[2, t+1] = yadd_dist.mean() |
| 178 | + |
| 179 | + madd_dist = norm(ymeans[nx+nm], np.sqrt(yvar[nx+nm, nx+nm])) |
| 180 | + mbounds[:2, t+1] = madd_dist.ppf([0.1, .9]) |
| 181 | + mbounds[2, t+1] = madd_dist.mean() |
| 182 | + |
| 183 | + sadd_dist = norm(ymeans[nx+2*nm], np.sqrt(yvar[nx+2*nm, nx+2*nm])) |
| 184 | + sbounds[:2, t+1] = sadd_dist.ppf([0.1, .9]) |
| 185 | + sbounds[2, t+1] = sadd_dist.mean() |
| 186 | + |
| 187 | + |
| 188 | + #Mdist = lognorm(np.asscalar(np.sqrt(yvar[nx+nm, nx+nm])), scale=np.asscalar(np.exp(ymeans[nx+nm]- \ |
| 189 | + # t*(.5)*np.expand_dims(np.diag(H @ H.T), 1)))) |
| 190 | + #Sdist = lognorm(np.asscalar(np.sqrt(yvar[nx+2*nm, nx+2*nm])), |
| 191 | + # scale = np.asscalar(np.exp(-ymeans[nx+2*nm]))) |
| 192 | + #mbounds_mult[:2, t+1] = Mdist.ppf([.01, .99]) |
| 193 | + #mbounds_mult[2, t+1] = Mdist.mean() |
| 194 | + |
| 195 | + #sbounds_mult[:2, t+1] = Sdist.ppf([.01, .99]) |
| 196 | + #sbounds_mult[2, t+1] = Sdist.mean() |
| 197 | + |
| 198 | + ybounds[:, 0] = amf_future.x_0[2+nx] |
| 199 | + mbounds[:, 0] = mbounds[-1, 1] |
| 200 | + sbounds[:, 0] = -g @ amf_future.x_0[2:2+nx] |
| 201 | + |
| 202 | + #mbounds_mult[:, 0] = mbounds_mult[-1, 1] |
| 203 | + #sbounds_mult[:, 0] = np.exp(-g @ amf_future.x_0[2:2+nx]) |
| 204 | + |
| 205 | + return mbounds, sbounds, ybounds #, mbounds_mult, sbounds_mult |
0 commit comments