状态空间模型:表示与卡尔曼滤波

如果这篇文章只记一件事，记这个：一大批时间序列模型——ARMA、ARIMA、结构时间序列、（线性化后的）随机波动率、时变系数回归、动态因子模型——都可以写成状态空间形式。状态空间不是又一个模型，是一个框架：只要你的模型能装进”观测方程加状态方程”这两行，卡尔曼滤波就用同一套递推，一趟扫过数据，同时给你隐状态的最优估计和精确的似然函数。学一个算法，通吃一族模型——这笔买卖是本文要兑现的承诺。

推导按这个顺序走：先给表示和滤波递推，看似然怎么从滤波里免费掉出来；再讲平滑和缺失数据这两个框架自带的赠品；然后逐个把具体模型装进框架里；最后处理工程细节和非高斯世界的逃生通道。参考 Hamilton (1994, 第 13 章)、Durbin & Koopman (2012)、Harvey (1989)；金融应用见 Tsay (2010, 第 11 章)。

基本状态空间模型

线性高斯状态空间模型就两个方程。观测方程，你测到的东西：

状态方程，藏在底下的动态：

$\mathbf{y}_t$ 是 $K \times 1$ 的观测向量；状态 $\boldsymbol{\alpha}_t$ 是 $m \times 1$ 向量，不可直接观测。矩阵 $\mathbf{Z}_t, \mathbf{T}_t, \mathbf{H}_t, \mathbf{R}_t, \mathbf{Q}_t$ 和向量 $\mathbf{d}_t, \mathbf{c}_t$ 是参数的已知函数，可以时变。初始状态 $\boldsymbol{\alpha}_1 \sim N(\mathbf{a}_1, \mathbf{P}_1)$，$\boldsymbol{\varepsilon}_t$ 和 $\boldsymbol{\eta}_t$ 在时间上独立、彼此也独立。

两行方程分工明确：状态方程说世界怎么演化，观测方程说你隔着多厚的噪声看世界。这个表示的适用范围比看上去宽得多——后面会把一个个熟悉的模型塞进去给你看。但先把回报讲清楚：滤波。

卡尔曼滤波

