第 11a 章 局部水平模型与 Kalman 滤波

本章对应 Tsay 原书第 11 章 11.1 节（局部趋势模型）。金融数据里有很多量看不见：真实波动率、公允价格、时变 β、经济周期所处的阶段。能看到的只是带噪声的观测。状态空间模型（state-space model）把"看不见的状态怎样演化"和"观测怎样由状态产生"分开写，Kalman 滤波（Kalman filter）则给出在线更新状态估计的递推公式。本章用最简单的局部水平模型把全部思想讲清楚：滤波、预测、平滑、缺失值、初始化和极大似然估计。第 11b 章把这些结果推广到一般线性状态空间模型，并用时变 β 和动态对冲比率做完整的量化例子。

学习目标

1. 写出局部水平模型的观测方程与状态方程，说明它与 ARIMA(0,1,1)（简单指数平滑）的等价关系，会在两组参数之间换算。
2. 区分滤波、预测、平滑三类推断，知道哪一类可以用于回测、哪一类会引入未来信息。
3. 能用多元正态条件分布公式（定理 11.1）独立推导标量 Kalman 滤波的预测步与更新步，理解 Kalman 增益是"状态对新息的回归系数"。
4. 掌握后向平滑递推 \(q_{t-1},M_{t-1}\) 的来历，会处理缺失值和扩散初始化。
5. 会用预测误差分解写出精确对数似然，手写 Kalman 滤波完成参数的极大似然估计，并用标准化新息做模型诊断。

读前导读

这一章在解决什么问题。 很多你真正关心的量是看不见的：一只股票"真实"的波动率、公允价值、此刻的 β。你看到的只是带噪声的读数：今天的已实现波动率、最新成交价、一段窗口的回归斜率。本章回答的问题是：每来一个新读数，该把对真实值的估计调整多少？调太多，估计会被噪声带着跑；调太少，又跟不上真实值的变化。Kalman 滤波给出的答案是：按"我对当前估计有多不确定"与"新读数有多吵"的比例来调，这个比例就是 Kalman 增益。

和你熟悉的东西对应：EWMA（RiskMetrics 的 \(\lambda\)、技术分析的指数均线）就是本章模型在稳态下的特例，本章会告诉你平滑常数应该怎么从数据里估出来，而不是拍脑袋取 0.94。Kalman 增益的公式本质上是一个回归 β——"用新息去回归真实状态"，与 CFA 里"条件期望 = 回归预测"是同一件事。贝叶斯的视角也可以帮助理解：先验（昨天的估计）+ 新证据（今天的读数）= 后验，权重由两者的精度决定。

需要先想起来的数学。

• 方差与协方差的运算规则。 \(\mathrm{Var}(X+Y)=\mathrm{Var}(X)+\mathrm{Var}(Y)+2\mathrm{Cov}(X,Y)\)，独立时交叉项为 0；常数可以提出协方差：\(\mathrm{Cov}(aX,Y)=a\mathrm{Cov}(X,Y)\)。本章 11.2 节所有推导都只用到这两条。见 第 00 册第 07 章 概率中的分析工具。
• 正态分布的条件期望就是回归。 若 \((X,Y)\) 二元正态，则 \(E(X|Y)=\mu_X+\beta(Y-\mu_Y)\)，\(\beta=\mathrm{Cov}(X,Y)/\mathrm{Var}(Y)\)，条件方差 \(=\mathrm{Var}(X)(1-\rho^2)\)。例：\(\mathrm{Var}(X)=1\)、\(\mathrm{Var}(Y)=4\)、\(\mathrm{Cov}=1\)，观测到 \(Y\) 比均值高 2，则 \(X\) 的条件期望上调 \(0.25\times2=0.5\)，条件方差 \(1-1/4=0.75\)。这是定理 11.1 的核心。见第 00 册第 07 章。
• 递推与几何衰减。 形如 \(x_{t+1}=(1-K)x_t+Ky_t\) 的递推把历史观测按 \(K(1-K)^j\) 加权，就是 EWMA；半衰期 \(\ln0.5/\ln(1-K)\)。见 第 00 册第 04 章 级数与收敛。
• 一元二次方程。 稳态增益来自 \(p^2-qp-q=0\)，用求根公式取正根即可。两根之积等于常数项/首项系数，这解释了 11.2 节"两根互为倒数"。
• 变量替换与雅可比行列式。 把随机向量 \(\boldsymbol y\) 线性变换为 \(\boldsymbol v=\boldsymbol K\boldsymbol y\) 时，密度要乘以 \(|\det\boldsymbol K|^{-1}\)；行列式为 1 时密度不变。只在 11.5 节推论一用到。见 第 00 册第 05 章 多元微积分与优化。

另外，记号 \(\mu_{t|j}\) 读作"用截至 \(j\) 期的信息对 \(\mu_t\) 的估计"，竖线右边是"已知什么"，左边是"估计哪一期"。\(\mu_{t|t-1}\) 是预测，\(\mu_{t|t}\) 是滤波，\(\mu_{t|T}\) 是平滑。

怎么读这一章。 核心必读：11.1（模型）、11.3（三类推断与定理 11.1）、11.4（Kalman 滤波全部四小节）、11.9（似然与诊断）以及 Python 示例。11.2 的 ARIMA 等价性要看结论（"局部水平模型 = 指数平滑"），推导可以配合下面的讲解框读。11.5 和 11.6 的平滑推导较长，第一次只需掌握"先前向滤波、再后向平滑"的算法和"平滑不能用于回测"的结论。11.7 缺失值、11.8 初始化读结论即可。建议顺序：11.1 → 11.3 → 11.4 → 示例前半（滤波与估计）→ 11.2 → 11.6 → 其余各节。

11.1 一个动机：已实现波动率里的噪声

第 03b 章 3.15 节介绍过已实现波动率（realized volatility）：用日内高频收益的平方和来度量当天的波动。它比日收益平方精确得多，但仍不是"真实"波动——买卖价反弹、价格离散化、非同步交易等市场微观结构噪声（market microstructure noise）都会混进去。一个自然的建模方式是：

• 真实的对数波动率 \(\mu_t\) 缓慢变化，看不见；
• 我们观测到的对数已实现波动率 \(y_t\) 等于 \(\mu_t\) 加上一个测量误差。

把这两句话写成公式，就是局部趋势模型（local trend model），Durbin & Koopman 称之为局部水平模型（local-level model）：

