From ca9d286e82e0875000669dcaa529e836a7ac80b6 Mon Sep 17 00:00:00 2001 From: Matt McKay Date: Thu, 22 Oct 2026 14:44:49 +1100 Subject: [PATCH 1/9] Update translation: lectures/measurement_models.md --- lectures/measurement_models.md | 392 +++++++++++++++++++++++++++------ 1 file changed, 329 insertions(+), 63 deletions(-) diff --git a/lectures/measurement_models.md b/lectures/measurement_models.md index 270f9981..82038f46 100644 --- a/lectures/measurement_models.md +++ b/lectures/measurement_models.md @@ -27,6 +27,7 @@ translation: A classical model of measurements initially collected by an agency::Computing coefficients in a Wold moving average representation: 计算 Wold 移动平均表示中的系数 A classical model of measurements initially collected by an agency::Gaussian likelihood: 高斯似然 A classical model of measurements initially collected by an agency::Forecast-error-variance decomposition: 预测误差方差分解 + A classical model of measurements initially collected by an agency::Granger causality without an ordering: 不依赖排序的格兰杰因果关系 A classical model of measurements initially collected by an agency::Wold impulse responses: Wold 脉冲响应 A model of optimal estimates reported by an agency: 一个由机构报告最优估计的模型 A model of optimal estimates reported by an agency::Innovations representation for filtered data: 滤波数据的新息表示 @@ -36,6 +37,7 @@ translation: A model of optimal estimates reported by an agency::Wold impulse responses: Wold 脉冲响应 Simulation: 模拟 Summary: 总结 + Exercises: 练习 --- (measurement_models)= @@ -163,7 +165,7 @@ c_t + k_{t+1} = f k_t + \theta_t, \qquad \beta f^2 > 1. :label: shock_process a(L)\,\theta_t = \varepsilon_t, ``` -其中 $L$ 是后移（或"滞后"）算子，$a(z) = 1 - a_1 z - a_2 z^2 - \cdots - a_r z^r$ 的所有零点都在单位圆之外。 +其中 $L$ 是后移（或"滞后"）算子，$\varepsilon_t$ 是标量白噪声，而 $a(z) = 1 - a_1 z - a_2 z^2 - \cdots - a_r z^r$ 的所有零点都在单位圆之外。 ### 最优决策规则 @@ -173,15 +175,17 @@ $c_t$ 的最优决策规则是 :label: opt_decision c_t = \frac{-\alpha}{f-1} + \left(1 - \frac{1}{\beta f^2}\right) + \left[ \frac{L - f^{-1} a(f^{-1})^{-1} a(L)}{L - f^{-1}}\,\theta_t - + f k_t, + + f k_t + \right], \qquad k_{t+1} = f k_t + \theta_t - c_t, ``` 其中 $\alpha = u_1[1-(\beta f)^{-1}]/u_2$。 -方程 {eq}`shock_process` 和 {eq}`opt_decision` 展现了理性预期模型所特有的跨方程约束。 +方程 {eq}`shock_process` 和 {eq}`opt_decision` 展现了理性预期模型所特有的跨方程约束：冲击过程的参数 $a_1, \ldots, a_r$ 也出现在消费的决策规则中。 ### 净收入与加速数 @@ -207,13 +211,15 @@ y_{nt} = (f-1)k_t + \theta_t. 假设2不那么重要，只影响各种常数。 -在这两个假设下，{eq}`opt_decision` 简化为 +在这两个假设下，$\alpha = 0$，{eq}`opt_decision` 简化为 ```{math} :label: simple_crule -c_t = (1-f^{-1})\,\theta_t + (f-1)\,k_t. +c_t = (1-f^{-1})\,\theta_t + (f-1)\,k_t, ``` +从而资本存量遵循 $k_{t+1} = k_t + f^{-1} \theta_t$。 + 当 {eq}`simple_crule`、{eq}`net_income` 和 {eq}`tech_constraint` 结合起来时，最优计划满足 ```{math} @@ -237,13 +243,17 @@ y_{nt} = \theta_t + (1-\beta)(\theta_{t-1} + \theta_{t-2} + \cdots). 这与 {cite:t}`Chow1968` 实证记录的机制相同（见 {doc}`chow_business_cycles`）。 -方程 {eq}`income_process` 表明，可支配收入的一阶差分是一个一阶移动平均过程，其新息等于禀赋冲击 $\theta_t$ 的新息。 +方程 {eq}`income_process` 表明，净收入的一阶差分 $y_{nt} - y_{n,t-1} = \theta_t - \beta \theta_{t-1}$ 是一个一阶移动平均过程，其新息等于禀赋冲击 $\theta_t$ 的新息。 正如 {cite:t}`Muth1960` 所展示的，这样的过程可以通过几何分布滞后或"适应性预期"方案最优地预测。 +消费本身是一个鞅：{eq}`simple_crule` 与 $k_t - k_{t-1} = f^{-1}\theta_{t-1}$ 蕴含 $c_t - c_{t-1} = (1 - \beta)\theta_t$，这正是 {cite:t}`Hall1978` 所强调的消费随机游走性质。 + +因此，消费、资本存量和净收入都具有单位根，而投资 $k_{t+1} - k_t = f^{-1}\theta_t$ 则是序列不相关的。 + ### 加速数之谜 -当所有变量都被准确测量并由单一冲击 $\theta_t$ 驱动时，$(c_t,\, k_{t+1}-k_t,\, y_{nt})$ 的谱密度矩阵在所有频率上的秩为1。 +当所有变量都被准确测量并由单一冲击 $\theta_t$ 驱动时，协方差平稳向量 $(c_t - c_{t-1},\, k_{t+1}-k_t,\, y_{nt} - y_{n,t-1})$ 的谱密度矩阵在所有频率上的秩为1。 每个变量都是同一白噪声的一个可逆单边分布滞后，因此没有任何变量格兰杰引起任何其他变量。 @@ -267,10 +277,7 @@ y_{nt} = \theta_t + (1-\beta)(\theta_{t-1} + \theta_{t-2} + \cdots). + (v_{k,t+1} - v_{kt}), ``` -```{math} -:label: income_process_ma -y_{nt} = \theta_t + (1-\beta)(\theta_{t-1} + \theta_{t-2} + \cdots). -``` +其中 $y_{nt}$ 由 {eq}`income_process` 给出。 在这种情况下，收入格兰杰引起消费和投资，但不被它们格兰杰引起。 @@ -280,7 +287,7 @@ y_{nt} = \theta_t + (1-\beta)(\theta_{t-1} + \theta_{t-2} + \cdots). 在这种情况下，每个观测到的序列都将共同信号 $\theta_t$ 与特异性测量噪声混合在一起。 -测量误差方差较低的序列更紧密地跟踪 $\theta_t$，因此它的新息包含更多关于其他序列未来值的信息。 +**相对于真实序列的变动而言**测量误差较小的序列更紧密地跟踪 $\theta_t$，因此它的新息包含更多关于其他序列未来值的信息。 因此，在预测误差方差分解中，对测量较好的序列的冲击在其他变量的 $j$ 步向前预测误差中占更大份额。 @@ -302,17 +309,19 @@ x_t = \begin{bmatrix} k_t \\ \theta_t \end{bmatrix}, z_t = \begin{bmatrix} y_{nt} \\ c_t \\ \Delta k_t \end{bmatrix}, ``` -使得无误差数据由状态空间系统描述 +其中 $\Delta k_t \equiv k_{t+1} - k_t$ 表示 $t$ 时刻的投资。 + +无误差数据由状态空间系统描述 ```{math} :label: true_ss \begin{aligned} -x_{t+1} &= A x_t + \varepsilon_t, \\ -z_t &= C x_t. +x_{t+1} &= A x_t + w_t, \\ +z_t &= C x_t, \end{aligned} ``` -其中 $\varepsilon_t = \begin{bmatrix} 0 \\ \theta_t \end{bmatrix}$ 具有协方差 $E \varepsilon_t \varepsilon_t^\top = Q$，而各矩阵为 +其中 $w_t = \begin{bmatrix} 0 \\ \theta_{t+1} \end{bmatrix}$ 具有协方差 $E w_t w_t^\top = Q$，而各矩阵为 ```{math} A = \begin{bmatrix} @@ -332,8 +341,12 @@ Q = \begin{bmatrix} \end{bmatrix}. ``` +这里我们用 $t$ 来标注在 $t$ 与 $t+1$ 之间到达的冲击 $w_t$；这与下文测量误差的时间约定相匹配。 + $Q$ 是奇异的，因为只有一个随机性来源 $\theta_t$；给定 $\theta_t$，资本存量 $k_t$ 确定性地演化。 +注意 $C$ 的各行满足 $C_c + C_{\Delta k} = C_y$，这正是国民收入恒等式 $c_t + \Delta k_t = y_{nt}$。 + ```{code-cell} ipython3 # 真实经济的基线结构矩阵 f = 1.05 @@ -365,7 +378,7 @@ Q = np.array([ 响应清楚地展示了投资加速数：对净收入 $y_n$ 的全部影响发生在滞后0处，而消费仅调整 $1 - f^{-1} \approx 0.048$，投资吸收余下部分。 -从滞后1开始，经济处于新的稳态 +从滞后1开始，经济处于新的稳态，拥有永久更高的资本、收入和消费 ```{code-cell} ipython3 def table2_irf(A, C, n_lags=6): @@ -395,14 +408,14 @@ display(Latex(df_to_latex_array(table2))) v_{t+1} = D v_t + \eta_t, ``` -其中 $\eta_t$ 是一个向量白噪声，满足 $E \eta_t \eta_t^\top = \Sigma_\eta$ 且对所有 $t, s$ 有 $E \varepsilon_t v_s^\top = 0$。 +其中 $\eta_t$ 是一个向量白噪声，满足 $E \eta_t \eta_t^\top = \Sigma_\eta$ 且对所有 $t, s$ 有 $E w_t v_s^\top = 0$。 参数为 ```{math} -D = \operatorname{diag}(0.6, 0.7, 0.3), +D = \operatorname{diag}(\rho_y, \rho_c, \rho_{\Delta k}) = \operatorname{diag}(0.6, 0.7, 0.3), \qquad -\sigma_\eta = (0.05, 0.035, 0.65), +\Sigma_\eta = \operatorname{diag}(\sigma_{\eta}^2), \quad \sigma_\eta = (0.05, 0.035, 0.65), ``` 因此 $v_t$ 的无条件协方差为 @@ -411,9 +424,13 @@ D = \operatorname{diag}(0.6, 0.7, 0.3), R = \operatorname{diag}\!\left(\frac{\sigma_{\eta,i}^2}{1 - \rho_i^2}\right). ``` -新息方差对消费最小（$\sigma_\eta = 0.035$），其次是收入（$\sigma_\eta = 0.05$），投资最大（$\sigma_\eta = 0.65$）。 +新息标准差 $\sigma_\eta$ 对消费最小（$0.035$），其次是收入（$0.05$），投资最大（$0.65$）。 + +然而，对于格兰杰因果关系不对称性而言，重要的是每个测量误差**相对于**相应真实序列中信号的大小。 -正如 {cite:t}`Sargent1989` 以及我们上面的讨论，对于格兰杰因果关系不对称性而言重要的是整个系统中的总体测量质量：产出测量相对较好，而投资测量相对较差。 +真实序列的新息为：收入是 $\theta_t$（标准差为 $1$），消费是 $(1-\beta)\theta_t$（标准差为 $0.048$），投资是 $\beta \theta_t$（标准差为 $0.95$）。 + +相对于这些新息，收入测量得非常好（比率约为 $0.05$），而消费和投资的测量都较差（比率约为 $0.7$）。 ```{code-cell} ipython3 ρ = np.array([0.6, 0.7, 0.3]) @@ -441,7 +458,7 @@ display(Latex(df_to_latex_matrix(pd.DataFrame(D), 'D'))) def steady_state_kalman(A, C_obs, Q, R, W=None, tol=1e-13, max_iter=200_000): """ Solve steady-state Kalman equations for - x_{t+1} = A x_t + w_{t+1} + x_{t+1} = A x_t + w_t y_t = C_obs x_t + v_t with cov(w)=Q, cov(v)=R, cov(w,v)=W. """ @@ -485,7 +502,7 @@ def steady_state_kalman(A, C_obs, Q, R, W=None, tol=1e-13, max_iter=200_000): 总体向量自回归及其如何依赖于状态空间系统和测量误差过程的参数，为如何解释 $\bar z_t$ 的估计向量自回归提供了洞见。 -构造向量自回归作为将 $\bar z_t$ 样本的似然作为自由参数 $\{A, C, D, Q, R\}$ 的函数进行计算的中间步骤也是有用的。 +构造向量自回归作为将 $\bar z_t$ 样本的似然作为自由参数 $\{A, C, D, Q, \Sigma_\eta\}$ 的函数进行计算的中间步骤也是有用的。 用于构造向量自回归表示的特定方法也被证明在构造最优报告机构模型时作为中间步骤是有用的。 @@ -504,7 +521,7 @@ def steady_state_kalman(A, C_obs, Q, R, W=None, tol=1e-13, max_iter=200_000): ```{math} :label: model1_qd \tilde z_t = \bar z_{t+1} - D \bar z_t, \qquad -\bar\nu_t = C \varepsilon_t + \eta_t, \qquad +\bar\nu_t = C w_t + \eta_t, \qquad \bar C = CA - DC. ``` @@ -513,22 +530,24 @@ def steady_state_kalman(A, C_obs, Q, R, W=None, tol=1e-13, max_iter=200_000): ```{math} :label: model1_transformed \begin{aligned} -x_{t+1} &= A x_t + \varepsilon_t, \\ +x_{t+1} &= A x_t + w_t, \\ \tilde z_t &= \bar C\, x_t + \bar\nu_t, \end{aligned} ``` -其中 $(\varepsilon_t, \bar\nu_t)$ 是一个白噪声过程，满足 +其中 $(w_t, \bar\nu_t)$ 是一个白噪声过程，满足 ```{math} :label: model1_covs -E \begin{bmatrix} \varepsilon_t \end{bmatrix} -\begin{bmatrix} \varepsilon_t^\top & \bar\nu_t^\top \end{bmatrix} +E \begin{bmatrix} w_t \\ \bar\nu_t \end{bmatrix} +\begin{bmatrix} w_t \\ \bar\nu_t \end{bmatrix}^\top = \begin{bmatrix} Q & W_1 \\ W_1^\top & R_1 \end{bmatrix}, \qquad -R_1 = C Q C^\top + R, \quad W_1 = Q C^\top. +R_1 = C Q C^\top + \Sigma_\eta, \quad W_1 = Q C^\top. ``` +注意，$R_1$ 涉及测量误差的**新息** $\eta_t$ 的协方差矩阵 $\Sigma_\eta$，而不是测量误差本身的无条件协方差矩阵 $R$：拟差分消除了 $v_t$ 中序列相关的部分。 + 具有协方差 {eq}`model1_covs` 的系统 {eq}`model1_transformed` 由五个矩阵 $[A, \bar C, Q, R_1, W_1]$ 刻画。 ### 新息表示 @@ -586,13 +605,12 @@ V_1 = E\, u_t u_t^\top = \bar C\, S_1\, \bar C^\top + R_1. ```{code-cell} ipython3 C_bar = C @ A - D @ C -R1 = C @ Q @ C.T + R +R1 = C @ Q @ C.T + Σ_η W1 = Q @ C.T K1, S1, V1 = steady_state_kalman(A, C_bar, Q, R1, W1) ``` - ### 计算 Wold 移动平均表示中的系数 为了数值计算 {eq}`model1_wold` 中的移动平均系数，定义增广状态 @@ -604,9 +622,9 @@ r_t = \begin{bmatrix} \hat x_{t-1} \\ \bar z_{t-1} \end{bmatrix}, 具有动态 ```{math} -r_{t+1} = F_1 r_t + G_1 u_t, +r_{t+1} = F_1 r_t + G_1 u_{t-1}, \qquad -\bar z_t = H_1 r_t + u_t, +\bar z_t = H_1 r_t + u_{t-1}, ``` 其中 @@ -627,6 +645,8 @@ I H_1 = [\bar C \;\; D]. ``` +（回想一下，$u_{t-1}$ 是 $\bar z_t$ 的新息。） + 移动平均系数则为 $\psi_0 = I$ 和 $\psi_j = H_1 F_1^{j-1} G_1$（对 $j \geq 1$）。 ```{code-cell} ipython3 @@ -674,11 +694,11 @@ decomp1 = fev_contributions(psi1, V1, n_horizons=20) ```{math} :label: model1_loglik -\mathcal{L}^* = -T\ln 2\pi - \tfrac{1}{2}T\ln|V_1| +\mathcal{L}^* = -\tfrac{3}{2}T\ln 2\pi - \tfrac{1}{2}T\ln|V_1| - \tfrac{1}{2}\sum_{t=0}^{T-1} u_t^\top V_1^{-1} u_t, ``` -其中 $u_t$ 是由下面 {eq}`model1_recursion` 定义的 $\{\bar z_t\}$ 的函数。 +其中 $u_t$ 是由下面 {eq}`model1_recursion` 定义的 $\{\bar z_t\}$ 的函数，第一项中的 $3$ 是被测量序列的数量。 为了使用 {eq}`model1_innov` 计算 $\{u_t\}$，有用的做法是将其表示为 @@ -700,12 +720,16 @@ u_t &= -\bar C\,\hat x_t + \tilde z_t, 为了度量每个新息的相对重要性，我们分解每个测量变量的 $j$ 步向前预测误差方差。 -写 $\bar z_{t+j} - E_t \bar z_{t+j} = \sum_{i=0}^{j-1} \psi_i u_{t+j-i}$。 +写 $\bar z_{t+j} - E_t \bar z_{t+j} = \sum_{i=0}^{j-1} \psi_i u_{t+j-i}$，其中为简化记号，此处 $u_{t+j-i}$ 表示 $\bar z_{t+j-i}$ 的新息。 设 $P$ 是 $V_1$ 的下三角 Cholesky 因子，使得正交化的新息为 $e_t = P^{-1} u_t$。 那么正交化新息 $k$ 对变量 $m$ 的 $j$ 步向前方差的贡献为 $\sum_{i=0}^{j-1} (\psi_i P)_{mk}^2$。 +我们按 $y_n$、$c$、$\Delta k$ 的顺序排列变量。 + +在此排序下，第一个正交化新息就是 $y_n$ 的新息，第二个是 $c$ 的新息中与 $y_n$ 新息正交的部分，第三个是 $\Delta k$ 的新息中与前两者都正交的部分。 + 下表显示了在视界1到20下每个正交化新息对 $y_n$、$c$ 和 $\Delta k$ 预测误差方差的累积贡献。 每个面板固定一个正交化新息并报告其对每个变量预测误差方差的累积贡献。 @@ -744,6 +768,10 @@ display(Latex('$' + r' \quad '.join(parts) + '$')) 这与论文的信息相符，即在单一共同指数模型中，测量相对最好的序列具有最强的预测能力。 +由于收入被排在第一位，该分解把收入新息与其他新息共享的全部方差都归因于收入新息。 + +下面我们通过一个不依赖于排序的计算来确认这种格兰杰因果关系模式。 + 让我们看看新息的协方差矩阵 ```{code-cell} ipython3 @@ -762,6 +790,54 @@ print(np.sort(np.linalg.eigvalsh(V1))[::-1].round(4)) 第一个特征值远大于其他特征值，这与存在一个占主导地位的共同冲击 $\theta_t$ 相一致 +### 不依赖排序的格兰杰因果关系 + +$\bar z_t$ 的总体向量自回归直接展示了格兰杰因果关系。 + +对 {eq}`model1_recursion` 进行迭代给出 + +$$ +\tilde z_t = \sum_{j=1}^\infty \bar C (A - K_1 \bar C)^{j-1} K_1 \tilde z_{t-j} + u_t . +$$ + +令 $B_j = \bar C (A - K_1 \bar C)^{j-1} K_1$。 + +代入 $\tilde z_t = \bar z_{t+1} - D \bar z_t$ 并整理各项，得到向量自回归 + +$$ +\bar z_{t+1} = \sum_{j=1}^\infty \Pi_j \bar z_{t+1-j} + u_t, +\qquad +\Pi_1 = D + B_1, \quad \Pi_j = B_j - B_{j-1} D \ \ (j \geq 2). +$$ + +变量 $i$ 不格兰杰引起变量 $m$，当且仅当所有 $\Pi_j$ 的 $(m, i)$ 元素都为零。 + +以下是前两个系数矩阵；行是方程，列是滞后变量。 + +```{code-cell} ipython3 +M_filter = A - K1 @ C_bar +print("eigenvalues of A - K1 C_bar:", np.linalg.eigvals(M_filter).round(4)) + +B_coef = [C_bar @ np.linalg.matrix_power(M_filter, j) @ K1 for j in range(3)] +Π = [D + B_coef[0]] + [B_coef[j] - B_coef[j-1] @ D for j in range(1, 3)] + +eq_labels = [r'\text{eq. } ' + lab for lab in labels] +for j in range(2): + df = pd.DataFrame(Π[j], index=eq_labels, columns=labels) + display(Latex(r'$\Pi_' + str(j + 1) + ' = ' + + df_to_latex_array(df).strip('$') + '$')) +``` + +$A - K_1 \bar C$ 的特征值都在单位圆内，所以该向量自回归是收敛的。 + +在投资方程中，滞后收入的系数很大（约为 $-0.28$）。 + +在收入方程中，滞后投资的系数极小（约为 $-0.004$）。 + +滞后收入有助于预测投资，但滞后投资几乎无助于预测收入，这正是 {cite:t}`Sargent1989` 试图解释的格兰杰因果关系模式。 + +{ref}`meas_ex2` 计算了 {cite:t}`Geweke1982` 提出的一种不依赖排序的格兰杰因果关系强度度量。 + ### Wold 脉冲响应 Wold 表示中的脉冲响应使用正交化新息报告（$V_1$ 的 Cholesky 分解，排序为 $y_n$、$c$、$\Delta k$）。 @@ -801,7 +877,7 @@ display(Latex('$' + r' \quad '.join(parts) + '$')) 在后续滞后处，收入新息在所有三个变量中产生持续的响应，因为作为测量最好的序列，它的新息由真实的永久冲击 $\theta_t$ 主导。 -消费和投资新息产生的响应根据它们各自测量误差的 AR(1) 结构衰减（$\rho_c = 0.7$，$\rho_{\Delta k} = 0.3$），对其他变量几乎没有溢出。 +消费和投资新息产生的响应以接近它们各自测量误差自回归系数的速率衰减（$\rho_c = 0.7$，$\rho_{\Delta k} = 0.3$），对其他变量几乎没有溢出。 ## 一个由机构报告最优估计的模型 @@ -819,9 +895,9 @@ display(Latex('$' + r' \quad '.join(parts) + '$')) 为了准备其估计，报告机构自己计算卡尔曼滤波器以获得新息表示 {eq}`model1_innov`。 -报告机构不报告含误差数据 $\bar z_t$，而是报告 $\tilde z_t = G \hat x_t$，其中 $G$ 是一个"选择矩阵"，对于机构所报告的数据可能等于 $C$。 +报告机构不报告含误差数据 $\bar z_t$，而是报告 $\hat z_t = G \hat x_t$，其中 $G$ 是一个"选择矩阵"，对于机构所报告的数据可能等于 $C$。 -数据 $G \hat x_t = E[G x_t \mid \bar z_t, \bar z_{t-1}, \ldots, \hat x_0]$。 +报告数据为 $\hat z_t = G \hat x_t = E[G x_t \mid \bar z_t, \bar z_{t-1}, \ldots, \hat x_0]$。 报告数据的状态空间表示则为 @@ -829,15 +905,15 @@ display(Latex('$' + r' \quad '.join(parts) + '$')) :label: model2_state \begin{aligned} \hat x_{t+1} &= A \hat x_t + K_1 u_t, \\ -\tilde z_t &= G \hat x_t, +\hat z_t &= G \hat x_t, \end{aligned} ``` 其中 {eq}`model2_state` 的第一行来自新息表示 {eq}`model1_innov`。 -注意 $u_t$ 是 $\bar z_{t+1}$ 的新息，而*不是* $\tilde z_t$ 的新息。 +注意 $u_t$ 是 $\bar z_{t+1}$ 的新息，而*不是* $\hat z_t$ 的新息。 -为了得到 $\tilde z_t$ 的 Wold 表示以及 $\tilde z_t$ 样本的似然函数，我们需要为 {eq}`model2_state` 得到一个新息表示。 +为了得到 $\hat z_t$ 的 Wold 表示以及 $\hat z_t$ 样本的似然函数，我们需要为 {eq}`model2_state` 得到一个新息表示。 ### 滤波数据的新息表示 @@ -847,18 +923,18 @@ display(Latex('$' + r' \quad '.join(parts) + '$')) :label: model2_obs \begin{aligned} \hat x_{t+1} &= A \hat x_t + K_1 u_t, \\ -\tilde z_t &= G \hat x_t + \eta_t, +\hat z_t &= G \hat x_t + \zeta_t, \end{aligned} ``` -其中 $\eta_t$ 是一个第二类白噪声测量误差过程（"打字错误"），假定其协方差矩阵 $R_2$ 非常小。 +其中 $\zeta_t$ 是一个第二类白噪声测量误差过程（"打字错误"），假定其协方差矩阵 $R_2$ 非常小。 联合噪声的协方差矩阵为 ```{math} :label: model2_Q -E \begin{bmatrix} K_1 u_t \\ \eta_t \end{bmatrix} - \begin{bmatrix} K_1 u_t \\ \eta_t \end{bmatrix}^\top +E \begin{bmatrix} K_1 u_t \\ \zeta_t \end{bmatrix} + \begin{bmatrix} K_1 u_t \\ \zeta_t \end{bmatrix}^\top = \begin{bmatrix} Q_2 & 0 \\ 0 & R_2 \end{bmatrix}, ``` @@ -874,7 +950,7 @@ E \begin{bmatrix} K_1 u_t \\ \eta_t \end{bmatrix} :label: model2_innov \begin{aligned} \check{x}_{t+1} &= A \check{x}_t + K_2 a_t, \\ -\tilde z_t &= G \check{x}_t + a_t, +\hat z_t &= G \check{x}_t + a_t, \end{aligned} ``` @@ -883,14 +959,14 @@ E \begin{bmatrix} K_1 u_t \\ \eta_t \end{bmatrix} ```{math} :label: model2_innov_defs \begin{aligned} -a_t &= \tilde z_t - E[\tilde z_t \mid \tilde z_{t-1}, \tilde z_{t-2}, \ldots], \\ -\check{x}_t &= E[\hat x_t \mid \tilde z_{t-1}, \tilde z_{t-2}, \ldots, \check{x}_0], \\ +a_t &= \hat z_t - E[\hat z_t \mid \hat z_{t-1}, \hat z_{t-2}, \ldots], \\ +\check{x}_t &= E[\hat x_t \mid \hat z_{t-1}, \hat z_{t-2}, \ldots, \check{x}_0], \\ S_2 &= E[(\hat x_t - \check{x}_t)(\hat x_t - \check{x}_t)^\top], \\ [K_2, S_2] &= \text{kalmanfilter}(A, G, Q_2, R_2, 0). \end{aligned} ``` -因此 $\{a_t\}$ 是报告数据 $\tilde z_t$ 的新息过程，其新息协方差为 +因此 $\{a_t\}$ 是报告数据 $\hat z_t$ 的新息过程，其新息协方差为 ```{math} :label: model2_V2 @@ -899,11 +975,11 @@ V_2 = E\, a_t a_t^\top = G\, S_2\, G^\top + R_2. ### Wold 表示 -由 {eq}`model2_innov` 得到 $\tilde z_t$ 的 Wold 移动平均表示为 +由 {eq}`model2_innov` 得到 $\hat z_t$ 的 Wold 移动平均表示为 ```{math} :label: model2_wold -\tilde z_t = \bigl[G(I - AL)^{-1} K_2 L + I\bigr] a_t, +\hat z_t = \bigl[G(I - AL)^{-1} K_2 L + I\bigr] a_t, ``` 系数为 $\psi_0 = I$ 和 $\psi_j = G A^{j-1} K_2$（对 $j \geq 1$）。 @@ -912,21 +988,21 @@ V_2 = E\, a_t a_t^\top = G\, S_2\, G^\top + R_2. ### 高斯似然 -当使用类似于模型1的方法时，可以通过首先利用下式从 $\tilde z_t$ 的观测中计算一个 $\{a_t\}$ 序列来计算 $\tilde z_t$ 的高斯对数似然 +当使用类似于模型1的方法时，可以通过首先利用下式从 $\hat z_t$ 的观测中计算一个 $\{a_t\}$ 序列来计算 $\hat z_t$ 的高斯对数似然 ```{math} :label: model2_recursion \begin{aligned} -\check{x}_{t+1} &= (A - K_2 G)\,\check{x}_t + K_2 \tilde z_t, \\ -a_t &= -G\,\check{x}_t + \tilde z_t. +\check{x}_{t+1} &= (A - K_2 G)\,\check{x}_t + K_2 \hat z_t, \\ +a_t &= -G\,\check{x}_t + \hat z_t. \end{aligned} ``` -$T$ 个观测样本 $\{\tilde z_t\}$ 的似然函数则为 +$T$ 个观测样本 $\{\hat z_t\}$ 的似然函数则为 ```{math} :label: model2_loglik -\mathcal{L}^{**} = -T\ln 2\pi - \tfrac{1}{2}T\ln|V_2| +\mathcal{L}^{**} = -\tfrac{3}{2}T\ln 2\pi - \tfrac{1}{2}T\ln|V_2| - \tfrac{1}{2}\sum_{t=0}^{T-1} a_t^\top V_2^{-1} a_t. ``` @@ -942,6 +1018,10 @@ $T$ 个观测样本 $\{\tilde z_t\}$ 的似然函数则为 在我们的参数化中 $G = C$。 +由于 $C$ 的各行满足国民收入恒等式，报告数据精确满足 $\hat c_t + \Delta \hat k_t = \hat y_{nt}$。 + +因此当 $R_2 = 0$ 时，新息协方差矩阵 $V_2$ 是奇异的，全部三个报告序列的似然 {eq}`model2_loglik` 是退化的；{ref}`meas_ex3` 探讨了这种"随机奇异性"。 + ```{code-cell} ipython3 Q2 = K1 @ V1 @ K1.T ε = 1e-6 @@ -996,9 +1076,16 @@ display(Latex(df_to_latex_matrix(df_v2))) ```{code-cell} ipython3 print('Eigenvalues of covariance matrix:') -print(np.sort(np.linalg.eigvalsh(V2))[::-1].round(4)) +print(np.array2string(np.sort(np.linalg.eigvalsh(V2))[::-1], + formatter={'float': '{:.2e}'.format})) ``` +一个特征值很大。 + +第二个特征值为正但极小：它反映了机构对资本存量估计的修正，这些修正仅通过微小系数 $f - 1$ 进入报告数据。 + +第三个特征值等于 $\epsilon$，如果 $R_2 = 0$，由于国民收入恒等式，它将恰好为零。 + 正如 {cite:t}`Sargent1989` 所强调的，尽管两个测量模型共享相同的底层参数，它们却对经济动态产生了截然不同的推断。 ### Wold 脉冲响应 @@ -1074,7 +1161,7 @@ def simulate_series(seed=7909, T=80, k0=10.0): # 通过模型1变换滤波器的滤波数据 xhat_prev = np.array([k0, 0.0]) - z_prev = np.zeros(3) + z_prev = C @ xhat_prev # 不可观测的 z̄_{-1} 的模型隐含值 z_filt = np.empty((T, 3)) k_filt = np.empty(T) @@ -1219,11 +1306,13 @@ plot_true_vs_other(t, sim["k_true"], sim["k_filt"], 事实上，来自模型1的卡尔曼滤波估计消除了大部分测量噪声，并紧密地跟踪真实值。 +改善效果最明显的是投资，它的测量误差最大；而收入本来就测量得很好。 + 在真实模型中，国民收入恒等式 $c_t + \Delta k_t = y_{n,t}$ 精确成立。 独立的测量误差在测量数据中打破了这个会计恒等式。 -卡尔曼滤波器近似地恢复了它。 +由于滤波序列都具有 $C \hat x_t$ 的形式，而 $C$ 的各行满足该恒等式，因此滤波数据精确满足它。 下图通过显示测量和滤波数据两者的残差 $c_t + \Delta k_t - y_{n,t}$ 来证实这一点 @@ -1250,7 +1339,7 @@ plt.tight_layout() plt.show() ``` -正如我们所预测的，测量数据的残差大而波动，而滤波数据的残差在数值上为0。 +正如我们所预测的，测量数据的残差大而波动，而滤波数据的残差在舍入误差范围内为0。 ## 总结 @@ -1262,8 +1351,185 @@ plt.show() 模型1重现了实证加速数文献中记录的**格兰杰因果关系**模式：收入似乎格兰杰引起消费和投资，{cite:t}`Sargent1989` 将这一结果归因于原始报告数据中的测量误差和信号提取。 +对这一模式而言，重要的是每个序列的测量误差相对于真实序列变化的大小。 + 模型2使用滤波数据，将几乎所有方差归因于单个结构冲击 $\theta_t$，并且*无法*重现格兰杰因果关系模式。 {doc}`卡尔曼滤波器 ` 有效地从数据中剥离了测量噪声，因此滤波序列紧密地跟踪真实值。 -原始测量误差打破了国民收入会计恒等式，但接近零的残差表明滤波器近似地恢复了它。 \ No newline at end of file +原始测量误差打破了国民收入会计恒等式，而滤波数据则精确地满足该恒等式，这使得这三个报告序列在随机意义上是退化的（stochastically singular）。 + +## 练习 + +```{exercise} +:label: meas_ex1 + +在假设 $a(L) = 1$ 且 $\beta f = 1$ 的条件下，使用决策规则 {eq}`simple_crule` 和运动规律 $k_{t+1} = k_t + f^{-1} \theta_t$ 证明： + +1. $c_t - c_{t-1} = (1 - \beta)\theta_t$，因此消费是一个鞅 +1. $c_t - \beta c_{t-1} = (1-\beta) y_{nt}$，这就是弗里德曼的消费函数 {eq}`friedman_consumption` +1. $y_{nt} - y_{n,t-1} = \theta_t - \beta \theta_{t-1}$，这就是方程 {eq}`income_process` +1. $(k_{t+1} - k_t) - \beta (k_t - k_{t-1}) = \beta (y_{nt} - y_{n,t-1})$，这就是加速数 {eq}`mm_accelerator` + +然后使用 `sim` 中的真实序列，通过 $\theta_t = y_{nt} - (f-1)k_t$ 恢复 $\theta_t$，并对这四个恒等式进行数值验证。 +``` + +```{solution-start} meas_ex1 +:class: dropdown +``` + +由 $k_t - k_{t-1} = \beta \theta_{t-1}$ 和 {eq}`simple_crule` 可得： + +$$ +c_t - c_{t-1} = (f-1)\beta\theta_{t-1} + (1-\beta)(\theta_t - \theta_{t-1}) = (1-\beta)\theta_t +$$ + +这是因为 $(f-1)\beta = 1 - \beta$。 + +类似地，$y_{nt} - y_{n,t-1} = (f-1)\beta \theta_{t-1} + \theta_t - \theta_{t-1} = \theta_t - \beta\theta_{t-1}$。 + +对于第二个命题，将 $k_t = k_{t-1} + \beta\theta_{t-1}$ 代入 {eq}`simple_crule` 和 {eq}`net_income`： + +$$ +\begin{aligned} +c_t - \beta c_{t-1} & = (f-1)(1-\beta) k_{t-1} + (1-\beta)\theta_t + \left[(f-1)\beta - \beta(1-\beta)\right]\theta_{t-1} \\ +(1-\beta) y_{nt} & = (1-\beta)(f-1) k_{t-1} + (1-\beta)\theta_t + (1-\beta)(f-1)\beta\,\theta_{t-1} +\end{aligned} +$$ + +利用 $(f-1)\beta = 1 - \beta$ 可知两边相等。 + +最后，$k_{t+1} - k_t = \beta\theta_t$ 意味着 $(1 - \beta L)(k_{t+1} - k_t) = \beta(\theta_t - \beta\theta_{t-1}) = \beta (y_{nt} - y_{n,t-1})$。 + +```{code-cell} ipython3 +y, c, dk, k = sim["y_true"], sim["c_true"], sim["dk_true"], sim["k_true"] +θ = y - (f - 1) * k + +print("martingale: ", np.allclose(np.diff(c), (1 - β) * θ[1:])) +print("Friedman: ", np.allclose(c[1:] - β * c[:-1], (1 - β) * y[1:])) +print("MA(1): ", np.allclose(np.diff(y), θ[1:] - β * θ[:-1])) +print("accelerator: ", np.allclose(dk[1:] - β * dk[:-1], β * np.diff(y))) +``` + +```{solution-end} +``` + +```{exercise} +:label: meas_ex2 + +基于 Cholesky 分解的预测误差方差分解取决于变量的排序方式。 + +{cite:t}`Geweke1982` 提出了一种不依赖于排序的格兰杰因果关系度量方法。 + +对于被测量的序列 $i$ 和 $j$，令 $\sigma^2_{j}$ 为仅基于其自身历史得到的 $\bar z_{j,t+1}$ 的一步向前预测误差方差，令 $\sigma^2_{j \mid i}$ 为基于 $\bar z_{i}$ 和 $\bar z_{j}$ 两者历史得到的预测误差方差。 + +那么 $F_{i \to j} = \log(\sigma^2_j / \sigma^2_{j \mid i}) \geq 0$，当且仅当 $\bar z_i$ 不是 $\bar z_j$ 的格兰杰原因时取等号。 + +由于测量误差是相互独立的 AR(1) 过程，两个方差均可以通过将 `steady_state_kalman` 应用于准差分系统 {eq}`model1_transformed`、仅使用与所包含序列对应的 $\bar C$ 和 $R_1$ 的行以及 $W_1$ 的列来计算。 + +1. 在基准测量误差参数下，计算 $F_{y \to \Delta k}$、$F_{\Delta k \to y}$、$F_{y \to c}$ 和 $F_{c \to y}$。 +1. 交换收入和投资的测量误差新息标准差后重新计算，使得 $\sigma_\eta = (0.65, 0.035, 0.05)$。 +1. 使用 $\sigma_\eta = (0.3, 0.3, 0.3)$ 重新计算，并解释为什么相等的测量误差方差**不会**产生对称的格兰杰因果关系。 +``` + +```{solution-start} meas_ex2 +:class: dropdown +``` + +以下是一种解法。 + +```{code-cell} ipython3 +def geweke(σ_η, ρ=np.array([0.6, 0.7, 0.3]), pair=(0, 2)): + """ + Return (F_{i->j}, F_{j->i}) for the measured series pair = (i, j). + """ + D = np.diag(ρ) + Σ = np.diag(np.asarray(σ_η)**2) + C_bar = C @ A - D @ C + R1 = C @ Q @ C.T + Σ + W1 = Q @ C.T + + def pev(idx): + idx = list(idx) + return steady_state_kalman(A, C_bar[idx], Q, + R1[np.ix_(idx, idx)], W1[:, idx])[2] + + i, j = pair + V_i, V_j, V_ij = pev([i])[0, 0], pev([j])[0, 0], pev([i, j]) + return np.log(V_j / V_ij[1, 1]), np.log(V_i / V_ij[0, 0]) + + +cases = {"baseline": (0.05, 0.035, 0.65), + "swapped y and Δk": (0.65, 0.035, 0.05), + "equal": (0.3, 0.3, 0.3)} + +for name, σ in cases.items(): + F_yk, F_ky = geweke(σ, pair=(0, 2)) + F_yc, F_cy = geweke(σ, pair=(0, 1)) + print(f"{name:18s} y→Δk {F_yk:.4f} Δk→y {F_ky:.4f} " + f"y→c {F_yc:.4f} c→y {F_cy:.4f}") +``` + +在基准情形下，收入是投资和消费的格兰杰原因，而两者都不是收入的格兰杰原因。 + +交换收入和投资的测量误差会颠倒两者之间格兰杰因果关系的方向。 + +在测量误差方差相等的情况下，投资作为收入的格兰杰原因的程度超过收入作为投资的格兰杰原因，并且收入强烈地是消费的格兰杰原因。 + +真正重要的不是测量误差的方差本身，而是其相对于每个序列中**信号**的大小。 + +真实投资的新息 $\beta \theta_t$ 的标准差为 $0.95$，而真实消费的新息 $(1-\beta)\theta_t$ 的标准差仅为 $0.048$。 + +因此，新息标准差为 $0.3$ 的测量误差对投资的干扰远小于其对消费的淹没程度。 + +在基准情形下，消费的测量误差新息（$0.035$）在绝对数值上虽小，但相对于真实消费中极小的新息而言并不小，这正是收入成为消费的格兰杰原因的原因。 + +```{solution-end} +``` + +```{exercise} +:label: meas_ex3 + +在模型 2 中，该机构报告 $\hat z_t = C \hat x_t$。 + +1. 证明对于任意估计值 $\hat x_t$，报告的数据都**精确**满足国民收入恒等式 $\hat c_t + \Delta \hat k_t - \hat y_{nt} = 0$。 +1. 解释为什么这意味着当 $R_2 = 0$ 时，报告数据的新息协方差矩阵 $V_2$ 是奇异的。 +1. 对于 $\epsilon = 10^{-2}, 10^{-4}, 10^{-6}$ 的 $R_2 = \epsilon I$，计算 $V_2$ 的特征值并作出解释。 +1. 这种"随机奇异性"对于想要从全部三个报告序列中计算似然 {eq}`model2_loglik` 的计量经济学家意味着什么？ +``` + +```{solution-start} meas_ex3 +:class: dropdown +``` + +$C$ 的各行满足 $C_{c} + C_{\Delta k} - C_{y} = [(f-1) + 0 - (f-1), \; (1 - f^{-1}) + f^{-1} - 1] = [0, 0]$。 + +因此，对于每个 $\hat x_t$ 都有 $\hat c_t + \Delta \hat k_t - \hat y_{nt} = (C_c + C_{\Delta k} - C_y)\hat x_t = 0$。 + +由于向量 $b = [-1, 1, 1]$ 对所有 $t$ 满足 $b\, C \hat x_t = 0$，当 $R_2 = 0$ 时它同样满足 $b\, a_t = 0$，因此 $b V_2 b^\top = 0$，$V_2$ 是奇异的。 + +```{code-cell} ipython3 +print("C_c + C_Δk - C_y =", C[1] + C[2] - C[0]) +print("max |identity residual| in reported data:", + np.max(np.abs(sim['c_filt'] + sim['dk_filt'] - sim['y_filt']))) + +for ε_R in (1e-2, 1e-4, 1e-6): + _, _, V2_ε = steady_state_kalman(A, C, Q2, ε_R * np.eye(3)) + print(f"ε = {ε_R:.0e}: eigenvalues of V2 = " + f"{np.array2string(np.sort(np.linalg.eigvalsh(V2_ε))[::-1], formatter={'float': '{:.2e}'.format})}") +``` + +$V_2$ 的一个特征值较大：它反映了有关永久性冲击 $\theta_t$ 的新信息。 + +第二个特征值很小但为正：它反映了该机构对资本存量估计的修正，这一修正仅通过较小的系数 $f - 1$ 进入报告数据。 + +第三个特征值近似等于 $\epsilon$，并且由于该恒等式的存在，随着 $\epsilon \to 0$ 而趋于零。 + +试图在 $R_2 = 0$ 的情形下计算全部三个报告序列的高斯似然的计量经济学家将面临一个退化分布：$\ln |V_2| = -\infty$。 + +在实践中，必须要么舍弃三个序列中的一个，要么如 {eq}`model2_obs` 所示加入 "第二类" 报告误差，使 $R_2 > 0$。 + +相比之下，在模型 1 的原始数据中，独立的测量误差打破了该恒等式，因此这三个被测量序列具有非奇异的新息协方差矩阵。 + +```{solution-end} +``` From a84b1c9cca24103523071bbf2ed1bee9e300855d Mon Sep 17 00:00:00 2001 From: Matt McKay Date: Thu, 22 Oct 2026 14:44:50 +1100 Subject: [PATCH 2/9] Update translation: .translate/state/measurement_models.md.yml --- .translate/state/measurement_models.md.yml | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/.translate/state/measurement_models.md.yml b/.translate/state/measurement_models.md.yml index b9ef0cf8..191a0de3 100644 --- a/.translate/state/measurement_models.md.yml +++ b/.translate/state/measurement_models.md.yml @@ -1,6 +1,6 @@ -source-sha: bb3310a3ed9e1346211f7daee56187883670b6f1 -synced-at: "2026-07-22" -model: claude-opus-4-8 -mode: NEW -section-count: 7 -tool-version: 0.20.0 +source-sha: 63880c0644afdcad5599295767c212cd34ac0280 +synced-at: "2026-10-22" +model: claude-sonnet-5 +mode: UPDATE +section-count: 8 +tool-version: 0.30.0 From 9357be5ba81dc021cfb6b9300fbac16d2c885652 Mon Sep 17 00:00:00 2001 From: Matt McKay Date: Thu, 22 Oct 2026 14:44:52 +1100 Subject: [PATCH 3/9] Update translation: lectures/svd_intro.md --- lectures/svd_intro.md | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/lectures/svd_intro.md b/lectures/svd_intro.md index b447e5f1..5fb62ae0 100644 --- a/lectures/svd_intro.md +++ b/lectures/svd_intro.md @@ -41,7 +41,7 @@ translation: 在这一讲的后续内容 {doc}`VAR 与 DMD ` 中，我们将介绍如何利用 SVD，借助**动态模态分解**（DMD）来快速计算一阶向量自回归（VAR）模型的低秩近似。 -与主成分分析类似，动态模态分解（DMD）也可以被看作是一种数据降维方法，它通过将数据投影到一组有限的因子上来表示数据中的显著模式。 +与主成分分析类似，DMD 也可以被看作是一种数据降维方法，它通过将数据投影到一组有限的因子上来表示数据中的显著模式。 ## 基本设定 @@ -69,7 +69,7 @@ translation: 在 $m \gg n$ 的情况下，即属性数量 $m$ 远大于个体数量 $n$，且在时间序列环境中 $n$ 等于数据集 $X$ 中所覆盖的时间段数量时，我们将采用不同的方法。 -在后续的 {doc}`VARs 和 DMDs ` 中，我们仍然使用奇异值分解，但目的是构建动态模态分解（DMD）。 +在后续的 {doc}`VARs 和 DMDs ` 中，我们仍然使用奇异值分解，但目的是构建**动态模态分解**（DMD）。 ## 奇异值分解 @@ -133,7 +133,7 @@ $m \times n$ 矩阵 $X$ 的这种结构为构建系统开启了大门 稍后我们还将描述**经济型**或**简化** SVD。 -在研究**简化** SVD之前，我们将进一步讨论**完全** SVD的性质。 +在研究**简化** SVD 之前，我们将进一步讨论**完全** SVD 的性质。 ## 四个基本子空间 @@ -688,7 +688,7 @@ $$ V^\top = \begin{bmatrix}V_1^\top \\V_2^\top \\\vdots\\V_n^\top \end{bmatrix} $$ -在方程{eq}`eq:PCA1`中，每个$m \times n$矩阵$U_{j}V_{j}^\top $显然是秩1的。 +在方程 {eq}`eq:PCA1` 中，每个$m \times n$矩阵$U_{j}V_{j}^\top $显然是秩1的。 因此，我们有 @@ -757,7 +757,7 @@ $$ XX^\top =P\Lambda P^\top . $$ (eq:XXo) -因此数据矩阵$X$可以写成: +因此数据矩阵$X$可以写成： $$ X =\begin{bmatrix}P_1|P_2|\ldots|P_m\end{bmatrix} @@ -783,17 +783,17 @@ X&=\sqrt{\lambda_1}P_1\tilde{\epsilon}_1+\sqrt{\lambda_2}P_2\tilde{\epsilon}_2+\ \end{aligned} $$ -这与SVD分解 +这与 $$ X=\sigma_1U_1{V_1}^{\top}+\sigma_2 U_2{V_2}^{\top}+\ldots+\sigma_{p} U_{p}{V_{p}}^{\top} $$ -是等价的，只要我们令: +是等价的，只要我们令： -* $U_j=P_j$ (变量在第j个主成分上的载荷向量) +* $U_j=P_j$（变量在第$j$个主成分上的载荷向量） -* ${V_k}^{\top}=\tilde{\epsilon}_k$ (第k个主成分) +* ${V_k}^{\top}=\tilde{\epsilon}_k$（第$k$个主成分） 由于对于给定的数据矩阵$X$，计算$P$和$U$有不同的算法，根据所使用的算法，我们可能会得到符号差异或特征向量顺序的不同。 From 8f7d6023db282f49751ed6534899f9ffd8f4b1e1 Mon Sep 17 00:00:00 2001 From: Matt McKay Date: Thu, 22 Oct 2026 14:44:53 +1100 Subject: [PATCH 4/9] Update translation: .translate/state/svd_intro.md.yml --- .translate/state/svd_intro.md.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.translate/state/svd_intro.md.yml b/.translate/state/svd_intro.md.yml index 936cdd37..8b42d8e1 100644 --- a/.translate/state/svd_intro.md.yml +++ b/.translate/state/svd_intro.md.yml @@ -3,4 +3,4 @@ synced-at: "2026-10-22" model: claude-sonnet-5 mode: UPDATE section-count: 12 -tool-version: 0.29.6 +tool-version: 0.30.0 From e41eb939e9fd2d7dc5d1b0c241268eaeb99b3138 Mon Sep 17 00:00:00 2001 From: Matt McKay Date: Thu, 22 Oct 2026 14:44:54 +1100 Subject: [PATCH 5/9] Update translation: lectures/var_dmd.md --- lectures/var_dmd.md | 1057 ++++++++++++++++++++++++++----------------- 1 file changed, 647 insertions(+), 410 deletions(-) diff --git a/lectures/var_dmd.md b/lectures/var_dmd.md index 85ff63c4..4248f6e0 100644 --- a/lectures/var_dmd.md +++ b/lectures/var_dmd.md @@ -4,153 +4,183 @@ jupytext: extension: .md format_name: myst format_version: 0.13 - jupytext_version: 1.16.7 + jupytext_version: 1.10.3 kernelspec: - display_name: Python 3 (ipykernel) + display_name: Python 3 language: python name: python3 translation: title: 向量自回归和动态模态分解 headings: - First-Order Vector Autoregressions: 一阶向量自回归 - Dynamic Mode Decomposition (DMD): 动态模态分解(DMD) - Representation 1: 表示法1 - Representation 2: 表示法 2 - Representation 3: 表示法3 - Representation 3::Decoder of $\check b$ as a linear projection: $\check b$ 的解码器作为线性投影 - Representation 3::An Approximation: 一种近似方法 - Representation 3::Using Fewer Modes: 使用更少的模态 - Source for Some Python Code: Python代码来源 + Overview: 概述 + First-order vector autoregressions: 一阶向量自回归 + First-order vector autoregressions::Computing $X^+$ with an SVD: 用 SVD 计算 $X^+$ + A reduced-rank VAR: 降秩向量自回归 + A reduced-rank VAR::Interpretation as a principal components regression: 作为主成分回归的解释 + The DMD algorithm: DMD 算法 + 'Right eigenvectors: DMD modes': 右特征向量：DMD 模态 + 'Right eigenvectors: DMD modes::Projected modes': 投影模态 + Left eigenvectors and a modal representation: 左特征向量与模态表示 + Left eigenvectors and a modal representation::A tempting alternative that does not work: 一个看似可行但实际不可行的替代方案 + Left eigenvectors and a modal representation::A state-space representation: 状态空间表示 + Long-run responses: 长期响应 + A simulated example: 模拟示例 + Source for some Python code: Python代码来源 + Exercises: 练习 --- # 向量自回归和动态模态分解 -本节讲座中，我们将应用在 {doc}`奇异值分解 ` 中学到的计算方法来研究： +## 概述 -* 一阶向量自回归 (VARs) -* 动态模态分解 (DMDs) -* 一阶向量自回归和动态模态分解之间的联系 +本讲座将我们在 {doc}`奇异值分解 ` 讲座中学到的计算方法应用于 + +* 一阶向量自回归（VAR） +* 动态模态分解（DMD） +* DMD 与一阶 VAR 之间的联系 + +{cite}`sargent2026dynamic` 详细研究了这些联系。 + +我们尤其关注**高瘦型**数据集，即变量数量 $m$ 超过时间周期数量 $n$ 的数据集。 + +这类数据集很常见。 + +例如，{cite}`SSY_CEX_2026` 研究了季度消费者支出调查数据，涉及三个横截面——私人收入、税后转移支付收入和消费——各自的100个分位数，以及总收入增长。 + +这使得他们拥有 $m = 301$ 个变量，但只有 $n = 133$ 个季度观测值。 + +他们使用 DMD 为这301个变量估计了一个秩为 $3$ 的一阶 VAR。 + +在此过程中，我们将学到 + +* DMD 计算的是一阶 VAR 系数矩阵的降秩估计量 +* 该估计量是一种**主成分回归**：将横截面压缩为少数主成分，然后用它们对下一期数据进行回归 +* DMD **模态**是估计系数矩阵的右特征向量 +* 与之匹配的**左**特征向量可以廉价地计算出来，它们与模态一起给出了估计 VAR 的精确模态表示，这对预测和计算长期响应非常有用 + +我们将使用以下导入。 + +```{code-cell} ipython3 +import numpy as np +import matplotlib.pyplot as plt +``` ## 一阶向量自回归 -我们想要拟合一个**一阶向量自回归** +我们希望拟合一个**一阶向量自回归** $$ -X_{t+1} = A X_t + C \epsilon_{t+1}, \quad \epsilon_{t+1} \perp X_t +X_{t+1} = A X_t + C \epsilon_{t+1}, \quad \epsilon_{t+1} \perp X_t $$ (eq:VARfirstorder) -其中，$\epsilon_{t+1}$ 是一个独立同分布的随机向量序列 $m \times 1$ 在时间$t+1$ 的分量，且该序列有零均值向量和单位协方差矩阵；而 $ m \times 1 $ 的向量 $ X_t $ 是： +其中 $\epsilon_{t+1}$ 是均值向量为零、协方差矩阵为单位矩阵的独立同分布 $m \times 1$ 随机向量序列在时间 $t+1$ 处的分量，而 $m \times 1$ 向量 $X_t$ 为 $$ -X_t = \begin{bmatrix} X_{1,t} & X_{2,t} & \cdots & X_{m,t} \end{bmatrix}^\top +X_t = \begin{bmatrix} X_{1,t} & X_{2,t} & \cdots & X_{m,t} \end{bmatrix}^\top $$ (eq:Xvector) -其中 $\cdot ^\top $ 表示共轭转置，$ X_{i,t} $ 是时间 $ t $ 时的变量 $ i $。 +这里 $\cdot^\top$ 表示矩阵转置，$X_{i,t}$ 是时间 $t$ 时的变量 $i$。 + +方程 {eq}`eq:VARfirstorder` 不含常数项。 + +这是因为我们假设每个变量都已经经过变换——例如，通过取差分、减去一个共同趋势、以及减去其时间序列均值——以使 $\{X_t\}$ 可以被视为一个均值为零的协方差平稳过程。 -我们想要拟合方程 {eq}`eq:VARfirstorder`。 +例如，{cite}`SSY_CEX_2026` 从每个对数分位数中减去对数私人收入的横截面均值，使用该均值的增长率作为一个附加变量，然后从所有变量中减去时间序列均值。 -我们的数据则存储在一个 $ m \times (n+1) $ 的矩阵 $ \tilde X $ 中 +我们的数据被组织成一个 $m \times (n+1)$ 矩阵 $\tilde X$ $$ \tilde X = \begin{bmatrix} X_1 \mid X_2 \mid \cdots \mid X_n \mid X_{n+1} \end{bmatrix} $$ -其中对于 $ t = 1, \ldots, n +1 $时，$ m \times 1 $ 的向量 $ X_t $ 由 {eq}`eq:Xvector` 给出。 - -因此，我们想要估计一个系统 {eq}`eq:VARfirstorder`，它由 $ m $ 个最小二乘回归组成，将**所有变量**对**所有变量**的一阶滞后值进行回归。 - -{eq}`eq:VARfirstorder` 的第 $i$ 个方程是将 $X_{i,t+1}$ 对向量 $X_t$ 进行回归。 - -我们按如下步骤进行。 +其中对于 $t = 1, \ldots, n+1$，$m \times 1$ 向量 $X_t$ 由 {eq}`eq:Xvector` 给出。 -从 $ \tilde X $ 中，我们构造以下两个 $m \times n$ 矩阵 +由 $\tilde X$，我们构造两个 $m \times n$ 矩阵 $$ X = \begin{bmatrix} X_1 \mid X_2 \mid \cdots \mid X_{n}\end{bmatrix} $$ -和 +以及 $$ X' = \begin{bmatrix} X_2 \mid X_3 \mid \cdots \mid X_{n+1}\end{bmatrix} $$ -这里的 $ ' $ 是矩阵 $ X' $ 的名称的一部分，并不表示矩阵转置。 +这里 $'$ 是矩阵 $X'$ 名称的一部分，并不表示矩阵转置。 -我们使用 $\cdot^\top $ 来表示矩阵转置或其在复矩阵中的扩展。 +在构造 $X$ 和 $X'$ 时，我们在每种情况下都从 $\tilde X$ 中删除了一列：对于 $X$ 是删除最后一列，对于 $X'$ 是删除第一列。 -在构造 $ X $ 和 $ X' $ 的过程中，我们都从 $ \tilde X $ 中删除了某一列，$ X $ 是删除最后一列，$ X' $ 则是删除第一列。 +我们希望估计一个由 {eq}`eq:VARfirstorder` 构成的系统，该系统由 $m$ 个最小二乘回归组成，即**所有变量**对**所有变量**的一期滞后值的回归。 -显然，$ X $ 和 $ X' $ 都是 $ m \times n $ 的矩阵。 +{eq}`eq:VARfirstorder` 的第 $i$ 个方程是 $X_{i,t+1}$ 对向量 $X_t$ 的回归。 -我们用 $ p \leq \min(m, n) $ 表示 $ X $ 的秩。 +因此，估计量 $\hat A$ 的第 $i$ 行是 $X_{i,t+1}$ 对 $X_{j,t}, j = 1, \ldots, m$ 的回归系数构成的 $1 \times m$ 向量。 -我们感兴趣的两种情况是： +最小二乘估计量 $\hat A$ 求解 -* $ n > > m $，即时间序列观测值的数量 $n$ 远大于变量的数量 $m$ -* $ m > > n $，即变量的数量 $m$ 远大于时间序列观测值的数量 $n$ +$$ +\min_{\check A} \| X' - \check A X \|_F +$$ (eq:ALSeqn) + +其中 $\| \cdot \|_F$ 表示矩阵的弗罗贝尼乌斯范数 -在考虑了这两种特殊情况的一般情况中，有一个通用的公式描述了 $A$ 的最小二乘估计量 $\hat A$。 +$$ + \|B\|_F = \sqrt{ \sum_{i} \sum_{j} |B_{ij}|^2 } . +$$ -但重要的细节有所不同。 +我们用 $p \leq \min(m, n)$ 表示 $X$ 的秩。 -这个通用的公式是： +{eq}`eq:ALSeqn` 的**最小范数**解为 -$$ -\hat A = X' X^+ +$$ +\hat A = X' X^+ $$ (eq:commonA) -其中 $X^+$ 是 $X$ 的广义逆矩阵，或伪逆。 +其中 $X^+$ 是 $X$ 的 [摩尔-彭若斯伪逆](https://en.wikipedia.org/wiki/Moore%E2%80%93Penrose_inverse)。 -关于**穆尔-彭罗斯广义逆矩阵**的详细信息，请参见[穆尔-彭罗斯广义逆矩阵](https://baike.baidu.com/item/穆尔-彭罗斯广义逆矩阵/22770999?fr=aladdin) +我们关注两种情形 -在我们的两种情况下，广义逆矩阵适用的公式有所不同。 +* $n \gg m$，即时间序列观测值 $n$ 远多于变量数 $m$ +* $m \gg n$，即变量数 $m$ 远多于时间序列观测值 $n$ -**短胖型情况：** +**矮胖情形：** -当$n >> m$时，即时间序列的观测值$n$远多于变量的数量$m$，且当$X$具有线性独立的**行**时，$X X^\top$的逆矩阵存在，且广义逆矩阵$X^+$为 +当 $n \gg m$ 且 $X$ 的**行**线性无关时，$X X^\top$ 具有逆矩阵，伪逆 $X^+$ 为 $$ -X^+ = X^\top (X X^\top )^{-1} +X^+ = X^\top (X X^\top )^{-1} $$ -这里$X^+$是一个满足$X X^+ = I_{m \times m}$的**右逆**。 +这里 $X^+$ 是满足 $X X^+ = I_{m \times m}$ 的**右逆**。 -在这种情况下，我们用于估计总体回归系数矩阵$A$的最小二乘估计量的公式{eq}`eq:commonA`就变为 +在这种情形下，问题 {eq}`eq:ALSeqn` 有唯一解，公式 {eq}`eq:commonA` 变为 -$$ +$$ \hat A = X' X^\top (X X^\top )^{-1} $$ (eq:Ahatform101) -这个计算最小二乘回归系数的公式在计量经济学中被广泛地使用。 - -它也被用于估计向量自回归。 +这一最小二乘回归系数公式在计量经济学中被广泛用于估计向量自回归。 -公式{eq}`eq:Ahatform101`的右边，正比于$X_{t+1}$和$X_t$的经验交叉二阶矩矩阵，乘以$X_t$二阶矩阵的逆矩阵。 +公式 {eq}`eq:Ahatform101` 的右边正比于 $X_{t+1}$ 和 $X_t$ 的经验交叉二阶矩矩阵，乘以 $X_t$ 的二阶矩矩阵的逆。 -**高瘦型情况：** +**瘦高情形：** -当$m > > n$时，即属性数量$m$远大于时间序列的观测值$n$，且当$X$的**列**线性独立时，$X^\top X$的逆矩阵存在，且广义逆矩阵$X^+$为 +当 $m \gg n$ 且 $X$ 的**列**线性无关时，$X^\top X$ 具有逆矩阵，伪逆 $X^+$ 为 $$ -X^+ = (X^\top X)^{-1} X^\top +X^+ = (X^\top X)^{-1} X^\top $$ -这里$X^+$是一个满足$X^+ X = I_{n \times n}$的**左逆**。 +这里 $X^+$ 是满足 $X^+ X = I_{n \times n}$ 的**左逆**。 -在这种情况下，我们用于估计$A$的最小二乘估计公式{eq}`eq:commonA`变为 +在这种情形下，公式 {eq}`eq:commonA` 变为 $$ -\hat A = X' (X^\top X)^{-1} X^\top +\hat A = X' (X^\top X)^{-1} X^\top $$ (eq:hatAversion0) -请比较{eq}`eq:Ahatform101`和{eq}`eq:hatAversion0`中$\hat A$的表达式。 - -这里我们特别关注公式{eq}`eq:hatAversion0`。 - -$\hat A$的第$i$行是一个$m \times 1$的向量，其中包含了$X_{i,t+1}$对$X_{j,t}, j = 1, \ldots, m$回归的系数。 - -如果我们使用公式{eq}`eq:hatAversion0`来计算$\hat A X$，我们发现 +如果我们使用公式 {eq}`eq:hatAversion0` 来计算 $\hat A X$，会发现 $$ \hat A X = X' @@ -158,598 +188,805 @@ $$ 因此回归方程**完全拟合**。 -这是**欠定最小二乘**模型中典型的结果。 +这是**欠定最小二乘**模型中的典型结果。 -再次重申，**高瘦**情况(见{doc}`奇异值分解`)指观测的数量$n$相对于向量$X_t$属性的数量$m$较小时，我们想要拟合方程{eq}`eq:VARfirstorder`。 +事实上，每一个满足 $N X = 0$ 的矩阵 $\check A = \hat A + N$ 也能完全拟合，因此问题 {eq}`eq:ALSeqn` 有无穷多个解。 -我们面临着最小二乘估计量是欠定的，且回归方程完全拟合的事实。 +公式 {eq}`eq:hatAversion0` 选取的是弗罗贝尼乌斯范数最小的那个解。 -接下来，我们想要更加高效地计算广义逆矩阵$X^+$。 +完美的样本内拟合是一个警示信号：当系数数量多于观测值数量时，估计量 $\hat A$ 不仅在追踪信号，也在追踪噪声。 -广义逆矩阵$X^+$将是我们$A$估计量的一个组成部分。 +我们将通过构造一个不再完全拟合的**降秩**估计量来应对这一问题。 -作为对$A$的估计量$\hat A$，我们想要形成一个$m \times m$的矩阵，来解决最小二乘最佳拟合问题 +### 用 SVD 计算 $X^+$ -$$ -\hat A = \textrm{argmin}_{\check A} || X' - \check A X ||_F -$$ (eq:ALSeqn) +计算伪逆 $X^+$ 的一种高效方法是从**约化**奇异值分解开始 + +$$ +X = U \Sigma V^\top +$$ (eq:SVDDMD) + +其中 $U$ 是 $m \times p$ 矩阵，$\Sigma$ 是由正奇异值 $\sigma_1 \geq \sigma_2 \geq \cdots \geq \sigma_p > 0$ 构成的 $p \times p$ 对角矩阵，$V$ 是 $n \times p$ 矩阵。 -其中 $|| \cdot ||_F$ 表示矩阵的Frobenius（或欧几里得）范数。 +正如我们在 {doc}`奇异值分解 ` 中所看到的，对于约化 SVD -Frobenius范数定义为 +* $U^\top U = I_{p \times p}$ 且 $V^\top V = I_{p \times p}$，但是 +* 除非 $p = m$，否则 $U U^\top \neq I_{m \times m}$；除非 $p = n$，否则 $V V^\top \neq I_{n \times n}$。 + +对于任意秩 $p$，$X$ 的摩尔-彭若斯伪逆为 + +$$ +X^{+} = V \Sigma^{-1} U^\top +$$ (eq:Xplusformula) + +当 $X$ 的列线性无关，即 $p = n$ 且 $V$ 为方正交矩阵时，我们可以确认公式 {eq}`eq:Xplusformula` 与我们先前的公式一致： $$ - ||A||_F = \sqrt{ \sum_{i=1}^m \sum_{j=1}^m |A_{ij}|^2 } +\begin{aligned} +(X^\top X)^{-1} X^\top + & = (V \Sigma U^\top U \Sigma V^\top )^{-1} V \Sigma U^\top \\ + & = (V \Sigma^2 V^\top )^{-1} V \Sigma U^\top \\ + & = V \Sigma^{-2} V^\top V \Sigma U^\top \\ + & = V \Sigma^{-1} U^\top +\end{aligned} $$ -方程{eq}`eq:ALSeqn`右侧的最小值解为 +其中第三行使用了 $V^{-1} = V^\top$。 + +将公式 {eq}`eq:Xplusformula` 代入公式 {eq}`eq:commonA` 得到 $$ -\hat A = X' X^{+} -$$ (eq:hatAform) +\hat A = X' V \Sigma^{-1} U^\top +$$ (eq:AhatSVDformula) -其中（可能是巨大的）$ n \times m $ 的矩阵 $ X^{+} = (X^\top X)^{-1} X^\top $ 同样是 $ X $ 的广义逆矩阵。 +公式 {eq}`eq:AhatSVDformula` 显示出 $\hat A$ 的一个弱点。 -对于我们感兴趣的一些情况，$X^\top X $ 可能接近奇异，这种情况会使某些数值算法变得不准确。 +较小的奇异值 $\sigma_j$ 通过 $1/\sigma_j$ 的形式进入公式。 -为了应对这种可能性，我们将使用高效的算法来构建公式{eq}`eq:hatAversion0`中 $\hat A$ 的**降秩近似**。 +因此，数据 $X$ 变化非常小的方向会获得非常大的权重，这会放大噪声。 -这种近似方式，让我们的向量自回归估计不再完全拟合。 +## 降秩向量自回归 -$ \hat A $ 的第 $ i $ 行是一个 $ m \times 1 $ 的回归系数向量，表示 $ X_{i,t+1} $ 对 $ X_{j,t}, j = 1, \ldots, m $ 的回归。 +一种补救方法是只保留 $r < p$ 个最大的奇异值。 -一种高效计算广义逆矩阵$X^+$的方式是从奇异值分解开始 +设 $U_r$ 为由 $U$ 的前 $r$ 列组成的 $m \times r$ 矩阵，设 $\Sigma_r$ 为 $\Sigma$ 左上角的 $r \times r$ 分块，设 $V_r$ 为由 $V$ 的前 $r$ 列组成的 $n \times r$ 矩阵。 + +我们对 $A$ 的 **秩-$r$ 估计量** 为 $$ -X = U \Sigma V^\top -$$ (eq:SVDDMD) +\hat A_r = X' V_r \Sigma_r^{-1} U_r^\top +$$ (eq:Ahat_r) + +$m \times m$ 矩阵 $\hat A_r$ 的秩至多为 $r$。 + +注意，$V_r \Sigma_r^{-1} U_r^\top$ 是 $\hat X_r = U_r \Sigma_r V_r^\top$ 的伪逆，而 $\hat X_r$ 正是 {doc}`奇异值分解 ` 中埃卡特-杨定理所描述的 $X$ 的最佳秩-$r$ 近似。 + +于是 $\hat A_r = X' \hat X_r^+$：我们通过在最小二乘公式 $\hat A = X' X^+$ 中用 $X$ 的埃卡特-杨近似代替 $X$，得到了秩-$r$ 估计量。 + +这正是 **动态模态分解** 所计算的估计量。 -这里我们提醒自己，这个**简化**SVD中，$X$是一个$m \times n$的数据矩阵，$U$是一个$m \times p$的矩阵，$\Sigma$是一个$p \times p$的矩阵，而$V$是一个$n \times p$的矩阵。 +动态模态分解由 {cite}`schmid2010` 引入。 -通过以下一系列等式，我们可以有效地构造相关的广义逆矩阵$X^+$。 +你可以在 {cite}`DMD_book` 和 {cite}`Brunton_Kutz_2019`（第7.2节）中阅读相关内容。 + +与 $\hat A_r$ 相关的拟合残差及其协方差矩阵为 $$ -\begin{aligned} -X^{+} & = (X^\top X)^{-1} X^\top \\ - & = (V \Sigma U^\top U \Sigma V^\top )^{-1} V \Sigma U^\top \\ - & = (V \Sigma \Sigma V^\top )^{-1} V \Sigma U^\top \\ - & = V \Sigma^{-1} \Sigma^{-1} V^\top V \Sigma U^\top \\ - & = V \Sigma^{-1} U^\top -\end{aligned} -$$ (eq:efficientpseudoinverse) +\hat a_{t+1} = X_{t+1} - \hat A_r X_t, \qquad +\hat \Omega = \frac{1}{n} \sum_{t=1}^{n} \hat a_{t+1} \hat a_{t+1}^\top += \frac{1}{n} (X' - \hat A_r X)(X' - \hat A_r X)^\top +$$ (eq:residualsOmega) -（由于$m > > n$，在简化SVD中$V^\top V = I_{p \times p}$，因此我们可以将前面的一系列等式同时用于简化SVD和完整SVD。） +因此，我们估计得到的秩-$r$ 向量自回归为 -因此，我们将使用方程{eq}`eq:SVDDMD`中$X$的奇异值分解来构造$X$的广义逆矩阵$X^+$，计算方法为： +$$ +X_{t+1} = \hat A_r X_t + \hat a_{t+1} +$$ (eq:reducedVAR) + +当 $r < p$ 时，回归不再完全拟合，$\hat \Omega$ 也不再为零。 + +在应用中，$r$ 通常非常小，例如三个或更少。 + +选择 $r$ 的一种常见方法是绘制奇异值 $\sigma_1, \sigma_2, \ldots$ 的图像，并在它们不再快速下降而趋于平稳的位置进行截断——即所谓的"肘部法则"。 + +{cite}`SSY_CEX_2026` 使用该法则选择了 $r = 3$。 + +### 作为主成分回归的解释 + +我们可以从统计学角度解释公式 {eq}`eq:Ahat_r`。 + +定义 $r \times 1$ 向量 $$ -X^{+} = V \Sigma^{-1} U^\top -$$ (eq:Xplusformula) +f_t = U_r^\top X_t +$$ -其中矩阵$\Sigma^{-1}$是通过将$\Sigma$中的每个非零元素替换为$\sigma_j^{-1}$构造而成。 +其各分量是截面数据 $X_t$ 的前 $r$ 个 **主成分**（参见 {doc}`奇异值分解 ` 中关于主成分分析的讨论）。 -我们可以将公式{eq}`eq:Xplusformula`与公式{eq}`eq:hatAform`结合使用来计算回归系数矩阵$\hat A$。 +将它们堆叠起来可得到 $r \times n$ 矩阵 $F = U_r^\top X = \Sigma_r V_r^\top$。 -因此，我们对$m \times m$的系数矩阵$A$的估计量$\hat A = X' X^+$为： +$X_{t+1}$ 对 $f_t$ 的最小二乘回归具有系数矩阵 $$ -\hat A = X' V \Sigma^{-1} U^\top -$$ (eq:AhatSVDformula) +X' F^\top (F F^\top)^{-1} = X' V_r \Sigma_r (\Sigma_r^2)^{-1} = X' V_r \Sigma_r^{-1} +$$ -## 动态模态分解(DMD) +因此 -我们转向与**动态模态分解**相关的$m >>n$ **高瘦型**情况。 +$$ +\hat A_r X_t = \left( X' V_r \Sigma_r^{-1} \right) f_t . +$$ -这里一个$ m \times n+1 $ 的数据矩阵 $ \tilde X $ 包含了比时间周期 $ n+1 $ 多得多的属性（或变量）$ m $。 +类似地，$f_{t+1} = U_r^\top X_{t+1}$ 对 $f_t$ 的最小二乘回归具有 $r \times r$ 系数矩阵 -动态模态分解由{cite}`schmid2010`首次提出， +$$ +\tilde A = U_r^\top X' V_r \Sigma_r^{-1} +$$ (eq:Atilde) -你可以在 {cite}`DMD_book` 和 {cite}`Brunton_Kutz_2019`（第7.2节）中阅读有关动态模态分解的内容。 +因此 $\tilde A$ 是 $r$ 个主成分的一阶向量自回归系数矩阵。 -**动态模态分解**（DMD）计算公式{eq}`eq:AhatSVDformula`所描述的最小二乘回归系数$\hat A$的一个秩$r < p$的近似。 +注意 $\tilde A = U_r^\top \hat A_r U_r$。 -我们将逐步构建一种适合应用的表示。 +因此，动态模态分解分两步进行： -我们将通过三种不同的表示方式来描述一阶线性动态系统（即我们的向量自回归），从而实现这一点。 +1. 利用 $X$ 本身的变异，将截面数据 $X_t$ 压缩为 $r$ 个主成分 $f_t$ +1. 将下一期数据 $X_{t+1}$ 对这些主成分进行回归 -**三种表示的指南：** 在实践中，我们主要关注表示3。 +统计学家将此称为 **主成分回归**。 -我们使用前两种表示来呈现一些有用的中间推导，这些步骤有助于我们理解表示3的内在原理。 +它类似于 {cite}`stock_watson2002` 的"扩散指数"预测方法。 -在应用中，我们将只使用**DMD模态**的一小部分子集来近似动态。 +```{note} +估计量 $\hat A_r$ **并非** 使 $\| X' - \check A X \|_F$ 最小化的秩-$r$ 矩阵。 -我们使用这样一个小的DMD模态子集来构建对$A$的降秩近似。 +那个最小化矩阵是 {cite}`anderson1951` 所提出的 **降秩回归** 估计量。 -为此，我们需要使用与表示法3相关的**简化**SVD，而不是与表示法1和2相关的**完整**SVD。 +它选取能够最好地*预测* $X_{t+1}$ 的 $r$ 个方向，而 $\hat A_r$ 则使用能够最好地*描述* $X_t$ 的 $r$ 个方向。 -**给不耐烦的读者的指南：** 在我们的应用中，我们将使用表示法3。 +{ref}`var_dmd_ex2` 比较了这两种估计量。 +``` -第一次阅读时，您可以跳过铺垫性的表示法1和2。 +## DMD 算法 -+++ +下面是 DMD 算法。 -## 表示法1 +1. 计算约简奇异值分解 $X = U \Sigma V^\top$。 +1. 保留 $r$ 个最大的奇异值，从而构成 $U_r$、$\Sigma_r$ 和 $V_r$。 +1. 构成 $r \times r$ 矩阵 $\tilde A = U_r^\top X' V_r \Sigma_r^{-1}$ 并计算其特征分解 -在这个表示法中，我们将使用$X$的**完整**SVD。 + $$ + \tilde A W = W \Lambda + $$ (eq:tildeAeigenred) -我们使用$U$的$m$个**列**，即$U^\top$的$m$个**行**，来定义一个$m \times 1$的向量$\tilde b_t$： + 其中 $\Lambda$ 是由特征值 $\lambda_1, \ldots, \lambda_r$ 组成的对角矩阵，$W$ 的列是对应的特征向量。 +1. 构成 $m \times r$ 矩阵 $\Phi$ 和 $r \times m$ 矩阵 $\Psi$ -$$ -\tilde b_t = U^\top X_t . -$$ (eq:tildeXdef2) + $$ + \Phi = X' V_r \Sigma_r^{-1} W, \qquad \Psi = (W \Lambda)^{-1} U_r^\top + $$ (eq:PhiPsi) -原始数据$X_t$可以表示为： +$\Lambda$ 对角线上的特征值 $\lambda_i$ 称为 **DMD 特征值**，$\Phi$ 的列 $\phi_i$ 称为 **DMD 模态**。 -$$ -X_t = U \tilde b_t -$$ (eq:Xdecoder) +我们将在下文看到，$\Psi$ 的行 $\psi_i$ 是与之匹配的**左**特征向量。 -（这里我们使用$b$来提醒自己我们正在创建一个**基**向量。） +请注意，其中的主要计算工作在于对 $m \times n$ 矩阵 $X$ 进行奇异值分解，以及对较小的 $r \times r$ 矩阵 $\tilde A$ 进行特征分解。 -由于我们现在使用的是**完全**SVD，$U U^\top = I_{m \times m}$。 +我们完全不需要构造或分解 $m \times m$ 矩阵。 -因此从方程{eq}`eq:tildeXdef2`可以得出，我们可以用$\tilde b_t$重新构造$X_t$。 +在整个过程中，我们假设 $\tilde A$ 具有 $r$ 个线性无关的特征向量和非零特征值，这是一般情形。 -特别地， +矩阵 $\tilde A$ 是实矩阵但不是对称矩阵，因此某些特征值可能以共轭复数对的形式出现。 -* 方程 {eq}`eq:tildeXdef2` 作为一个**编码器**，将 $m \times 1$ 向量 $X_t$ **旋转**成一个 $m \times 1$ 的向量 $\tilde b_t$ - -* 方程 {eq}`eq:Xdecoder` 作为一个**解码器**，通过旋转 $m \times 1$ 向量 $\tilde b_t$ 来**重新构造** $m \times 1$ 的向量 $X_t$ +在这种情况下，$\Phi$ 的相应列和 $\Psi$ 的相应行也会以共轭复数对的形式出现，而诸如 $\Phi \Lambda^j \Psi$ 之类的乘积则是实数。 -为 $m \times 1$ 的基向量 $\tilde b_t$ 定义一个转移矩阵： +对于复矩阵，我们用 $\cdot^*$ 表示共轭转置。 -$$ -\tilde A = U^\top \hat A U -$$ (eq:Atilde0) +```{note} +{cite}`tu_Rowley` 以不同的方式对模态进行缩放，将其定义为 $\Phi \Lambda^{-1}$ 的列。 -我们可以通过以下方式表示 $\hat A$： +在这种缩放下，与之匹配的左特征向量是 $\Lambda \Psi = W^{-1} U_r^\top$ 的行。 -$$ -\hat A = U \tilde A U^\top -$$ +对某一模态及其左特征向量按相反方向进行重新缩放，不会改变下文中的任何公式。 +``` + +## 右特征向量：DMD 模态 -$m \times 1$ 的基向量 $\tilde b_t$ 的动态由以下方程支配： +{cite}`tu_Rowley` 中的以下结果解释了为什么 $\Phi$ 的列被称为模态。 + +```{prf:proposition} +:label: prop-dmd-eigenvectors + +$\Phi$ 的列是 $\hat A_r$ 的特征向量： $$ -\tilde b_{t+1} = \tilde A \tilde b_t +\hat A_r \Phi = \Phi \Lambda $$ +``` -为了构建基于 $X_1$ 的 $X_t$ 未来值的预测 $\overline X_t$，我们可以对这个方程的两边应用解码器（即旋转器），从而推导出： +```{prf:proof} +利用公式 {eq}`eq:Ahat_r`、{eq}`eq:Atilde` 和 {eq}`eq:PhiPsi`， $$ -\overline X_{t+1} = U \tilde A^t U^\top X_1 +\begin{aligned} + \hat A_r \Phi & = (X' V_r \Sigma_r^{-1} U_r^\top) (X' V_r \Sigma_r^{-1} W) \\ + & = X' V_r \Sigma_r^{-1} \tilde A W \\ + & = X' V_r \Sigma_r^{-1} W \Lambda \\ + & = \Phi \Lambda +\end{aligned} $$ +``` + +对各列取相等关系，可得 $\hat A_r \phi_i = \lambda_i \phi_i$（$i = 1, \ldots, r$）。 -这里我们用 $\overline X_{t+1}, t \geq 1$ 表示预测值。 +我们还可以进一步说明。 -+++ +将 $\hat A_r$ 写为 $B U_r^\top$，其中 $B = X' V_r \Sigma_r^{-1}$ 是 $m \times r$ 矩阵。 -## 表示法 2 +因为 $B U_r^\top$ 的非零特征值与 $U_r^\top B = \tilde A$ 的特征值一致，所以 $\hat A_r$ 的特征值为 $\lambda_1, \ldots, \lambda_r$，再加上 $m - r$ 个零特征值。 -这种表示方法与{cite}`schmid2010`最初提出的方法有关。 +因此，估计的 VAR {eq}`eq:reducedVAR` 稳定当且仅当 $\max_i |\lambda_i| < 1$，该条件可以通过一个 $r \times r$ 矩阵来验证。 -它可以被视为推导后面将要介绍的表示3的一个中间步骤。 +（另见 {cite}`DDSE_book`，第 238 页。） -与表示1一样，我们继续： +### 投影模态 -* 使用**完整**SVD而**不是**简化SVD +{cite}`schmid2010` 最初使用的是**投影模态**（projected modes）$U_r W$，而非 $\Phi$。 -正如我们在{doc}`奇异值分解 ` 讲座中观察和阐释的那样 +因为 $U_r^\top \Phi = \tilde A W = W \Lambda$，所以有 - * (a) 对于完整SVD，$U U^\top = I_{m \times m}$和$U^\top U = I_{p \times p}$都是单位矩阵 - - * (b) 对于$X$的简化SVD，$U^\top U$不是单位矩阵。 +$$ +U_r U_r^\top \Phi = U_r W \Lambda +$$ -正如我们后面将会看到的，完整SVD对于我们最终想要做的事情来说过于局限，也就是说，我们需要应对由于使用$X$的简化SVD而导致$U^\top U$**不是**单位矩阵的情形。 +因此，在尺度意义下，投影模态是模态 $\Phi$ 向 $U_r$ 的列空间的正交投影。 -但现在，让我们假设我们使用的是完整SVD，这样条件(a)和(b)都得到满足。 +只有当 $\Phi = U_r W \Lambda$ 时，即 $\Phi$ 的列位于 $U_r$ 的列空间中时，投影模态才是 $\hat A_r$ 的特征向量。 -对方程{eq}`eq:Atilde0`中定义的$m \times m$矩阵$\tilde A = U^\top \hat A U$进行特征分解： +但这一点未必成立，因为 $\hat A_r U_r W = \Phi$。 -$$ -\tilde A = W \Lambda W^{-1} -$$ (eq:tildeAeigen) +{cite}`tu_Rowley` 所提出的构造 $\Phi = \hat A_r U_r W$ 在所有情况下都能得到 $\hat A_r$ 的特征向量。 -其中$\Lambda$是特征值的对角矩阵，$W$是一个$m \times m$的矩阵，其每一列都对应于$\Lambda$中行(特征值)的特征向量。 +## 左特征向量与模态表示 -当$U U^\top = I_{m \times m}$时（这在$X$的完全SVD中是成立的），可得： +现在我们来看 $\Psi$ 的各行。 -$$ -\hat A = U \tilde A U^\top = U W \Lambda W^{-1} U^\top -$$ (eq:eqeigAhat) +```{prf:proposition} +:label: prop-dmd-left -根据方程{eq}`eq:eqeigAhat`，对角矩阵$\Lambda$包含$\hat A$的特征值，而$\hat A$对应的特征向量是矩阵$UW$的列。 +设 $\Psi = (W \Lambda)^{-1} U_r^\top$。则 -因此，我们的一阶向量自回归所捕获的$X_t$动态的系统部分（即非随机部分）可以描述为： +1. $\Psi \Phi = I_{r \times r}$ +1. $\Psi \hat A_r = \Lambda \Psi$，因此 $\Psi$ 的各行是 $\hat A_r$ 的左特征向量 +1. $\hat A_r = \Phi \Lambda \Psi = \sum_{i=1}^r \lambda_i \phi_i \psi_i$ +``` + +```{prf:proof} +对于第一个命题， $$ -X_{t+1} = U W \Lambda W^{-1} U^\top X_t +\Psi \Phi = \Lambda^{-1} W^{-1} U_r^\top X' V_r \Sigma_r^{-1} W += \Lambda^{-1} W^{-1} \tilde A W = \Lambda^{-1} \Lambda = I 。 $$ -将上述方程两边同时乘以$W^{-1} U^\top $，得到： +对于第二个命题，由于 $U_r^\top \hat A_r = U_r^\top X' V_r \Sigma_r^{-1} U_r^\top = \tilde A U_r^\top$， -$$ -W^{-1} U^\top X_{t+1} = \Lambda W^{-1} U^\top X_t +$$ +\Psi \hat A_r = \Lambda^{-1} W^{-1} \tilde A U_r^\top += \Lambda^{-1} \Lambda W^{-1} U_r^\top = W^{-1} U_r^\top = \Lambda \Psi 。 $$ -或 +对于第三个命题， $$ -\hat b_{t+1} = \Lambda \hat b_t +\Phi \Lambda \Psi = X' V_r \Sigma_r^{-1} W \Lambda \Lambda^{-1} W^{-1} U_r^\top += X' V_r \Sigma_r^{-1} U_r^\top = \hat A_r 。 +$$ +``` + +定义 $r \times 1$ 的**模态坐标**向量 + $$ +b_t = \Psi X_t = (W \Lambda)^{-1} U_r^\top X_t +$$ (eq:modalcoords) + +计算 $b_t$ 只需要主成分 $f_t = U_r^\top X_t$ 以及一个 $r \times r$ 矩阵。 -其中，我们的**编码器**是 +由于 $\Psi \Phi = I$，对 $\hat A_r = \Phi \Lambda \Psi$ 进行迭代可得 -$$ -\hat b_t = W^{-1} U^\top X_t $$ +\hat A_r^j = \Phi \Lambda^j \Psi, \quad j \geq 1 +$$ (eq:Ahatpower) -我们的**解码器**是 +因此，所估计的向量自回归模型 {eq}`eq:reducedVAR` 蕴含 $$ -X_t = U W \hat b_t +b_{t+1} = \Lambda b_t + \Psi \hat a_{t+1} $$ -我们可以使用这种表示来构建一个基于$X_1$的$X_{t+1}$的预测器$\overline X_{t+1}$： +且在给定 $X_t$ 的条件下，对 $X_{t+j}$ 的 $j$ 步超前预测为 $$ -\overline X_{t+1} = U W \Lambda^t W^{-1} U^\top X_1 -$$ (eq:DSSEbookrepr) +\hat A_r^j X_t = \Phi \Lambda^j b_t = \sum_{i=1}^r \phi_i \lambda_i^j b_{i,t} +$$ (eq:modalforecast) -实际上， -{cite}`schmid2010`定义了一个$m \times m$的矩阵$\Phi_s$为 +每个模态 $\phi_i$ 贡献的分量，其振幅以速率 $\lambda_i$ 按几何级数衰减。 -$$ -\Phi_s = UW -$$ (eq:Phisfull) +矩阵 $\Phi \Psi$ 是向 $\Phi$ 的列空间上的斜投影。 -和一个广义逆矩阵 +它把 $X_t$ 分解为 $\Phi b_t$ 加上一个余项 $(I - \Phi \Psi) X_t$，该余项会被 $\hat A_r$ 映射为零，因为 $\hat A_r (I - \Phi \Psi) = \Phi \Lambda \Psi - \Phi \Lambda \Psi \Phi \Psi = 0$。 -$$ -\Phi_s^+ = W^{-1}U^\top -$$ (eq:Phisfullinv) +### 一个看似可行但实际不可行的替代方案 -{cite}`schmid2010`随后将方程{eq}`eq:DSSEbookrepr`表示为 +摩尔-彭若斯伪逆 $\Phi^+ = (\Phi^* \Phi)^{-1} \Phi^*$ 也满足 $\Phi^+ \Phi = I$。 -$$ -\overline X_{t+1} = \Phi_s \Lambda^t \Phi_s^+ X_1 -$$ (eq:schmidrep) +向量 $\check b_t = \Phi^+ X_t$ 是 $X_t$ 对 $\Phi$ 各列进行最小二乘回归所得的系数向量，而 $\Phi \check b_t$ 则是 $X_t$ 在 $\Phi$ 的列空间上的正交投影。 -基向量的分量$ \hat b_t = W^{-1} U^\top X_t \equiv \Phi_s^+ X_t$是 -DMD**投影模态**。 +这使得 $\check b_t$ 成为描述每个模态在 $X_t$ 中所占比重的一种自然的**描述性**概括。 -要理解为什么它们被称为**投影模态**，注意到 +但一般而言 -$$ -\Phi_s^+ = ( \Phi_s^\top \Phi_s)^{-1} \Phi_s^\top +$$ +\Phi \Lambda \Phi^+ \neq \hat A_r $$ -所以 $m \times p$ 的矩阵 +因此在所估计的向量自回归模型下，$\check b_t$ 的演化**并不**遵循 $\check b_{t+1} = \Lambda \check b_t$，并且 $\Phi \Lambda^j \Phi^+ X_t$ 也不是向量自回归模型的预测值 $\hat A_r^j X_t$。 -$$ -\hat b = \Phi_s^+ X -$$ +原因在于 $\Phi^+$ 和 $\Psi$ 都是 $\Phi$ 的左逆，但它们各自的行所张成的空间不同： -是 $m \times n$ 矩阵 $X$ 在 $m \times p$ 矩阵 $\Phi_s$ 上的回归系数矩阵。 +* $\Phi^+$ 的各行张成的空间与 $\Phi^*$ 的各行张成的空间相同 +* $\Psi$ 的各行张成的空间与 $U_r^\top$ 的各行张成的空间相同 -在讨论表示法3时，我们将在相关背景下进一步讨论这个解释，表示法3是由 Tu 等人 {cite}`tu_Rowley` 提出的。 +这两个左逆重合，当且仅当 $\Phi$ 与 $U_r$ 的列空间相同。 -当我们想要使用简化SVD时（这在实践中经常出现），使用表示法3更为合适。 +$\Phi$ 的列空间是 $X' V_r$ 的列空间，它描述了数据在**下一期**会走向何方。 -## 表示法3 +$U_r$ 的列空间则由**当前**数据 $X$ 变化最剧烈的 $r$ 个方向张成。 -与构建表示法1和表示法2的程序不同（它们都使用了**完整**SVD），我们现在使用**简化**SVD。 +这两个空间在特殊情形下会重合，例如当无噪声数据恰好由 $r$ 个因子驱动时（见 {ref}`var_dmd_ex1`），但一般情况下并不重合。 -同样，令 $p \leq \textrm{min}(m,n)$ 为 $X$ 的秩。 +```{note} +公式 $(W \Lambda)^{-1} U_r^\top X_1$ 出现在 {cite}`DDSE_book`（第240页）中，在该文献中它被作为 $\Phi^+ X_1$ 的一种计算成本较低的近似而给出。 -构造一个**简化**SVD +{prf:ref}`prop-dmd-left` 表明，它恰恰是与所估计的向量自回归模型的动态一致的坐标向量。 +``` -$$ -X = \tilde U \tilde \Sigma \tilde V^\top , -$$ +### 状态空间表示 -其中现在 $\tilde U$ 是 $m \times p$ 的矩阵，$\tilde \Sigma$ 是 $p \times p$ 的矩阵，而 $\tilde V^\top$ 是 $p \times n$ 的矩阵。 +借助模态表示，我们可以将所估计的向量自回归模型写成一个具有 $r \times 1$ 状态向量的线性状态空间系统。 -我们的 $A$ 的最小范数最小二乘近似器现在的表示为 +定义 $$ -\hat A = X' \tilde V \tilde \Sigma^{-1} \tilde U^\top -$$ (eq:Ahatwithtildes) +\hat x_t = \Lambda \Psi X_{t-1} = W^{-1} U_r^\top X_{t-1} = W^{-1} f_{t-1} +$$ +那么由 $X_t = \hat A_r X_{t-1} + \hat a_t = \Phi \Lambda \Psi X_{t-1} + \hat a_t$ 以及 $\hat x_{t+1} = \Lambda \Psi X_t$，可得 -**计算$\hat A$的主要特征向量** +$$ +\begin{aligned} +\hat x_{t+1} & = \Lambda \hat x_t + \Lambda \Psi \hat a_t \\ +X_t & = \Phi \hat x_t + \hat a_t +\end{aligned} +$$ (eq:dmd_statespace) -我们首先参照构建表示法1时使用的步骤，通过以下方式为旋转的$p \times 1$状态$\tilde b_t$定义一个转移矩阵： +其中我们利用了 $\Psi \Phi = I$。 -$$ -\tilde A =\tilde U^\top \hat A \tilde U -$$ (eq:Atildered) +在表示 {eq}`eq:dmd_statespace` 中 +* 状态 $\hat x_t$ 是以主成分向量自回归模型 $\tilde A$ 的特征向量基表示的滞后主成分向量 +* 转移矩阵 $\Lambda$ 是对角矩阵，因此状态的 $r$ 个分量各自独立演化，只是冲击之间存在相关性 +* $\Phi$ 的各列是 $m$ 个观测变量在状态的 $r$ 个分量上的**载荷** +* 一个 $m \times 1$ 的残差向量 $\hat a_t$ 同时驱动两个方程；根据最小二乘正规方程，它在样本内与 $f_{t-1}$ 正交，因而也与 $\hat x_t$ 正交 -**作为投影系数的解释** +{cite}`sargent2026dynamic` 将动态模态分解与此类线性状态空间模型联系了起来。 +## 长期响应 -{cite}`DDSE_book`指出$\tilde A$可以被理解为$\hat A$在$\tilde U$中$p$个模态上的投影。 +模态表示 {eq}`eq:Ahatpower` 也简化了长期响应的计算。 -要验证这一点，首先注意到，由于$ \tilde U^\top \tilde U = I$，因此： +如果对所有 $i$ 都有 $|\lambda_i| < 1$，则 $$ -\tilde A = \tilde U^\top \hat A \tilde U = \tilde U^\top X' \tilde V \tilde \Sigma^{-1} \tilde U^\top \tilde U -= \tilde U^\top X' \tilde V \tilde \Sigma^{-1} \tilde U^\top -$$ (eq:tildeAverify) - +\sum_{j=1}^\infty \hat A_r^j = \Phi \Lambda (I - \Lambda)^{-1} \Psi +$$ -接下来，我们将使用标准最小二乘公式计算$\hat A$在$\tilde U$上的投影的回归系数 +以及 $$ -(\tilde U^\top \tilde U)^{-1} \tilde U^\top \hat A = (\tilde U^\top \tilde U)^{-1} \tilde U^\top X' \tilde V \tilde \Sigma^{-1} \tilde U^\top = -\tilde U^\top X' \tilde V \tilde \Sigma^{-1} \tilde U^\top = \tilde A . -$$ +(I - \hat A_r)^{-1} = I + \Phi \Lambda (I - \Lambda)^{-1} \Psi +$$ (eq:longrun) + +矩阵 $\Lambda (I - \Lambda)^{-1}$ 是对角矩阵，其元素为 $\lambda_i / (1 - \lambda_i)$。 -至此，我们验证了$\tilde A$是$\hat A$在$\tilde U$上的最小二乘投影。 +因此长期响应由最持久的模态主导。 -**一个逆运算的挑战** +例如，{cite}`SSY_CEX_2026` 中的秩为3的 VAR 具有 DMD 特征值 $0.963$、$0.865$ 和 $0.547$，对应的长期权重大致为 $26$、$6.4$ 和 $1.2$。 -因为我们使用的是简化SVD，所以$\tilde U \tilde U^\top \neq I$。 +类似 {eq}`eq:longrun` 的公式出现在**加性泛函**中，这在 {doc}`advanced:additive_functionals` 中有所研究。 -因此， +假设 $X_t$ 的第一个分量是某个总量对数 $Y_t$ 的增长率，并以其均值 $\nu$ 的偏差来衡量，即 $Y_{t+1} - Y_t - \nu = e_1 X_{t+1}$，其中 $e_1 = \begin{bmatrix} 1 & 0 & \cdots & 0 \end{bmatrix}$。 + +如果 $X_t$ 遵循 VAR {eq}`eq:reducedVAR`，那么 $Y_t$ 可以分解为确定性趋势 $t \nu$、鞅 $\sum_{j=1}^t H \hat a_j$ 和平稳成分 $-g X_t$，再加上一个常数，其中 $$ -\hat A \neq \tilde U \tilde A \tilde U^\top , +g = e_1 \hat A_r (I - \hat A_r)^{-1}, \qquad H = e_1 (I - \hat A_r)^{-1} = e_1 + g $$ -所以我们不能简单地用$\tilde A$和$\tilde U$计算$\hat A$。 - -**死胡同** +{cite}`SSY_CEX_2026` 将此分解应用于总收入以及所有300个 CEX 分位数。 -我们可以先抱着最好的希望，继续构造$p \times p$矩阵$\tilde A$的特征分解： +公式 {eq}`eq:longrun` 意味着 $$ - \tilde A = \tilde W \Lambda \tilde W^{-1} -$$ (eq:tildeAeigenred) +g = e_1 \Phi \Lambda (I - \Lambda)^{-1} \Psi +$$ -其中$\Lambda$是包含$p$个特征值的对角矩阵，$\tilde W$的列是对应的特征向量。 +因此计算 $g$ 和 $H$ 不需要对 $m \times m$ 矩阵求逆。 -仿照表示法2中的步骤，我们可以轻松计算出一个$m \times p$的矩阵 +## 模拟示例 -$$ -\tilde \Phi_s = \tilde U \tilde W -$$ (eq:Phisred) +为了观察这些对象的实际作用，我们模拟了一个瘦高型数据集，该数据集与 {cite}`SSY_CEX_2026` 中的 CEX 数据具有一些共同特征。 -该矩阵对应到完整SVD中的{eq}`eq:Phisfull`。 +该数据集有 $m = 300$ 个横截面变量，我们可以将其视为由 $q \in [0, 1]$ 索引的分位数，以及 $n + 1 = 134$ 个季度观测值。 -此时，当$\hat A$由公式{eq}`eq:Ahatwithtildes`给出时，计算$\hat A \tilde \Phi_s$很有意思： +横截面由三个潜在因子 $z_t$ 驱动，它们遵循具有对角转移矩阵的一阶向量自回归（VAR）过程 $$ -\begin{aligned} -\hat A \tilde \Phi_s & = (X' \tilde V \tilde \Sigma^{-1} \tilde U^\top ) (\tilde U \tilde W) \\ - & = X' \tilde V \tilde \Sigma^{-1} \tilde W \\ - & \neq (\tilde U \tilde W) \Lambda \\ - & = \tilde \Phi_s \Lambda - \end{aligned} +z_{t+1} = D z_t + w_{t+1}, \qquad D = \textrm{diag}(0.95, 0.80, 0.50) $$ -$\hat A \tilde \Phi_s \neq \tilde \Phi_s \Lambda$意味着，与表示法2中的相应情况不同，$\tilde \Phi_s = \tilde U \tilde W$的列**不是**$\hat A$对应于矩阵$\Lambda$对角线上特征值的特征向量。 - -**一种可行的方法** - -我们继续寻找**能够**通过简化SVD计算的$\hat A$的特征向量，这里不妨定义一个$m \times p$的矩阵$\Phi$： +而观测数据为 $$ -\Phi \equiv \hat A \tilde \Phi_s = X' \tilde V \tilde \Sigma^{-1} \tilde W -$$ (eq:Phiformula) +X_t = G z_t + \sigma_v v_t +$$ -不难发现，$\Phi$的列**确实是**$\hat A$的特征向量。 +其中 $w_{t+1}$ 和 $v_t$ 为标准正态随机向量。 -这是Tu等人{cite}`tu_Rowley`证明的一个结果，我们下面来介绍。 +$G$ 的三列分别呈现出**水平**因子（使所有分位数一起变动）、**斜率**因子（使低分位数和高分位数朝相反方向变动）以及**曲率**因子（使尾部相对于中间部分变动）的形态。 +以下代码用于模拟数据。 -```{prf:proposition} -:label: prop-dmd-eigenvectors +```{code-cell} ipython3 +def simulate_data(m=300, T=134, D_diag=(0.95, 0.80, 0.50), σ_v=0.5, seed=1234): + """ + Simulate an m x T data matrix driven by three latent AR(1) factors. + """ + rng = np.random.default_rng(seed) + q = np.linspace(0, 1, m) + G = np.column_stack([np.ones(m), # level + 2 * (q - 0.5), # slope + 6 * (q - 0.5)**2 - 0.5]) # curvature + d = np.array(D_diag) + k = len(d) + z = np.zeros((k, T)) + z[:, 0] = rng.normal(size=k) / np.sqrt(1 - d**2) + for t in range(T - 1): + z[:, t+1] = d * z[:, t] + rng.normal(size=k) + X_tilde = G @ z + σ_v * rng.normal(size=(m, T)) + return X_tilde, G, z, q -$\Phi$的$p$列是$\hat A$的特征向量。 +X_tilde, G, z, q = simulate_data() +X_tilde = X_tilde - X_tilde.mean(axis=1, keepdims=True) # demean each variable +X_tilde.shape ``` -```{prf:proof} -根据公式{eq}`eq:Phiformula`，我们有 - -$$ -\begin{aligned} - \hat A \Phi & = (X' \tilde V \tilde \Sigma^{-1} \tilde U^\top ) (X' \tilde V \Sigma^{-1} \tilde W) \cr - & = X' \tilde V \tilde \Sigma^{-1} \tilde A \tilde W \cr - & = X' \tilde V \tilde \Sigma^{-1}\tilde W \Lambda \cr - & = \Phi \Lambda - \end{aligned} -$$ +下一个函数实现了 DMD 算法。 + +```{code-cell} ipython3 +def dmd(X_tilde, r): + """ + Rank-r DMD of the m x (n+1) data matrix X_tilde. + """ + X, X_prime = X_tilde[:, :-1], X_tilde[:, 1:] + U, σ, Vt = np.linalg.svd(X, full_matrices=False) + U_r, σ_r, V_r = U[:, :r], σ[:r], Vt[:r, :].T + + B = X_prime @ V_r / σ_r # X' V_r Σ_r^{-1} + A_tilde = U_r.T @ B # r x r VAR for principal components + λ, W = np.linalg.eig(A_tilde) + order = np.argsort(-np.abs(λ)) # most persistent modes first + λ, W = λ[order], W[:, order] + + Φ = B @ W # DMD modes (right eigenvectors) + Ψ = np.linalg.solve(W * λ, U_r.T) # (W Λ)^{-1} U_r^T (left eigenvectors) + A_hat = B @ U_r.T # rank-r estimator of A + a_hat = X_prime - A_hat @ X # residuals + Ω_hat = a_hat @ a_hat.T / X.shape[1] + + return dict(λ=λ, Φ=Φ, Ψ=Ψ, W=W, A_tilde=A_tilde, A_hat=A_hat, + Ω_hat=Ω_hat, σ=σ, U_r=U_r) +``` -因此 +为了选择 $r$，我们绘制 $X$ 的最大奇异值。 -$$ -\hat A \Phi = \Phi \Lambda . -$$ (eq:APhiLambda) +```{code-cell} ipython3 +σ = np.linalg.svd(X_tilde[:, :-1], compute_uv=False) -令 $\phi_i$ 为 $\Phi$ 的第 $i$ 列，$\lambda_i$ 为分解式 {eq}`eq:tildeAeigenred` 中 $\tilde A$ 对应的第 $i$ 个特征值。 +fig, ax = plt.subplots() +ax.plot(np.arange(1, 21), σ[:20], 'o') +ax.set_xlabel('index $j$') +ax.set_ylabel(r'singular value $\sigma_j$') +plt.show() +``` -将等式 {eq}`eq:APhiLambda` 两边的 $m \times 1$ 向量对应项相等得到： +奇异值在 $j = 1, 2, 3$ 时迅速下降，随后趋于平稳，因此按照拐点法则，我们将 $r$ 设为 3。 -$$ -\hat A \phi_i = \lambda_i \phi_i . -$$ +```{code-cell} ipython3 +r = 3 +res = dmd(X_tilde, r) +λ, Φ, Ψ, A_hat = res['λ'], res['Φ'], res['Ψ'], res['A_hat'] -这个等式证实了 $\phi_i$ 是 $\hat A$ 的特征向量，对应于 $\tilde A$ 和 $\hat A$ 的特征值 $\lambda_i$。 +print("DMD eigenvalues: ", np.round(λ, 3)) +print("true eigenvalues: ", np.array([0.95, 0.80, 0.50])) ``` -另见 {cite}`DDSE_book` (第238页)。 +两个持续性最强的特征值估计得很好。 +持续性最弱的那个则向零方向产生偏差。 -### $\check b$ 的解码器作为线性投影 +这是因为数据是对因子的含噪观测，所以 $\{X_t\}$ 并不严格是一阶 VAR 过程；主成分 $f_t$ 中的误差会削弱所估计的系数，这在持续性最弱的模态中表现得最为明显。 -根据特征分解 {eq}`eq:APhiLambda` ，我们可以将 $\hat A$ 表示为： +接下来我们验证 {prf:ref}`prop-dmd-eigenvectors` 和 {prf:ref}`prop-dmd-left`，并检验 $\Phi \Lambda \Phi^+$ 是否等于 $\hat A_r$。 -$$ -\hat A = \Phi \Lambda \Phi^+ . -$$ (eq:Aform12) +```{code-cell} ipython3 +Φ_pinv = np.linalg.pinv(Φ) +rel = lambda M: np.linalg.norm(M) / np.linalg.norm(A_hat) -从公式 {eq}`eq:Aform12` 我们可以推导出 $p \times 1$ 向量 $\check b_t$ 的动态： +print("max |A_hat Φ - Φ Λ| =", np.abs(A_hat @ Φ - Φ * λ).max()) +print("max |Ψ Φ - I| =", np.abs(Ψ @ Φ - np.eye(r)).max()) +print("max |Ψ A_hat - Λ Ψ| =", np.abs(Ψ @ A_hat - λ[:, None] * Ψ).max()) +print("relative error, Φ Λ Ψ =", rel(A_hat - (Φ * λ) @ Ψ)) +print("relative error, Φ Λ Φ^+ =", rel(A_hat - (Φ * λ) @ Φ_pinv)) +``` -$$ -\check b_{t+1} = \Lambda \check b_t -$$ +前四个数值在舍入误差范围内均为零，而 $\Phi \Lambda \Phi^+$ 与 $\hat A_r$ 的差距达到了百分之几。 -其中 +现在让我们来看一下各个模态。 -$$ -\check b_t = \Phi^+ X_t -$$ (eq:decoder102) +模态只能确定到某一比例因子，因此我们将每个模态归一化，使其最大绝对值为 1，并选择其符号以与 $G$ 相应列（我们以同样方式进行归一化）保持一致。 -由于 $m \times p$ 矩阵 $\Phi$ 有 $p$ 个线性独立的列，$\Phi$ 的广义逆矩阵为 +```{code-cell} ipython3 +def normalize(v): + return v / np.abs(v).max() -$$ -\Phi^{+} = (\Phi^\top \Phi)^{-1} \Phi^\top -$$ +labels = ['level', 'slope', 'curvature'] +Φ_real = np.real_if_close(Φ) -因此 +fig, axes = plt.subplots(1, 3, figsize=(12, 3.5)) +for i, ax in enumerate(axes): + g_i = normalize(G[:, i]) + φ_i = normalize(Φ_real[:, i]) + φ_i *= np.sign(φ_i @ g_i) + ax.plot(q, φ_i, lw=2, label=rf'mode $\phi_{i+1}$') + ax.plot(q, g_i, 'k--', label=f'true {labels[i]} loading') + ax.set_ylim(-1.1, 1.1) + ax.set_xlabel('quantile $q$') + ax.legend() +plt.tight_layout() +plt.show() +``` -$$ -\check b = (\Phi^\top \Phi)^{-1} \Phi^\top X -$$ (eq:checkbform) +每个模态都接近于某一真实载荷模式，但并不完全相同。 -$p \times n$ 矩阵 $\check b$ 可以被视为是 $m \times n$ 的矩阵 $X$ 在 $m \times p$ 的矩阵 $\Phi$ 上的最小二乘回归系数矩阵，因此 +由于抽样误差，所估计的特征向量会在一定程度上混合真实因子，这在第二个模态中表现得最为明显。 -$$ -\check X = \Phi \check b -$$ (eq:Xcheck_) +接下来我们将模态坐标 $b_t = \Psi X_t$ 与真实因子 $z_t$ 在标准化之后进行比较。 -是 $X$ 在 $\Phi$ 上的最小二乘投影的 $m \times n$ 矩阵。 +```{code-cell} ipython3 +b = np.real_if_close(Ψ @ X_tilde) +standardize = lambda x: (x - x.mean()) / x.std() -**$X$ 的方差分解** +fig, axes = plt.subplots(3, 1, figsize=(10, 7), sharex=True) +for i, ax in enumerate(axes): + corr = np.corrcoef(b[i], z[i])[0, 1] + ax.plot(np.sign(corr) * standardize(b[i]), lw=2, label=rf'$b_{{{i+1},t}}$') + ax.plot(standardize(z[i]), 'k--', label=rf'$z_{{{i+1},t}}$') + ax.set_ylabel(f'corr = {abs(corr):.2f}') + ax.legend(loc='upper right') +axes[-1].set_xlabel('$t$') +plt.tight_layout() +plt.show() +``` -根据这个 QuantEcon 讲座 中讨论的最小二乘的投影理论，我们可以将 $X$ 表示为 $X$ 在 $\Phi$ 上的投影 $\check X$ 和误差矩阵的和。 +DMD 能够很好地还原潜在因子。 -要验证这一点，注意到最小二乘投影 $\check X$ 与 $X$ 的关系是 +最后，我们比较以三种方式计算得到的向前 $j$ 步预测： -$$ -X = \check X + \epsilon -$$ +* 直接基于 $\hat A_r^j X_t$ 计算 +* 基于模态表示 {eq}`eq:modalforecast`，其中 $b_t = \Psi X_t$ +* 基于看似诱人的替代方案 $\Phi \Lambda^j \Phi^+ X_t$ -或 +```{code-cell} ipython3 +X = X_tilde[:, :-1] +for j in (1, 4, 12): + direct = np.linalg.matrix_power(A_hat, j) @ X + via_Ψ = np.real_if_close((Φ * λ**j) @ Ψ @ X) + via_pinv = np.real_if_close((Φ * λ**j) @ Φ_pinv @ X) + err = lambda F: np.linalg.norm(F - direct) / np.linalg.norm(direct) + print(f"j = {j:2d}: error with Ψ = {err(via_Ψ):.1e}, " + f"error with Φ^+ = {err(via_pinv):.3f}") +``` -$$ -X = \Phi \check b + \epsilon -$$ (eq:Xbcheck) +由 $\Psi$ 构建的预测精确地重现了 VAR 预测，而由 $\Phi^+$ 构建的预测则不然。 -其中 $\epsilon$ 是一个 $m \times n$ 的最小二乘误差矩阵，满足最小二乘正交条件 $\epsilon^\top \Phi =0$ 或 +在本例中，使用 $\Phi^+$ 所产生的预测差异较小，但正如 {ref}`var_dmd_ex1` 所示，这一差异会随着数据中噪声量的增加而增大。 -$$ -(X - \Phi \check b)^\top \Phi = 0_{m \times p} -$$ (eq:orthls) +## Python代码来源 -重新整理正交条件 {eq}`eq:orthls` 得到 $X^\top \Phi = \check b \Phi^\top \Phi$，这就推导出公式 {eq}`eq:checkbform`。 +你可以在 [PyDMD](https://github.com/PyDMD/PyDMD) 包中找到DMD的Python实现。 +## 练习 -### 一种近似方法 +```{exercise} +:label: var_dmd_ex1 -我们现在描述一种不使用公式 {eq}`eq:decoder102`的近似计算 $p \times 1$ 的向量 $\check b_t$ 的方法。 +在本讲座中我们看到，当且仅当 $\Phi$ 和 $U_r$ 的列空间相同时，$\Phi^+ = \Psi$。 -具体来说，以下论述改编自 {cite}`DDSE_book`（第240页）提供的一种高效计算方法，从而近似 $\check b_t$。 +1. 使用 `simulate_data(σ_v=0.0)` 模拟数据，使数据成为三个因子的精确线性组合，并计算一个秩为 3 的 DMD。通过数值方法验证 $\Phi^+ = \Psi$ 以及 $\Phi \Lambda \Phi^+ = \hat A_r$。 +1. 解释为什么在这种情况下 $\Phi$ 和 $U_r$ 的列空间重合。 +1. 对 $\sigma_v = 0.1, 0.5, 1.0$ 重复第 1 部分，并报告 $\Phi \Lambda \Phi^+$ 与 $\hat A_r$ 相差多远。 +``` -为方便起见，我们将在时间 $t=1$ 应用该方法。 +```{solution-start} var_dmd_ex1 +:class: dropdown +``` -对于 $t=1$，根据方程 {eq}`eq:Xbcheck`，我们有 +以下是一种解法。 + +```{code-cell} ipython3 +for σ_v in (0.0, 0.1, 0.5, 1.0): + X_sim = simulate_data(σ_v=σ_v)[0] + X_sim = X_sim - X_sim.mean(axis=1, keepdims=True) + out = dmd(X_sim, 3) + Φ_s, Ψ_s, λ_s, A_s = out['Φ'], out['Ψ'], out['λ'], out['A_hat'] + Φ_s_pinv = np.linalg.pinv(Φ_s) + gap_left = np.linalg.norm(Φ_s_pinv - Ψ_s) / np.linalg.norm(Ψ_s) + gap_A = np.linalg.norm(A_s - (Φ_s * λ_s) @ Φ_s_pinv) / np.linalg.norm(A_s) + print(f"σ_v = {σ_v:.1f}: |Φ^+ - Ψ|/|Ψ| = {gap_left:.1e}, " + f"|A_hat - Φ Λ Φ^+|/|A_hat| = {gap_A:.1e}") +``` -$$ - \check X_1 = \Phi \check b_1 -$$ (eq:X1proj) +当 $\sigma_v = 0$ 时，对每个 $t$ 都有 $X_t = G z_t$，因此 $X$ 和 $X'$ 的每一列都位于 $G$ 的三维列空间中。 -其中 $\check b_1$ 是一个 $p \times 1$ 的向量。 +由于 $X$ 的秩为三，$U_3$ 的列空间等于 $G$ 的列空间。 -回顾上面表示法1中的 $X_1 = U \tilde b_1$，其中 $\tilde b_1$ 是表示法1的时间1基向量，而 $U$ 来自完整SVD分解 $X = U \Sigma V^\top$。 +$\Phi = X' V_3 \Sigma_3^{-1} W$ 的列是 $X'$ 的列的线性组合，所以它们也位于 $G$ 的列空间中。 -从方程 {eq}`eq:Xbcheck` 可以得出: +由于 $\Phi$ 有三个线性无关的列，其列空间即为 $G$ 的列空间，也就是 $U_3$ 的列空间。 -$$ - U \tilde b_1 = X' \tilde V \tilde \Sigma^{-1} \tilde W \check b_1 + \epsilon_1 -$$ +于是 $\Psi$ 的行位于由 $\Phi^*$ 的行张成的空间中。 -其中 $\epsilon_1$ 是方程 {eq}`eq:Xbcheck` 中的最小二乘误差向量。 +所以 $\Psi - \Phi^+$ 的行可以写成 $C \Phi^*$，其中 $C$ 为某个矩阵。 -因此可得: +由于 $(\Psi - \Phi^+) \Phi = I - I = 0$，我们得到 $C \Phi^* \Phi = 0$，又因为 $\Phi^* \Phi$ 可逆，所以 $C = 0$。 -$$ -\tilde b_1 = U^\top X' V \tilde \Sigma^{-1} \tilde W \check b_1 + U^\top \epsilon_1 -$$ +有测量噪声时，$X'$ 中存在超出 $U_3$ 列空间的成分，两个左逆不再相同，且 $\Phi \Lambda \Phi^+$ 与 $\hat A_r$ 之间的差距会随着 $\sigma_v$ 增大而扩大。 -将误差项 $U^\top \epsilon_1$ 替换为零，并将来自 $X$ 的**完整**SVD的 $U$ 替换为来自**简化**SVD的 $\tilde U$，我们得到 $\tilde b_1$ 的近似值 $\hat b_1$: +```{solution-end} +``` -$$ - \hat b_1 = \tilde U^\top X' \tilde V \tilde \Sigma^{-1} \tilde W \check b_1 -$$ +```{exercise} +:label: var_dmd_ex2 -回顾方程 {eq}`eq:tildeAverify` 中的 $ \tilde A = \tilde U^\top X' \tilde V \tilde \Sigma^{-1}$。 +本练习将 DMD 估计量 $\hat A_r$ 与 {cite}`anderson1951` 的降秩回归估计量进行比较。 -因此可得: +降秩回归估计量是使 $\| X' - \check A X \|_F$ 最小化的秩-$r$ 矩阵。 -$$ - \hat b_1 = \tilde A \tilde W \check b_1 -$$ +它可以计算为 $\hat A^{RRR}_r = P_r \hat A$，其中 $\hat A = X' X^+$ 是最小范数最小二乘估计量，$P_r = Q_r Q_r^\top$，$Q_r$ 的列是拟合值 $\hat A X$ 的前 $r$ 个左奇异向量。 -因此，根据 $\tilde A$ 的特征分解 {eq}`eq:tildeAeigenred`，我们有 +使用 `simulate_data()` 生成的数据，在前 100 个时期上估计两种秩为 3 的估计量，然后 -$$ - \hat b_1 = \tilde W \Lambda \check b_1 -$$ +1. 计算每个估计量样本内的相对残差范数 $\| X' - \check A X \|_F / \| X' \|_F$ +1. 使用剩余时期的一步超前预测，计算样本外的同一统计量 +1. 对多个随机种子重复上述过程，并评论结果 +``` -因此， - -$$ - \hat b_1 = ( \tilde W \Lambda)^{-1} \tilde b_1 -$$ +```{solution-start} var_dmd_ex2 +:class: dropdown +``` -或者 +以下是一种解法。 + +```{code-cell} ipython3 +def rrr(X, X_prime, r): + "Reduced-rank regression estimator of rank r." + A_full = X_prime @ np.linalg.pinv(X) + Q = np.linalg.svd(A_full @ X, full_matrices=False)[0][:, :r] + return Q @ Q.T @ A_full + +def rel_resid(A, X, X_prime): + return np.linalg.norm(X_prime - A @ X) / np.linalg.norm(X_prime) + +n_est = 100 +print("seed DMD in DMD out RRR in RRR out") +for seed in range(5): + X_sim = simulate_data(seed=seed)[0] + X_sim = X_sim - X_sim.mean(axis=1, keepdims=True) + X_in, Xp_in = X_sim[:, :n_est-1], X_sim[:, 1:n_est] + X_out, Xp_out = X_sim[:, n_est-1:-1], X_sim[:, n_est:] + + A_dmd = dmd(X_sim[:, :n_est], 3)['A_hat'] + A_rrr = rrr(X_in, Xp_in, 3) + print(f"{seed:4d} {rel_resid(A_dmd, X_in, Xp_in):.3f} " + f"{rel_resid(A_dmd, X_out, Xp_out):.3f} " + f"{rel_resid(A_rrr, X_in, Xp_in):.3f} " + f"{rel_resid(A_rrr, X_out, Xp_out):.3f}") +``` -$$ - \hat b_1 = ( \tilde W \Lambda)^{-1} \tilde U^\top X_1 , -$$ (eq:beqnsmall) +从构造上看，降秩回归的样本内残差更小。 -这是对以下方程 {eq}`eq:decoder102` 中初始向量 $\check b_1$ 的计算效率较高的近似： +但在这些样本中，它的样本外预测效果比 DMD 估计量更差。 -$$ - \check b_1= \Phi^{+} X_1 -$$ (eq:bphieqn) +降秩回归从 $\hat A = X' V \Sigma^{-1} U^\top$ 出发，这里用到了 $X$ 的**所有**奇异值，包括那些会放大噪声的微小奇异值，然后选择在样本内最适合 $X'$ 的三个方向。 -（为了强调 {eq}`eq:beqnsmall` 是一个近似值，DMD的使用者有时将基向量 $\check b_t = \Phi^+ X_t $ 的分量称为**精确**DMD模态，将 $\hat b_t = ( \tilde W \Lambda)^{-1} \tilde U^\top X_t$ 的分量称为**近似**模态。） +而 DMD 估计量在回归之前就舍弃了小的奇异值，这起到了一种正则化的作用。 -在给定 $X_t$ 的条件下，我们可以通过精确模态计算解码后的 $\check X_{t+j}, j = 1, 2, \ldots $ ： +```{solution-end} +``` -$$ -\check X_{t+j} = \Phi \Lambda^j \Phi^{+} X_t -$$ (eq:checkXevoln) +```{exercise} +:label: var_dmd_ex3 +使用讲座中模拟数据的秩为 3 的 DMD， -或者通过近似模态计算解码的 $\hat X_{t+j}$: +1. 证明公式 {eq}`eq:Ahatpower`，并用它推导公式 {eq}`eq:longrun` +1. 通过与 $(I - \hat A_r)^{-1}$ 的直接计算进行比较，用数值方法验证公式 {eq}`eq:longrun` +1. 将 $X_t$ 的第一个分量视为增长率变量，分别用直接方法和模态公式 $g = e_1 \Phi \Lambda (I - \Lambda)^{-1} \Psi$ 计算 $g = e_1 \hat A_r (I - \hat A_r)^{-1}$，并报告每个模态对 $g$ 的贡献 +1. 用 $\Phi^+$ 代替 $\Psi$ 重新计算模态公式，并报告答案变化了多少 +``` -$$ - \hat X_{t+j} = \Phi \Lambda^j (\tilde W \Lambda)^{-1} \tilde U^\top X_t . -$$ (eq:checkXevoln2) +```{solution-start} var_dmd_ex3 +:class: dropdown +``` -然后我们可以使用解码后的 $\check X_{t+j}$ 或 $\hat X_{t+j}$ 来预测 $X_{t+j}$。 +对于第 1 部分，{prf:ref}`prop-dmd-left` 给出 $\hat A_r = \Phi \Lambda \Psi$ 以及 $\Psi \Phi = I$。 +因此 $\hat A_r^2 = \Phi \Lambda (\Psi \Phi) \Lambda \Psi = \Phi \Lambda^2 \Psi$，由归纳法可得对所有 $j \geq 1$，$\hat A_r^j = \Phi \Lambda^j \Psi$。 -### 使用更少的模态 +当对所有 $i$ 都有 $|\lambda_i| < 1$ 时， -在实际应用中，我们通常只使用少数几个模态，通常不多于三个。 +$$ +\sum_{j=1}^\infty \hat A_r^j = \Phi \left( \sum_{j=1}^\infty \Lambda^j \right) \Psi = \Phi \Lambda (I - \Lambda)^{-1} \Psi +$$ -前面的一些公式假设中，我们保留了与 $X$ 的奇异值相关的所有 $p$ 个模态。 +而且由于 $\hat A_r$ 的所有特征值都在单位圆内，有 $(I - \hat A_r)^{-1} = I + \sum_{j=1}^\infty \hat A_r^j$。 -我们可以调整公式，描述只保留 $r < p$ 个最大奇异值的情况。 +以下是第 2--4 部分的代码。 -在这种情况下，我们只需将 $\tilde \Sigma$ 替换为相应的 $r\times r$ 奇异值矩阵，将 $\tilde U$ 替换为对应于 $r$ 个最大奇异值的 $m \times r$ 的矩阵，将 $\tilde V$ 替换为对应于 $r$ 个最大奇异值的 $n \times r$ 的矩阵。 +```{code-cell} ipython3 +m = A_hat.shape[0] +I = np.eye(m) +weights = λ / (1 - λ) -上述所有重要公式都有其对应的形式。 +LR_direct = np.linalg.inv(I - A_hat) +LR_modal = np.real_if_close(I + (Φ * weights) @ Ψ) +print("max |direct - modal| for (I - A_hat)^{-1}:", + np.abs(LR_direct - LR_modal).max()) +e1 = np.zeros(m) +e1[0] = 1.0 +g_direct = e1 @ A_hat @ LR_direct +g_modal = np.real_if_close(e1 @ (Φ * weights) @ Ψ) +print("max |g_direct - g_modal|:", np.abs(g_direct - g_modal).max()) -## Python代码来源 +print("\nlong-run weights λ/(1-λ):", np.round(weights.real, 2)) +for i in range(r): + g_i = np.real_if_close(Φ[0, i] * weights[i] * Ψ[i, :]) + share = np.linalg.norm(g_i) / np.linalg.norm(g_direct) + print(f"mode {i+1}: |contribution to g| / |g| = {share:.2f}") + +g_pinv = np.real_if_close(e1 @ (Φ * weights) @ Φ_pinv) +print("\nrelative error in g using Φ^+:", + np.linalg.norm(g_pinv - g_direct) / np.linalg.norm(g_direct)) +``` + +模态公式与直接计算的结果相符，误差仅来自舍入，并且除了构造 $\Phi$ 和 $\Psi$ 之外，只需要进行 $r \times r$ 规模的运算。 + +最持久的模态，即具有最大权重 $\lambda_i/(1-\lambda_i)$ 的模态，对 $g$ 的贡献最大。 -你可以在这里找到DMD的Python实现： +这些比例之和不必为一，也可能有某一项超过一，因为不同模态的贡献可能部分相互抵消。 -https://mathlab.sissa.it/pydmd +用 $\Phi^+$ 代替 $\Psi$ 会得到一个不同的、错误的答案。 + +```{solution-end} +``` From be6fc89e12bc94026fb271f5efb97f3471e0fb76 Mon Sep 17 00:00:00 2001 From: Matt McKay Date: Thu, 22 Oct 2026 14:44:55 +1100 Subject: [PATCH 6/9] Update translation: .translate/state/var_dmd.md.yml --- .translate/state/var_dmd.md.yml | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/.translate/state/var_dmd.md.yml b/.translate/state/var_dmd.md.yml index 5ddfada8..c393fb5a 100644 --- a/.translate/state/var_dmd.md.yml +++ b/.translate/state/var_dmd.md.yml @@ -1,6 +1,6 @@ -source-sha: 78030a3a27f6527046675bcd8a8d27995ca25af6 -synced-at: "2026-07-18" +source-sha: 63880c0644afdcad5599295767c212cd34ac0280 +synced-at: "2026-10-22" model: claude-sonnet-5 -mode: RESYNC -section-count: 6 -tool-version: 0.17.0 +mode: UPDATE +section-count: 10 +tool-version: 0.30.0 From 9f32b6368c49bed06f0c891bc4ce273b149649b2 Mon Sep 17 00:00:00 2001 From: Matt McKay Date: Thu, 22 Oct 2026 14:44:56 +1100 Subject: [PATCH 7/9] Update translation: lectures/_static/quant-econ.bib --- lectures/_static/quant-econ.bib | 49 +++++++++++++++++++++++++++++++++ 1 file changed, 49 insertions(+) diff --git a/lectures/_static/quant-econ.bib b/lectures/_static/quant-econ.bib index 5d2b1fb8..a77683af 100644 --- a/lectures/_static/quant-econ.bib +++ b/lectures/_static/quant-econ.bib @@ -5873,3 +5873,52 @@ @article{Leonard1983 pages = {461--479}, year = {1983} } + +@article{Geweke1982, + author = {Geweke, John}, + title = {Measurement of Linear Dependence and Feedback Between Multiple + Time Series}, + journal = {Journal of the American Statistical Association}, + volume = {77}, + number = {378}, + pages = {304--313}, + year = {1982} +} + +@article{sargent2026dynamic, + title = {Dynamic Mode Decompositions and Vector Autoregressions}, + author = {Sargent, Thomas J. and Selvakumar, Yatheesan J. and Yang, Ziyue}, + journal = {International Economic Review}, + year = {2026}, + publisher = {Wiley Online Library} +} + +@techreport{SSY_CEX_2026, + author = {Sargent, Thomas J. and Selvakumar, Yatheesan J. and Yang, Ziyue}, + title = {Aggregate Shocks and Cross-Section Dynamics: Quantifying + Redistribution and Insurance in {US} Household Data}, + institution = {New York University}, + year = {2026} +} + +@article{stock_watson2002, + author = {Stock, James H. and Watson, Mark W.}, + title = {Forecasting Using Principal Components from a Large Number of + Predictors}, + journal = {Journal of the American Statistical Association}, + volume = {97}, + number = {460}, + pages = {1167--1179}, + year = {2002} +} + +@article{anderson1951, + author = {Anderson, T. W.}, + title = {Estimating Linear Restrictions on Regression Coefficients for + Multivariate Normal Distributions}, + journal = {Annals of Mathematical Statistics}, + volume = {22}, + number = {3}, + pages = {327--351}, + year = {1951} +} From bb0a705e525a25904b6edcf55bdf03c150b47085 Mon Sep 17 00:00:00 2001 From: Matt McKay Date: Thu, 22 Oct 2026 15:04:25 +1100 Subject: [PATCH 8/9] Drop svd_intro from this resync: #303 already delivered it The resync of QuantEcon/lecture-python.myst#1055 regenerates all three of its lectures, but svd_intro.md was delivered by #303 from the same English. Applying #1055's changes again altered 9 lines of a translation that was already current, so svd_intro.md and its state file go back to main's versions. Co-Authored-By: Claude Opus 5.5 --- .translate/state/svd_intro.md.yml | 2 +- lectures/svd_intro.md | 18 +++++++++--------- 2 files changed, 10 insertions(+), 10 deletions(-) diff --git a/.translate/state/svd_intro.md.yml b/.translate/state/svd_intro.md.yml index 8b42d8e1..936cdd37 100644 --- a/.translate/state/svd_intro.md.yml +++ b/.translate/state/svd_intro.md.yml @@ -3,4 +3,4 @@ synced-at: "2026-10-22" model: claude-sonnet-5 mode: UPDATE section-count: 12 -tool-version: 0.30.0 +tool-version: 0.29.6 diff --git a/lectures/svd_intro.md b/lectures/svd_intro.md index 5fb62ae0..b447e5f1 100644 --- a/lectures/svd_intro.md +++ b/lectures/svd_intro.md @@ -41,7 +41,7 @@ translation: 在这一讲的后续内容 {doc}`VAR 与 DMD ` 中，我们将介绍如何利用 SVD，借助**动态模态分解**（DMD）来快速计算一阶向量自回归（VAR）模型的低秩近似。 -与主成分分析类似，DMD 也可以被看作是一种数据降维方法，它通过将数据投影到一组有限的因子上来表示数据中的显著模式。 +与主成分分析类似，动态模态分解（DMD）也可以被看作是一种数据降维方法，它通过将数据投影到一组有限的因子上来表示数据中的显著模式。 ## 基本设定 @@ -69,7 +69,7 @@ translation: 在 $m \gg n$ 的情况下，即属性数量 $m$ 远大于个体数量 $n$，且在时间序列环境中 $n$ 等于数据集 $X$ 中所覆盖的时间段数量时，我们将采用不同的方法。 -在后续的 {doc}`VARs 和 DMDs ` 中，我们仍然使用奇异值分解，但目的是构建**动态模态分解**（DMD）。 +在后续的 {doc}`VARs 和 DMDs ` 中，我们仍然使用奇异值分解，但目的是构建动态模态分解（DMD）。 ## 奇异值分解 @@ -133,7 +133,7 @@ $m \times n$ 矩阵 $X$ 的这种结构为构建系统开启了大门 稍后我们还将描述**经济型**或**简化** SVD。 -在研究**简化** SVD 之前，我们将进一步讨论**完全** SVD 的性质。 +在研究**简化** SVD之前，我们将进一步讨论**完全** SVD的性质。 ## 四个基本子空间 @@ -688,7 +688,7 @@ $$ V^\top = \begin{bmatrix}V_1^\top \\V_2^\top \\\vdots\\V_n^\top \end{bmatrix} $$ -在方程 {eq}`eq:PCA1` 中，每个$m \times n$矩阵$U_{j}V_{j}^\top $显然是秩1的。 +在方程{eq}`eq:PCA1`中，每个$m \times n$矩阵$U_{j}V_{j}^\top $显然是秩1的。 因此，我们有 @@ -757,7 +757,7 @@ $$ XX^\top =P\Lambda P^\top . $$ (eq:XXo) -因此数据矩阵$X$可以写成： +因此数据矩阵$X$可以写成: $$ X =\begin{bmatrix}P_1|P_2|\ldots|P_m\end{bmatrix} @@ -783,17 +783,17 @@ X&=\sqrt{\lambda_1}P_1\tilde{\epsilon}_1+\sqrt{\lambda_2}P_2\tilde{\epsilon}_2+\ \end{aligned} $$ -这与 +这与SVD分解 $$ X=\sigma_1U_1{V_1}^{\top}+\sigma_2 U_2{V_2}^{\top}+\ldots+\sigma_{p} U_{p}{V_{p}}^{\top} $$ -是等价的，只要我们令： +是等价的，只要我们令: -* $U_j=P_j$（变量在第$j$个主成分上的载荷向量） +* $U_j=P_j$ (变量在第j个主成分上的载荷向量) -* ${V_k}^{\top}=\tilde{\epsilon}_k$（第$k$个主成分） +* ${V_k}^{\top}=\tilde{\epsilon}_k$ (第k个主成分) 由于对于给定的数据矩阵$X$，计算$P$和$U$有不同的算法，根据所使用的算法，我们可能会得到符号差异或特征向量顺序的不同。 From 91f3106c2088b1dff0503f273e5edeed7d5aa592 Mon Sep 17 00:00:00 2001 From: Matt McKay Date: Thu, 22 Oct 2026 15:05:31 +1100 Subject: [PATCH 9/9] var_dmd: restore the 2025 reviewers' conventions lost in the resync MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The resync translated the rewritten sections of QuantEcon/lecture-python.myst#1055 afresh, so the conventions Songyi Wang and Humphrey Yang set in their 2025 review of this lecture did not survive. Restored, text only (structure unchanged: 992 lines, 16 headings, 28 directives, 88 math lines, 12 code cells, 2 links): - pseudo-inverse is 广义逆矩阵 and Moore-Penrose is 穆尔-彭罗斯 (7 lines), with the reviewers' Baidu Baike link and their note on the choice; - the case labels are 短胖型情况 and 高瘦型情况, and tall-skinny is 高瘦型; - reduced SVD is 简化 (the resync wrote 约化); - a dimension before 矩阵 takes 的: `$m \times n$ 的矩阵` (18 places). Co-Authored-By: Claude Opus 5.5 --- lectures/var_dmd.md | 50 ++++++++++++++++++++++----------------------- 1 file changed, 25 insertions(+), 25 deletions(-) diff --git a/lectures/var_dmd.md b/lectures/var_dmd.md index 4248f6e0..246d2a93 100644 --- a/lectures/var_dmd.md +++ b/lectures/var_dmd.md @@ -87,7 +87,7 @@ $$ (eq:Xvector) 例如，{cite}`SSY_CEX_2026` 从每个对数分位数中减去对数私人收入的横截面均值，使用该均值的增长率作为一个附加变量，然后从所有变量中减去时间序列均值。 -我们的数据被组织成一个 $m \times (n+1)$ 矩阵 $\tilde X$ +我们的数据被组织成一个 $m \times (n+1)$ 的矩阵 $\tilde X$ $$ \tilde X = \begin{bmatrix} X_1 \mid X_2 \mid \cdots \mid X_n \mid X_{n+1} \end{bmatrix} @@ -95,7 +95,7 @@ $$ 其中对于 $t = 1, \ldots, n+1$，$m \times 1$ 向量 $X_t$ 由 {eq}`eq:Xvector` 给出。 -由 $\tilde X$，我们构造两个 $m \times n$ 矩阵 +由 $\tilde X$，我们构造两个 $m \times n$ 的矩阵 $$ X = \begin{bmatrix} X_1 \mid X_2 \mid \cdots \mid X_{n}\end{bmatrix} @@ -137,16 +137,16 @@ $$ \hat A = X' X^+ $$ (eq:commonA) -其中 $X^+$ 是 $X$ 的 [摩尔-彭若斯伪逆](https://en.wikipedia.org/wiki/Moore%E2%80%93Penrose_inverse)。 +其中 $X^+$ 是 $X$ 的 [穆尔-彭罗斯广义逆矩阵](https://baike.baidu.com/item/穆尔-彭罗斯广义逆矩阵/22770999?fr=aladdin)。 我们关注两种情形 * $n \gg m$，即时间序列观测值 $n$ 远多于变量数 $m$ * $m \gg n$，即变量数 $m$ 远多于时间序列观测值 $n$ -**矮胖情形：** +**短胖型情况：** -当 $n \gg m$ 且 $X$ 的**行**线性无关时，$X X^\top$ 具有逆矩阵，伪逆 $X^+$ 为 +当 $n \gg m$ 且 $X$ 的**行**线性无关时，$X X^\top$ 具有逆矩阵，广义逆矩阵 $X^+$ 为 $$ X^+ = X^\top (X X^\top )^{-1} @@ -164,9 +164,9 @@ $$ (eq:Ahatform101) 公式 {eq}`eq:Ahatform101` 的右边正比于 $X_{t+1}$ 和 $X_t$ 的经验交叉二阶矩矩阵，乘以 $X_t$ 的二阶矩矩阵的逆。 -**瘦高情形：** +**高瘦型情况：** -当 $m \gg n$ 且 $X$ 的**列**线性无关时，$X^\top X$ 具有逆矩阵，伪逆 $X^+$ 为 +当 $m \gg n$ 且 $X$ 的**列**线性无关时，$X^\top X$ 具有逆矩阵，广义逆矩阵 $X^+$ 为 $$ X^+ = (X^\top X)^{-1} X^\top @@ -200,20 +200,20 @@ $$ ### 用 SVD 计算 $X^+$ -计算伪逆 $X^+$ 的一种高效方法是从**约化**奇异值分解开始 +计算广义逆矩阵 $X^+$ 的一种高效方法是从**简化**奇异值分解开始 $$ X = U \Sigma V^\top $$ (eq:SVDDMD) -其中 $U$ 是 $m \times p$ 矩阵，$\Sigma$ 是由正奇异值 $\sigma_1 \geq \sigma_2 \geq \cdots \geq \sigma_p > 0$ 构成的 $p \times p$ 对角矩阵，$V$ 是 $n \times p$ 矩阵。 +其中 $U$ 是 $m \times p$ 的矩阵，$\Sigma$ 是由正奇异值 $\sigma_1 \geq \sigma_2 \geq \cdots \geq \sigma_p > 0$ 构成的 $p \times p$ 对角矩阵，$V$ 是 $n \times p$ 的矩阵。 -正如我们在 {doc}`奇异值分解 ` 中所看到的，对于约化 SVD +正如我们在 {doc}`奇异值分解 ` 中所看到的，对于简化 SVD * $U^\top U = I_{p \times p}$ 且 $V^\top V = I_{p \times p}$，但是 * 除非 $p = m$，否则 $U U^\top \neq I_{m \times m}$；除非 $p = n$，否则 $V V^\top \neq I_{n \times n}$。 -对于任意秩 $p$，$X$ 的摩尔-彭若斯伪逆为 +对于任意秩 $p$，$X$ 的穆尔-彭罗斯广义逆矩阵为 $$ X^{+} = V \Sigma^{-1} U^\top @@ -249,7 +249,7 @@ $$ (eq:AhatSVDformula) 一种补救方法是只保留 $r < p$ 个最大的奇异值。 -设 $U_r$ 为由 $U$ 的前 $r$ 列组成的 $m \times r$ 矩阵，设 $\Sigma_r$ 为 $\Sigma$ 左上角的 $r \times r$ 分块，设 $V_r$ 为由 $V$ 的前 $r$ 列组成的 $n \times r$ 矩阵。 +设 $U_r$ 为由 $U$ 的前 $r$ 列组成的 $m \times r$ 的矩阵，设 $\Sigma_r$ 为 $\Sigma$ 左上角的 $r \times r$ 分块，设 $V_r$ 为由 $V$ 的前 $r$ 列组成的 $n \times r$ 的矩阵。 我们对 $A$ 的 **秩-$r$ 估计量** 为 @@ -257,9 +257,9 @@ $$ \hat A_r = X' V_r \Sigma_r^{-1} U_r^\top $$ (eq:Ahat_r) -$m \times m$ 矩阵 $\hat A_r$ 的秩至多为 $r$。 +$m \times m$ 的矩阵 $\hat A_r$ 的秩至多为 $r$。 -注意，$V_r \Sigma_r^{-1} U_r^\top$ 是 $\hat X_r = U_r \Sigma_r V_r^\top$ 的伪逆，而 $\hat X_r$ 正是 {doc}`奇异值分解 ` 中埃卡特-杨定理所描述的 $X$ 的最佳秩-$r$ 近似。 +注意，$V_r \Sigma_r^{-1} U_r^\top$ 是 $\hat X_r = U_r \Sigma_r V_r^\top$ 的广义逆矩阵，而 $\hat X_r$ 正是 {doc}`奇异值分解 ` 中埃卡特-杨定理所描述的 $X$ 的最佳秩-$r$ 近似。 于是 $\hat A_r = X' \hat X_r^+$：我们通过在最小二乘公式 $\hat A = X' X^+$ 中用 $X$ 的埃卡特-杨近似代替 $X$，得到了秩-$r$ 估计量。 @@ -303,7 +303,7 @@ $$ 其各分量是截面数据 $X_t$ 的前 $r$ 个 **主成分**（参见 {doc}`奇异值分解 ` 中关于主成分分析的讨论）。 -将它们堆叠起来可得到 $r \times n$ 矩阵 $F = U_r^\top X = \Sigma_r V_r^\top$。 +将它们堆叠起来可得到 $r \times n$ 的矩阵 $F = U_r^\top X = \Sigma_r V_r^\top$。 $X_{t+1}$ 对 $f_t$ 的最小二乘回归具有系数矩阵 @@ -352,14 +352,14 @@ $$ (eq:Atilde) 1. 计算约简奇异值分解 $X = U \Sigma V^\top$。 1. 保留 $r$ 个最大的奇异值，从而构成 $U_r$、$\Sigma_r$ 和 $V_r$。 -1. 构成 $r \times r$ 矩阵 $\tilde A = U_r^\top X' V_r \Sigma_r^{-1}$ 并计算其特征分解 +1. 构成 $r \times r$ 的矩阵 $\tilde A = U_r^\top X' V_r \Sigma_r^{-1}$ 并计算其特征分解 $$ \tilde A W = W \Lambda $$ (eq:tildeAeigenred) 其中 $\Lambda$ 是由特征值 $\lambda_1, \ldots, \lambda_r$ 组成的对角矩阵，$W$ 的列是对应的特征向量。 -1. 构成 $m \times r$ 矩阵 $\Phi$ 和 $r \times m$ 矩阵 $\Psi$ +1. 构成 $m \times r$ 的矩阵 $\Phi$ 和 $r \times m$ 的矩阵 $\Psi$ $$ \Phi = X' V_r \Sigma_r^{-1} W, \qquad \Psi = (W \Lambda)^{-1} U_r^\top @@ -369,9 +369,9 @@ $\Lambda$ 对角线上的特征值 $\lambda_i$ 称为 **DMD 特征值**，$\Phi$ 我们将在下文看到，$\Psi$ 的行 $\psi_i$ 是与之匹配的**左**特征向量。 -请注意，其中的主要计算工作在于对 $m \times n$ 矩阵 $X$ 进行奇异值分解，以及对较小的 $r \times r$ 矩阵 $\tilde A$ 进行特征分解。 +请注意，其中的主要计算工作在于对 $m \times n$ 的矩阵 $X$ 进行奇异值分解，以及对较小的 $r \times r$ 的矩阵 $\tilde A$ 进行特征分解。 -我们完全不需要构造或分解 $m \times m$ 矩阵。 +我们完全不需要构造或分解 $m \times m$ 的矩阵。 在整个过程中，我们假设 $\tilde A$ 具有 $r$ 个线性无关的特征向量和非零特征值，这是一般情形。 @@ -420,11 +420,11 @@ $$ 我们还可以进一步说明。 -将 $\hat A_r$ 写为 $B U_r^\top$，其中 $B = X' V_r \Sigma_r^{-1}$ 是 $m \times r$ 矩阵。 +将 $\hat A_r$ 写为 $B U_r^\top$，其中 $B = X' V_r \Sigma_r^{-1}$ 是 $m \times r$ 的矩阵。 因为 $B U_r^\top$ 的非零特征值与 $U_r^\top B = \tilde A$ 的特征值一致，所以 $\hat A_r$ 的特征值为 $\lambda_1, \ldots, \lambda_r$，再加上 $m - r$ 个零特征值。 -因此，估计的 VAR {eq}`eq:reducedVAR` 稳定当且仅当 $\max_i |\lambda_i| < 1$，该条件可以通过一个 $r \times r$ 矩阵来验证。 +因此，估计的 VAR {eq}`eq:reducedVAR` 稳定当且仅当 $\max_i |\lambda_i| < 1$，该条件可以通过一个 $r \times r$ 的矩阵来验证。 （另见 {cite}`DDSE_book`，第 238 页。） @@ -489,7 +489,7 @@ $$ b_t = \Psi X_t = (W \Lambda)^{-1} U_r^\top X_t $$ (eq:modalcoords) -计算 $b_t$ 只需要主成分 $f_t = U_r^\top X_t$ 以及一个 $r \times r$ 矩阵。 +计算 $b_t$ 只需要主成分 $f_t = U_r^\top X_t$ 以及一个 $r \times r$ 的矩阵。 由于 $\Psi \Phi = I$，对 $\hat A_r = \Phi \Lambda \Psi$ 进行迭代可得 @@ -517,7 +517,7 @@ $$ (eq:modalforecast) ### 一个看似可行但实际不可行的替代方案 -摩尔-彭若斯伪逆 $\Phi^+ = (\Phi^* \Phi)^{-1} \Phi^*$ 也满足 $\Phi^+ \Phi = I$。 +穆尔-彭罗斯广义逆矩阵 $\Phi^+ = (\Phi^* \Phi)^{-1} \Phi^*$ 也满足 $\Phi^+ \Phi = I$。 向量 $\check b_t = \Phi^+ X_t$ 是 $X_t$ 对 $\Phi$ 各列进行最小二乘回归所得的系数向量，而 $\Phi \check b_t$ 则是 $X_t$ 在 $\Phi$ 的列空间上的正交投影。 @@ -620,11 +620,11 @@ $$ g = e_1 \Phi \Lambda (I - \Lambda)^{-1} \Psi $$ -因此计算 $g$ 和 $H$ 不需要对 $m \times m$ 矩阵求逆。 +因此计算 $g$ 和 $H$ 不需要对 $m \times m$ 的矩阵求逆。 ## 模拟示例 -为了观察这些对象的实际作用，我们模拟了一个瘦高型数据集，该数据集与 {cite}`SSY_CEX_2026` 中的 CEX 数据具有一些共同特征。 +为了观察这些对象的实际作用，我们模拟了一个高瘦型数据集，该数据集与 {cite}`SSY_CEX_2026` 中的 CEX 数据具有一些共同特征。 该数据集有 $m = 300$ 个横截面变量，我们可以将其视为由 $q \in [0, 1]$ 索引的分位数，以及 $n + 1 = 134$ 个季度观测值。