给定状态空间模型，卡尔曼滤波递推地计算滤波状态估计 $\mathbf{a}_t = E[\boldsymbol{\alpha}_t \mid \mathbf{y}_1, \ldots, \mathbf{y}_t]$ 和它的误差协方差 $\mathbf{P}_t = E[(\boldsymbol{\alpha}_t - \mathbf{a}_t)(\boldsymbol{\alpha}_t - \mathbf{a}_t)' \mid \mathbf{y}_1, \ldots, \mathbf{y}_t]$。每个 $t$ 走两步。

预测步，用状态方程把 $t-1$ 的知识推到 $t$：

更新步，把新观测 $\mathbf{y}_t$ 吸收进来。先算预测误差（创新）：

和创新协方差：

然后是整个算法的灵魂，卡尔曼增益：

它衡量”这次观测意外了多少”应该换算成”状态估计修正多少”——状态不确定性大、观测噪声小，增益就大，新观测说了算；反过来增益就小，先验说了算。更新：

预测、更新，交替到 $T$，一趟结束。

现在看承诺里的另一半。创新 $\mathbf{v}_t$ 是模型没料到的部分，在正确设定下 $\mathbf{v}_t \sim N(\mathbf{0}, \mathbf{F}_t)$ 且相互独立——这叫预测误差分解。它直接给出数据的对数似然：

$\mathbf{v}_t$ 和 $\mathbf{F}_t$ 全是滤波的副产品。也就是说，跑一遍滤波，状态估计和精确似然同时到手；把 $\ell(\boldsymbol{\theta})$ 塞给数值优化器，就是极大似然估计。这就是框架的全部魔法：表示统一了，算法就统一了，估计也就统一了。

有一个不起眼但躲不掉的问题：递推从哪起步？状态平稳时用平稳分布——$\mathbf{a}_0$ 取无条件均值，$\mathbf{P}_0$ 取无条件协方差。状态里有非平稳成分（随机游走趋势、积分序列）时，无条件分布不存在，用 diffuse 初始化：$\mathbf{a}_0 = \mathbf{0}$、$\mathbf{P}_0 = \kappa \mathbf{I}$、$\kappa \to \infty$，表示”起点一无所知”。天真地取一个大 $\kappa$ 数值上很危险，前几期的预测误差分解会失真；正确做法是精确 diffuse 卡尔曼滤波（Koopman, 1997; Durbin & Koopman, 2012, 第 5 章），它把初始瞬态解析地处理掉，多数正经软件都实现了。

平滑:让未来说话

滤波估计 $\mathbf{a}_t$ 只用截至 $t$ 的数据，适合实时场景。事后分析时你有整个样本，该问的是 $\hat{\boldsymbol{\alpha}}_t = E[\boldsymbol{\alpha}_t \mid \mathbf{y}_1, \ldots, \mathbf{y}_T]$——平滑估计，让未来的观测也对过去的状态发言。Rauch-Tung-Striebel (RTS) 平滑器在滤波之后从 $T$ 倒着跑回 $1$：

平滑协方差一般小于滤波协方差——未来数据确实带来了关于 $\boldsymbol{\alpha}_t$ 的信息。用事后视角把趋势从噪声里剥出来、为经济分析重建未观测成分的完整路径，都是平滑的活。但用的时候记住它的身份：平滑用了未来数据，任何实时决策或样本外评估只能用滤波估计，这条纪律在所有隐状态模型里通用。

缺失数据:白送的功能

框架还有一个几乎白送的好处。某期 $\mathbf{y}_t$ 缺失？跳过那期的更新步就行：

滤波器靠状态方程继续传播，不确定性相应变大，下一个观测到来时自动收窄。不需要插值、不需要删样本，似然照常算。这是状态空间模型在不规则采样数据（临床测量、卫星观测、停牌的资产价格）上大行其道的原因。

把模型装进框架

框架的价值要靠具体模型来兑现。挨个装几个。

ARMA。以 ARMA(1,1) 为例：$y_t = \phi y_{t-1} + \theta \varepsilon_{t-1} + \varepsilon_t$。定义状态 $\boldsymbol{\alpha}_t = (y_t, \theta \varepsilon_t)'$，则

第一行展开验证：$y_{t+1} = \phi y_t + \theta \varepsilon_t + \varepsilon_{t+1}$，正是 ARMA(1,1)。对应到标准记号就是 $\mathbf{T} = \begin{pmatrix} \phi & 1 \\ 0 & 0 \end{pmatrix}$、$\mathbf{Z} = (1, 0)$、$\mathbf{R} = (1, \theta)'$、$\mathbf{H} = 0$（观测无噪声）、$\mathbf{Q} = \sigma^2$。一般的 ARMA($p,q$) 用同样的伴随形式套路，状态维度 $\max(p, q+1)$。这个表示不只是智力练习：主流软件对 ARMA 做精确 MLE，内部走的就是卡尔曼滤波——它把”条件在初值上”的近似似然换成了精确似然，还顺手解决了 MA 部分不可观测的问题。

结构时间序列。Harvey (1989) 的路线是不从 ARMA 出发，直接按可解释成分搭模型：

趋势 $\mu_t$、周期或季节成分 $\gamma_t$、不规则项 $\varepsilon_t$，每个成分给一套状态动态。趋势从简到繁有两档：局部水平（随机游走）$\mu_t = \mu_{t-1} + \eta_t$；局部线性趋势——水平带漂移、漂移自身也是随机游走：

周期成分用阻尼三角过程，频率 $\lambda_c$、阻尼 $\rho$，需要一对状态：

季节成分同样有两种写法：虚拟变量式（$s-1$ 个状态，约束 $\sum_{j=1}^{s-1} \gamma_{j,t} = 0$，各自走随机游走）或三角式（$\gamma_t = \sum_{j=1}^{[s/2]} a_j \cos(\lambda_j t) + b_j \sin(\lambda_j t)$，系数随机演化）。状态向量就是全部未观测成分摞起来——局部线性趋势加虚拟季节是 $(\mu_t, \nu_t, \gamma_{1,t}, \ldots, \gamma_{s-1,t})'$，维度 $2 + (s-1)$。状态空间矩阵从结构方程机械拼装，最小的例子是局部水平模型：

装好之后一切照旧：滤波给似然，平滑把趋势、季节、周期各自的路径提取出来，预测递推着往前推。

时变系数回归。$\mathbf{y}_t = \mathbf{X}_t \boldsymbol{\beta}_t + \varepsilon_t$，让系数走随机游走 $\boldsymbol{\beta}_{t+1} = \boldsymbol{\beta}_t + \boldsymbol{\eta}_t$——状态就是系数向量，$\mathbf{Z}_t = \mathbf{X}_t$ 时变。卡尔曼滤波输出的正是系数的时变估计，比滚动窗口回归干净得多。

动态因子模型。$K$ 条观测序列由 $r$ 个隐因子驱动（$r \ll K$），状态是因子向量（可带 VAR 动态），观测方程用 $K \times r$ 载荷矩阵 $\boldsymbol{\Lambda}$ 把因子映到序列上。估计用 EM 或两步法（PCA 初始化加卡尔曼精修）。

随机波动率。对数方差 $\log \sigma_t^2$ 作为隐状态走 AR(1)。观测方程严格说是非线性非高斯的，取对数平方收益后可以用卡尔曼做 QML 近似，认真做用粒子滤波——这是往下一节的引子。顺带一提本框架的出身：卡尔曼 (1960) 的原始应用是导航追踪，状态是位置和速度，观测是带噪的测量，飞机和卫星先用上，计量经济学是后来的移民。

有两点边界要说清。其一，GARCH 不在这份名单里：它的条件方差给定过去是确定的（所谓 observation-driven），不含真正的隐状态噪声，不是天然的状态空间模型——带隐藏噪声的对应物是上面的随机波动率。其二，“能装进框架”不等于”必须用框架的思维”：对普通 ARMA，软件已经把卡尔曼封装在 arima 接口后面，你不需要手写矩阵；框架思维真正的用武之地，是你要自己搭一个带隐成分的模型的时候。

估计的工程面

原理一句话——优化器套滤波——但工程细节决定成败。参数经常带约束（方差为正、相关系数在区间内），要么用有界优化器（多数状态空间软件用 L-BFGS-B），要么变换参数（方差取 log、相关取 logit）到无约束空间。初始条件按前述规则：平稳成分用无条件矩，非平稳成分用 diffuse，混合情形交给精确 diffuse 滤波统一处理。收敛判据看对数似然和参数的变化量双双降到容差以下（比如 $10^{-6}$）。标准误来自观测信息矩阵的逆：

数值上用有限差分或自动微分算 Hessian；走 EM 路线时可以用 Louis (1982) 恒等式直接得到信息矩阵，更省。

说到 EM：它是直接极大化的替代方案，对状态空间模型经常更稳定，因为两步都有闭式。E 步跑卡尔曼平滑器，算出平滑状态 $E[\boldsymbol{\alpha}_t \mid \mathbf{y}_{1:T}, \theta^{\text{old}}]$ 和平滑叉积 $E[\boldsymbol{\alpha}_t \boldsymbol{\alpha}_{t-1}' \mid \mathbf{y}_{1:T}, \theta^{\text{old}}]$；M 步拿这些充分统计量更新参数，比如状态转移矩阵和状态噪声协方差：

观测方程参数类似。EM 单调爬升但只保证局部极大，多组随机起点是必要的保险。

估计完还要对账，而对账的工具已经在手里：预测误差分解说正确设定下标准化创新 $\mathbf{v}_t / \sqrt{\mathbf{F}_t}$ 应该近似 iid 标准正态。画它的 ACF、跑 Ljung-Box 查独立性，Jarque-Bera（多元用 Mahalanobis 距离）查正态性，ARCH-LM 查异方差。任何一项挂了都说明模型误设，常见的修法是加成分（趋势、季节、外生变量）、放开参数让它时变、或者换非高斯观测分布——最后这条通向下一节。

离开线性高斯世界

卡尔曼滤波只对线性高斯模型精确。现实经常不配合，逃生通道按代价从低到高排。最便宜的是准极大似然（QML）：条件均值和方差设定正确时，照用高斯似然也能得到一致估计，配三明治标准误即可——高阶矩（偏度、峰度）错了不伤一致性，多数金融应用到这一步就够了。模型非线性但不太狂野时用扩展卡尔曼滤波：在当前状态估计附近线性化，对线性化系统跑标准滤波，近似但常常好使。非线性更强时无迹卡尔曼滤波更稳：选一小组 sigma 点穿过非线性变换去传播状态分布，不用求导。都不行时上粒子滤波：序贯蒙特卡洛，一群粒子按状态方程演化、按观测似然重加权，粒子数够多时一致，代价是计算量。选择的原则和这篇文章的主线一致：能用便宜工具解决的，别上贵的。

代码:AR(1) 加噪声

把最小的例子完整写一遍——隐状态是 AR(1)，观测带噪声——滤波和似然一共三十行：

import numpy as np
from scipy.optimize import minimize

def kalman_loglik(params, y):
    phi, sigma_state, sigma_obs = params
    T = len(y)
    
    # 初始化
    a = 0.0
    P = sigma_state**2 / (1 - phi**2) if abs(phi) < 1 else 100.0
    
    loglik = 0
    for t in range(T):
        # 预测
        a_pred = phi * a
        P_pred = phi**2 * P + sigma_state**2
        
        # 创新
        v = y[t] - a_pred
        F = P_pred + sigma_obs**2
        
        loglik += -0.5 * (np.log(2*np.pi*F) + v**2 / F)
        
        # 更新
        K = P_pred / F
        a = a_pred + K * v
        P = (1 - K) * P_pred
    
    return -loglik


result = minimize(kalman_loglik, x0=[0.9, 1.0, 1.0], args=(y,),
                  method='Nelder-Mead', options={'xatol': 1e-6})
phi_hat, sig_state, sig_obs = result.x

生产环境别自己写——statsmodels、R 的 KFAS、MATLAB 的 ssm 都有成熟实现，含 diffuse 初始化和数值稳定的协方差更新。手写一遍的价值在于把递推焊进脑子。

什么时候用它,代价是什么

回到开头的承诺做个收尾。状态空间是对的选择，当数据生成过程真有隐状态（随机波动率、动态因子、时变系数），当多条序列共享同一个隐动态，当数据有缺失或采样不规则，或者当你要的输出本身就是未观测成分的滤波或平滑路径（趋势、周期、时变系数）。反过来，没有隐成分、标准接口就能跑的模型不需要你亲自搭框架，高度非线性的动态则可能超出扩展和无迹滤波的舒适区、更适合模拟类方法。

计算代价心里要有数：滤波每步的矩阵运算由 $m \times m$ 乘法主导，一趟滤波约 $O(T m^3)$（状态矩阵有稀疏或块结构时可以压下来），平滑同量级，MLE 还要乘上优化器的迭代次数。小状态维度的日常应用快得无感；状态维度上到成百上千（大型动态因子模型、集成方法的场景）时，计算就成了要认真设计的事——稀疏结构、并行的粒子或集成方法、sigma 点近似都是那个世界的工具。对本文的读者，更常见的情形是 $m$ 个位数、$T$ 几千：框架的全部好处，几乎不花计算的钱。

参考文献

• Durbin, J., & Koopman, S. J. (2012). Time Series Analysis by State Space Methods (2nd ed.). Oxford University Press.
• Hamilton, J. D. (1994). Time Series Analysis. Princeton University Press. Ch. 13.
• Harvey, A. C. (1989). Forecasting, Structural Time Series Models and the Kalman Filter. Cambridge University Press.
• Kalman, R. E. (1960). A new approach to linear filtering and prediction problems. Journal of Basic Engineering, 82(1), 35-45.
• Koopman, S. J. (1997). Exact initial Kalman filter and smoother for non-stationary time series models. Journal of the American Statistical Association, 92(440), 1630-1638.
• Louis, T. A. (1982). Finding the observed information matrix when using the EM algorithm. Journal of the Royal Statistical Society B, 44(2), 226-233.
• Rauch, H. E., Tung, F., & Striebel, C. T. (1965). Maximum likelihood estimates of linear dynamic systems. AIAA Journal, 3(8), 1445-1450.
• Tsay, R. S. (2010). Analysis of Financial Time Series (3rd ed.). Wiley. Ch. 11.
• Shumway, R. H., & Stoffer, D. S. (2017). Time Series Analysis and Its Applications (4th ed.). Springer.