其中 \(\{e_t\}\) 与 \(\{\eta_t\}\) 是相互独立的高斯白噪声，初值 \(\mu_1\) 给定或服从已知分布且与噪声独立。

术语如下：

• \(\mu_t\) 叫状态（state），这里是一个随机游走，也称趋势（trend）。它不可观测。
• (11.1) 叫观测方程（observation equation），\(e_t\) 是测量误差（measurement error）。
• (11.2) 叫状态方程或状态转移方程（state transition equation），\(\eta_t\) 是状态新息。

注意 \(y_t\) 本身没有任何直接的动态，它的序列相关完全来自状态 \(\mu_t\)。参数 \(\sigma_e\) 衡量噪声大小，\(\sigma_\eta\) 衡量真实波动率变化有多快。两者之比 \(q=\sigma_\eta^2/\sigma_e^2\) 叫信噪比（signal-to-noise ratio），后面会看到它决定了滤波"相信新数据"的程度。

局部水平模型是 Harvey (1993) 结构时间序列模型（structural time series model）中最简单的一种，也是一般线性高斯状态空间模型的特例。原书推荐的进一步阅读有 Durbin & Koopman (2001)、Kim & Nelson (1999，经济应用与区制转换)、Harvey (1993)、Hamilton (1994)、Shumway & Stoffer (2000)、West & Harrison (1997，贝叶斯预测) 等。

11.2 局部水平模型就是 ARIMA(0,1,1)

11.2.1 推导

先看两个极端。若 \(\sigma_e=0\)，则 \(y_t=\mu_t\) 就是随机游走，即 ARIMA(0,1,0)。若 \(\sigma_e>0\)，对 \(y_t\) 做一阶差分：由 (11.2)，\(\mu_t=\mu_{t-1}+\eta_{t-1}\)，于是

\(w_t\) 是高斯的，且

自协方差在滞后 1 之后截尾，所以 \(w_t\) 是 MA(1)，\(y_t\) 是 ARIMA(0,1,1)：

让两边的方差和滞后 1 自协方差相等：

两式相除得到关于 \(\theta\) 的二次方程 \(\theta^2-(2+q)\theta+1=0\)（\(q=\sigma_\eta^2/\sigma_e^2\)），两根互为倒数，取 \(|\theta|<1\) 的那个（可逆解），再由 (11.5) 得 \(\sigma_a^2\)。因为 \(\sigma_e^2>0\)，(11.5) 说明这样得到的 \(\theta\) 一定为正。

推导拆解：\(w_t\) 的矩怎么算。\(w_t=\eta_{t-1}+e_t-e_{t-1}\)，三项两两独立。 方差：独立变量之和的方差等于方差之和，\(\mathrm{Var}(w_t)=\sigma_\eta^2+\sigma_e^2+\sigma_e^2\)。 滞后 1 协方差：\(w_{t-1}=\eta_{t-2}+e_{t-1}-e_{t-2}\)。把两式逐项配对，只有 \(w_t\) 中的 \(-e_{t-1}\) 与 \(w_{t-1}\) 中的 \(+e_{t-1}\) 是同一个随机变量，其余都独立，所以协方差 \(=-\sigma_e^2\)。 滞后 2 及以上：\(w_t\) 与 \(w_{t-2}\) 没有任何共同项，协方差为 0。 右边 MA(1) 的矩：\((1-\theta B)a_t=a_t-\theta a_{t-1}\)，方差 \((1+\theta^2)\sigma_a^2\)，滞后 1 协方差 \(-\theta\sigma_a^2\)。两边对齐即 (11.4)(11.5)。 二次方程：(11.4) 除以 (11.5) 得 \((1+\theta^2)/\theta=(2\sigma_e^2+\sigma_\eta^2)/\sigma_e^2=2+q\)，两边乘 \(\theta\) 整理即得。由韦达定理，两根之积为 1，所以一根在 \((0,1)\)、另一根大于 1。

白话解释：差分后出现负的滞后 1 自相关，是"噪声"的指纹。今天的读数被噪声推高，明天的差分就会显得偏低，因为噪声不会延续。买卖价反弹让成交价收益出现负自相关，也是同一道理。\(\theta\) 越接近 1，说明噪声越占主导。

这正是原书第 2 章提到的简单指数平滑（simple exponential smoothing）模型。

11.2.2 反方向

反过来，给定 \(\theta>0\) 的 ARIMA(0,1,1)，由 (11.5) 得 \(\sigma_e^2=\theta\sigma_a^2\)，代入 (11.4) 得 \(\sigma_\eta^2=(1-\theta)^2\sigma_a^2\)，两者都为正，可以写回局部水平模型。若 \(\theta<0\)，则 \(\sigma_e^2<0\) 不可能，此时模型仍能写成状态空间形式，只是不带观测误差的另一种形式（第 11b 章的 Akaike、Harvey 表示）。

一个 ARIMA 模型有多种状态空间写法。仅从数据出发，选 ARIMA 形式还是状态空间形式并不关键，取决于分析目的和对问题的理解：如果你关心的是"真实波动率是多少""噪声占多大比例"，状态空间形式的参数更有解释力。

11.2.3 例 11.1：Alcoa 已实现波动率

数据为 Alcoa 股票 2003-01-02 至 2004-05-07 的日已实现波动率，共 340 个观测。日已实现波动率是日内 10 分钟对数收益（%）的平方和，不含隔夜收益和开盘后第一个 10 分钟收益，数据来自 NYSE TAQ，分析其对数。ARIMA 拟合为

\(\hat\theta\) 的标准误为 0.028。残差 \(Q(12)=12.4\)（p 值 0.33），平方残差 \(Q(12)=8.2\)（p 值 0.77），没有 ARCH 效应。

由于 \(\hat\theta>0\)，可以换算成局部水平模型：\(\sigma_e^2=0.858\times0.5184^2=0.2306\)，\(\sigma_\eta^2=(1+0.858^2)\times0.2687-2\times0.2306=0.0054\)，即 \(\sigma_e=0.480\)，\(\sigma_\eta=0.0736\)。直接用极大似然估计局部水平模型得到 \(\hat\sigma_\eta=0.0735\)，\(\hat\sigma_e=0.4803\)，两者几乎一致。

