第 11b 章 线性状态空间模型与时变 β

本章对应 Tsay 原书第 11 章 11.2–11.7 节。第 11a 章用局部水平模型讲清了 Kalman 滤波的全部思想；本章把它推广到向量状态、时变系统矩阵的一般线性高斯状态空间模型。先学会把常见模型（时变系数 CAPM、回归、ARMA、带 ARMA 误差的回归、趋势–季节–周期分解）写成状态空间形式，再给出一般的滤波、平滑、扰动平滑、缺失值与预测公式，最后用原书的两个实证例子和三个量化例子收尾：时变 β 的动态对冲、配对交易的动态对冲比率。

学习目标

1. 写出一般线性高斯状态空间模型的状态方程和观测方程，认识各系统矩阵的含义。
2. 能把时变系数 CAPM、线性回归、ARMA（Akaike、Harvey、Aoki 三种表示）、带 ARMA 误差的回归、结构时间序列模型写成状态空间形式。
3. 掌握一般 Kalman 滤波 (11.64) 与固定区间平滑器 (11.74) 的推导要点，了解扰动平滑的用途与稳态 Riccati 方程。
4. 会用"把未来当作缺失值"的方法做多步预测，并处理部分分量缺失。
5. 能用 Kalman 滤波在线估计时变 β 和配对交易的动态对冲比率，明白模型设定（价差是白噪声还是 AR(1)）怎样决定交易信号的质量，避免平滑值带来的前视偏差。

读前导读

这一章在解决什么问题。 你在 CFA 里用 OLS 估过 β：拿过去 60 个月的数据回归一次，得到一个数，然后假设它不变。现实中 β 会漂移，业界的应对通常是"滚动窗口"，但窗口取 20 天还是 120 天全凭经验。状态空间模型换了个思路：把 β 当作一个看不见、会慢慢游走的"状态"，每天观察到的收益只是它带噪声的影子；Kalman 滤波则是一套记账规则，每来一个新数据，就按"这次的意外有多大、我原来有多不确定"来修正对 β 的估计。第 11a 章用一维的局部水平模型讲清了这个思想，本章把它推广到向量和矩阵：状态可以同时包含 α、β、价差、趋势、季节，系统矩阵可以随时间变化。

本章第二个要点是"写成状态空间形式"这门手艺。回归、ARMA、趋势–季节分解看起来各不相同，但都能装进同一对方程 (11.26)–(11.27)，然后用同一个滤波程序处理。这像会计里把各种业务都记成借贷分录：格式统一了，同一套总账就能处理所有业务。

需要先想起来的数学。