结论很有意思：测量误差的方差（0.23）是状态新息方差（0.0054）的四十多倍。日内高频收益确实受到相当大的测量误差影响，单日的已实现波动率读数很"毛"，需要平滑。

11.3 三类推断与一条关键定理

11.3.1 滤波、预测、平滑

设 \(F_t=\{y_1,\dots,y_t\}\) 是到 \(t\) 时为止的信息，并暂时假设参数 \(\sigma_e,\sigma_\eta\) 已知。对不可观测的 \(\mu_t\)，有三类推断：

• 滤波（filtering）：给定 \(F_t\) 估计 \(\mu_t\)，即从当前及过去的观测中去除测量误差。
• 预测（prediction）：给定 \(F_t\) 预测 \(\mu_{t+h}\) 或 \(y_{t+h}\)，\(h>0\)。
• 平滑（smoothing）：给定全部样本 \(F_T\)（\(T>t\)）估计 \(\mu_t\)。

原书有一个好记的比喻：读一张字迹潦草的便条，滤波是根据已经读过的内容辨认当前这个字，预测是猜下一个字，平滑是读完全文之后回头再辨认某个字。

对量化交易的含义：滤波和预测只用到 \(t\) 时已有的信息，可以在回测和实盘中使用；平滑用到了 \(t\) 之后的数据，适合事后分析（例如研究历史上 β 怎样变化），但如果把平滑值当作交易信号放进回测，就犯了前视偏差（look-ahead bias）。这是状态空间模型在量化研究中最常见的错误之一。

11.3.2 记号与新息

记

1 步预测误差（one-step-ahead forecast error）也叫新息（innovation）：

由 (11.1)，\(y_{t|t-1}=\mu_{t|t-1}\)，所以

新息有两个基本性质：\(E(v_t)=0\)；\(v_t\) 与过去观测不相关，\(\mathrm{Cov}(v_t,y_j)=0\)（\(j<t\)），在正态假设下就是独立。因此 \(V_t\) 不依赖 \(F_{t-1}\) 的具体取值，\(\mathrm{Var}(v_t|F_{t-1})=\mathrm{Var}(v_t)\)。

更重要的是，\(y_t\) 和 \(v_t\) 在给定 \(F_{t-1}\) 时一一对应（差一个已知常数 \(\mu_{t|t-1}\)），所以

新的观测带来的"新东西"全部装在 \(v_t\) 里，并且 \(v_t\) 与旧信息无关。这正是 Kalman 滤波能递推的根本原因。

11.3.3 定理 11.1（多元正态的条件分布）

设 \(x,y,z\) 联合正态，\(w=(y',z')'\) 的协方差阵 \(\Sigma_{ww}\) 非奇异，且 \(\Sigma_{yz}=0\)（\(y\) 与 \(z\) 不相关）。则

1. \(E(x|y)=\mu_x+\Sigma_{xy}\Sigma_{yy}^{-1}(y-\mu_y)\)；
2. \(\mathrm{Var}(x|y)=\Sigma_{xx}-\Sigma_{xy}\Sigma_{yy}^{-1}\Sigma_{yx}\)；
3. \(E(x|y,z)=E(x|y)+\Sigma_{xz}\Sigma_{zz}^{-1}(z-\mu_z)\)；
4. \(\mathrm{Var}(x|y,z)=\mathrm{Var}(x|y)-\Sigma_{xz}\Sigma_{zz}^{-1}\Sigma_{zx}\)。

第 1、2 条是多元正态条件分布的标准结果（第 02 册、第 03 册），也就是总体线性回归：条件均值是 \(x\) 对 \(y\) 的回归，条件方差是回归残差方差。第 3、4 条是关键：当新加入的变量 \(z\) 与已有信息 \(y\) 不相关时，\(\Sigma_{ww}\) 是分块对角阵，求逆可以分开做，于是条件均值和条件方差可以增量更新——旧的结论保留，只加上 \(z\) 带来的修正项。把 \(y\) 换成 \(F_{t-1}\)、\(z\) 换成 \(v_t\)，就是 Kalman 滤波。

白话解释：标量情形下，第 1 条就是 \(E(x|y)=\mu_x+\beta(y-\mu_y)\)，\(\beta=\mathrm{Cov}(x,y)/\mathrm{Var}(y)\)——你在 CFA 里熟悉的回归预测；第 2 条就是"残差方差 = 总方差 − 被解释的部分"。第 3、4 条说的是多元回归里的一个事实：如果两个解释变量互不相关，多元回归的系数就等于分别做两个一元回归的系数，可以先用 \(y\) 回归，再把 \(z\) 的贡献"加上去"，不用重新估计。

金融直觉：想象你在给一家公司估值。昨天你已经根据全部历史信息得出估计（对应 \(E(x|y)\)）。今天公布了一份财报，你只关心其中"超出市场预期"的部分（对应 \(z\)，它与已有信息不相关，就像盈利意外）。新的估值 = 昨天的估值 + 敏感度 × 盈利意外，敏感度就是回归系数。你不必把全部历史信息重新处理一遍。Kalman 滤波的新息 \(v_t\) 就是"读数意外"。

11.4 Kalman 滤波

11.4.1 更新步

目标：已知 \(\mu_t|F_{t-1}\sim N(\mu_{t|t-1},\Sigma_{t|t-1})\)，观测到新数据 \(y_t\) 后得到 \(\mu_t|F_t\)。由于 \(F_t=\{F_{t-1},v_t\}\)，只需知道 \((\mu_t,v_t)'\) 在给定 \(F_{t-1}\) 时的联合分布。

均值已知：\(E(\mu_t|F_{t-1})=\mu_{t|t-1}\)，\(E(v_t|F_{t-1})=0\)。方差也已知：\(\Sigma_{t|t-1}\) 与 \(V_t\)。只差协方差：

第二个等号用到 \(e_t\) 与 \(\mu_t\) 独立，以及 \(E[\mu_{t|t-1}(\mu_t-\mu_{t|t-1})|F_{t-1}]=0\)（给定 \(F_{t-1}\) 时 \(\mu_{t|t-1}\) 是常数）。于是

套用定理 11.1 的第 3、4 条：

其中

叫 Kalman 增益（Kalman gain）。它就是 \(\mu_t\) 对 \(v_t\) 的回归系数，决定新冲击 \(v_t\) 有多大比例被计入状态估计。\(K_t\) 介于 0 和 1 之间：

• 如果当前对状态很不确定（\(\Sigma_{t|t-1}\) 大）而观测噪声小（\(\sigma_e^2\) 小），\(K_t\to1\)，几乎完全相信新数据；
• 反过来，状态已经估得很准而观测很吵，\(K_t\to0\)，新数据只做小幅修正。

(11.11) 还说明更新后的方差一定变小：每看到一个新数据，不确定性都会下降。

推导拆解：把定理 11.1 第 3、4 条代入（\(x=\mu_t\)，\(z=v_t\)，\(\mu_z=0\)）。 均值：修正项 \(=\Sigma_{xz}\Sigma_{zz}^{-1}(z-0)=\mathrm{Cov}(\mu_t,v_t)\cdot V_t^{-1}\cdot v_t=(\Sigma_{t|t-1}/V_t)v_t\)，即 (11.10)。 方差：减掉 \(\Sigma_{xz}\Sigma_{zz}^{-1}\Sigma_{zx}=\Sigma_{t|t-1}^2/V_t\)，提出 \(\Sigma_{t|t-1}\) 得 \(\Sigma_{t|t-1}(1-K_t)\)，即 (11.11)。 数值例子：设 \(\mu_{t|t-1}=1.00\)，\(\Sigma_{t|t-1}=0.04\)，\(\sigma_e^2=0.16\)，今天读数 \(y_t=1.50\)。则 \(v_t=0.50\)，\(V_t=0.20\)，\(K_t=0.2\)，更新后 \(\mu_{t|t}=1.00+0.2\times0.50=1.10\)，\(\Sigma_{t|t}=0.04\times0.8=0.032\)。读数高出 0.5，但因为噪声方差是状态不确定性的 4 倍，只采信其中 20%。

金融直觉：贝叶斯的说法是"精度加权平均"。精度 = 方差的倒数：先验精度 \(1/0.04=25\)，读数精度 \(1/0.16=6.25\)，于是 \(\mu_{t|t}=(25\times1.00+6.25\times1.50)/(25+6.25)=1.10\)，后验精度 \(25+6.25=31.25\)，对应方差 0.032。两种算法结果一样。这和把两个分析师的盈利预测按各自历史准确度加权，是同一个逻辑。

11.4.2 预测步

由状态方程 (11.2)，\(\eta_t\) 与 \(F_t\) 独立，

状态往前走一步，不确定性增加 \(\sigma_\eta^2\)。观测到 \(y_{t+1}\) 后再做更新，如此交替，这就是 Kalman (1960) 的算法。

11.4.3 算法汇总

给定初值 \(\mu_1\sim N(\mu_{1|0},\Sigma_{1|0})\)，对 \(t=1,\dots,T\)：

每一步只做几次加减乘除，计算量是 \(O(T)\)，非常适合实盘的逐笔或逐 bar 更新。Kalman 滤波有很多推导方法（最小二乘、贝叶斯、投影），用定理 11.1 是最简短的一种。初值的选择见 11.8 节，未知参数 \(\sigma_e,\sigma_\eta\) 用极大似然估计，似然本身也由 Kalman 滤波算出（11.9 节）。

11.4.4 稳态：Kalman 滤波就是指数平滑

注意 (11.14) 中 \(K_t\) 和 \(\Sigma_{t|t-1}\) 的递推完全不依赖数据。对时不变模型，\(\Sigma_{t|t-1}\) 很快收敛到一个常数 \(P\)。记 \(p=P/\sigma_e^2\)，稳态满足

稳态增益

稳态时 (11.14) 的状态递推变成

这就是平滑常数为 \(K\) 的指数加权移动平均（EWMA）。与 11.2 节对照，\(K=1-\theta\)。换句话说，指数平滑是局部水平模型在稳态下的最优滤波，平滑常数由信噪比唯一决定：噪声越大（\(q\) 越小），\(K\) 越小，均线越"慢"。这为技术分析里"均线参数怎么选"提供了一个有统计依据的答案：先估计信噪比，再算 \(K\)。

推导拆解：稳态方程的由来与化简。 第一步，把 (11.14) 中 \(\Sigma\) 的递推合成一步：\(\Sigma_{t+1|t}=\Sigma_{t|t-1}(1-K_t)+\sigma_\eta^2=\frac{\Sigma_{t|t-1}\sigma_e^2}{\Sigma_{t|t-1}+\sigma_e^2}+\sigma_\eta^2\)（因为 \(1-K_t=\sigma_e^2/V_t\)）。稳态时左右两边都等于 \(P\)。 第二步，两边除以 \(\sigma_e^2\)，记 \(p=P/\sigma_e^2\)：\(p=\frac p{p+1}+q\)。两边乘 \((p+1)\)：\(p^2+p=p+q(p+1)\)，即 \(p^2-qp-q=0\)，取正根。 第三步，\(K=p/(1+p)\)。由 \(p^2=q(p+1)\) 得 \(p/(1+p)=q/p\)，再把 \(p\) 代入并对分母有理化，就得到正文的 \(K\)。 第四步，状态递推：\(\mu_{t+1|t}=\mu_{t|t-1}+K(y_t-\mu_{t|t-1})=(1-K)\mu_{t|t-1}+Ky_t\)。 数值感：例 11.1 中 \(q=0.0054/0.2306\approx0.023\)，\(K\approx0.14\)，对应 EWMA 的 \(\lambda=1-K\approx0.86\)，半衰期约 4.5 天（练习 2）。为什么 \(\Sigma_{t|t-1}\) 会收敛：它的递推里不含数据 \(y_t\)，是一个确定性的数列，从任何初值出发都会被拉到同一个不动点。

11.4.5 例 11.1 续

对 Alcoa 数据取 \(\Sigma_{1|0}=\infty\)（扩散初始化）、\(\mu_{1|0}=0\) 运行滤波。原书图 11.2(a) 显示滤波状态 \(\mu_{t|t}\) 比原序列平滑得多；图 11.2(b) 的 1 步预测误差 \(v_t\) 稳定地围绕 0 波动，这些都是样本外的 1 步预测误差。

11.5 预测误差的性质

11.5.1 新息是观测的线性变换

给定与数据独立的初值，\(v_t\) 是 \(y_1,\dots,y_t\) 的线性函数：

写成矩阵形式

\(K\) 是对角元为 1 的下三角阵，\(k_{i,i-1}=-K_{i-1}\)，\(k_{ij}=-(1-K_{i-1})(1-K_{i-2})\cdots(1-K_{j+1})K_j\)（\(j\le i-2\)）。由于增益序列只依赖 \(\Sigma_{1|0},\sigma_e^2,\sigma_\eta^2\)，矩阵 \(K\) 不依赖数据和 \(\mu_{1|0}\)。

由此得到两个重要推论。

推论一：正态下 \(\{v_t\}\) 相互独立。 变换 (11.15) 的雅可比行列式为 1，所以

推论二：Kalman 滤波实现了协方差阵的 Cholesky 分解。 记 \(\Omega=\mathrm{Cov}(y)\)，则

这和第 10 章用 Cholesky 分解把相关的收益正交化是同一回事。正是因为有这个分解，似然计算不需要对 \(T\times T\) 的 \(\Omega\) 求逆（11.9 节）。

白话解释：推论一说的是"换一种记账方式，概率不变"。\(y_1,\dots,y_T\) 彼此高度相关（都围绕同一个缓慢变化的 \(\mu_t\)），直接写联合密度需要一个 \(T\times T\) 协方差矩阵。新息 \(v_1,\dots,v_T\) 是从 \(y\) 中依次扣掉"可预测部分"后剩下的意外，彼此独立。\(\boldsymbol K\) 是对角元为 1 的下三角阵，行列式是对角元之积，等于 1，所以变换不拉伸也不压缩概率，\(p(\boldsymbol v)=p(\boldsymbol y)\)。 第一个等号到第二个等号用的是概率的乘法规则：联合密度 = 第一个的密度 × 给定第一个时第二个的条件密度 × ……；而给定 \(F_{j-1}\) 时，\(y_j\) 的条件分布就是 \(v_j\) 的分布平移一个已知常数。

金融直觉：这与第 10 章 Cholesky 的逻辑完全一样：把相关的收益依次对前面的收益回归、只保留残差，残差之间就不相关了。这里"前面的收益"换成了"过去的观测"。

11.5.2 状态估计误差的递推

记状态预测误差 \(x_t=\mu_t-\mu_{t|t-1}\)，则 \(\mathrm{Var}(x_t|F_{t-1})=\Sigma_{t|t-1}\)，且

其中 \(L_t=1-K_t=\sigma_e^2/V_t\)，\(x_1=\mu_1-\mu_{1|0}\)。第二式的推导：\(x_{t+1}=\mu_t+\eta_t-\mu_{t|t-1}-K_tv_t=x_t+\eta_t-K_t(x_t+e_t)\)。(11.16) 本身又是一个时变状态空间模型：状态是 \(x_t\)，观测是 \(v_t\)。下面推导平滑时会反复用到它。

11.6 状态平滑

11.6.1 平滑状态

目标是 \(\mu_t|F_T\sim N(\mu_{t|T},\Sigma_{t|T})\)。利用：

• \(\{v_t,\dots,v_T\}\) 互相独立，与 \(F_{t-1}\) 独立，均值 0，方差 \(V_j\)；
• \(F_T\) 等价于 \(\{F_{t-1},v_t,\dots,v_T\}\)。

反复使用定理 11.1 第 3 条：

需要的协方差可以由 (11.16) 算出。因为 \(v_j\) 与 \(\mu_{t|t-1}\) 不相关，\(\mathrm{Cov}(\mu_t,v_j)=\mathrm{Cov}(x_t,v_j)\)，而

（\(x_{t+1}\) 中只有 \(L_tx_t\) 一项与 \(x_t\) 相关，其余是未来的噪声。）代入得 \(\mu_{t|T}=\mu_{t|t-1}+\Sigma_{t|t-1}q_{t-1}\)，

它满足后向递推 \(q_{t-1}=v_t/V_t+L_tq_t\)，\(q_T=0\)（11.19）。于是后向平滑递推为

直观上，\(q_{t-1}\) 是"\(t\) 之后的全部新息对 \(\mu_t\) 的加权修正"，权重按 \(L_j<1\) 的乘积几何衰减——离 \(t\) 越远的未来观测，对 \(\mu_t\) 的修正越小。

推导拆解：(11.17) 是把定理 11.1 第 3 条连用多次。\(F_T=\{F_{t-1},v_t,v_{t+1},\dots,v_T\}\)，而 \(v_t,\dots,v_T\) 互相独立、又都与 \(F_{t-1}\) 独立，所以每加入一个 \(v_j\) 就加一项修正 \(\mathrm{Cov}(\mu_t,v_j)V_j^{-1}v_j\)，互不干扰。 协方差为什么是 \(\Sigma_{t|t-1}L_tL_{t+1}\cdots\)：由 (11.16)，\(v_j=x_j+e_j\)，而 \(x_j=L_{j-1}x_{j-1}+(\text{与 }x_t\text{ 无关的未来噪声})\)。从 \(x_t\) 走到 \(x_j\)，每走一步乘一个 \(L\)，所以 \(\mathrm{Cov}(x_t,x_j)=\mathrm{Var}(x_t)\prod_{i=t}^{j-1}L_i\)，再加上与 \(x_t\) 无关的 \(e_j\) 不改变协方差。 后向递推：把 (11.18) 的第一项单独写出，剩下的项都带一个公因子 \(L_t\)，提出后正好是 \(q_t\) 的定义，所以 \(q_{t-1}=v_t/V_t+L_tq_t\)。

白话解释：平滑就是"读完全文再回头改字"。滤波在 \(t\) 时刻只能看到过去，平滑还会参考 \(t\) 之后的读数：如果之后几天的读数都系统性偏高，说明 \(t\) 时刻的真实水平可能也被低估了，于是往上修正。正因为用了未来信息，它在回测里是禁用的。

11.6.2 平滑状态方差

由定理 11.1 第 4 条，

后向递推

实施时分两步：先前向运行 (11.14)，存下 \(v_t,V_t,K_t,\mu_{t|t-1},\Sigma_{t|t-1}\)；再从 \(t=T\) 往回跑 (11.20) 和 (11.23)。两遍都是 \(O(T)\)，不需要对 \(T\times T\) 矩阵求逆。

例 11.1 续：原书图 11.3 是滤波状态 \(\mu_{t|t}\) 及 95% 逐点区间，图 11.4 是平滑状态 \(\mu_{t|T}\) 及区间。平滑状态更平滑、区间更窄，因为它用了更多信息。\(\mu_{1|1}\) 的区间宽度取决于初值方差 \(\Sigma_{1|0}\)。

11.7 缺失值

设 \(y_{\ell+1},\dots,y_{\ell+h}\) 缺失。保持原来的时间刻度，对缺失段内的 \(t\)，

于是 \(E(\mu_t|F_{t-1})=\mu_{\ell+1|\ell}\)，\(\mathrm{Var}(\mu_t|F_{t-1})=\Sigma_{\ell+1|\ell}+(t-\ell-1)\sigma_\eta^2\)，也就是

这等价于在缺失时点令 \(v_t=0\)、\(K_t=0\) 继续运行 (11.14)：没有新数据，就没有新息，也没有增益；状态估计保持不变，不确定性每期增加 \(\sigma_\eta^2\)。缺失段一结束，第一笔新数据会以较大的增益修正状态。平滑时，缺失点上令 \(q_{t-1}=q_t\)、\(M_{t-1}=M_t\)（即 \(V_t^{-1}\) 项取 0、\(L_t=1\)）。

实用价值：股票停牌、不同交易所假日不一致、宏观数据与日频数据混合（混频）、数据源偶发断档，都可以不做任何插值，直接让 Kalman 滤波"空跑"过去。

11.8 初始化的影响

第一步有 \(v_1=y_1-\mu_{1|0}\)，\(V_1=\Sigma_{1|0}+\sigma_e^2\)，

令 \(\Sigma_{1|0}\to\infty\)，得 \(\mu_{2|1}=y_1\)，\(\Sigma_{2|1}=\sigma_e^2+\sigma_\eta^2\)。这相当于把 \(y_1\) 当作固定值、令 \(\mu_1\sim N(y_1,\sigma_e^2)\)。这种做法叫扩散初始化（diffuse initialization），表示对初始状态一无所知。

对平滑的影响：\(t=T,\dots,2\) 的平滑结果不受影响。对 \(\mu_1\)，\(L_1=\sigma_e^2/V_1\)，

建议：对 \(\mu_1\) 所知甚少时用扩散初始化。实现上取一个很大的有限数（如 \(10^7\)）即可，但计算似然时要去掉第 1 个观测的贡献（它的 \(V_1\) 近似无穷大）。如果不愿意引入方差无穷大的随机变量，也可以把 \(\mu_1\) 当作一个额外参数与 \(\sigma_e,\sigma_\eta\) 一起估计，这与第 02a 章、第 08a 章中的精确极大似然密切相关。

11.9 参数估计：预测误差分解

参数 \(\sigma_e,\sigma_\eta\) 未知时，用极大似然估计。由 11.5 节的推论一，预测误差分解（prediction error decomposition）给出

对任意一组参数，跑一遍 Kalman 滤波就得到 \(\{v_t,V_t\}\)，进而得到对数似然；再交给数值优化器（第 04 册）即可。有缺失值时，缺失点不计入求和。实务上常对方差做对数参数化（优化 \(\ln\sigma^2\)），保证估计值为正。

白话解释：(11.25) 的每一项 \(\ln V_t+v_t^2/V_t\) 和 GARCH 似然里的 \(\ln\sigma_t^2+a_t^2/\sigma_t^2\) 形式完全相同：\(v_t\) 是 1 步预测误差，\(V_t\) 是模型对这个误差方差的预测。好参数要同时做到两件事：预测误差小（\(v_t^2\) 小），并且对误差大小的判断诚实（\(V_t\) 既不夸大也不低估）。所以极大似然估计本质上是在找"1 步预测表现最好"的 \((\sigma_\eta,\sigma_e)\)——这和量化里用样本外预测误差选参数的直觉一致，只是这里每个预测都严格只用了过去的信息。 对数参数化：优化器在整条实数轴上搜索 \(\ln\sigma^2\)，取指数后自动为正，免去了"方差不能为负"的约束。

估计完成后要做模型诊断：标准化新息 \(\tilde v_t=v_t/\sqrt{V_t}\) 在模型正确时应近似为 iid 标准正态。原书对 Alcoa 数据得到 \(Q(25)=23.37\)（p 值 0.56），25 阶 ARCH LM 检验统计量 18.48（p 值 0.82），说明模型是充分的。

软件说明

原书用 Koopman, Shephard & Doornik (1999) 的 SsfPack（S-Plus/OX）。其记号与命令简述如下，读原书代码时可对照：系统矩阵 \(\delta,\Phi,\Omega,\Sigma\) 对应 mDelta、mPhi、mOmega、mSigma；SsfFit 做极大似然估计，KalmanFil 滤波，KalmanSmo 平滑，SsfMomentEst(task="STFIL") 与 (task="STSMO") 给出滤波/平滑状态及方差；mSigma 中初值方差填 \(-1\) 表示扩散初始化。对 Alcoa 数据，SsfFit 给出 \((\hat\sigma_\eta,\hat\sigma_e)=(0.07350827,0.48026284)\)，mOmega 对角元为 0.0054035 与 0.2306524。

在 Python 中，statsmodels.tsa.UnobservedComponents(y, level="llevel") 就是局部水平模型；更一般的模型可以继承 statsmodels.tsa.statespace.MLEModel 自行指定系统矩阵。不过手写一遍 (11.14) 是理解 Kalman 滤波最快的途径，下面的实战代码就是这样做的。

量化实战

应用场景

1. 波动率去噪：日已实现波动率、隐含波动率、成交量这类"读数很毛"的序列，用局部水平模型分离真实信号与测量噪声。滤波值可作为次日波动预测、仓位调整（波动率目标策略）或期权定价的输入。
2. 自适应均线：11.4.4 节说明，EWMA 的平滑常数不必拍脑袋，可由估计出的信噪比 \(q\) 算出 \(K\)。信噪比随品种、频率不同，用同一套代码估计即可。
3. 数据工程：停牌、跨市场日历不一致、数据断档，直接令 \(v_t=0,K_t=0\) 跑过去，不需要先插值。
4. 回测纪律：回测只能用 \(\mu_{t|t}\) 或 \(\mu_{t+1|t}\)；\(\mu_{t|T}\) 只能用于事后研究。参数 \(\sigma_e,\sigma_\eta\) 若用全样本估计，也会带来轻微的前视，严格的做法是滚动或扩展窗口重估。

Python 示例：手写 Kalman 滤波、平滑与极大似然

下面模拟一条与例 11.1 量级相近的"对数已实现波动率"序列（\(T=340\)，\(\sigma_\eta=0.07\)，\(\sigma_e=0.48\)），完成：手写滤波与似然、极大似然估计并与 statsmodels 对照、与 ARIMA(0,1,1) 参数互换、检验稳态增益 \(K=1-\theta\)、后向平滑、新息诊断和缺失值处理。

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

rng = np.random.default_rng(3)

# ---------- 1. 模拟"对数已实现波动率"：真实对数波动 mu_t 为随机游走，观测含微观结构噪声 ----------
T = 340
sig_eta, sig_e = 0.07, 0.48          # 与例 11.1 的估计值同量级
mu = 1.0 + np.cumsum(np.r_[0, rng.normal(0, sig_eta, T - 1)])
y = mu + rng.normal(0, sig_e, T)

# ---------- 2. 局部水平模型的 Kalman 滤波 (11.14) 与对数似然 (11.25) ----------
def kalman_local_level(y, s_eta2, s_e2, mu10=0.0, P10=1e7):
    n = len(y)
    m_pred, P_pred = np.empty(n), np.empty(n)   # mu_{t|t-1}, Sigma_{t|t-1}
    v, V, K = np.empty(n), np.empty(n), np.empty(n)
    m, P = mu10, P10
    for t in range(n):
        m_pred[t], P_pred[t] = m, P
        if np.isnan(y[t]):                       # 缺失值：v_t=0, K_t=0 (11.24)
            v[t], V[t], K[t] = 0.0, np.inf, 0.0
            P = P + s_eta2
            continue
        v[t] = y[t] - m
        V[t] = P + s_e2
        K[t] = P / V[t]
        m = m + K[t] * v[t]                      # mu_{t+1|t} = mu_{t|t}
        P = P * (1 - K[t]) + s_eta2
    return m_pred, P_pred, v, V, K

def loglik(params, y, skip=1):
    s_eta2, s_e2 = np.exp(params)                 # 对数参数化保证方差为正
    _, _, v, V, _ = kalman_local_level(y, s_eta2, s_e2)
    ok = np.isfinite(V)
    ok[:skip] = False                             # 扩散初始化：去掉第 1 个观测的似然贡献
    return -0.5 * np.sum(np.log(2 * np.pi * V[ok]) + v[ok] ** 2 / V[ok])

res = minimize(lambda p: -loglik(p, y), x0=np.log([0.01, 0.1]), method="Nelder-Mead")
s_eta2_hat, s_e2_hat = np.exp(res.x)
print(f"手写 Kalman MLE : sigma_eta={np.sqrt(s_eta2_hat):.4f}, sigma_e={np.sqrt(s_e2_hat):.4f}, logL={-res.fun:.2f}")

# ---------- 3. 与 statsmodels 的 UnobservedComponents 对照 ----------
uc = sm.tsa.UnobservedComponents(y, level="llevel").fit(disp=False)
print("statsmodels UC  : sigma_eta=%.4f, sigma_e=%.4f" % tuple(np.sqrt(uc.params[[1, 0]])))

# ---------- 4. ARIMA(0,1,1) 与局部水平模型的换算 (11.4)(11.5) ----------
ar = sm.tsa.ARIMA(y, order=(0, 1, 1)).fit()
theta = -ar.params[0]                            # statsmodels 记 (1+ma B)，Tsay 记 (1-theta B)
sa2 = ar.params[1]
se2_from_arima = theta * sa2
seta2_from_arima = (1 + theta ** 2) * sa2 - 2 * se2_from_arima
print(f"ARIMA(0,1,1)    : theta={theta:.3f}, sigma_a={np.sqrt(sa2):.4f} -> "
      f"sigma_e={np.sqrt(se2_from_arima):.4f}, sigma_eta={np.sqrt(seta2_from_arima):.4f}")

# ---------- 5. 稳态增益 K 与 1-theta ----------
m_pred, P_pred, v, V, K = kalman_local_level(y, s_eta2_hat, s_e2_hat)
q = s_eta2_hat / s_e2_hat                        # 信噪比
K_ss = (-q + np.sqrt(q ** 2 + 4 * q)) / 2         # 标量 Riccati 方程的稳态增益
print(f"K_1..K_5 = {np.round(K[:5], 3)},  K_T={K[-1]:.4f}, 解析稳态 K={K_ss:.4f}, 1-theta={1 - theta:.4f}")

# ---------- 6. 后向平滑 (11.20)(11.23) ----------
L = 1 - K
qv, M = 0.0, 0.0
mu_s, P_s = np.empty(T), np.empty(T)
for t in range(T - 1, -1, -1):
    qv = v[t] / V[t] + L[t] * qv
    M = 1 / V[t] + L[t] ** 2 * M
    mu_s[t] = m_pred[t] + P_pred[t] * qv
    P_s[t] = P_pred[t] - P_pred[t] ** 2 * M
mu_f = m_pred + P_pred / V * v                    # 滤波值 mu_{t|t}
P_f = P_pred * (1 - K)
rmse = lambda a: np.sqrt(np.mean((a[5:] - mu[5:]) ** 2))
print(f"相对真实状态的 RMSE: 原始观测 {rmse(y):.3f}, 滤波 {rmse(mu_f):.3f}, 平滑 {rmse(mu_s):.3f}")
print(f"平均标准差: 滤波 {np.sqrt(P_f[5:]).mean():.3f}, 平滑 {np.sqrt(P_s[5:]).mean():.3f}")

# ---------- 7. 标准化新息诊断 ----------
from statsmodels.stats.diagnostic import acorr_ljungbox
z = (v / np.sqrt(V))[1:]
lb = acorr_ljungbox(z, lags=[12]); lb2 = acorr_ljungbox(z ** 2, lags=[12])
print(f"标准化新息 Q(12)={lb.lb_stat.iloc[0]:.2f} (p={lb.lb_pvalue.iloc[0]:.2f}), "
      f"平方 Q(12)={lb2.lb_stat.iloc[0]:.2f} (p={lb2.lb_pvalue.iloc[0]:.2f})")

# ---------- 8. 缺失值：挖掉第 200-209 天 ----------
y_miss = y.copy(); y_miss[200:210] = np.nan
mp, Pp, *_ = kalman_local_level(y_miss, s_eta2_hat, s_e2_hat)
print("缺失段预测方差 Sigma_{t|t-1}:", np.round(Pp[199:212], 4))

关键输出：

手写 Kalman MLE : sigma_eta=0.0700, sigma_e=0.4641, logL=-246.98
statsmodels UC  : sigma_eta=0.0700, sigma_e=0.4641
ARIMA(0,1,1)    : theta=0.860, sigma_a=0.5004 -> sigma_e=0.4641, sigma_eta=0.0700
K_1..K_5 = [1.    0.506 0.346 0.269 0.226],  K_T=0.1399, 解析稳态 K=0.1399, 1-theta=0.1399
相对真实状态的 RMSE: 原始观测 0.465, 滤波 0.175, 平滑 0.144
平均标准差: 滤波 0.174, 平滑 0.128
标准化新息 Q(12)=11.87 (p=0.46), 平方 Q(12)=7.11 (p=0.85)
缺失段预测方差 Sigma_{t|t-1}: [0.035  0.035  0.0399 0.0448 0.0497 0.0546 0.0595 0.0644 0.0693 0.0742
 0.0791 0.084  0.0653]

读输出：

• 手写滤波、statsmodels 和 ARIMA 换算三条路径给出完全相同的 \((\hat\sigma_\eta,\hat\sigma_e)\)，验证了 11.2 节的等价性。模拟数据的 \(\hat\theta=0.860\)，与例 11.1 的 0.858 很接近。
• 增益从扩散初始化的 \(K_1=1\) 迅速收敛到稳态 0.1399，与解析公式和 \(1-\hat\theta\) 一致：这条序列的最优 EWMA 平滑常数约为 0.14。
• 原始观测离真实状态的 RMSE 约 0.47（就是噪声 \(\sigma_e\)）；滤波把它降到 0.175，平滑进一步降到 0.144，标准差也相应缩小。回测时只能用滤波值那一档的精度。
• 标准化新息的 Ljung–Box 检验不显著，模型充分。
• 缺失段内 \(\Sigma_{t|t-1}\) 每天增加 \(\hat\sigma_\eta^2\approx0.0049\)，数据恢复后的第一笔观测立即把方差拉回来（最后一个数 0.0653）。

本章小结

局部水平模型 \(y_t=\mu_t+e_t\)、\(\mu_{t+1}=\mu_t+\eta_t\) 是最简单的状态空间模型，它与 \(\theta>0\) 的 ARIMA(0,1,1) 一一对应，稳态 Kalman 滤波就是平滑常数 \(K=1-\theta\) 的指数平滑。Kalman 滤波来自多元正态的条件分布公式：因为新息 \(v_t\) 与过去信息独立，条件均值和方差可以增量更新，增益 \(K_t=\Sigma_{t|t-1}/V_t\) 是状态对新息的回归系数。新息序列相互独立，给出观测协方差阵的 Cholesky 分解，因而精确似然可写成 \(-\frac12\sum(\ln V_t+v_t^2/V_t)\)，跑一遍滤波就能算出。平滑用后向递推 \(q_{t-1},M_{t-1}\)，结果更平滑、区间更窄，但用到了未来数据，不能用于回测信号。缺失值只需令 \(v_t=0,K_t=0\)；初值未知时用扩散初始化。

练习

基础

1. 由 (11.4)(11.5) 推出 \(\theta\) 满足 \(\theta^2-(2+q)\theta+1=0\)，并证明两根互为倒数、可逆根在 \((0,1)\) 内。 提示：两式相除，\((1+\theta^2)/\theta=2+q\)。
2. 用例 11.1 的 \(\hat\theta=0.858\)、\(\hat\sigma_a=0.5184\) 算出 \(\sigma_e\)、\(\sigma_\eta\) 和稳态增益 \(K\)，并解释 \(K\) 作为 EWMA 平滑常数的含义（半衰期约多少天？）。 提示：\(\sigma_e=0.480\)，\(\sigma_\eta=0.0736\)，\(K=0.142\)；半衰期 \(\ln0.5/\ln(1-K)\approx4.5\) 天。
3. 证明 (11.9) 中 \(E[\mu_{t|t-1}(\mu_t-\mu_{t|t-1})|F_{t-1}]=0\)，并说明为什么需要 \(e_t\) 与 \(\mu_t\) 独立。
4. 写出局部水平模型在 \(t=1,2\) 的 \(v_t\)，验证 (11.15) 中 \(K\) 的第 3 行元素。
5. 某数据在第 \(\ell+1\) 到 \(\ell+5\) 期缺失，已知 \(\Sigma_{\ell+1|\ell}=0.03\)，\(\sigma_\eta^2=0.005\)。求 \(\Sigma_{\ell+6|\ell+5}\)，并说明第 \(\ell+6\) 期增益为何比稳态大。 提示：\(0.03+5\times0.005=0.055\)。

进阶

6. 证明 (11.16)，并由它推出 \(\mathrm{Cov}(x_t,v_{t+2})=\Sigma_{t|t-1}L_tL_{t+1}\)。
7. 推导扩散初始化下 \(\mu_{1|T}\to y_1+\sigma_e^2q_1\)，\(\Sigma_{1|T}\to\sigma_e^2-\sigma_e^4M_1\)。
8. 修改本章代码：把 \(\sigma_e\) 改为 0.1（噪声很小），重新估计并观察稳态增益和滤波 RMSE 的变化；再把 \(\sigma_\eta\) 改为 0（状态恒定），解释此时 Kalman 滤波退化为什么估计量（样本均值的递推形式）。
9. 用本章代码的极大似然改为滚动 250 日窗口估计 \((\sigma_\eta,\sigma_e)\)，每个窗口末端只取滤波值，构造一个无前视的波动率预测序列，并与全样本参数的版本比较差异。
10. （原书习题 11.2）用 20 分钟收益构造的 Alcoa 已实现波动率：拟合 ARIMA(0,1,1)，估计局部趋势模型，画滤波与平滑状态及 95% 区间。没有原始数据时，可以用第 05 章方法模拟日内价格（含买卖价反弹）来构造。

原书推荐习题：11.2（已实现波动率的局部趋势模型，理解微观结构噪声与 ARIMA 等价性）。其余习题见第 11b 章。

原书对照

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