• 矩阵乘法、转置与逆。\(T_ts_t\) 是矩阵乘以列向量；\(A'\) 是转置；\(V_t^{-1}\) 是逆矩阵，标量情形就是倒数。一个规则要记牢：\((AB)'=B'A'\)。若 \(s=(\alpha,\beta)'\)，\(Z=(1,r_M)\)，则 \(Zs=\alpha+\beta r_M\)，这正是 CAPM 的右边。见 第 00 册第 06 章 线性代数速成。
• 协方差矩阵与线性变换。若 \(x\) 的协方差矩阵是 \(\Sigma\)，则 \(Ax\) 的协方差矩阵是 \(A\Sigma A'\)。这就是组合方差 \(w'\Sigma w\) 的矩阵版本（\(w'\) 是 \(1\times n\) 的 \(A\)）。本章 (11.56) \(\Sigma_{t+1|t}=T_t\Sigma_{t|t}T_t'+R_tQ_tR_t'\) 就是这条规则加上"独立变量方差相加"。见 第 00 册第 06 章。
• 多元正态的条件分布。若 \((x,y)\) 联合正态，则 \(E(x|y)=E(x)+\mathrm{Cov}(x,y)\mathrm{Var}(y)^{-1}(y-E(y))\)，\(\mathrm{Var}(x|y)=\mathrm{Var}(x)-\mathrm{Cov}(x,y)\mathrm{Var}(y)^{-1}\mathrm{Cov}(y,x)\)。标量时第一式就是"回归预测 = 均值 + β×偏离"，β 恰为 \(\mathrm{Cov}/\mathrm{Var}\)。这就是第 11a 章定理 11.1，本章所有更新公式都由它推出。见 第 00 册第 07 章 概率中的分析工具。
• 行列式。似然里的 \(|V_t|\) 是 \(V_t\) 的行列式，可理解为多维"方差的大小"；\(k=1\) 时就是 \(V_t\) 本身。例如 \(\mathrm{diag}(2,3)\) 的行列式是 6。见 第 00 册第 06 章。
• 递推与求和交换。平滑公式把一个很长的求和写成"从后往前一步步累加"的递推，和用逆向递推算年金现值是同一个思路。见 第 00 册第 04 章 级数与收敛。

记号提示：\(s_{t|t-1}\) 读作"用 \(t-1\) 及以前的信息对 \(s_t\) 的估计"，竖线右边是"信息截止到哪天"；\(F_t\) 表示截至 \(t\) 的全部观测（信息集）；\(\mathrm{diag}(a,b)\) 是对角线为 \(a,b\)、其余为 0 的矩阵；\(I_p\) 是 \(p\) 阶单位阵。

怎么读这一章。 必读：11.10.1（模型）、11.11.1–11.11.2（时变 CAPM 与回归）、11.12.2 的算法汇总、11.15（缺失值与预测）以及量化实战全部，特别是例三。11.11.3 的三种 ARMA 表示第一次只需看懂 Harvey 形式，Akaike 与 Aoki 可以跳过。11.12.1 和 11.13–11.14 的推导记号很重，原书也建议偏应用的读者只记结论 (11.64)、(11.74)、(11.82)；建议先读讲解框弄懂"每个式子在干什么"，以后需要自己写代码或改模型时再逐行推。SsfPack 的 mJPhi 等细节只在用 R/S-PLUS 复现原书时才需要。

11.10 一般线性高斯状态空间模型

11.10.1 模型

许多经济金融模型都能写成状态空间形式：ARIMA、含不可观测成分的动态线性模型、时变系数回归、随机波动率模型（第 12b 章）等。一般的高斯线性状态空间模型为

• \(s_t\)：\(m\) 维状态向量（state vector），不可观测；
• \(y_t\)：\(k\) 维观测向量；
• \(d_t\)、\(c_t\)：已知的确定性向量（截距）；
• \(T_t\)（\(m\times m\)）：转移矩阵；\(Z_t\)（\(k\times m\)）：观测矩阵；
• \(R_t\)（\(m\times n\)）：通常由单位阵的部分列组成，用来指定哪些状态分量有随机冲击；
• \(\eta_t\sim N(0,Q_t)\)，\(e_t\sim N(0,H_t)\)，二者相互独立（这一条可以放宽）；
• 初始状态 \(s_1\sim N(\mu_{1|0},\Sigma_{1|0})\)，与所有噪声独立。

(11.27) 是观测方程（measurement equation），(11.26) 是状态方程或转移方程，它是一阶马尔可夫的。\(T_t,R_t,Q_t,Z_t,H_t\) 统称系统矩阵（system matrices），通常很稀疏，可以是未知参数 \(\theta\) 的函数，用极大似然估计。注意记号冲突：\(T_t\) 是转移矩阵，而样本量也记为 \(T\)，靠下标区分。

局部水平模型就是 \(m=k=1\)，\(d_t=c_t=0\)，\(T_t=Z_t=R_t=1\)，\(Q_t=\sigma_\eta^2\)，\(H_t=\sigma_e^2\) 的特例。

白话解释：两个方程分工明确。状态方程说"看不见的东西怎么演化"：明天的状态 = 今天的状态经过 \(T_t\) 变换 + 一点随机冲击。观测方程说"看得见的东西怎么由看不见的东西生成"：今天的观测 = 状态经过 \(Z_t\) 映射 + 测量噪声。\(Q_t\) 大表示状态变得快，\(H_t\) 大表示观测噪声多；两者之比（信噪比）决定了滤波对新数据有多"敏感"。 "一阶马尔可夫"的意思是：只要知道 \(s_t\)，再多知道 \(s_{t-1},s_{t-2},\dots\) 对预测 \(s_{t+1}\) 没有帮助。这不是限制，因为可以把需要的滞后项都塞进状态向量（11.11.3 节的 ARMA 就是这么做的）。 \(R_t\) 的作用举例：状态有 3 个分量但只有前 2 个受冲击，就取 \(R_t=\begin{bmatrix}1&0\\0&1\\0&0\end{bmatrix}\)，\(\eta_t\) 是 2 维的；\(R_tQ_tR_t'\) 是 \(3\times3\) 的状态冲击协方差，第 3 行第 3 列为 0。

11.10.2 紧凑形式与 SsfPack 记号

把两个方程叠起来：

\(\Omega_t=\mathrm{Cov}(u_t)=\mathrm{diag}(R_tQ_tR_t',H_t)\)。原书使用的 SsfPack 正是按这种紧凑形式输入：mDelta、mPhi、mOmega 分别对应 \(\delta,\Phi,\Omega\)，初值放在 \((m+1)\times m\) 的矩阵 \(\Sigma=\begin{bmatrix}\Sigma_{1|0}\\\mu_{1|0}'\end{bmatrix}\)（mSigma）中。

扩散初始化：\(\Sigma_{1|0}=\Sigma_*+\lambda\Sigma_\infty\)，\(\lambda\) 很大（可趋于无穷），SsfPack 中把对应对角元写成 \(-1\) 表示扩散。系统矩阵可以时不变，也可以时变；时变时 SsfPack 用数据矩阵 mX 存放时变数值，再用与 \(\Phi_t\) 同维的索引矩阵 mJPhi（以及 mJDelta、mJOmega）指出哪个元素取 mX 的哪一列，其余元素填 \(-1\)。

白话解释：扩散初始化就是"开头我对状态一无所知"。给初始方差一个极大的数（代码里用 \(10^6\) 或 \(10^8\)），等于说先验几乎没有分量，第一批观测一进来就几乎完全由数据决定估计值。这和贝叶斯里的"无信息先验"是同一个意思。代价是前一两期的新息方差 \(V_t\) 极大，算似然时通常要丢掉这几期（本章代码 loglik 的 burn=2 就是这个用途）。

11.11 把常见模型写成状态空间形式

11.11.1 时变系数 CAPM

\(r_t\)、\(r_{M,t}\) 是资产和市场的超额收益，\(e_t\sim N(0,\sigma_e^2)\)，\(\eta_t\sim N(0,\sigma_\eta^2)\)，\(\epsilon_t\sim N(0,\sigma_\epsilon^2)\) 互相独立。α 与 β 都是随机游走。对应的状态空间形式：

这里时变的是观测矩阵 \(Z_t\)：它包含了外生变量 \(r_{M,t}\)。紧凑形式为 \(\Phi_t=\begin{bmatrix}1&0\\0&1\\1&r_{M,t}\end{bmatrix}\)，\(\Omega=\mathrm{diag}(\sigma_\eta^2,\sigma_\epsilon^2,\sigma_e^2)\)，扩散初始化 \(\Sigma=\begin{bmatrix}-1&0\\0&-1\\0&0\end{bmatrix}\)。原书示例：GM 1990-01 至 2003-12 月度超额收益，S&P 500 为市场，取 \((\sigma_\eta,\sigma_\epsilon,\sigma_e)=(0.02,0.04,0.1)\)，mX 为 \(168\times2\) 的矩阵 cbind(1, sp)，mJPhi 的第 3 行设为 \((1,2)\)，表示 \(\Phi_t\) 第 3 行两个元素分别取 mX 的第 1、2 列。

这个模型是量化交易中最有用的状态空间模型之一：Kalman 滤波给出的 \(\beta_{t|t-1}\) 就是"今天开盘前对今天 β 的最优估计"，可以直接作为对冲比率。

11.11.2 线性回归与递推最小二乘

常系数回归 \(y_t=x_t'\beta+e_t\) 也能写成状态空间形式：令状态 \(s_t=\beta\) 不随时间变化，

即 \(T_t=I_p\)，\(Z_t=x_t'\)，\(Q_t=0\)，\(H_t=\sigma_e^2\)，用扩散初始化。此时 Kalman 滤波就是递推最小二乘（recursive least squares, RLS）：\(s_{t|t}\) 是用前 \(t\) 个观测做 OLS 的系数，\(s_{T|T}\) 就是全样本 OLS 估计，\(\Sigma_{T|T}=\sigma_e^2(X'X)^{-1}\)。本章实战例一会数值验证这一点。

把状态方程放松为 \(\beta_{t+1}=\beta_t+R_t\eta_t\)，\(\eta_t\sim N(0,1)\)，\(R_t=(\sigma_1,\dots,\sigma_p)'\)，就得到随机系数回归；某个 \(\sigma_i=0\) 表示该系数不变。时变 CAPM 就是它的特例。SsfPack 中 GetSsfReg(cbind(1,sp)) 直接生成市场模型的状态空间形式：mPhi\(=[I_2;\,0\ 0]\)，mOmega\(=\mathrm{diag}(0,0,1)\)，mJPhi 第 3 行为 \((1,2)\)。

11.11.3 ARMA 模型的三种表示

零均值 ARMA(\(p,q\)) 模型 \(\phi(B)y_t=\theta(B)a_t\)（11.30）。令 \(m=\max(p,q+1)\)，补零后写成

其中 \(\phi_i=0\)（\(i>p\)），\(\theta_j=0\)（\(j>q\)）。

(1) Akaike 方法（1975）。状态取"产生所有未来预测所需的最小变量集合"：

观测方程 \(y_t=Zs_t\)，\(Z=(1,0,\dots,0)\)（11.32）。推导靠两件事：

• 第一个分量：\(s_{1,t+1}=y_{t+1}=y_{t+1|t}+a_{t+1}=s_{2t}+a_{t+1}\)（11.33）。
• 预测更新公式：用 MA 表示 \(y_t=\sum_i\psi_ia_{t-i}\)（\(\psi_0=1\)，\(\psi_1=\phi_1-\theta_1\)，\(\psi_2=\phi_1\psi_1+\phi_2-\theta_2\)，…，一般地 \(\psi_j=\sum_{i=1}^{j}\phi_i\psi_{j-i}-\theta_j\)）（11.34），可得

即新信息全部体现在新息 \(a_{t+1}\) 上，以权重 \(\psi_{j-1}\) 修正旧预测。

推导拆解：为什么是 \(\psi_{j-1}\)？

1. 用 MA(\(\infty\)) 表示写出 \(y_{t+j}=a_{t+j}+\psi_1a_{t+j-1}+\cdots+\psi_{j-1}a_{t+1}+\psi_ja_t+\cdots\)。
2. 站在 \(t\) 时刻，未来冲击 \(a_{t+1},\dots,a_{t+j}\) 的条件期望都是 0，所以 \(y_{t+j|t}=\psi_ja_t+\psi_{j+1}a_{t-1}+\cdots\)。
3. 站在 \(t+1\) 时刻，\(a_{t+1}\) 已经知道了，于是 \(y_{t+j|t+1}=\psi_{j-1}a_{t+1}+\psi_ja_t+\cdots\)。
4. 两式相减：\(y_{t+j|t+1}-y_{t+j|t}=\psi_{j-1}a_{t+1}\)。 金融上看，这就像分析师修正盈利预测：新一季的意外 \(a_{t+1}\) 出来后，对未来第 \(j\) 期的预测调整 \(\psi_{j-1}\) 倍，\(\psi\) 衰减得越慢，意外的影响越持久。最后一个分量还要用 AR 递推：\(y_{t+m|t+1}=\sum_{i=1}^m\phi_iy_{t+m-i|t}+\psi_{m-1}a_{t+1}\)（11.36）。合起来：

即 \(s_{t+1}=Ts_t+R\eta_t\)，\(\eta_t=a_{t+1}\sim N(0,\sigma_a^2)\)（11.38）。

(2) Harvey 方法（1993, §4.4）。取 \(s_{1t}=y_t\)，其余分量递推定义：

\(s_{2,t+1}=\phi_2s_{1t}+s_{3t}-\theta_1\eta_t\)，…，\(s_{m,t+1}=\phi_ms_{1t}-\theta_{m-1}\eta_t\)，其中 \(\eta_t=a_{t+1}\)。于是

\(Z=(1,0,\dots,0)\)，没有测量误差。它的好处是 AR、MA 系数直接出现在系统矩阵里，估计时最方便。

(3) Aoki 方法（1987, 第 4 章）。先看 MA(\(q\))：\(y_t=\theta(B)a_t\)，取 \(s_t=(a_{t-q},\dots,a_{t-1})'\)，转移矩阵是上移位矩阵：

此时 \(a_t\) 同时出现在状态方程和观测方程中（两种噪声相关）。AR(\(p\)) 模型 \(\phi(B)z_t=a_t\) 有两种 Aoki 形式：(i) \(s_t=(z_{t-p+1},\dots,z_t)'\)，用伴随矩阵（末行 \((\phi_p,\dots,\phi_1)\)）转移，冲击向量 \((0,\dots,0,1)'a_{t+1}\)，\(z_t=(0,\dots,0,1)s_t\)（11.42）；(ii) 把最后一个分量换成 \(z_t-a_t\)，冲击向量 \((0,\dots,1,\phi_1)'a_t\)，\(z_t=(0,\dots,0,1)s_t+a_t\)（11.43）。对 ARMA(\(p,q\))（设 \(q<p\)），引入辅助变量 \(z_t=a_t/\phi(B)\)，则 \(\phi(B)z_t=a_t\)，\(y_t=\theta(B)z_t\)；用 (11.42) 作转移方程时观测方程为 \(y_t=(-\theta_{p-1},\dots,-\theta_1,1)s_t\)（11.44），用 (11.43) 时再加上 \(a_t\)（11.45）。

小结：ARMA 有许多状态空间表示，估计和预测时任选其一，结果相同。反过来，对时不变状态空间模型，由 Cayley–Hamilton 定理（第 01 册）可证观测 \(y_t\) 服从 ARMA(\(m,m\))，\(m\) 为状态维数。

SsfPack 示例与符号陷阱。GetSsfArma 采用 Harvey 方法。

• AR(1) \(y_t=0.6y_{t-1}+a_t\)，\(a_t\sim N(0,0.4^2)\)：mPhi\(=(0.6,1)'\)，mOmega\(=\mathrm{diag}(0.16,0)\)，mSigma\(=(0.25,0)'\)。平稳 AR(1) 的初值取无条件分布：\(\Sigma_{1|0}=0.16/(1-0.36)=0.25\)，\(\mu_{1|0}=0\)。
• ARMA(2,1) \(y_t=1.2y_{t-1}-0.35y_{t-2}+a_t-0.25a_{t-1}\)，\(a_t\sim N(0,1.1^2)\)：\(T=\begin{bmatrix}1.2&1\\-0.35&0\end{bmatrix}\)，\(Z=(1,0)\)，\(R=(1,-0.25)'\)，mOmega 左上块 \(\sigma_a^2RR'=\begin{bmatrix}1.21&-0.3025\\-0.3025&0.075625\end{bmatrix}\)。状态 \(s_{1t}=y_t\)，\(s_{2t}=-0.35y_{t-1}-0.25a_t\)，mSigma 是其无条件协方差 \(\begin{bmatrix}4.0607&-1.4874\\-1.4874&0.5731\end{bmatrix}\)（解离散 Lyapunov 方程 \(\Sigma=T\Sigma T'+\sigma_a^2RR'\) 可验证）。

推导拆解：Lyapunov 方程从哪来？平稳时状态的协方差不随时间变，即 \(\mathrm{Var}(s_{t+1})=\mathrm{Var}(s_t)=\Sigma\)。对 \(s_{t+1}=Ts_t+R\eta_t\) 两边取协方差：\(Ts_t\) 的协方差是 \(T\Sigma T'\)（线性变换规则），\(R\eta_t\) 的协方差是 \(\sigma_a^2RR'\)，两者独立所以相加，得到 \(\Sigma=T\Sigma T'+\sigma_a^2RR'\)。 标量情形就是 AR(1) 的无条件方差：\(\sigma^2=\phi^2\sigma^2+\sigma_a^2\Rightarrow\sigma^2=\sigma_a^2/(1-\phi^2)\)，上一条的 \(0.16/(1-0.36)=0.25\) 正是这样算的。矩阵情形是 3 个未知数（\(\Sigma\) 对称，只有 \(\Sigma_{11},\Sigma_{12},\Sigma_{22}\)）的线性方程组，可以手解，也可以用 scipy.linalg.solve_discrete_lyapunov。

注意 SsfPack 的 MA 多项式写成 \(\theta(B)=1+\theta_1B+\cdots\)，与本书 \(1-\theta_1B-\cdots\) 符号相反，所以上例要输入 ma=-0.25。statsmodels 的 ARIMA 也采用 \(1+\theta_1B\) 的约定（第 11a 章代码里取了负号）。

11.11.4 带 ARMA 误差的线性回归

\(x_t=1\) 时就是非零均值的 ARMA。设 \(s_t\) 是 \(z_t\) 的（如 Harvey 形式）状态向量，把回归系数作为常数状态 \(\beta_t=\beta\)（11.48）拼进去：\(s_t^*=(s_t',\beta_t')'\)，

SsfPack GetSsfRegArma(X, ar=c(1.2,-0.35), ma=c(-0.25)) 对 \(y_t=\beta_0+\beta_1x_t+z_t\)、\(z_t=1.2z_{t-1}-0.35z_{t-2}+a_t-0.25a_{t-1}\) 生成 \(5\times4\) 的 mPhi（前两列 ARMA 部分，后两列单位阵，末行 \((1,0,0,0)\) 并通过 mJPhi 第 5 行 \((-1,-1,1,2)\) 引用 mX 的两列）、mOmega 左上块 \(\begin{bmatrix}1&-0.25\\-0.25&0.0625\end{bmatrix}\)、mSigma 中 ARMA 部分为无条件协方差 \(\begin{bmatrix}3.356&-1.229\\-1.229&0.4736\end{bmatrix}\)，回归系数部分为扩散初始化。

这个构造的价值在于：状态向量可以随意拼接。实战例三正是这样把"缓慢漂移的对冲比率"和"均值回复的 AR(1) 价差"拼成一个三维状态。

11.11.5 结构时间序列模型

结构时间序列模型（structural time series model, STSM），也叫不可观测成分模型（unobserved component model）：

分别是趋势、季节、周期和不规则成分。

• 趋势（可能含两个单位根）：
  \[\mu_{t+1}=\mu_t+\beta_t+\eta_t,\quad \eta_t\sim N(0,\sigma_\eta^2);\qquad \beta_t=\beta_{t-1}+\varsigma_t,\quad \varsigma_t\sim N(0,\sigma_\varsigma^2),\tag{11.52}\]
  \(\mu_1,\beta_1\sim N(0,\xi)\)，\(\xi\) 很大（如 \(10^8\)）。\(\sigma_\varsigma=0\) 时为带漂移的随机游走；\(\sigma_\varsigma=\sigma_\eta=0\) 时为确定性线性趋势。这叫局部线性趋势模型。
• 季节：\((1+B+\cdots+B^{s-1})\gamma_t=\omega_t\)，\(\omega_t\sim N(0,\sigma_\omega^2)\)（11.53），\(s\) 为季节周期，即一个周期内季节效应之和在噪声意义下为零；\(\sigma_\omega=0\) 时季节确定。
• 周期：
  \[\begin{bmatrix}\omega_{t+1}\\\omega_{t+1}^*\end{bmatrix}=\delta\begin{bmatrix}\cos\lambda_c&\sin\lambda_c\\-\sin\lambda_c&\cos\lambda_c\end{bmatrix}\begin{bmatrix}\omega_t\\\omega_t^*\end{bmatrix}+\begin{bmatrix}\varepsilon_t\\\varepsilon_t^*\end{bmatrix},\tag{11.54}\]
  \((\varepsilon_t,\varepsilon_t^*)'\sim N(0,\sigma_\varepsilon^2(1-\delta^2)I_2)\)，\(\omega_0,\omega_0^*\sim N(0,\sigma_\varepsilon^2)\)，\(\delta\in(0,1]\) 是阻尼因子，\(\lambda_c=2\pi/q\)（\(q\) 为周期长度）；\(\delta=1\) 时为确定性正余弦波。

SsfPack 的 GetSsfStsm 最多支持 10 个周期成分，参数对应关系为：irregular→\(\sigma_e\)，level→\(\sigma_\eta\)，slope→\(\sigma_\varsigma\)，seasonalDummy/seasonalTrig/seasonalHS→\((\sigma_\omega,s)\)，Cycle0…Cycle9→\((\sigma_\varepsilon,\lambda_c,\delta)\)。例如局部水平模型 \(\sigma_e=0.4,\sigma_\eta=0.2\)：GetSsfStsm(irregular=0.4, level=0.2) 得 mPhi\(=(1,1)'\)，mOmega\(=\mathrm{diag}(0.04,0.16)\)，mSigma\(=(-1,0)'\)。Python 中对应的是 statsmodels.tsa.UnobservedComponents（参数 level、trend、seasonal、cycle 等）。

11.12 一般 Kalman 滤波

本节和下两节沿用第 11a 章的思路推导，参考 Durbin & Koopman (2001, 第 4 章)。原书说明这部分记号繁重，偏应用的读者可以只记结论 (11.64)、(11.74)、(11.82)。

11.12.1 推导

记 \(s_j|F_i\sim N(s_{j|i},\Sigma_{j|i})\)。

预测步。由 (11.26)，\(\eta_t\) 与 \(F_t\) 独立：

新息。\(y_{t|t-1}=c_t+Z_ts_{t|t-1}\)，

与标量情形一样，\(E(v_t|F_{t-1})=0\)，\(v_t\) 与 \(F_{t-1}\) 独立，\(\{v_t\}\) 是独立的正态向量序列。

更新步。\(C_t=\mathrm{Cov}(s_t,v_t|F_{t-1})=\Sigma_{t|t-1}Z_t'\)，由定理 11.1（第 11a 章）：

（假设 \(H_t\) 可逆，从而 \(V_t\) 可逆。）

推导拆解：更新步怎么套定理 11.1。

1. 在信息 \(F_{t-1}\) 下，\(s_t\) 与 \(v_t\) 联合正态：\(s_t\) 均值 \(s_{t|t-1}\)、方差 \(\Sigma_{t|t-1}\)；\(v_t\) 均值 0、方差 \(V_t\)。
2. 协方差：\(\mathrm{Cov}(s_t,v_t)=\mathrm{Cov}(s_t,\ Z_t(s_t-s_{t|t-1})+e_t)=\Sigma_{t|t-1}Z_t'\)。这里 \(s_{t|t-1}\) 在 \(F_{t-1}\) 下是已知常数，\(e_t\) 与 \(s_t\) 独立，所以只剩 \(\mathrm{Var}(s_t)Z_t'\)。
3. 知道 \(F_t\) 等于知道 \(F_{t-1}\) 加上 \(v_t\)（\(y_t\) 和 \(v_t\) 可以互推），于是套条件正态公式：均值 \(=s_{t|t-1}+C_tV_t^{-1}v_t\)，方差 \(=\Sigma_{t|t-1}-C_tV_t^{-1}C_t'\)。 直观上，\(C_tV_t^{-1}\) 就是"状态对新息的回归系数"：新息每意外 1 个单位，状态估计调整多少。

金融直觉：以时变 β 为例，\(v_t=r_t-\alpha_{t|t-1}-\beta_{t|t-1}r_{M,t}\) 是"按昨天估的 α、β 预测今天个股收益，实际差了多少"。β 的调整量约为 \(\frac{\mathrm{Var}(\beta)\,r_{M,t}}{V_t}v_t\)：市场今天涨跌幅越大（\(r_{M,t}\) 大），这次意外越能归因到 β 上；市场几乎不动时，意外多半是特质噪声，β 基本不调。这正是人工判断 β 时会用的逻辑。

合并。把更新步代入预测步：

\(K_t\) 是一般情形的 Kalman 增益。方差递推：

验证：\(T_t\Sigma_{t|t-1}L_t'=T_t\Sigma_{t|t-1}T_t'-T_t\Sigma_{t|t-1}Z_t'K_t'=T_t\Sigma_{t|t-1}T_t'-K_tV_tK_t'\)，这正是把 (11.62) 代入 (11.56) 的结果。

11.12.2 算法汇总

给定 \(s_{1|0},\Sigma_{1|0}\)，对 \(t=1,\dots,T\)：

需要同期滤波值 \(s_{t|t}\) 时（比如想知道"收盘后对今天 β 的估计"），用等价的写法：

本章实战代码就是按这个版本写的。

11.12.3 稳态与计算量

对时不变模型，\(\Sigma_{t|t-1}\) 收敛到 \(\Sigma^*\)，它满足离散代数 Riccati 方程（discrete algebraic Riccati equation）

达到稳态后 \(V_t\)、\(K_t\)、\(\Sigma_{t+1|t}\) 都是常数，不必再更新，可以大大节省计算（第 11a 章的局部水平模型就是它的标量版本）。一般情形每步主要是 \(m\times m\) 矩阵乘法和 \(k\times k\) 求逆，总计算量 \(O(T(m^3+k^3))\)，关于样本长度是线性的。

白话解释：方差递推 (11.63) 不依赖观测值 \(y_t\)，只依赖系统矩阵。所以对时不变模型，\(\Sigma_{t|t-1}\) 会走向一个"不确定性的均衡点"：每期状态冲击 \(RQR'\) 增加的不确定性，恰好被一次观测更新减少的不确定性抵消。Riccati 方程就是"增加 = 减少"这个平衡条件。到达稳态后 Kalman 滤波退化为固定权重的指数平滑，和第 11a 章局部水平模型的结论一致。 \(O(\cdot)\) 记号（"大 O"）表示"量级不超过"，\(O(T)\) 就是"随样本长度线性增长"，样本翻倍、耗时约翻倍。

11.12.4 状态估计误差

记 \(x_t=s_t-s_{t|t-1}\)，\(\mathrm{Var}(x_t|F_{t-1})=\Sigma_{t|t-1}\)，则

\(x_1=s_1-s_{1|0}\)。推导：\(x_{t+1}=T_tx_t+R_t\eta_t-K_tv_t\)，再代入 \(v_t\)。与标量情形一样可证 \(\{v_t\}\) 相互独立，且 \(\{v_t,\dots,v_T\}\) 与 \(F_{t-1}\) 独立。

11.13 状态平滑

由定理 11.1 第 3 条，

由 (11.65) 递推得

所以 \(s_{t|T}=s_{t|t-1}+\Sigma_{t|t-1}q_{t-1}\)（11.69），

满足 \(q_{t-1}=Z_t'V_t^{-1}v_t+L_t'q_t\)，\(q_T=0\)（11.70）。同理由定理 11.1 第 4 条得平滑方差，\(M_t=\mathrm{Var}(q_t)\)（11.72–11.73）。合起来就是 de Jong (1989) 的固定区间平滑器（fixed interval smoother）：

实施两步：先前向跑 (11.64)，存下 \(v_t,V_t,K_t,s_{t|t-1},\Sigma_{t|t-1}\)；再从 \(q_T=0,M_T=0\) 往回跑 (11.74)。

白话解释：平滑 = 滤波值 + "未来所有意外对今天状态的修正"。(11.66) 说，第 \(t\) 天之后每一天的新息 \(v_j\) 都含有一点关于 \(s_t\) 的信息，按协方差加权后加回来。\(q_{t-1}\) 把这些加权的未来新息打包成一个数（向量），\(L_t'\) 是"信息每往回传一天要打的折扣"。所以它天然要从 \(T\) 往回算：\(q_T=0\) 是因为最后一天之后没有未来；往前每一步，把当天的新息加进去，再把已有的累计值打折。这和从到期日往回折现债券现金流的逆向递推完全同构。 \(M_t\) 是 \(q_t\) 的方差，它衡量"未来信息让我们对 \(s_t\) 多确定了多少"，所以平滑方差 \(\Sigma_{t|T}\) 总是不大于滤波方差 \(\Sigma_{t|t-1}\)。

金融直觉：平滑值是"事后诸葛亮"。例如复盘 2020 年 3 月某股的 β，平滑值会参考 4–6 月的数据，所以比当时实时能得到的估计准得多。做研究、归因分析时可以用；回测交易信号时绝不能用，例二的输出会量化这个"作弊"带来的虚假改进。

11.14 扰动平滑

有时我们关心的不是状态，而是噪声本身的平滑估计 \(e_{t|T}=E(e_t|F_T)\) 和 \(\eta_{t|T}=E(\eta_t|F_T)\)。它们是模型诊断的工具：观测噪声的平滑值异常大提示异常值（第 12a 章会从贝叶斯角度再讨论），状态噪声的平滑值异常大提示结构突变（例如某天公司并购导致 β 跳变）。

推导与状态平滑相同：\(e_{t|T}=\sum_{j\ge t}E(e_tv_j')V_j^{-1}v_j\)（11.75），其中 \(E(e_tv_t')=H_t\)，\(E(e_tv_j')=E(e_tx_j')Z_j'\)（\(j>t\)）（11.76），而由 (11.65)，\(E(e_tx_{t+1}')=-H_tK_t'\)，\(E(e_tx_{t+2}')=-H_tK_t'L_{t+1}'\)，…（11.77）。整理得

\(o_t=V_t^{-1}v_t-K_t'q_t\) 叫平滑测量误差（smoothing measurement error）。状态噪声只影响 \(t+1\) 以后的观测，\(E(\eta_tv_{t+1}')=Q_tR_t'Z_{t+1}'\)，\(E(\eta_tx_{t+2}')=Q_tR_t'L_{t+1}'\)，…（11.79），得

Koopman (1993) 由此给出平滑状态的另一种前向递推：

只需存 \(q_t\)，省内存。

扰动平滑算法汇总（\(t=T,\dots,1\)，\(q_T=0\)，\(M_T=0\)）：

用平滑扰动除以其标准差，就得到"辅助残差"（auxiliary residuals），超过 \(\pm2\) 或 \(\pm3\) 的点值得逐一检查。

11.15 缺失值与预测

11.15.1 缺失值

情形一：\(y_t\) 整体缺失（\(t=\ell+1,\dots,\ell+h\)）。令 \(v_t=0\)、\(K_t=0\)，滤波照常：

平滑时 \(q_{t-1}=T_t'q_t\)，\(M_{t-1}=T_t'M_tT_t\)。

情形二：\(y_t\) 部分分量缺失。设能观测到的部分为 \(y_t^*=Jy_t\)，\(J\) 由 \(I_k\) 的部分行组成（选择矩阵）。在该时点改用

滤波与平滑公式不变。比如同时跟踪 A 股和港股的一对双重上市股票，某天港股休市，就只用 A 股的观测更新。易于处理缺失值是状态空间模型的一大优点。

11.15.2 预测

站在预测原点 \(t\)，最小均方误差预测 \(y_t(j)=E(y_{t+j}|F_t)\)。最简单的做法是把 \(y_{t+1},\dots,y_{t+h}\) 当作缺失值，继续跑 Kalman 滤波。1 步预测为 \(y_t(1)=c_{t+1}+Z_{t+1}s_{t+1|t}\)，误差 \(e_t(1)=Z_{t+1}(s_{t+1}-s_{t+1|t})+e_{t+1}\)，方差就是 \(V_{t+1}\)。\(j\) 步：

这正是 \(v_{t+j}=0\)、\(K_{t+j}=0\) 时的 (11.64)。注意时变 CAPM 这类模型的 \(Z_{t+j}\) 含未来的 \(r_{M,t+j}\)，预测时要么给定情景值，要么另建模型预测它。

预测误差 \(\{v_t\}\) 用于构造似然 \(\ln L=-\frac12\sum_t(k\ln2\pi+\ln|V_t|+v_t'V_t^{-1}v_t)\)；标准化预测误差 \(D_t^{-1/2}v_t\)（\(D_t=\mathrm{diag}\{V_t(1,1),\dots,V_t(k,k)\}\)）用于模型检验。

推导拆解：这个似然叫"预测误差分解"。

1. 联合密度可以拆成条件密度的连乘：\(p(y_1,\dots,y_T)=\prod_tp(y_t|F_{t-1})\)（概率的乘法规则反复使用）。
2. 在 \(F_{t-1}\) 下，\(y_t\) 服从正态，均值 \(c_t+Z_ts_{t|t-1}\)、方差 \(V_t\)；减去均值就是 \(v_t\)，所以 \(p(y_t|F_{t-1})\) 就是 \(N(0,V_t)\) 在 \(v_t\) 处的密度。
3. \(k\) 维正态密度取对数为 \(-\frac12(k\ln2\pi+\ln|V_t|+v_t'V_t^{-1}v_t)\)，对 \(t\) 求和即得上式。\(k=1\) 时就是 \(-\frac12(\ln2\pi+\ln V_t+v_t^2/V_t)\)，代码 loglik 正是这个。 好处是：滤波一遍就得到似然，不必处理 \(T\times T\) 的大协方差矩阵；再用数值优化器对 \(Q,H\) 等参数求最大化即可。

11.16 原书应用

例 11.2：GM 的市场模型与时变 CAPM

数据为 GM 月度简单超额收益（%），1990-01 至 2003-12，市场为 S&P 500 超额收益。

固定系数。\(r_t=\alpha+\beta r_{M,t}+e_t\)（11.86）的 OLS 结果：截距 0.1982（标准误 0.6302，t=0.31），斜率 1.0457（标准误 0.1453，t=7.20）；\(R^2=0.2378\)，调整 \(R^2=0.2332\)，DW=2.029，Jarque–Bera 2.53（p=0.28），Ljung–Box 24.21（p=0.34），残差标准误 8.13（166 自由度）。模型充分。

用 SsfPack 估计同一模型：GetSsfReg(mX) 建立状态空间形式，把 mOmega[3,3] 设为 \(\exp(\text{parm})\)（对数参数化保证为正），SsfFit(c.start=10, ...) 得 \(\sqrt{\exp(\hat\theta)}=8.130114\)；SsfMomentEst(task="STSMO") 给出平滑状态 \((0.1982025,1.045702)\)，标准差 \((0.6302091,0.1453139)\)——与 OLS 完全一致，印证了 11.11.2 节"常系数回归的 Kalman 滤波就是 RLS"。

时变 CAPM。对 (11.29) 用对数参数化估计三个方差，得

α 与 β 的新息标准差都接近 0。原书图 11.5 显示：平滑的 \(\alpha_t\) 约为 0.206，几乎不变；\(\beta_t\) 在 1.02–1.06 之间小幅变动。纵轴刻度很窄，说明对 GM 来说固定系数模型已经足够。

这个"否定"的结论同样有价值：时变 β 模型应当让数据来决定 β 是否真的在变。如果 \(\hat\sigma_\epsilon\approx0\)，时变模型自动退化为常系数模型，不会凭空制造出变化。

例 11.3：强生季度 EPS 的趋势–季节分解

数据为 Johnson & Johnson 1960–1980 季度每股收益的对数（第 02b 章 2.10 节用过的数据）。模型

即 \(\gamma_t=-\sum_{j=1}^3\gamma_{t-j}+\omega_t\)，有三个参数 \(\sigma_e,\sigma_\eta,\sigma_\omega\)。状态空间形式：

\(\mathrm{Cov}(\eta_t,\omega_t)=\mathrm{diag}(\sigma_\eta^2,\sigma_\omega^2)\)。用 GetSsfStsm(irregular, level, seasonalDummy=c(σω,4)) 并对数参数化估计，得

不规则成分几乎为零，EPS 的变化几乎完全由趋势与季节解释。平滑后的趋势 \(\mu_{t|T}\) 与季节 \(\gamma_{t|T}\)（\(T=84\)）及 95% 区间见原书图 11.6：季节模式随时间演变，这是确定性季节虚拟变量做不到的。图 11.7 给出 1 步预测误差和平滑响应残差。

注意：成分分解不唯一，依赖模型设定和约束。改用 seasonalTrig（三角函数季节）会得到另一套分解，所以解释"趋势"和"季节"时要谨慎；但只要是有效的分解，对预测没有影响。

量化实战

应用场景

1. 时变 β 与动态对冲：用 (11.29) 在线估计个股、基金或策略对市场（或期货）的暴露，\(\beta_{t|t-1}\) 作为当天的对冲比率。与滚动窗口 OLS 相比，它没有"窗口长度"这个硬参数，用似然从数据中学出 β 该变多快；也不会在大事件滑出窗口时产生跳变。
2. 配对交易的动态对冲比率：把第 08b 章的静态协整回归 \(y_t=\alpha+\gamma x_t+\epsilon_t\) 改为随机游走系数，即得动态对冲比率。本节例三会说明：价差本身该怎样进入模型，直接决定了信号能不能交易。
3. 因子暴露与风格漂移：把 \(Z_t\) 换成多个因子收益（第 09 章），就是时变因子暴露，用于基金风格漂移监控与收益归因。
4. 结构时间序列：盈利、销量、宏观数据的趋势–季节分解；动态 Nelson–Siegel 利率曲线模型也是状态空间模型（三个因子作为状态，各期限收益率作为观测）。
5. 回测纪律：信号只能来自 \(s_{t|t-1}\) 或 \(s_{t|t}\)；超参数（方差）也应只用训练期估计。平滑值只用于事后研究。

Python 示例：时变 β、动态对冲与配对交易

代码先实现通用的 Kalman 滤波（同期滤波版本）和固定区间平滑器，支持时变 \(Z_t\) 和缺失值；然后做三个例子：

• 例一：验证常系数回归的 Kalman 滤波 = OLS（对应例 11.2 的第一部分）。
• 例二：日度数据，个股 β 是随机游走。只用前 500 天估计 \(\sigma_\beta,\sigma_e\)，在后 1000 天比较 Kalman 预测 β、滚动 OLS 和静态 OLS 的跟踪误差和对冲效果。
• 例三：配对交易。\(Y\) 与 \(X\) 的对数价格满足 \(y_t=0.3+\gamma_tx_t+s_t\)，对冲比率 \(\gamma_t\) 缓慢漂移，价差 \(s_t\) 是 \(\phi=0.9\) 的 AR(1)。比较三种 Kalman 设定得到的交易信号。

import numpy as np
from scipy.optimize import minimize
import statsmodels.api as sm

# ================= 通用 Kalman 滤波与平滑 (11.64)(11.74)，观测为标量 =================
def kalman_filter(y, Z, T, RQR, H, s10, P10):
    """y:(n,) 可含 nan；Z:(n,m) 时变观测向量；T:(m,m)；RQR:(m,m)；H:标量"""
    n, m = Z.shape
    sp, Pp = np.zeros((n, m)), np.zeros((n, m, m))      # s_{t|t-1}, Sigma_{t|t-1}
    sf, Pf = np.zeros((n, m)), np.zeros((n, m, m))      # s_{t|t},   Sigma_{t|t}
    v, V, Kg = np.zeros(n), np.full(n, np.inf), np.zeros((n, m))
    s, P = s10.astype(float), P10.astype(float)
    for t in range(n):
        sp[t], Pp[t] = s, P
        if np.isnan(y[t]):                               # 缺失：v_t=0, K_t=0
            sf[t], Pf[t] = s, P
        else:
            z = Z[t]
            v[t] = y[t] - z @ s
            C = P @ z                                    # C_t = Sigma_{t|t-1} Z_t'
            V[t] = z @ C + H
            sf[t] = s + C * v[t] / V[t]
            Pf[t] = P - np.outer(C, C) / V[t]
            Kg[t] = T @ C / V[t]                         # K_t = T Sigma Z' V^{-1}
        s = T @ sf[t]
        P = T @ Pf[t] @ T.T + RQR
    return dict(sp=sp, Pp=Pp, sf=sf, Pf=Pf, v=v, V=V, K=Kg)

def kalman_smoother(kf, Z, T):
    n, m = Z.shape
    q, M = np.zeros(m), np.zeros((m, m))
    ss, Ps = np.zeros((n, m)), np.zeros((n, m, m))
    for t in range(n - 1, -1, -1):
        if np.isfinite(kf["V"][t]):
            z = Z[t]; L = T - np.outer(kf["K"][t], z)
            q = z * kf["v"][t] / kf["V"][t] + L.T @ q
            M = np.outer(z, z) / kf["V"][t] + L.T @ M @ L
        else:
            q = T.T @ q; M = T.T @ M @ T
        ss[t] = kf["sp"][t] + kf["Pp"][t] @ q
        Ps[t] = kf["Pp"][t] - kf["Pp"][t] @ M @ kf["Pp"][t]
    return ss, Ps

def loglik(kf, burn=2):
    v, V = kf["v"][burn:], kf["V"][burn:]
    ok = np.isfinite(V)
    return -0.5 * np.sum(np.log(2 * np.pi * V[ok]) + v[ok] ** 2 / V[ok])

# ================= 例一：常系数回归的 Kalman 滤波 = 递推最小二乘 (11.46) =================
rng = np.random.default_rng(42)
n = 168                                               # 与例 11.2 相同：14 年月度数据
rm = rng.normal(0.5, 4.5, n)                          # 市场超额收益(%)
r = 0.2 + 1.05 * rm + rng.normal(0, 8.13, n)          # 个股超额收益(%)
X = np.column_stack([np.ones(n), rm])
ols = sm.OLS(r, X).fit()
kf = kalman_filter(r, X, np.eye(2), np.zeros((2, 2)), ols.scale, np.zeros(2), 1e8 * np.eye(2))
print("[例一] OLS 系数:", np.round(ols.params, 4), "标准误:", np.round(ols.bse, 4))
print("[例一] KF s_T|T:", np.round(kf["sf"][-1], 4), "标准差:", np.round(np.sqrt(np.diag(kf["Pf"][-1])), 4))

# ================= 例二：日度时变 beta 与动态对冲 =================
rng = np.random.default_rng(2)
n, ntr = 1500, 500                                     # 前 500 天估计超参数，后 1000 天样本外
rm = rng.normal(0.0003, 0.01, n)                       # 指数(期货)日收益
beta = 1.0 + np.cumsum(rng.normal(0, 0.02, n))         # 真实 beta：随机游走
r = beta * rm + rng.normal(0, 0.01, n)                 # 个股日收益
Z = rm[:, None]
def kf_beta(p, end=n):
    return kalman_filter(r[:end], Z[:end], np.eye(1), np.diag([np.exp(p[0])]), np.exp(p[1]),
                         np.zeros(1), 1e6 * np.eye(1))
res = minimize(lambda p: -loglik(kf_beta(p, ntr)), np.log([1e-4, 1e-4]), method="Nelder-Mead")
print(f"\n[例二] 训练期 MLE: sigma_beta={np.sqrt(np.exp(res.x[0])):.4f} (真值 0.02), "
      f"sigma_e={np.sqrt(np.exp(res.x[1])):.4f} (真值 0.01)")
kf = kf_beta(res.x)
ss, _ = kalman_smoother(kf, Z, np.eye(1))
te = slice(ntr, n)
def report(name, b):
    rmse = np.sqrt(np.mean((b[te] - beta[te]) ** 2))
    extra = np.mean(((b[te] - beta[te]) * rm[te]) ** 2) / 0.01 ** 2   # 对冲误差额外方差 / 特质方差
    print(f"  {name:22s} beta RMSE={rmse:.3f}   额外对冲误差方差/特质方差={extra:6.1%}")
for w in (20, 60, 120):
    br = np.full(n, np.nan)
    for t in range(w, n):                              # 只用 t-1 及以前的数据
        br[t] = np.sum(rm[t - w:t] * r[t - w:t]) / np.sum(rm[t - w:t] ** 2)
    report(f"滚动 {w} 日 OLS", br)
report("训练期静态 OLS", np.full(n, np.sum(rm[:ntr] * r[:ntr]) / np.sum(rm[:ntr] ** 2)))
report("Kalman 预测 beta_t|t-1", kf["sp"][:, 0])
report("Kalman 平滑 beta_t|T(前视)", ss[:, 0])

# ================= 例三：配对交易的动态对冲比率 =================
rng = np.random.default_rng(1)
x = np.log(50) + np.cumsum(rng.normal(0, 0.015, n))   # X 的对数价格
bt_true = 1.0 + np.cumsum(rng.normal(0, 0.002, n))    # 真实对冲比率，缓慢漂移
s = np.zeros(n)
for t in range(1, n):
    s[t] = 0.9 * s[t - 1] + rng.normal(0, 0.008)       # 可交易的均值回复价差
y = 0.3 + bt_true * x + s                              # Y 的对数价格

def backtest(z, b, entry=1.0):
    """z>entry 做空价差(空 Y 多 b 份 X)，z<-entry 做多；z 回到 0 平仓。b 用 t 时已知值"""
    pos = np.zeros(n)
    for t in range(ntr, n - 1):
        if pos[t - 1] == 0:
            pos[t] = -1 if z[t] > entry else (1 if z[t] < -entry else 0)
        else:
            pos[t] = 0 if pos[t - 1] * z[t] >= 0 else pos[t - 1]
    pnl = (pos[:-1] * (np.diff(y) - b[:-1] * np.diff(x)))[ntr:]
    return pnl.mean() / pnl.std() * np.sqrt(252), int((np.diff(pos[ntr:]) != 0).sum())

Z2 = np.column_stack([np.ones(n), x])
def kf_pair(sa2, sb2, se2):                            # 观测噪声当作白噪声的"朴素"设定
    return kalman_filter(y, Z2, np.eye(2), np.diag([sa2, sb2]), se2, np.zeros(2), 1e6 * np.eye(2))
def nl_naive(p):
    sa2, sb2, se2 = np.exp(p)
    return -loglik(kalman_filter(y[:ntr], Z2[:ntr], np.eye(2), np.diag([sa2, sb2]), se2,
                                 np.zeros(2), 1e6 * np.eye(2)))
p_naive = np.exp(minimize(nl_naive, np.log([1e-6, 1e-6, 1e-4]), method="Nelder-Mead",
                          options=dict(maxiter=3000)).x)
print(f"\n[例三] 朴素模型 MLE 标准差: sigma_alpha={np.sqrt(p_naive[0]):.4f}, "
      f"sigma_beta={np.sqrt(p_naive[1]):.4f}, sigma_e={np.sqrt(p_naive[2]):.4f}")

def show(name, kf_, z):
    ac = np.corrcoef(z[ntr:-1], z[ntr + 1:])[0, 1]
    sh, ntrade = backtest(z, kf_["sp"][:, 1])
    print(f"  {name:30s} 信号一阶自相关={ac:5.2f}  年化Sharpe={sh:5.2f}  开平仓次数={ntrade}")

kf_n = kf_pair(*p_naive)
show("朴素 MLE(alpha,beta 都游走)", kf_n, kf_n["v"] / np.sqrt(kf_n["V"]))
kf_d = kf_pair(0.0, 1e-6, np.var(s[:ntr]))            # 经验做法：beta 只允许缓慢变化
show("固定 delta=1e-6(alpha 固定)", kf_d, kf_d["v"] / np.sqrt(kf_d["V"]))

# 正确设定：把价差写成 AR(1) 状态 (11.11.4 节的思路)，状态 = (alpha, beta_t, s_t)
Z3 = np.column_stack([np.ones(n), x, np.ones(n)])
def kf_ar(p, end=n):
    phi = np.tanh(p[0]); sb2, ss2, h = np.exp(p[1:])
    P0 = np.diag([1e6, 1e6, ss2 / (1 - phi ** 2)])
    return kalman_filter(y[:end], Z3[:end], np.diag([1, 1, phi]), np.diag([0, sb2, ss2]), h,
                         np.zeros(3), P0), phi, ss2
r3 = minimize(lambda p: -loglik(kf_ar(p, ntr)[0]), [1.0, np.log(1e-6), np.log(1e-4), np.log(1e-6)],
              method="Nelder-Mead", options=dict(maxiter=5000, xatol=1e-6, fatol=1e-6))
kf3, phi, ss2 = kf_ar(r3.x)
print(f"  AR(1) 价差模型 MLE: phi={phi:.3f} (真值 0.9), sigma_beta={np.sqrt(np.exp(r3.x[1])):.4f} (真值 0.002), "
      f"sigma_s={np.sqrt(ss2):.4f} (真值 0.008)")
show("AR(1) 价差状态模型", kf3, kf3["sf"][:, 2] / np.sqrt(ss2 / (1 - phi ** 2)))
sh_o, _ = backtest(s / np.sqrt(0.008 ** 2 / 0.19), bt_true)
print(f"  {'上帝视角(真实价差与真实 beta)':30s} 年化Sharpe={sh_o:5.2f}")

关键输出：

[例一] OLS 系数: [0.3775 1.1121] 标准误: [0.6315 0.1604]
[例一] KF s_T|T: [0.3775 1.1121] 标准差: [0.6315 0.1604]

[例二] 训练期 MLE: sigma_beta=0.0171 (真值 0.02), sigma_e=0.0103 (真值 0.01)
  滚动 20 日 OLS            beta RMSE=0.244   额外对冲误差方差/特质方差=  6.2%
  滚动 60 日 OLS            beta RMSE=0.191   额外对冲误差方差/特质方差=  3.6%
  滚动 120 日 OLS           beta RMSE=0.192   额外对冲误差方差/特质方差=  3.7%
  训练期静态 OLS              beta RMSE=0.181   额外对冲误差方差/特质方差=  3.3%
  Kalman 预测 beta_t|t-1   beta RMSE=0.155   额外对冲误差方差/特质方差=  2.3%
  Kalman 平滑 beta_t|T(前视) beta RMSE=0.105   额外对冲误差方差/特质方差=  1.1%

[例三] 朴素模型 MLE 标准差: sigma_alpha=0.0089, sigma_beta=0.0018, sigma_e=0.0021
  朴素 MLE(alpha,beta 都游走)         信号一阶自相关=-0.01  年化Sharpe= 0.51  开平仓次数=374
  固定 delta=1e-6(alpha 固定)        信号一阶自相关= 0.77  年化Sharpe= 1.69  开平仓次数=145
  AR(1) 价差模型 MLE: phi=0.887 (真值 0.9), sigma_beta=0.0019 (真值 0.002), sigma_s=0.0088 (真值 0.008)
  AR(1) 价差状态模型                   信号一阶自相关= 0.89  年化Sharpe= 1.08  开平仓次数=32
  上帝视角(真实价差与真实 beta)             年化Sharpe= 2.35

读输出：

• 例一：Kalman 滤波最后一步的状态与方差和 OLS 的系数与标准误逐位相同，这就是例 11.2 中 SsfPack 与 OLS 一致的原因。
• 例二：只用 500 天训练数据，MLE 就把 \(\sigma_\beta\)、\(\sigma_e\) 估到真值附近。样本外，Kalman 预测 β 的跟踪误差（0.155）低于所有滚动窗口和静态 OLS；对冲后多出来的残差方差从滚动 20 日的 6.2% 降到 2.3%。滚动窗口面临两难：窗口短则噪声大（20 日），窗口长则跟不上（120 日），Kalman 滤波用信噪比自动找到折中。平滑值看起来最好（1.1%），但它用了未来数据，在回测里用它就是前视偏差。
• 例三是本章最值得体会的地方。
  • 朴素设定把价差当白噪声观测误差，并让 α 也做随机游走。MLE 发现"让 α 跟着价差走"能让似然最大，于是 \(\hat\sigma_\alpha\) 比 \(\hat\sigma_e\) 还大，价差被吸进了状态：标准化新息的自相关约为 0，已经没有均值回复可以交易，频繁开平仓（374 次）。似然最优的滤波，未必产生可交易的信号，原因是模型把 AR(1) 价差错设成了白噪声。
  • 经验做法：固定 α、把 β 的新息方差压到很小（\(\delta=10^{-6}\)），强迫 β 只能缓慢变化，价差留在新息里，信号自相关 0.77。业界常这样调 \(\delta\)，但它是拍脑袋的超参数，需要样本外检验。
  • 正确设定：按 11.11.4 节的思路，把 AR(1) 价差作为第三个状态分量。MLE 只用训练期就把 \(\phi\)、\(\sigma_\beta\)、\(\sigma_s\) 都估得接近真值，滤波价差 \(s_{t|t}\) 的自相关 0.89，与真实价差的性质一致。

  金融直觉：为什么朴素模型会"吃掉"价差？似然只奖励"预测误差小"。如果允许 α 随机游走，滤波就会让 α 每天追着价差跑，新息变得很小、似然很高，但价差这个可交易的偏离也被当成"均衡水平变了"吸收掉了。好比把公允价值每天都调到市价，估值偏差永远为零，也就永远没有买卖信号。所以设定状态空间模型时，先想清楚哪部分是"会回归的偏离"、哪部分是"真的在变的参数"，再让似然去估参数。

  • Sharpe 只是一次模拟、不计交易成本的结果，三者的排序会随随机种子变化，不要把它们当作结论；稳健的结论是信号自相关和参数恢复。"上帝视角" 2.35 是用真实价差和真实 β 能达到的上限，与它的差距反映了从噪声中识别 \(\gamma_t\) 与 \(s_t\) 的固有难度：\(\gamma_t\) 每天漂移 0.002，乘上 \(x_t\approx4\) 后和价差新息（0.008）同量级，两者本来就难以分开。
• 实盘时还要加上：交易成本、对冲比率变动带来的调仓、\(\phi\) 与 \(\sigma_s\) 的滚动重估、价差半衰期 \(\ln0.5/\ln\phi\)（这里约 6.6 天）对持仓周期的约束。配对交易的协整基础见第 08b 章。

本章小结

一般线性高斯状态空间模型由状态方程 \(s_{t+1}=d_t+T_ts_t+R_t\eta_t\) 和观测方程 \(y_t=c_t+Z_ts_t+e_t\) 组成。时变系数 CAPM、线性回归（此时 Kalman 滤波即 RLS）、ARMA（Akaike、Harvey、Aoki 等多种表示）、带 ARMA 误差的回归、趋势–季节–周期分解都可以纳入，状态向量可以按需拼接。一般 Kalman 滤波 (11.64) 与标量版本结构相同，增益 \(K_t=T_t\Sigma_{t|t-1}Z_t'V_t^{-1}\)；固定区间平滑器和扰动平滑都用后向递推 \(q_{t-1},M_{t-1}\)；缺失值令 \(v_t=0,K_t=0\)，部分缺失用选择矩阵；多步预测就是把未来当作缺失值。实证上，GM 的 β 基本不变，强生 EPS 的季节模式随时间演变。在量化中，Kalman 滤波是估计时变 β 和动态对冲比率的标准工具，但模型设定（尤其是价差的动态）决定了信号是否可交易，而平滑值不能进入回测。

练习

基础

1. （原书习题 11.1）ARMA(1,1) \(y_t-0.8y_{t-1}=a_t+0.4a_{t-1}\)，\(a_t\sim N(0,0.49)\)，分别用 Akaike、Harvey、Aoki 方法写成状态空间形式。 提示：按本书记号 \(\theta_1=-0.4\)，\(m=2\)。Harvey：\(T=\begin{bmatrix}0.8&1\\0&0\end{bmatrix}\)，\(R=(1,0.4)'\)；Akaike：\(\psi_1=\phi_1-\theta_1=1.2\)，\(T=\begin{bmatrix}0&1\\0&0.8\end{bmatrix}\)，\(R=(1,1.2)'\)。
2. 把例 11.3 的模型改成"局部线性趋势 + 季节"，写出 5 维状态向量的 \(T\)、\(R\)、\(Z\)。
3. 验证 (11.63) 与 (11.56)+(11.62) 等价。
4. 对 11.11.3 节的 ARMA(2,1) 例子，解 Lyapunov 方程 \(\Sigma=T\Sigma T'+\sigma_a^2RR'\)，验证 mSigma 的数值 4.0607、\(-1.4874\)、0.5731。
5. 对时变 CAPM，写出某天市场数据缺失（\(r_{M,t}\) 未知）与个股数据缺失（\(r_t\) 未知）时滤波该怎样处理，二者有何不同。 提示：两种情况都没有观测方程可用，令 \(v_t=0,K_t=0\)；但若只有 \(r_{M,t}\) 缺失而 \(r_t\) 已知，也可以把 \(r_{M,t}\) 也作为状态建模——那是另一个模型。

进阶

6. 证明：对 (11.46)，扩散初始化下 \(s_{T|T}\) 等于 OLS 估计，\(\Sigma_{T|T}=\sigma_e^2(X'X)^{-1}\)。 提示：用信息形式 \(\Sigma_{t|t}^{-1}=\Sigma_{t|t-1}^{-1}+x_tx_t'/\sigma_e^2\)。
7. 推导 (11.78)：\(e_{t|T}=H_t(V_t^{-1}v_t-K_t'q_t)\)。
8. 在例二代码中加入"一天 β 跳变 +0.5"的结构突变，用扰动平滑 \(\eta_{t|T}\) 的标准化值定位突变日，并比较 Kalman 与滚动 60 日 OLS 适应突变的速度。
9. 在例三代码中加入每次开平仓 5 个基点的成本，再比较三种设定；并把入场阈值从 1.0 调到 1.5、2.0，观察 Sharpe 与换手的变化。
10. （原书习题 11.3）Pfizer 与 S&P 500 1990–2003 月度超额收益：拟合固定系数市场模型；拟合时变 CAPM，估计 \(\alpha_t,\beta_t\) 新息的标准差，画平滑估计。没有数据时可以用本章例二的模拟框架代替，取月度参数。
11. （原书习题 11.4、11.5）AR(3) 加测量误差 \(y_t=x_t+e_t\) 写成状态空间形式；若 \(E(e_t)=c\ne0\) 如何修改？对美国 PPI（1947-01 至 2009-11）对数差分去均值后的 AR(3)，假设有独立测量误差，估计状态新息方差与 \(\sigma_e^2\)，画平滑 \(x_t\) 与滤波响应残差。 提示：\(E(e_t)=c\) 时令 \(c_t=c\) 作为观测方程截距（或把 \(c\) 作为常数状态）。

原书推荐习题：11.3（时变 CAPM，量化中最常用的应用）、11.1（三种 ARMA 状态空间表示）、11.4、11.5（带测量误差的 AR 模型，练习构造系统矩阵与处理非零均值误差）。

原书对照

（原书印刷页码约等于 PDF 页码减 20。）