30 线性高斯状态空间模型

参考：

30.1 卡尔曼滤波算法

30.1.1 模型和问题

考虑如下的线性高斯状态空间模型：

\[\begin{align} \boldsymbol y_t =& Z_t \boldsymbol \alpha_t + \boldsymbol \varepsilon_t, \ \boldsymbol \varepsilon_t \sim \text{N}(0, H_t), \tag{30.1} \\ \boldsymbol \alpha_{t+1} =& T_t \boldsymbol \alpha_t + R_t \boldsymbol \eta_t, \ \boldsymbol \eta_t \sim \text{N}(0, Q_t), \tag{30.2} \\ \boldsymbol \alpha_1 \sim& \text{N}(\boldsymbol a_1, P_1) . \end{align}\]

其中\(\boldsymbol y_t\)是\(t\)时刻的观测值， 为\(p \times 1\)向量，设有\(t=1,2,\dots,n\)时刻的观测值； \(\boldsymbol \alpha_t\)是\(t\)时刻系统的状态， 是不可观测的\(m\times 1\)随机向量。 (30.1)称为观测方程， (30.2)称为状态方程。 \(\{\boldsymbol \varepsilon_t\}\)和\(\{\boldsymbol \eta_t \}\)相互独立， 都是独立同分布向量白噪声列， \(\boldsymbol \varepsilon_t\)为\(p\times 1\)随机向量， \(\boldsymbol \eta_t\)为\(r \times 1\)随机向量，\(r \leq m\)。 设各矩阵\(Z_t, T_t, R_t, H_t, Q_t\)已知， \(Z_t\)和\(T_{t-1}\)允许依赖于\(\boldsymbol y_1, \dots, \boldsymbol y_{t-1}\)， 初始状态\(\boldsymbol\alpha_1\)服从\(N(\boldsymbol a_1, P_1)\)， 设\(\boldsymbol a_1, P_1\)已知， \(\boldsymbol\alpha_1\)与\(\{\boldsymbol \varepsilon_t\}\)和\(\{\boldsymbol \eta_t \}\)独立。 当参数未知时， 设\(\boldsymbol\psi\)为未知参数， 矩阵\(Z_t, T_t, R_t, H_t, Q_t\)可以依赖于未知参数\(\boldsymbol\psi\)。

各个矩阵的维数如下：

向量 矩阵
\(\boldsymbol y_t\) \(p \times 1\) \(Z_t\) \(p \times m\)
\(\boldsymbol \alpha_t\) \(m \times 1\) \(T_t\) \(m \times m\)
\(\boldsymbol \varepsilon_t\) \(p \times 1\) \(H_t\) \(p \times p\)
\(\boldsymbol \eta_t\) \(r \times 1\) \(R_t\) \(m \times r\)
\(Q_t\) \(m \times r\)
\(\boldsymbol a_1\) \(p \times 1\) \(P_1\) \(m \times m\)

记 \[ \boldsymbol Y_t = \begin{pmatrix} \boldsymbol y_1 \\ \boldsymbol y_2 \\ \vdots \\ \boldsymbol y_t \end{pmatrix} . \]

设\(\boldsymbol a_1, P_1\)已知， \(\boldsymbol \alpha_1 \sim \text{N}(\boldsymbol a_1, P_1)\)。

归纳地， 设\(\boldsymbol \alpha_t | \boldsymbol Y_{t-1} \sim \text{N}(\boldsymbol a_t, P_t)\)。 在获得\(\boldsymbol y_t\)后， 递推计算\(\boldsymbol \alpha_t | \boldsymbol Y_{t}\)的条件分布， 和\(\boldsymbol \alpha_{t+1} | \boldsymbol Y_{t}\)的条件分布。 由定理28.1可知这两个条件分布正态。 记\(\boldsymbol a_{t|t} = E(\boldsymbol \alpha_t | \boldsymbol Y_t)\), \(P_{t|t} = \text{Var}(\boldsymbol \alpha_t | \boldsymbol Y_t)\), \(\boldsymbol a_{t+1} = E(\boldsymbol \alpha_{t+1} | \boldsymbol Y_t)\), \(P_{t+1} = \text{Var}(\boldsymbol \alpha_{t+1} | \boldsymbol Y_t)\)。

30.1.2 更新和预报步

记 \[\begin{align} \boldsymbol v_t =& \boldsymbol y_t - E(\boldsymbol y_t | \boldsymbol Y_{t-1}) = \boldsymbol y_t - E(Z_t \boldsymbol \alpha_t + \boldsymbol\varepsilon_t | \boldsymbol Y_{t-1}) = \boldsymbol y_t - Z_t \boldsymbol a_t . \tag{30.3} \end{align}\] 由定理28.1可知\(\boldsymbol v_t\)与\(\boldsymbol y_1, \dots, \boldsymbol y_{t-1}\)独立，\(E \boldsymbol v_t = \boldsymbol 0\)。 \(\sigma(\boldsymbol Y_{t-1}, \boldsymbol v_t) = \sigma(\boldsymbol Y_t)\)， 所以给定\(\boldsymbol Y_t\)的条件分布， 等于给定\(\boldsymbol Y_{t-1}\)和\(\boldsymbol v_t\)的条件分布。

令\(F_t = \text{Var}(\boldsymbol v_t)\)， 因\(\boldsymbol v_t\)与\(\boldsymbol Y_{t-1}\)独立， 所以有 \[\begin{align} F_t =& \text{Var}(\boldsymbol v_t | \boldsymbol Y_{t-1}) = \text{Var}(Z_t \boldsymbol\alpha_t + \boldsymbol\varepsilon_t - Z_t \boldsymbol a_t | \boldsymbol Y_{t-1}) \\ =& \text{Var}(Z_t (\boldsymbol\alpha_t - \boldsymbol a_t) + \boldsymbol\varepsilon_t | \boldsymbol Y_{t-1}) \\ =& Z_t P_t Z_t^T + H_t . \tag{30.4} \end{align}\]

由\(\sigma(\boldsymbol Y_{t-1}, \boldsymbol v_t) = \sigma(\boldsymbol Y_t)\)可得 \[\begin{aligned} \boldsymbol a_{t|t} =& E(\boldsymbol\alpha_t | \boldsymbol Y_{t-1}, \boldsymbol v_t) \\ =& E(\boldsymbol\alpha_t | \boldsymbol Y_{t-1}) + \text{Cov}(\boldsymbol\alpha_t, \boldsymbol v_t) [\text{Var}(\boldsymbol v_t)]^{-1} \boldsymbol v_t \end{aligned}\] 其中 \[\begin{align} \text{Cov}(\boldsymbol\alpha_t, \boldsymbol v_t) =& E \left\{ E\left[ \boldsymbol\alpha_t (Z_t \boldsymbol\alpha_t + \boldsymbol\varepsilon_t - Z_t \boldsymbol a_t)^T | \boldsymbol Y_{t-1} \right] \right\} \\ =& E \left\{ E\left[ \boldsymbol\alpha_t (\boldsymbol\alpha_t - \boldsymbol a_t)^T Z_t^T | \boldsymbol Y_{t-1} \right] \right\} \\ =& P_t Z_t^T . \tag{30.5} \end{align}\] 于是得 \[\begin{align} \boldsymbol a_{t|t} = \boldsymbol a_t + P_t Z_t^T F_t^{-1} \boldsymbol v_t . \tag{30.6} \end{align}\]

仍利用定理28.1可知 \[\begin{align} P_{t|t} =& \text{Var}(\boldsymbol\alpha_t | \boldsymbol Y_{t-1}, \boldsymbol v_t) \\ =& \text{Var}(\boldsymbol\alpha_t | \boldsymbol Y_{t-1}) - \text{Cov}(\boldsymbol\alpha_t, \boldsymbol v_t) [\text{Var}(\boldsymbol v_t)]^{-1} \text{Cov}(\boldsymbol\alpha_t, \boldsymbol v_t)^T \\ =& P_t - P_t Z_t^T F_t^{-1} Z_t P_t . \tag{30.7} \end{align}\]

以上设\(\boldsymbol v_t\)的方差阵\(F_t\)可逆， 对于设计正确的模型这一般是能满足的。 称(30.6)和(30.7)为卡尔曼滤波的“更新步”。

下面递推计算\(\boldsymbol a_{t+1}\)和\(P_{t+1}\)。 \[\begin{aligned} \boldsymbol a_{t+1} =& T_t \boldsymbol a_{t|t} \\ =& T_t \boldsymbol a_t + T_t P_t Z_t^T F_t^{-1} \boldsymbol v_t . \\ P_{t+1} =& \text{Var}(\boldsymbol\alpha_{t+1} | \boldsymbol Y_t) \\ =& \text{Var}(T_t \boldsymbol\alpha_t + R_t \boldsymbol\eta_t | \boldsymbol Y_t) \\ =& T_t P_{t|t} T_t^T + R_t Q_t R_t^T . \end{aligned}\] 记 \[\begin{align} K_t =& T_t P_t Z_t^T F_t^{-1} . \tag{30.8} \end{align}\] 则 \[\begin{align} \boldsymbol a_{t+1} =& T_t \boldsymbol a_t + K_t \boldsymbol v_t . \tag{30.9} \end{align}\] 而 \[\begin{align} P_{t+1} =& T_t(P_t - P_t Z_t^T F_t^{-1} Z_t P_t)T_t^T + R_t Q_t R_t^T \\ =& T_t P_t T_t^T - T_t P_t Z_t^T F_t^{-1} Z_t P_t T_t^T + R_t Q_t R_t^T \\ =& T_t P_t T_t^T - T_t P_t Z_t^T K_t^T + R_t Q_t R_t^T \\ =& T_t P_t (T_t - K_t Z_t)^T + R_t Q_t R_t^T . \tag{30.10} \end{align}\] 称(30.9)和(30.10)为卡尔曼滤波的“预测步”。

递推公式(30.6),(30.7),(30.9)和(30.10)就构成了卡尔曼滤波的递推算法。 在已知\(\boldsymbol Y_{t-1}\)的基础上， 新增了观测\(\boldsymbol y_t\)后， 可以更新关于\(\boldsymbol\alpha_t\)和\(\boldsymbol\alpha_{t+1}\)的条件分布。

如果直接计算模型(30.1)和(30.2)的似然函数， 因为\(\boldsymbol Y_n\)是\(pn\)维正态分布， 所以需要求一个\(pn \times pn\)阶矩阵的逆矩阵， 计算量是\(O(p^3 n^3)\)阶的。 如果使用卡尔曼滤波， 则只需要计算\(n\)步， 每一步仅需要对\(p\times p\)矩阵\(F_t\)求逆， 计算量是\(O(p^3 n)\)阶的。

即使\(\boldsymbol y_t\)是\(p \times 1\)向量， 在实际计算时， 还是将其转化为\(p\)个一维观测值更有利于提高计算效率。

30.1.3 卡尔曼滤波公式

设\(\boldsymbol\alpha_1 \sim \text{N}(\boldsymbol a_1, P_1)\)， \(\boldsymbol a_1\), \(P_1\)已知， 观测值为\(\boldsymbol y_1, \dots, \boldsymbol y_n\)。 对\(t=1, 2, \dots, n\)，计算 \[\begin{equation} \begin{aligned} \boldsymbol v_t =& \boldsymbol y_t - Z_t \boldsymbol a_t, & F_t =& Z_t P_t Z_t^T + H_t, \\ \boldsymbol a_{t|t} =& \boldsymbol a_t + P_t Z_t^T F_t^{-1} \boldsymbol v_t, & P_{t|t} =& P_t - P_t Z_t^T F_t^{-1} Z_t P_t, \\ K_t =& T_t P_t Z_t^T F_t^{-1}, & L_t =& T_t - K_t Z_t, \\ \boldsymbol a_{t+1} =& T_t \boldsymbol a_t + K_t \boldsymbol v_t, & P_{t+1} =& T_t P_t L_t^T + R_t Q_t R_t^T . \end{aligned} \tag{30.11} \end{equation}\]

(30.11)称为卡尔曼滤波。 如果不需要\(\boldsymbol\alpha_t | \boldsymbol Y_t\)的分布参数也可以略过\(\boldsymbol a_{t|t}\)和\(P_{t|t}\)的计算； 如果计算了\(\boldsymbol a_{t|t}\)和\(P_{t|t}\)， 则可以如下计算\(\boldsymbol a_{t+1}\)和\(P_{t+1}\)： \[ \boldsymbol a_{t+1} = T_t \boldsymbol a_{t|t}, \quad P_{t+1} = T_t P_{t|t} T_t^T + R_t Q_t R_t^T . \]

30.1.4 带有偏移项的模型

可以在观测方程和状态方程中分别加入一个均值偏移项： \[\begin{aligned} \boldsymbol y_t =& Z_t \boldsymbol \alpha_t + \boldsymbol d_t + \boldsymbol \varepsilon_t, \ \boldsymbol \varepsilon_t \sim \text{N}(0, H_t), \\ \boldsymbol \alpha_{t+1} =& T_t \boldsymbol \alpha_t + \boldsymbol c_t + R_t \boldsymbol \eta_t, \ \boldsymbol \eta_t \sim \text{N}(0, Q_t), \\ \boldsymbol \alpha_1 \sim& \text{N}(\boldsymbol a_1, P_1) . \end{aligned}\]

其中\(p \times 1\)向量\(\boldsymbol d_t\)已知， \(m \times 1\)向量\(\boldsymbol c_t\)已知。 这时， 卡尔曼滤波公式变成： \[\begin{aligned} \boldsymbol v_t =& \boldsymbol y_t - Z_t \boldsymbol a_t - \boldsymbol d_t, & F_t =& Z_t P_t Z_t^T + H_t, \\ \boldsymbol a_{t|t} =& \boldsymbol a_t + P_t Z_t^T F_t^{-1} \boldsymbol v_t, & P_{t|t} =& P_t - P_t Z_t^T F_t^{-1} Z_t P_t, \\ && K_t =& T_t P_t Z_t^T F_t^{-1}, \\ \boldsymbol a_{t+1} =& T_t \boldsymbol a_t + K_t \boldsymbol v_t + \boldsymbol c_t, & P_{t+1} =& T_t P_t (T_t - K_t Z_t)^T + R_t Q_t R_t^T . \end{aligned}\] 即仅有\(\boldsymbol v_t\)和\(\boldsymbol a_{t+1}\)公式改变了。

30.1.5 稳定状态

若模型(30.1)、(30.2)中的均值\(Z_t, H_t, T_t, R_t, Q_t\)都是不随时间变化的， 则\(n \to \infty\)时\(P_t\), \(K_t\), \(F_t\)极限存在。 记极限为\(\bar P, \bar K, \bar F\)， 应满足方程 \[\begin{aligned} \bar P =& T \bar P T^T - T \bar P Z^T \bar F^{-1} Z \bar P T^T + R Q R^T , \\ \bar K =& T \bar P Z^T \bar F^{-1} , \\ \bar L =& T - \bar K Z, \\ \bar F =& Z \bar P Z^T + H . \end{aligned}\] 当\(t\)充分大时， 可以恒定取\(P_t = P_{t|t} = \bar P\), \(K_t = \bar K\), \(L_t = \bar L\), \(F_t = \bar F\)， 不需要再计算\(P_{t|t}\)、\(K_t\)、\(L_t\)、\(F_t\)和\(P_{t+1}\)， 可以节省计算量。

30.1.6 状态和观测值的一步预测误差

观测值的一步预报误差是\(\boldsymbol v_t = \boldsymbol y_t - E(\boldsymbol y_t | \boldsymbol Y_{t-1}) = \boldsymbol y_t - Z_t \boldsymbol a_t\)。 由定理28.1可知\(\boldsymbol v_t\)与\(\boldsymbol Y_{t-1}\)独立， 从而各个\(\boldsymbol v_t\)相互独立。 有 \[ p(\boldsymbol y_1, \dots, \boldsymbol y_n) = p(\boldsymbol y_1) \prod_{t=2}^n p(\boldsymbol y_t | \boldsymbol Y_{t-1}) = \prod_{t=1}^n p(\boldsymbol v_t) \]

令 \[ \boldsymbol x_t = \boldsymbol\alpha_t - E(\boldsymbol\alpha_t | \boldsymbol Y_{t-1}) = \boldsymbol\alpha_t - \boldsymbol a_t, \] 称\(\boldsymbol x_t\)为状态的一步预测误差。 下面给出\(\boldsymbol x_t, \boldsymbol v_t\)的模型。 \[\begin{aligned} \boldsymbol v_t =& \boldsymbol y_t - Z_t \boldsymbol a_t \\ =& Z_t \boldsymbol\alpha_t + \boldsymbol\varepsilon_t - Z_t \boldsymbol a_t \\ =& Z_t \boldsymbol x_t + \boldsymbol\varepsilon_t . \end{aligned}\] 而 \[\begin{aligned} \boldsymbol x_{t+1} =& \boldsymbol\alpha_{t+1} - \boldsymbol a_{t+1} \\ =& T_t \boldsymbol\alpha_t + R_t \boldsymbol\eta_t - T_t \boldsymbol a_t - K_t \boldsymbol v_t \\ =& T_t \boldsymbol x_t + R_t \boldsymbol\eta_t - K_t Z_t \boldsymbol x_t - K_t \boldsymbol\varepsilon_t \\ =& L_t \boldsymbol x_t + (R_t \boldsymbol\eta_t - K_t \boldsymbol\varepsilon_t) . \end{aligned}\] 统一写成 \[\begin{align} \boldsymbol v_t =& Z_t \boldsymbol x_t + \boldsymbol\varepsilon_t, \tag{30.12} \\ \boldsymbol x_{t+1} =& L_t \boldsymbol x_t + (R_t \boldsymbol\eta_t - K_t \boldsymbol\varepsilon_t) . \tag{30.13} \end{align}\] 这是类似于SSM的模型， \(\boldsymbol x_t\)作为状态向量， \(\boldsymbol v_t\)作为观测值向量。 与模型(30.1)、(30.2)不完全相符的是这两个方程的误差项不是相互独立的。

30.2 状态平滑算法

30.2.1 问题和记号

考虑给定所有观测\(\boldsymbol Y_n\)条件下状态\(\boldsymbol\alpha_t\)的条件分布， 包括联合分布。 \((\boldsymbol\alpha_1, \dots, \boldsymbol\alpha_t)\)在\(\boldsymbol Y_n\)给定条件下的联合分布也是多元正态分布， 所以只需要计算\(\hat{\boldsymbol\alpha}_t = E(\boldsymbol\alpha_t | \boldsymbol Y_n)\)， \(V_t = \text{Var}(\boldsymbol\alpha_t | \boldsymbol Y_n)\)， \(\text{Cov}(\boldsymbol\alpha_t, \boldsymbol\alpha_s | \boldsymbol Y_n)\)。 称\(\hat{\boldsymbol\alpha}_t\)为状态平滑， 简称平滑。

仍设\(\boldsymbol\alpha_1 \sim \text{N}(\boldsymbol a_1, P_1)\), \(\boldsymbol a_1, P_1\)已知。

若考虑\(E(\boldsymbol\alpha_t | \boldsymbol y_s, \dots, \boldsymbol y_n)\)， 其中\(s, n\)给定， \(t \in [s, n]\)，则称这样的问题为固定区间平滑。

若考虑\(E(\boldsymbol\alpha_t | \boldsymbol y_1, \dots, \boldsymbol y_n)\)， 但\(t\)固定， \(n=t+1, t+2, \dots\)可变， 则称为“定点”平滑。

若考虑\(E(\boldsymbol\alpha_{n-j} | \boldsymbol y_1, \dots, \boldsymbol y_n)\)， 其中\(j\)固定， \(n=j+1, j+2, \dots\)可变， 则称为“定滞后”平滑。

本节考虑的是给定所有样本\(\boldsymbol y_1, \dots, \boldsymbol y_n\)条件下\(\boldsymbol\alpha_t\)的估计， 是定区间平滑问题。

30.2.2 状态平滑

记 \[ \boldsymbol v_{t:n} = \begin{pmatrix} \boldsymbol v_t \\ \vdots \\ \boldsymbol v_n \end{pmatrix} . \] \(\boldsymbol Y_n\)与\((\boldsymbol Y_{t-1}, \boldsymbol v_{t:n})\)可以互相决定。 所以条件分布\(p(\boldsymbol\alpha_t | \boldsymbol Y_n)\)， 等于\(p(\boldsymbol\alpha_t | \boldsymbol Y_{t-1}, \boldsymbol v_{t:n})\)， 其中\(\boldsymbol Y_{t-1}\)与\(\boldsymbol v_{t:n}\)相互独立。 利用定理28.1来计算状态平滑： \[\begin{aligned} \hat{\boldsymbol\alpha}_t =& E(\boldsymbol\alpha_t | \boldsymbol Y_n) = E(\boldsymbol\alpha_t | \boldsymbol Y_{t-1}, \boldsymbol v_{t:n}) \\ =& \boldsymbol a_t + \text{Cov}(\boldsymbol\alpha_t, \boldsymbol v_{t:n} | \boldsymbol Y_{t-1}) [\text{Var}(\boldsymbol v_{t:n} | \boldsymbol Y_{t-1})]^{-1} \boldsymbol v_{t:n} \\ =& \boldsymbol a_t + \sum_{j=t}^n \text{Cov}(\boldsymbol\alpha_t, \boldsymbol v_j | \boldsymbol Y_{t-1}) F_j^{-1} \boldsymbol v_j . \end{aligned}\] 这里以\(\boldsymbol Y_{t-1}\)条件下的条件分布为基础计算条件期望。 注意\(\boldsymbol v_j\)与\(\boldsymbol Y_{t-1}\)相互独立， 所以\(\text{Var}(\boldsymbol v_t | \boldsymbol Y_{t-1}) = \text{Var}(\boldsymbol v_t) = F_t\)。

需要计算\(\text{Cov}(\boldsymbol\alpha_t, \boldsymbol v_j | \boldsymbol Y_{t-1})\)， \(j=t,t+1,\dots,n\)。 \[\begin{aligned} \text{Cov}(\boldsymbol\alpha_t, \boldsymbol v_j | \boldsymbol Y_{t-1}) =& E(\boldsymbol\alpha_t \boldsymbol v_j^T | \boldsymbol Y_{t-1}) \\ =& E[\boldsymbol\alpha_t (Z_j \boldsymbol x_j + \boldsymbol\varepsilon_j)^T | \boldsymbol Y_{t-1}] \\ =& E(\boldsymbol\alpha_t \boldsymbol x_j^T | \boldsymbol Y_{t-1}) Z_j^T . \end{aligned}\] 利用(30.13)可得 \[\begin{aligned} E(\boldsymbol\alpha_t \boldsymbol x_t^T | \boldsymbol Y_{t-1}) =& E[\boldsymbol\alpha_t (\boldsymbol\alpha_t - \boldsymbol a_t)^T | \boldsymbol Y_{t-1}] = \text{Var}(\boldsymbol\alpha_t | \boldsymbol Y_{t-1}) = P_t, \\ E(\boldsymbol\alpha_t \boldsymbol x_{t+1}^T | \boldsymbol Y_{t-1}) =& E[\boldsymbol\alpha_t (L_t \boldsymbol x_t + R_t \boldsymbol\eta_t - K_t \boldsymbol\varepsilon_t)^T | \boldsymbol Y_{t-1}] \\ =& E(\boldsymbol\alpha_t \boldsymbol x_t^T | \boldsymbol Y_{t-1}) L_t^T = P_t L_t^T, \\ E(\boldsymbol\alpha_t \boldsymbol x_{t+2}^T | \boldsymbol Y_{t-1}) =& E[\boldsymbol\alpha_t (L_{t+1} \boldsymbol x_{t+1} + R_{t+1} \boldsymbol\eta_{t+1} - K_{t+1} \boldsymbol\varepsilon_{t+1})^T | \boldsymbol Y_{t-1}] \\ =& E(\boldsymbol\alpha_t \boldsymbol x_{t+1}^T | \boldsymbol Y_{t-1}) L_{t+1}^T = P_t L_t^T L_{t+1}^T, \\ & \cdots\cdots \\ E(\boldsymbol\alpha_t \boldsymbol x_{n}^T | \boldsymbol Y_{t-1}) =& P_t L_t^T L_{t+1}^T \cdots L_{n-1}^T . \end{aligned}\]

当\(t=n\)时， \(L_t^T L_{t+1}^T \cdots L_{n-1}^T\)应看作\(I_m\)。 当\(t=n-1\)时， \(L_t^T L_{t+1}^T \cdots L_{n-1}^T\)应看作\(L_{n-1}^T\)。 有了\(\text{Cov}(\boldsymbol\alpha_t, \boldsymbol v_j | \boldsymbol Y_{t-1})\)公式后， 就可得 \[\begin{aligned} \hat{\boldsymbol\alpha}_n =& \boldsymbol a_n + P_n Z_n^T F_n^{-1} \boldsymbol v_n, \\ \hat{\boldsymbol\alpha}_{n-1} =& \boldsymbol a_{n-1} + P_{n-1} Z_{n-1}^T F_{n-1}^{-1} \boldsymbol v_{n-1} + P_{n-1} L_n^T Z_{n-1}^T F_{n-1}^{-1} \boldsymbol v_{n}, \\ \hat{\boldsymbol\alpha}_t =& E(\boldsymbol\alpha_t | \boldsymbol Y_n) \\ =& \boldsymbol a_t + P_t Z_t^T F_t^{-1} \boldsymbol v_t + P_t L_t^T Z_{t+1}^T F_{t+1}^{-1} \boldsymbol v_{t+1} + P_t L_t^T L_{t+1}^T Z_{t+2}^T F_{t+2}^{-1} \boldsymbol v_{t+2} \\ &\ + \cdots + P_t L_t^T L_{t+1}^T \cdots L_{n-1}^T Z_n^T F_n^{-1} \boldsymbol v_n . \ t=n-2,n-3, \dots, 2, 1 . \end{aligned}\]

记 \[\begin{aligned} \boldsymbol r_n =& \boldsymbol 0, \\ \boldsymbol r_{n-1} =& Z_n^T F_n^{-1} \boldsymbol v_n, \end{aligned}\] 并令 \[\begin{align} \boldsymbol r_{t-1} =& Z_t^T F_t^{-1} \boldsymbol v_t + L_t^T Z_{t+1}^T F_{t+1}^{-1} \boldsymbol v_{t+1} + L_t^T L_{t+1}^T Z_{t+2}^T F_{t+2}^{-1} \boldsymbol v_{t+2} \\ &\ + \cdots + L_t^T L_{t+1}^T \cdots L_{n-1}^T Z_n^T F_n^{-1} \boldsymbol v_n . \ t=n-2,n-3, \dots, 2, 1. \tag{30.14} \end{align}\] 注意 \[\begin{aligned} \boldsymbol r_{t} =& Z_{t+1}^T F_{t+1}^{-1} \boldsymbol v_{t+1} + L_{t+1}^T Z_{t+2}^T F_{t+2}^{-1} \boldsymbol v_{t+2} \\ &\ + \cdots + L_{t+1}^T L_{t+2}^T \cdots L_{n-1}^T Z_n^T F_n^{-1} \boldsymbol v_n . \end{aligned}\] 比较可得 \[\begin{aligned} \boldsymbol r_{t-1} =& Z_t^T F_t^{-1} \boldsymbol v_t + L^T \boldsymbol r_t, \ t=n, n-1, \dots, 1 . \end{aligned}\]

于是，状态平滑的算法为： \[\begin{equation} \begin{aligned} \boldsymbol r_n =& \boldsymbol 0, \\ \boldsymbol r_{t-1} =& Z_t^T F_t^{-1} \boldsymbol v_t + L^T \boldsymbol r_t, \ t=n, n-1, \dots, 1 . \\ \hat{\boldsymbol\alpha}_t =& \boldsymbol a_t + P_t \boldsymbol r_{t-1} . \end{aligned} \tag{30.15} \end{equation}\]

\(\boldsymbol r_{t-1}\)是\(\boldsymbol v_t, \boldsymbol v_{t+1}, \dots, \boldsymbol v_n\)的线性组合。 \(\boldsymbol a_t\)是\(\boldsymbol Y_{t-1}\)的线性组合。 于是\(\boldsymbol a_t + P_t \boldsymbol r_{t-1}\)是\(1:(t-1)\)的信息与\(t:n\)的信息的组合。

30.2.3 状态平滑方差阵

下面给出\(\text{Var}(\boldsymbol\alpha_t | \boldsymbol Y_n)\)的算法。 由定理28.1可得 \[\begin{aligned} V_t =& \text{Var}(\boldsymbol \alpha_t | \boldsymbol Y_{t-1}, \boldsymbol v_{t:n}) \\ =& \text{Var}(\boldsymbol \alpha_t | \boldsymbol Y_{t-1}) - \text{Cov}(\boldsymbol \alpha_t, \boldsymbol v_{t:n} | \boldsymbol Y_{t-1}) [\text{Var}(\boldsymbol v_{t:n} | \boldsymbol Y_{t-1})]^{-1} [\text{Cov}(\boldsymbol \alpha_t, \boldsymbol v_{t:n} | \boldsymbol Y_{t-1})]^T \\ =& P_t - \sum_{j=t}^n \text{Cov}(\boldsymbol \alpha_t, \boldsymbol v_{j} | \boldsymbol Y_{t-1}) [\text{Var}(\boldsymbol v_{j} | \boldsymbol Y_{t-1})]^{-1} [\text{Cov}(\boldsymbol \alpha_t, \boldsymbol v_{j} | \boldsymbol Y_{t-1})]^T \\ =& P_t - \sum_{j=t}^n \text{Cov}(\boldsymbol \alpha_t, \boldsymbol v_{j} | \boldsymbol Y_{t-1}) F_j^{-1} [\text{Cov}(\boldsymbol \alpha_t, \boldsymbol v_{j} | \boldsymbol Y_{t-1})]^T . \end{aligned}\] 由上一小节可知 \[\begin{aligned} \text{Cov}(\boldsymbol \alpha_t, \boldsymbol v_{j} | \boldsymbol Y_{t-1}) = P_t L_t^T L_{t+1}^T \cdots L_{j-1}^T Z_j^T, \end{aligned}\] 可得 \[\begin{aligned} V_t =& P_t - P_t Z_t^T F_t^{-1} Z_t P_t - P_t L_t^T Z_{t+1}^T F_{t+1}^{-1} Z_{t+1} L_t P_t \\ & - P_t L_t^T L_{t+1}^T Z_{t+2}^T F_{t+2}^{-1} Z_{t+2} L_{t+1} L_t P_t - \cdots \\ & - P_t L_t^T L_{t+1}^T \cdots L_{n-1}^T Z_n^T F_n^{-1} Z_n L_{n-1} \cdots L_{t+1} L_t P_t \\ =& P_t - P_t N_{t-1} P_t, \end{aligned}\] 其中 \[\begin{equation} \begin{aligned} N_{t-1} =& Z_t^T F_t^{-1} Z_t + L_t^T Z_{t+1}^T F_{t+1}^{-1} Z_{t+1} L_t \\ & + L_t^T L_{t+1}^T Z_{t+2}^T F_{t+2}^{-1} Z_{t+2} L_{t+1} L_t + \cdots \\ & + L_t^T L_{t+1}^T \cdots L_{n-1}^T Z_n^T F_n^{-1} Z_n L_{n-1} \cdots L_{t+1} L_t . \end{aligned} \tag{30.16} \end{equation}\]

当\(t=n\)时\(L_t^T L_{t+1}^T \cdots L_{n-1}^T\)应看作\(I_m\)， 当\(t=n-1\)时应看作\(L_{n-1}^T\)。 对比 \[\begin{aligned} N_t =& Z_{t+1}^T F_{t+1}^{-1} Z_{t+1} + L_{t+1}^T Z_{t+2}^T F_{t+2}^{-1} Z_{t+2} L_{t+1} \\ & + \cdots + L_{t+1}^T \cdots L_{n-1}^T Z_n^T F_{n}^{-1} Z_n L_{n-1} \cdots L_{t+1}, \end{aligned}\] 可见 \[\begin{aligned} N_{t-1} =& Z_t F_t^{-1} Z_t + L_t^T N_t L_t, \ t=n,n-1,\dots, 2, 1. \end{aligned}\] 并应取\(N_n = 0\)。

于是， \(\text{Var}(\boldsymbol\alpha_t | \boldsymbol Y_n)\)的递推算法为 \[\begin{equation} \begin{aligned} N_n =& 0_{m \times m}, \\ N_{t-1} =& Z_t^T F_t^{-1} Z_t + L_t^T N_t L_t, \ t=n, n-1, \dots, 2, 1; \\ V_t =& P_t - P_t N_{t-1} P_t . \end{aligned} \tag{30.17} \end{equation}\]

由(30.14)和(30.16)， 以及\(\boldsymbol v_t, \boldsymbol v_{t+1}, \dots, \boldsymbol v_n\)相互独立， 可知 \[\begin{aligned} \text{Var}(\boldsymbol r_{t-1}) = N_{t-1}, \ t=n,n-1,\dots,1 . \end{aligned}\]

30.2.4 状态平滑算法汇总

将状态平滑和方差阵计算的算法罗列如下： \[\begin{equation} \begin{aligned} \boldsymbol r_n =& \boldsymbol 0_{m \times 1}, & N_n =& 0_{m \times m} , \\ \boldsymbol r_{t-1} =& Z_t^T F_t^{-1} \boldsymbol v_t + L_t^T \boldsymbol r_t, & N_{t-1} =& Z_t^T F_t^{-1} Z_t + L_t^T N_t L_t, \\ &&& t=n,n-1,\dots,1 . \\ \hat{\boldsymbol\alpha}_t =& \boldsymbol a_t + P_t \boldsymbol r_{t-1}, & V_t =& P_t - P_t N_{t-1} P_t . \end{aligned} \tag{30.18} \end{equation}\] 这称为“状态平滑递推公式”。 其中 \[\begin{aligned} \hat{\boldsymbol\alpha}_t =& E(\boldsymbol\alpha_t | \boldsymbol Y_n), & V_t =& \text{Var}(\boldsymbol\alpha_t | \boldsymbol Y_n) . \end{aligned}\] \(\boldsymbol r_{t-1}\)与\(N_{t-1}\)满足 \[\begin{aligned} \text{Var}(\boldsymbol r_{t-1}) = N_{t-1} . \end{aligned}\]

结合滤波递推公式(30.11)和状态平滑递推公式(30.18)， 可以统称为“卡尔曼滤波平滑算法”。 (30.11)是向前递推， 可以获得\(t=1,2,\dots,n\)时刻的\(\boldsymbol a_t\), \(P_t\), \(\boldsymbol v_t\)，\(K_t\), \(L_t\)， 并可以计算似然函数值; (30.18)是向后递推， 可以获得\(t=n,n-1,\dots,1\)时刻的平滑分布参数\(\hat{\boldsymbol\alpha}_t\), \(V_t\)， 以及\(\boldsymbol r_{t-1}\), \(N_{t-1}\)。

30.2.5 平滑结果的在线更新

基于\(\boldsymbol Y_n\)得到平滑结果后， 又获得了新观测\(\boldsymbol y_{n+1}\)， 要更新平滑结果。 仍由定理28.1可知 \[\begin{aligned} \hat{\boldsymbol\alpha}_{t|n} =& E(\boldsymbol\alpha_t | \boldsymbol Y_n), \\ \hat{\boldsymbol\alpha}_{t|n+1} =& E(\boldsymbol\alpha_t | \boldsymbol Y_{n+1}) = E(\boldsymbol\alpha_t | \boldsymbol Y_n, \boldsymbol v_{n+1}) \\ =& E(\boldsymbol\alpha_t | \boldsymbol Y_n) - \text{Cov}(\boldsymbol\alpha_t, \boldsymbol v_{n+1} | \boldsymbol Y_n) [\text{Var}(\boldsymbol v_{n+1} | \boldsymbol Y_n)]^{-1} \boldsymbol v_{n+1} \\ =& \hat{\boldsymbol\alpha}_{t|n} + P_t L_t^T \cdots L_n^T Z_{t+1}^T F_{n+1}^{-1} \boldsymbol v_{n+1}, \ t=1,2,\dots,n; \\ \hat{\boldsymbol\alpha}_{n+1|n+1} =& \boldsymbol a_{n+1} + P_{n+1} Z_{n+1}^T F_{n+1}^{-1} \boldsymbol v_{n+1} . \\ V_{t|n} =& \text{Var}(\boldsymbol\alpha_t | \boldsymbol Y_{n}), \\ V_{t|n+1} =& \text{Var}(\boldsymbol\alpha_t | \boldsymbol Y_{n+1}) = \text{Var}(\boldsymbol\alpha_t | \boldsymbol Y_{n}, \boldsymbol v_{n+1}) \\ =& \text{Var}(\boldsymbol\alpha_t | \boldsymbol Y_{n}) - \text{Cov}(\boldsymbol\alpha_t, \boldsymbol v_{n+1} | \boldsymbol Y_{n}) [\text{Var}(\boldsymbol v_{n+1} | \boldsymbol Y_{n})]^{-1} [\text{Cov}(\boldsymbol\alpha_t, \boldsymbol v_{n+1} | \boldsymbol Y_{n})]^T \\ =& V_{t|n} - P_t L_t^T \cdots L_n^T Z_{n+1}^T F_{n+1}^{-1} Z_{n+1} L_n \cdots L_t P_t, \ t=1,2,\dots,n ; \\ V_{n+1|n+1} =& P_{n+1} - P_{n+1} Z_{n+1}^T F_{n+1}^{-1} Z_{n+1} P_{n+1} . \end{aligned}\]

30.2.6 定点平滑与定滞后平滑

若对给定\(t\)， 令\(n=t+1, t+2, \dots\)要计算\(\boldsymbol\alpha_t | \boldsymbol Y_{n}\)的条件分布， 可以先针对\(1:t\)时间点进行滤波平滑， 然后利用上一小节结果获得\(n=t+1\)时刻的平滑结果， 反复利用上一小节结果即可递推得到\(n=t+2,t+3,\dots\)的平滑结果。

如果对固定的滞后值\(j\)， 要对\(n=j+1,j+2,\dots\)， 计算\(\hat{\boldsymbol\alpha}_{n-j|n}\)， 需要对固定的\(n\)反向递推计算 \[\begin{aligned} \boldsymbol r_n^{(n)} = \boldsymbol 0, \\ \boldsymbol r_{t-1}^{(n)} =& Z_t^T F_t^{-1} \boldsymbol v_t + L_t^T \boldsymbol r_{t}^{(n)}, \ t=n, n-1, \dots, n-j. \end{aligned}\] 然后得 \[\begin{aligned} \hat{\boldsymbol\alpha}_{n-j|n} =& \boldsymbol a_{n-j} + P_{n-j} \boldsymbol r_{n-j-1}^{(n)} . \end{aligned}\] 对每个\(n\)，反向递推计算 \[\begin{aligned} N_n^{(n)} =& 0, \\ N_{t-1}^{(n)} =& Z_t^T F_t^{-1} Z_t + L_t N_{t}^{(n)} L_t, \ t=n,n-1,\dots,n-j, \end{aligned}\] 然后得 \[\begin{aligned} V_{n-j|n} =& \text{Var}(\boldsymbol\alpha_{n-j} | \boldsymbol Y_n) = P_{n-j} - P_{n-j} N_{n-j-1}^{(n)} P_{n-j} . \end{aligned}\]

30.3 扰动项的平滑算法

考虑给定\(\boldsymbol Y_n\)后， 观测扰动项（误差）\(\boldsymbol\varepsilon_t\)和系统扰动项（误差）\(\boldsymbol\eta_t\)的估计问题。 这可以用于参数估计和模型诊断。

30.3.1 条件均值

记\(\hat{\boldsymbol\varepsilon}_t = E(\boldsymbol\varepsilon | \boldsymbol Y_n)\)。 由定理28.1有 \[\begin{aligned} \hat{\boldsymbol\varepsilon}_t =& E(\boldsymbol\varepsilon_t | \boldsymbol Y_t, \boldsymbol v_t, \dots, \boldsymbol v_n) \\ =& E(\boldsymbol\varepsilon_t | \boldsymbol Y_t) + \sum_{j=t}^n \text{Cov}(\boldsymbol\varepsilon_t, \boldsymbol v_j) F_j^{-1} \boldsymbol v_j \\ =& \sum_{j=t}^n E(\boldsymbol\varepsilon_t \boldsymbol v_j^T) F_j^{-1} \boldsymbol v_j . \end{aligned}\] 由(30.12)， \[\begin{aligned} E(\boldsymbol\varepsilon_t \boldsymbol v_j^T) =& E(\boldsymbol\varepsilon_t \boldsymbol x_j^T) Z_j^T + E(\boldsymbol\varepsilon_t \boldsymbol\varepsilon_j^T) \\ =& \begin{cases} H_t, & j = t, \\ E(\boldsymbol\varepsilon_t \boldsymbol x_j^T) Z_j^T, & j=t+1, \dots, n . \end{cases} \end{aligned}\] 由(30.13)， \[\begin{aligned} E(\boldsymbol\varepsilon_t \boldsymbol x_{t+1}^T) =& E \{ \boldsymbol\varepsilon_t [L_t \boldsymbol x_t + (R_t \boldsymbol\eta_t - K_t \boldsymbol\varepsilon_t)]^T \} \\ =& E(\boldsymbol\varepsilon_t \boldsymbol x_t^T) L_t^T + E(\boldsymbol\varepsilon_t \boldsymbol\eta_t^T) R_t^T - E(\boldsymbol\varepsilon_t \boldsymbol\varepsilon_t^T) K_t^T \\ =& -H_t K_t^T; \\ E(\boldsymbol\varepsilon_t \boldsymbol x_{t+2}^T) =& E \{ \boldsymbol\varepsilon_t [L_{t+1} \boldsymbol x_{t+1} + (R_{t+1} \boldsymbol\eta_{t+1} - K_{t+1} \boldsymbol\varepsilon_{t+1})]^T \} \\ =& E (\boldsymbol\varepsilon_t \boldsymbol x_{t+1}^T) L_{t+1}^T = -H_t K_t^T L_{t+1}^T; \\ & \cdots\cdots \\ E(\boldsymbol\varepsilon_t \boldsymbol x_{n}^T) =& -H_t K_t^T L_{t+1}^T \cdots L_{n-1}^T , \ t=1,2,\dots,n-1 . \end{aligned}\] 其中\(L_{t+1}^T \cdots L_{n-1}^T\)当\(t=n-1\)时应视为\(I_m\)， 当\(t=n-2\)时应视为\(L_{n-1}^T\)。

于是， \[\begin{equation} \begin{aligned} \hat{\boldsymbol\varepsilon}_t =& H_t \left( F_t^{-1} \boldsymbol v_t - K_t^T Z_{t+1}^T F_{t+1}^{-1} \boldsymbol v_{t+1} - K_t^T L_{t+1}^T Z_{t+2}^T F_{t+2}^{-1} \boldsymbol v_{t+2} \right. \\ & \left. - \cdots - K_t^T L_{t+1}^T \cdots L_{n-1}^T Z_{n}^T F_{n}^{-1} \boldsymbol v_{n} \right) \\ =& H_t (F_t^{-1} \boldsymbol v_t - K_t^T \boldsymbol r_t) \\ =& H_t \boldsymbol u_t, \ t=n,n-1,\dots,1 . \end{aligned} \tag{30.19} \end{equation}\] 其中 \[\begin{equation} \begin{aligned} \boldsymbol u_t =& F_t^{-1} \boldsymbol v_t - K_t^T \boldsymbol r_t, \end{aligned} \tag{30.20} \end{equation}\] 称为“平滑误差”。

记\(\hat{\boldsymbol\eta}_t = E(\boldsymbol\eta_t | \boldsymbol Y_n)\)。 由定理28.1有 \[\begin{aligned} \hat{\boldsymbol\eta}_t =& E(\boldsymbol\eta_t | \boldsymbol Y_{t-1}, \boldsymbol v_t, \dots, \boldsymbol v_n) \\ =& E(\boldsymbol\eta_t | \boldsymbol Y_{t-1}) + \sum_{j=t}^n \text{Cov}(\boldsymbol\eta_t, \boldsymbol v_j) F_j^{-1} \boldsymbol v_j \\ =& \sum_{j=t}^n \text{Cov}(\boldsymbol\eta_t, \boldsymbol v_j) F_j^{-1} \boldsymbol v_j . \end{aligned}\] 注意\(\boldsymbol v_t\)仅依赖于截止到\(\boldsymbol\alpha_t\)和\(\boldsymbol Y_{t-1}\)， \(\boldsymbol\eta_t\)与这两项独立， 所以上式\(j=t\)的项等于0。

对\(j=t+1,t+2,\dots,n\)， 利用(30.12)， 以及\(\boldsymbol x_t\)仅依赖于\(\boldsymbol\alpha_t\)和\(\boldsymbol Y_{t-1}\)， 可知\(\boldsymbol x_t\)与\(\boldsymbol\eta_t\)独立， \[\begin{aligned} \text{Cov}(\boldsymbol\eta_t, \boldsymbol v_j) =& E(\boldsymbol\eta_t \boldsymbol v_j^T) = E[\boldsymbol\eta_t (Z_j \boldsymbol x_j + \boldsymbol\varepsilon_j)^T] \\ =& E(\boldsymbol\eta_t \boldsymbol x_j^T) Z_j^T . \end{aligned}\] 其中\(E(\boldsymbol\eta_t \boldsymbol x_t^T) = 0\)， \[\begin{aligned} E(\boldsymbol\eta_t \boldsymbol x_{t+1}^T) =& E \{ \boldsymbol\eta_t [L_t \boldsymbol x_t + (R_t \boldsymbol\eta_t -K_t \boldsymbol\varepsilon_t)]^T\} \\ =& Q_t R_t^T; \\ E(\boldsymbol\eta_t \boldsymbol x_{t+2}^T) =& E \{ \boldsymbol\eta_t [L_{t+1} \boldsymbol x_{t+1} + (R_{t+1} \boldsymbol\eta_{t+1} - K_{t+1} \boldsymbol\varepsilon_{t+1})]^T\} \\ =& E (\boldsymbol\eta_t \boldsymbol x_{t+1}^T ) L_{t+1}^T Q_t R_t^T L_{t+1}^T; \\ E(\boldsymbol\eta_t \boldsymbol x_{t+2}^T) =& Q_t R_t^T L_{t+1}^T L_{t+2}^T; \\ & \cdots\cdots \\ E(\boldsymbol\eta_t \boldsymbol x_{n}^T) =& Q_t R_t^T L_{t+1}^T \cdots L_{n-1}^T . \end{aligned}\] 于是 \[\begin{equation} \begin{aligned} \hat{\boldsymbol\eta}_t =& Q_t R_t^T \left( Z_{t+1}^T F_{t+1}^{-1} \boldsymbol v_{t+1} + L_{t+1}^T Z_{t+2}^T F_{t+2}^{-1} \boldsymbol v_{t+2} \right. \\ & \left. + \cdots + L_{t+1}^T \cdots L_{n-1}^T Z_{n}^T F_{n}^{-1} \boldsymbol v_{n} \right) \\ =& Q_t R_t^T \boldsymbol r_t, \ t=n,n-1,\dots,1 . \end{aligned} \tag{30.21} \end{equation}\]

这也给出了\(\boldsymbol r_t\)的一个解释。 \(\boldsymbol r_t\)包含了\(\hat{\boldsymbol\eta}_t\)的信息。

30.3.2 条件方差阵

由定理28.1有 \[\begin{aligned} \text{Var}(\boldsymbol\varepsilon_t | \boldsymbol Y_n) =& \text{Var}(\boldsymbol\varepsilon_t | \boldsymbol Y_{t-1}, \boldsymbol v_t, \dots, \boldsymbol v_n) \\ =& \text{Var}(\boldsymbol\varepsilon_t | \boldsymbol Y_{t-1}) - \sum_{j=t}^n \text{Cov}(\boldsymbol\varepsilon_t, \boldsymbol v_j) F_j^{-1} [\text{Cov}(\boldsymbol\varepsilon_t, \boldsymbol v_j)]^T \\ =& H_t - \sum_{j=t}^n E(\boldsymbol\varepsilon_t \boldsymbol v_j^T) F_j^{-1} [E(\boldsymbol\varepsilon_t \boldsymbol v_j^T)]^T . \end{aligned}\] 由前一小节可知其中 \[\begin{aligned} E(\boldsymbol\varepsilon_t \boldsymbol v_t^T) =& H_t, \\ E(\boldsymbol\varepsilon_t \boldsymbol v_j^T) =& -H_t K_t^T L_{t+1}^T \cdots L_{j-1}^T Z_j^T, \ j=t+1, \dots, n . \end{aligned}\] 从而 \[\begin{aligned} \text{Var}(\boldsymbol\varepsilon_t | \boldsymbol Y_n) =& H_t - H_t F_t^{-1} H_t - H_t K_t^T Z_{t+1}^T F_{t+1}^{-1} Z_{t+1} K_t H_t \\ & - H_t K_t^T L_{t+1}^T Z_{t+2}^T F_{t+2}^{-1} Z_{t+2} L_{t+1} K_t H_t \\ & - H_t K_t^T L_{t+1}^T L_{t+2}^T Z_{t+3}^T F_{t+3}^{-1} Z_{t+3} L_{t+2} L_{t+1} K_t H_t \\ & - \cdots \\ & - H_t K_t^T L_{t+1}^T \cdots L_{n-1}^T Z_{n}^T F_{n}^{-1} Z_{n} L_{n-1} \cdots L_{t+1} K_t H_t \\ =& H_t - H_t \left[ F_t^{-1} + K_t^T Z_{t+1}^T F_{t+1}^{-1} Z_{t+1} K_t \right. \\ & + K_t L_{t+1}^T Z_{t+2}^T F_{t+2}^{-1} Z_{t+2} L_{t+1} K_t \\ & + K_t L_{t+1}^T L_{t+2}^T Z_{t+3}^T F_{t+3}^{-1} Z_{t+3} L_{t+2} L_{t+1} K_t \\ & + \cdots \\ & \left. + K_t L_{t+1}^T \cdots L_{n-1}^T Z_{n}^T F_{n}^{-1} Z_{n} L_{n-1} \cdots L_{t+1} K_t \right] \\ =& H_t - H_t \left[ F_t^{-1} + K_t^T \left( Z_{t+1}^T F_{t+1}^{-1} Z_{t+1} \right. \right. \\ & + L_{t+1}^T Z_{t+2}^T F_{t+2}^{-1} Z_{t+2} L_{t+1} \\ & + L_{t+1}^T L_{t+2}^T Z_{t+3}^T F_{t+3}^{-1} Z_{t+3} L_{t+2} L_{t+1} \\ & + \cdots \\ & \left. \left. + L_{t+1}^T \cdots L_{n-1}^T Z_{n}^T F_{n}^{-1} Z_{n} L_{n-1} \cdots L_{t+1} \right) K_t \right] H_t \\ =& H_t - H_t (F_t^{-1} - K_t^T N_t K_t) H_t \\ =& H_t - H_t D_t H_t, \ t=n,n-1,\dots,1 . \end{aligned}\] 其中 \[\begin{aligned} D_t =& F_t^{-1} - K_t^T N_t K_t . \end{aligned}\]

类似地， \[\begin{aligned} \text{Var}(\boldsymbol\eta_t | \boldsymbol Y_n) =& \text{Var}(\boldsymbol\eta_t | \boldsymbol Y_{t-1}, \boldsymbol v_t, \dots, \boldsymbol v_n) \\ =& \text{Var}(\boldsymbol\eta_t | \boldsymbol Y_{t-1}) - \sum_{j=t}^n \text{Cov}(\boldsymbol\eta_t, \boldsymbol v_j) F_j^{-1} [\text{Cov}(\boldsymbol\eta_t, \boldsymbol v_j)]^T \\ =& Q_t - \sum_{j=t}^n E(\boldsymbol\eta_t \boldsymbol v_j^T) F_j^{-1} [E(\boldsymbol\eta_t \boldsymbol v_j^T)]^T, \end{aligned}\] 其中 \[\begin{aligned} E(\boldsymbol\eta_t \boldsymbol v_t^T) =& 0, \\ E(\boldsymbol\eta_t \boldsymbol v_{t+1}^T) =& Q_t R_t^T Z_{t+1}^T, \\ & \cdots\cdots \\ E(\boldsymbol\eta_t \boldsymbol v_{n}^T) =& Q_t R_t^T L_{t+1}^T \cdots L_{n-1}^T Z_n^T, \end{aligned}\] 于是 \[\begin{aligned} \text{Var}(\boldsymbol\eta_t | \boldsymbol Y_n) =& Q_t - Q_t R_t^T \left( Z_{t+1}^T F_{t+1}^{-1} Z_{t+1} \right. \\ & + L_{t+1}^T Z_{t+2}^T F_{t+2}^{-1} Z_{t+2} L_{t+1} \\ & + \cdots \\ & \left. + L_{t+1}^T \cdots L_{n-1}^T Z_n^T F_n^{-1} Z_n L_{n-1} \cdots L_{t+1} \right) R_t Q_t \\ =& Q_t - Q_t R_t^T N_t R_t Q_t, \ t=n,n-1,\dots,1 . \end{aligned}\]

30.3.3 扰动项平滑的递推计算公式汇总

为了获得\(\boldsymbol Y_n\)下\(\boldsymbol\varepsilon_t\)和\(\boldsymbol\eta_t\)的条件期望和条件方差， 应反向递推如下： \[\begin{equation} \begin{aligned} \boldsymbol r_n =& \boldsymbol 0, & N_n =& 0 , \\ \hat{\boldsymbol\varepsilon}_t =& H_t(F_t^{-1} \boldsymbol v_t - K_t^T \boldsymbol r_t), & \text{Var}(\boldsymbol\varepsilon_t | \boldsymbol Y_n) =& H_t - H_t(F_t^{-1} + K_t^T N_t K_t) H_t, \\ \hat{\boldsymbol\eta}_t =& Q_t R_t^T \boldsymbol r_t, & \text{Var}(\boldsymbol\eta_t | \boldsymbol Y_n) =& Q_t - Q_t R_t^T N_t R_t Q_t, \\ \boldsymbol r_{t-1} =& Z_t^T F_t^{-1} \boldsymbol v_t + L_t^T \boldsymbol r_t, & N_{t-1} =& Z_t^T F_t^{-1} Z_t + L_t^T N_t L_t, \\ & t=n,n-1, \dots, 1 . \end{aligned} \tag{30.22} \end{equation}\]

引入\(\boldsymbol u_t\)与\(D_t\)， 算法可以写成 \[\begin{equation} \begin{aligned} \boldsymbol r_n =& \boldsymbol 0, & N_n =& 0 , \\ \boldsymbol u_t =& F_t^{-1} \boldsymbol v_t - K_t^T \boldsymbol r_t, & D_t =& F_t^{-1} + K_t^T N_t K_t, \\ \hat{\boldsymbol\varepsilon}_t =& H_t \boldsymbol u_t, & \text{Var}(\boldsymbol\varepsilon_t | \boldsymbol Y_n) =& H_t - H_t D_t H_t, \\ \hat{\boldsymbol\eta}_t =& Q_t R_t^T \boldsymbol r_t, & \text{Var}(\boldsymbol\eta_t | \boldsymbol Y_n) =& Q_t - Q_t R_t^T N_t R_t Q_t, \\ \boldsymbol r_{t-1} =& Z_t^T \boldsymbol u_t + T_t^T \boldsymbol r_t, & N_{t-1} =& Z_t^T D_t Z_t + T_t^T N_t T_t \\ &&& - Z_t^T K_t^T N_t T_t - T_t^T N_t K_t Z_t, \\ & t=n,n-1, \dots, 1 . \end{aligned} \tag{30.23} \end{equation}\]

算法(30.23)直接利用了原始的\(T_t\)和\(Z_t\)， 这两个矩阵常常是稀疏矩阵， 所以后面的算法通常效率更高。

(30.23)中\(\boldsymbol r_{t-1}\)递推公式推导如下： \(\boldsymbol u_t = F_t^{-1} \boldsymbol v_t - K_t^T \boldsymbol r_t\)， 所以\(F_t^{-1} \boldsymbol v_t = \boldsymbol u_t + K_t^T \boldsymbol r_t\)， 代入(30.22)中\(\boldsymbol r_{t-1}\)递推公式可得 \[\begin{aligned} \boldsymbol r_{t-1} =& Z_t^T (\boldsymbol u_t + K_t^T \boldsymbol r_t) + L_t^T \boldsymbol r_t \\ =& Z_t^T \boldsymbol u_t + (K_t Z_t + L_t)^T \boldsymbol r_t \\ =& Z_t^T \boldsymbol u_t + T_t^T \boldsymbol r_t, \end{aligned}\] 这里用了\(L_t = T_t - K_t Z_t\)。

(30.23)中\(N_{t-1}\)的递推公式推导如下： 在(30.22)关于\(N_{t-1}\)的递推公式中， 代入\(F_t^{-1} = D_t - K_t^T N_t K_t\)， \(L_t = T_t - K_t Z_t\)， 化简后即得该递推公式。

30.4 关于平滑

由 \[\begin{aligned} \hat{\boldsymbol\alpha}_t =& \boldsymbol a_t + P_t \boldsymbol r_{t-1} , \end{aligned}\] 可见 \[\begin{aligned} \boldsymbol r_{t-1} =& P_t^{-1}(\hat{\boldsymbol\alpha}_t - \boldsymbol a_t). \end{aligned}\]

还有许多其它的平滑算法， 这里给出的算法是比较通用且高效的。

30.5 平滑分布的协方差

考虑\(\text{Cov}(\boldsymbol\alpha_t, \boldsymbol\alpha_s | \boldsymbol Y_n)\)的计算， 配合前面已有的条件均值、条件方差阵公式， 可以确定\((\boldsymbol\alpha_1, \dots, \boldsymbol\alpha_n)\)在\(\boldsymbol Y_n\)条件下的联合分布。 \(t=s\)时已有递推计算公式， 见30.2.3。

内容待完成。

30.6 从平滑分布抽样

通过对\(\boldsymbol Y_n\)条件下从\(\boldsymbol\alpha_t\), \(\boldsymbol\varepsilon_t\), \(\boldsymbol\eta_t\)条件分布抽样， 可以生成多条类似的样本轨道， 包括状态的轨道， 借以研究模型的表现。 更重要的是， 研究这样的抽样算法， 可以用于解决非高斯、非线性的状态空间模型推断问题。 这种抽样方法称为“抽样平滑”。

这一节考虑在\(\boldsymbol Y_n\)条件下从\(\boldsymbol\alpha_t\), \(\boldsymbol\varepsilon_t\), \(\boldsymbol\eta_t\)条件分布抽样， 算法可称为“向前滤波、向后抽样”算法。

30.6.1 抽样平滑的均值校正方法

考虑\(\boldsymbol Y_n\)条件下从\(\boldsymbol\varepsilon_t\), \(\boldsymbol\eta_t\)的条件分布抽样。 记 \[\begin{aligned} \boldsymbol w = \begin{pmatrix} \boldsymbol\varepsilon_1 \\ \boldsymbol\eta_1 \\ \vdots \\ \boldsymbol\varepsilon_n \\ \boldsymbol\eta_n \end{pmatrix} . \end{aligned}\] 令 \[\begin{aligned} \hat{\boldsymbol w} =& E(\boldsymbol w | \boldsymbol Y_n), & W =& \text{Var}(\boldsymbol w | \boldsymbol Y_n) . \end{aligned}\] 由线性高斯状态空间模型的高斯过程性质可知\(\boldsymbol w\)在\(\boldsymbol Y_n\)条件下的分布是\(\text{N}(\hat{\boldsymbol w}, W)\)。 30.3.3中已经给出了\(\hat{\boldsymbol w}\)的算法。 而计算\(W\)是应尽可能避免的， 其中的协方差阵会涉及到复杂与繁重的计算， 直接从这样的高维正态分布抽样涉及到很大的矩阵求逆， 算法不稳定。 这里要介绍的均值校正方法可以避免这样的高维抽样困难。

\(\boldsymbol w\)的无条件分布为 \[ p(\boldsymbol w) \sim \text{N}(\boldsymbol 0, \Phi), \quad \Phi = \text{diag}(H_1, Q_1, \dots, H_n, Q_n) . \]

设\(\boldsymbol w^+\)为从无条件分布\(p(\boldsymbol w)\)中抽取的随机向量。 这只要分别抽取\(\boldsymbol\varepsilon_t^+\), \(\boldsymbol\eta_t^+\), \(t=1,\dots,n\)。

设\(\boldsymbol\alpha_1 \sim p(\boldsymbol\alpha_1)\)已知， 从此无条件分布抽取随机向量\(\boldsymbol\alpha_1^+\)， 然后令\(\boldsymbol y_1^+ = Z_1 \boldsymbol\alpha_1^+ + \boldsymbol\varepsilon_1^+\)， \(\boldsymbol\alpha_2^+ = T_1 \boldsymbol\alpha_1^+ + R_1 \boldsymbol\eta_1^+\)， 如此递推可以得到\((\boldsymbol\alpha_1^+, \dots, \boldsymbol\alpha_n^+)\)和\((\boldsymbol y_1^+, \dots, \boldsymbol y_n^+)\)。 这样， 我们得到了一组新的观测轨道\((\boldsymbol y_1^+, \dots, \boldsymbol y_n^+)\)， 这组观测值符合原始的(30.1)-(30.1)模型分布。

记 \[\begin{aligned} \boldsymbol y^+ = \begin{pmatrix} \boldsymbol y_1^+ \\ \vdots \\ \boldsymbol y_n^+ \end{pmatrix} . \end{aligned}\] 则\(\boldsymbol y^+\)与\(\boldsymbol Y_n\)同分布， 且\((\boldsymbol w^+, \boldsymbol y^+)\)与\((\boldsymbol w, \boldsymbol Y_n)\)同分布。

将\(\boldsymbol y^+\)看成是一组观测数据， 重新进行向前滤波、向后平滑， 得到\(\boldsymbol w\)基于\(\boldsymbol y^+\)的平滑结果 \[ \hat{\boldsymbol w}^+ = E(\boldsymbol w | \boldsymbol y^+) . \] \(\hat{\boldsymbol w}^+\)与\(\hat{\boldsymbol w}\)同分布。

由正态分布性质， 条件方差不依赖于作为条件的随机变量， 所以 \[\begin{aligned} W =& \text{Var}(\boldsymbol w | \boldsymbol Y_n) = \text{Var}(\boldsymbol w^+ | \boldsymbol y^+) \\ =& E \{ [\boldsymbol w^+ - \hat{\boldsymbol w}^+] [\boldsymbol w^+ - \hat{\boldsymbol w}^+]^T \}, \end{aligned}\] 即 \[ \boldsymbol w^+ - \hat{\boldsymbol w}^+ \sim \text{N}(0, W) . \] 于是 \[ \tilde{\boldsymbol w} = \boldsymbol w^+ - \hat{\boldsymbol w}^+ + \hat{\boldsymbol w} \sim \text{N}(\hat{\boldsymbol w}, W) . \] \(\tilde{\boldsymbol w}\)是\(\boldsymbol w\)在\(\boldsymbol Y_n\)条件下的条件分布的抽样。 这避免了计算很大的方差阵\(W\)以及直接从高维的正态分布抽样的问题。

注意这种方法假定初始分布\(\boldsymbol\alpha_1 \sim \text{N}(\boldsymbol a_1, P_1)\)已知； 如果初始分布有部分未知参数， 或者是发散先验， 算法需要进行修改。

30.6.2 状态向量的抽样

记\(\boldsymbol\alpha = (\boldsymbol\alpha_1^T, \dots, \boldsymbol\alpha_n^T)^T\)。 记\(\hat{\boldsymbol\alpha} = E(\boldsymbol\alpha | \boldsymbol Y_n)\)， 这只要使用向前滤波、向后平滑即可计算。

考虑从条件分布\(p(\boldsymbol\alpha | \boldsymbol Y_n)\)抽样的问题， 设这样的样本为\(\tilde{\boldsymbol\alpha}\)。

类似上一小节， 可以从无条件分布产生\(\boldsymbol w^+\), \(\boldsymbol\alpha^+\)以及\(\boldsymbol y^+\)， 然后从\(\boldsymbol y^+\)平滑得到\(\hat{\boldsymbol\alpha}^+ = E(\boldsymbol\alpha | \boldsymbol y^+)\)。 然后只要令 \[ \tilde{\boldsymbol\alpha} = \boldsymbol\alpha^+ - \hat{\boldsymbol\alpha}^+ + \hat{\boldsymbol\alpha}, \] 则\(\tilde{\boldsymbol\alpha}\)就是条件分布\(p(\boldsymbol\alpha | \boldsymbol Y_n)\)的抽样。

30.7 缺失值处理

线性高斯状态空间模型对于有缺失值的情形很容易处理。 设\(\boldsymbol y_t\)在\(t=\tau, \dots, \tau^*\)区间缺失， 其中\(1 < \tau \leq \tau* < n\)。 一种办法是将时间轴重新定义， 令\(\tau^*+1\)为\(\tau-1\)后面的点。 这只要适当调整系统方程中误差项方差阵和转移矩阵即可。

另一种方法则更方便， 不需要调整时间轴。 对\(t=\tau, \dots, \tau^*\)， 有 \[\begin{aligned} \boldsymbol a_{t|t} =& E(\boldsymbol\alpha_t | \boldsymbol Y_t) = E(\boldsymbol\alpha_t | \boldsymbol Y_{t-1}) = \boldsymbol a_t, \\ P_{t|t} =& \text{Var}(\boldsymbol\alpha_t | \boldsymbol Y_t) = \text{Var}(\boldsymbol\alpha_t | \boldsymbol Y_{t-1}) = P_t, \\ \boldsymbol a_{t+1} =& E(\boldsymbol\alpha_{t+1} | \boldsymbol Y_t) = E(T_t \boldsymbol\alpha_t + R_t \boldsymbol\eta_t | \boldsymbol Y_{t-1}) = T_t \boldsymbol a_t, \\ P_{t+1} =& \text{Var}(\boldsymbol\alpha_{t+1} | \boldsymbol Y_t) = \text{Var}(T_t \boldsymbol\alpha_t + R_t \boldsymbol\eta_t | \boldsymbol Y_{t-1}) \\ =& T_t P_t T_t^T + R_t Q_t R_t^T . \end{aligned}\] 对比原来的向前滤波公式， 这相当于取\(Z_t=0\), \(K_t=0\), \(L_t=T^T\)对\(t=\tau, \dots, \tau^*\)。 反向平滑的公式变成 \[\begin{aligned} \boldsymbol r_{t-1} =& T_t^T \boldsymbol r_t, & N_{t-1} =& T_t^T N_t T_t, \ t=\tau^*, \dots, \tau . \end{aligned}\]

算法仍是取消高斯分布限制下的最小方差线性无偏估计， 以及贝叶斯估计。 对缺失值的这样的简单处理能力是状态空间模型的优势之一。

如果\(\boldsymbol y_t\)仅有部分分量缺失怎么办？ 设\(\boldsymbol y_t\)是没有缺失的， \(\boldsymbol y_t^*\)有部分分量缺失， 则\(\boldsymbol y_t^* = W_t \boldsymbol y_t\)， 其中\(W_t\)是\(I_p\)的部分行组成。 这时观测方程变成 \[\begin{aligned} \boldsymbol y_t^* = Z_t^* \boldsymbol\alpha_t + \boldsymbol\varepsilon_t^*, \ \boldsymbol\varepsilon_t^* \sim \text{N}(0, H_t^*), \end{aligned}\] 其中\(Z_t^* = W_t Z_t\)， \(\boldsymbol\varepsilon_t^* = W_t \boldsymbol\varepsilon_t\)， \(H_t^* = W_t H_t W_t^T\)。 只要将算法推广到允许\(\boldsymbol y_t\)的维数可变， 就可以按正常的卡尔曼滤波算法进行计算。 还可以按\(\hat{\boldsymbol y}_t = Z_t \hat{\boldsymbol\alpha}_t\)进行缺失值估计。

从滤波、平滑公式可以看出， 实际上\(\boldsymbol y_t\)的维数可以推广到依赖于\(t\)， 这不会改变滤波、平滑公式。

可以将\(\boldsymbol y_t\)转化成每一个分量单独作为一步， 这样就可以统一地处理有部分分量缺失的问题。

30.8 向前预测

只要将未来值作为缺失值向前滤波， 就可以产生最小均方误差的预测。 条件期望是最小均方误差的， 所以若\(\boldsymbol y_1, \dots, \boldsymbol y_n\)已知， 要预测\(\boldsymbol y_{n+j}\), \(j=1,\dots,J\)， 应计算 \[\begin{aligned} E(\boldsymbol y_{n+j} | \boldsymbol Y_n) =& E(Z_{n+j} \boldsymbol\alpha_{n+j} + \boldsymbol\varepsilon_{n+j} | \boldsymbol Y_n) \\ =& Z_{n+j} E(\boldsymbol\alpha_{n+j} | \boldsymbol Y_n) . \end{aligned}\] 易见 \[\begin{aligned} E(\boldsymbol\alpha_{n+1} | \boldsymbol Y_n) =& \boldsymbol a_{n+1} = T_n \boldsymbol a_n, \\ \text{Var}(\boldsymbol\alpha_{n+1} | \boldsymbol Y_n) =& P_{n+1}, \\ E(\boldsymbol\alpha_{n+j+1} | \boldsymbol Y_n) =& E(T_{n+j} \boldsymbol\alpha_{n+j} + R_{n+j} \boldsymbol\eta_{n+j} | \boldsymbol Y_n) \\ =& T_{n+j} E(\boldsymbol\alpha_{n+j} | \boldsymbol Y_n), \ j=1,2,\dots,J-1 . \\ \text{Var}(\boldsymbol\alpha_{n+j+1} | \boldsymbol Y_n) =& T_{n+j} \text{Var}(\boldsymbol\alpha_{n+j} | \boldsymbol Y_n) T_{n+j}^T + R_{n+j} Q_{n+j} R_{n+j}^T . \end{aligned}\] 如此递推可得\(E(\boldsymbol y_{n+j} | \boldsymbol Y_n)\)的值， 以及其均方误差为 \[\begin{aligned} \text{Var}(\boldsymbol y_{n+j} | \boldsymbol Y_n) =& Z_{n+j} \text{Var}(\boldsymbol\alpha_{n+j} | \boldsymbol Y_n) Z_{n+j}^T + H_{n+j} . \end{aligned}\]

这相当于将\(\boldsymbol y_{n+j}\)看成缺失值然后向前进行卡尔曼滤波。 对\(j=1,2,\dots,J\)， 令\(Z_{n+j}=0\), \(K_{n+j}=0\), \(L_{n+j} = T_{n+j}\)。

30.9 初始分布问题

30.9.1 问题介绍

在实际的线性高斯状态空间模型中， \(\boldsymbol\alpha_1\)的初始分布\(\text{N}(\boldsymbol a_1, P_1)\)通常是未知的， 或者部分参数未知。 所以\(\boldsymbol\alpha_1\)的一般模型为 \[\begin{equation} \boldsymbol\alpha_1 = \boldsymbol a + A \boldsymbol\delta + R_0 \boldsymbol\eta_0 , \tag{30.24} \end{equation}\] 其中\(\boldsymbol a\)为\(m \times 1\)已知值， 一般是\(\boldsymbol 0\)。 \(\boldsymbol\eta_0\)是\(m-q\)维随机向量， \(\boldsymbol\eta_0 \sim \text{N}(\boldsymbol 0, Q_0)\)， \(Q_0\)已知， 代表了初始分布中方差已知的成分， \(R_0\)是\(m \times (m-q)\)矩阵， 其列向量来自于\(I_m\)矩阵中的列向量。 \(\boldsymbol\delta\)是\(q \times 1\)未知向量， 可能是非随机的， 可以参与到最大似然估计中； 也可能看成是方差无穷大的正态随机向量。 \(A\)是\(m \times q\)矩阵， 其列向量来自\(I_m\)中的列向量。 要求\(A^T R_0 = 0\)， 且\(A\)的列向量与\(R_0\)的列向量的集合包含\(I_m\)的\(g\)个列向量(\(g \leq m\))。

先考虑\(\boldsymbol\delta\)是方差无穷大的正态随机向量的处理方法。 设\(\boldsymbol\delta \sim \text{N}(\boldsymbol 0, \kappa I_q)\)， 其中\(\kappa \to \infty\)， 称为发散的分布。 这时\(\boldsymbol\alpha_1 \sim \text{N}(\boldsymbol a, P_1)\)， \[\begin{equation} P_1 = \kappa P_\infty + P_*, \tag{30.25} \end{equation}\] 其中\(P_\infty = A A^T\)是对角元素仅取1和0的对角阵， 共有\(q\)个1和\(m-q\)个0， \(P_* = R_0 Q_0 R_0^T\)。 (30.25)这样的初始分布设定称为“发散初始化”。 当\(P_\infty\)的\((i,i)\)元素非零时不妨设\(\boldsymbol a\)的第\(i\)元素为零。

发散先验时进行卡尔曼滤波的一个简化做法是取足够大的\(\kappa\)值， 然后从已知的\(\boldsymbol a\)和\(P_1\)出发进行递推， 但这样会造成较大的舍入误差。 更好的做法是推导当\(\kappa \to \infty\)时按照(30.25)分布的精确递推算法。

推导的思想是将矩阵乘法表示为\(\kappa^{-1}\)的幂级数形式， 令\(\kappa\to\infty\)得到主要部分。 算法中涉及\(F_t^{-1}\)， 对于\(p=1\)情形， 与\(P_\infty\)对应的\(F_t\)部分或者是0， 或者是一个取正值的标量； 对于\(p>1\)的情形， \(F_t\)中相应于\(P_\infty\)的部分有可能不满秩， 我们先假设\(F_t\)中相应于\(P_\infty\)的部分或者满秩， 或者等于0。 将算法转化为\(p=1\)可以避免\(F_t\)不满秩问题。

30.9.2 精确初始化方法

这里用\(O(\kappa^{-j})\)表示函数\(f(\kappa)\)满足\(\lim_{\kappa\to\infty} \kappa^j f(\kappa)\)存在有限。

当\(P_1\)有(30.25)那样的分解时， 可以证明 \[\begin{equation} P_t = \kappa P_{\infty,t} + P_{*,t} + O(\kappa^{-1}), \ t=2,\dots,n, \tag{30.26} \end{equation}\] 其中\(P_{\infty,t}\)和\(P_{*,t}\)不依赖于\(\kappa\)。 一般情况下存在\(d\)，\(d\)远小于\(n\)， 当\(t>d\)时\(P_{\infty,t}=0\), 对\(t=d+1,\dots,n\)可以取\(P_t = P_{*,t}\)， 使用正常的向前滤波算法即可， 这样就消除了发散先验的影响。 如果初始分布参数\(\boldsymbol a_1\)和\(P_1\)完全已知， 则\(P_{\infty,t}=0\), \(t=1,\dots, n\)， 即\(d=0\)， 不需要使用发散先验。

30.9.2.1 初始化阶段

对\(t=1,\dots,d\)， \(P_{\infty,t} \neq 0\)， \(P_t\)的方差值中存在\(\infty\)， 但我们可以正常地进行卡尔曼滤波递推， 关键是在求\(F_t^{-1}\)时， 可以将\(\kappa^{-1}\)项令\(\kappa\to\infty\)从而忽略掉。

取\(P_{\infty,1} = P_\infty = A A^T\)， 这是一个主对角元素仅取1和0的对角阵。 取\(P_{*,1} = P_* = R_0 Q_0 R_0^T\)。

由(30.26)以及\(F_t = Z_t P_t Z_t^T + H_t\)可知\(F_t\)也有分解 \[\begin{equation} F_t = \kappa F_{\infty,t} + F_{*,t} + O(\kappa^{-1}) . \tag{30.27} \end{equation}\] 对比(30.26)可得 \[\begin{aligned} F_{\infty,t} =& Z_t P_{\infty,t} Z_t^T, & F_{*,t} =& Z_t P_{*,t} Z_t^T + H_t . \end{aligned}\] 令\(M_t = P_t Z_t^T\)， 则 \[\begin{aligned} M_t =& \kappa M_{\infty,t} + M_{*,t} + O(\kappa^{-1}), \\ M_{\infty,t} =& P_{\infty,t} Z_t^T, & M_{*,t} =& P_{*,t} Z_t^T . \end{aligned}\]

注意\(M_{\infty,t} = 0 \Longleftrightarrow F_{\infty,t} = 0\)， \(P_{\infty,t} = 0 \Longrightarrow M_{\infty,t} = 0\)。 一旦\(P_{\infty,t} = 0\)， 则令\(d=t-1\)，对\(t=d+1, \dots, n\)， 就不需要再考虑发散先验部分， 令\(P_t = P_{*,t}\)按照正常的卡尔曼滤波向前递推即可。

在\(t=1,\dots,d\)步， 需要计算\(F_t^{-1}\)， 要考虑三种情况：

  • \(F_t = 0\);
  • \(F_t\)正定；
  • \(F_t \neq 0\)但不满秩。

我们仅处理前两种情况。 当观测为一元时间序列(\(p=1\))时， 只有前两种情况； 当观测为多元时， 一般也较少遇到第三种情况， 一旦遇到， 可以采取将多元观测转化为一元观测的办法。

考虑\(F_t^{-1}\)求解。 对\(F_t^{-1}\)的(30.27)关于\(\kappa^{-1}\)做泰勒展开， 得 \[\begin{aligned} F_t^{-1} =& [ \kappa F_{\infty,t} + F_{*,t} + O(\kappa^{-1})]^{-1} \\ =& F_t^{(0)} + \kappa^{-1} F_t^{(1)} + \kappa^{-2} F_t^{(2)} + O(\kappa^{-3}). \end{aligned}\] 用待定系数法求解， 当\(F_{\infty,t}\)正定时，可解得 \[\begin{aligned} F_{\infty,t} =& Z_t P_{\infty,t} Z_t^T, & F_{*,t} =& Z_t P_{*,t} Z_t^T + H_t , \\ M_{\infty,t} =& P_{\infty,t} Z_t^T, & M_{*,t} =& P_{*,t} Z_t^T ,\\ F_t^{(0)} =& 0, \qquad F_t^{(1)} = F_{\infty,t}^{-1}, & F_t^{(2)} =& -F_{\infty,t}^{-1} F_{*,t} F_{\infty,t}^{-1} . \\ K_t =& K_t^{(0)} + \kappa^{-1} K_t^{(1)} + O(\kappa^{-2}), \\ K_t^{(0)} =& T_t M_{\infty,t} F_t^{(1)}, & K_t^{(1)} =& T_t M_{*,t} F_t^{(1)} + T_t M_{\infty,t} F_t^{(2)}, \\ L_t =& L_t^{(0)} + \kappa^{-1} L_t^{(1)} + O(\kappa^{-2}), \\ L_t^{(0)} =& T_t - K_t^{(0)} Z_t, & L_t^{(1)} =& -K_t^{(1)} Z_t . \end{aligned}\]

如果\(F_{\infty,t}=0\)， 则\(M_{\infty,t}=0\)， \[\begin{aligned} F_{\infty,t} =& 0, & F_{*,t} =& Z_t P_{*,t} Z_t^T + H_t , \\ M_{\infty,t} =& 0, & M_{*,t} =& P_{*,t} Z_t^T ,\\ F_t =& F_{*,t} + O(\kappa^{-1}), & M_t =& M_{*,t} + O(\kappa^{-1}), \\ F_t^{-1} =& F_{*,t}^{-1} + O(\kappa^{-1}), \\ F_t^{(0)} =& F_{*,t}^{-1}, & F_t^{(1)} =& F_t^{(2)} = 0, \\ K_t =& T_t M_{*,t} F_{*,t}^{-1} + O(\kappa^{-1}), \\ K_t^{(0)} =& T_t M_{*,t} F_{*,t}^{-1}, & K_t^{(1)} =& 0, \\ L_t =& T_t - K_t Z_t, \\ L_t^{(0)} =& T_t - K_t^{(0)} Z_t , & L_t^{(1)} =& 0 . \end{aligned}\]

关于\(\boldsymbol a_t\)的递推， 易见 \[ \boldsymbol a_t = \boldsymbol a_t^{(0)} + \kappa^{-1} \boldsymbol a_t^{(1)} + O(\kappa^{-2}) , \] 递推时令\(\kappa\to\infty\)， 所以不必计算\(\boldsymbol a_t^{(1)}\)。 取初值\(\boldsymbol a_1^{(0)} = \boldsymbol a\), \(\boldsymbol a_1^{(1)} = \boldsymbol 0\)。 类似有 \[ \boldsymbol v_t = \boldsymbol v_t^{(0)} + \kappa^{-1} \boldsymbol v_t^{(1)} + O(\kappa^{-2}) , \] 递推时不必计算\(\boldsymbol v_t^{(1)}\)， 取\(\boldsymbol v_t^{(0)} = \boldsymbol y_t - Z_t \boldsymbol a_t^{(0)}\)， \(\boldsymbol v_t^{(1)} = -Z_t \boldsymbol a_t^{(1)}\)。 对\(t=1,\dots,d\)只要计算 \[\begin{aligned} \boldsymbol v_t^{(0)} =& \boldsymbol y_t - Z_t \boldsymbol a_t^{(0)}, \\ \boldsymbol a_{t+1}^{(0)} =& T_t \boldsymbol a_t^{(0)} + K_t^{(0)} \boldsymbol v_t^{(0)} . \end{aligned}\]

关于\(P_{\infty,t}\)和\(P_{*,t}\), 可得 \[\begin{aligned} P_{\infty,t+1} =& T_t P_{\infty,t} [L_t^{(0)}]^T, \\ P_{*,t+1} =& T_t P_{\infty,t} [L_t^{(1)}]^T + T_t P_{*,t} [L_t^{(0)}]^T + R_t Q_t R_t^T . \end{aligned}\] 当\(F_t=0\)时上面的\(L_t^{(1)}=0\)。

30.9.2.2 初始化以后的衔接

考虑\(t=d+1, \dots, n\)的向前递推计算。 可以证明， 在模型设置合理的情况下， 存在\(d \geq 0\)， 使得\(P_{\infty,t} \neq 0\), \(t \leq d\); \(P_{\infty,t} = 0\), \(t > d\)。 证明见(Durbin and Koopman 2012) P.129节5.2.2。

考虑(30.24)中\(\boldsymbol\delta\)， 这是\(q\)维的发散先验部分， 当\(t > d\)时\(P(\boldsymbol\delta | \boldsymbol Y_t)\)为有限值， 从而\(P_t\)有限。 对\(t=d+1, \dots, n\)， 只要取\(\boldsymbol a_{d+1} = \boldsymbol a_{d+1}^{(0)}\), \(P_{d+1} = P_{*,d+1}\)， 就可以按已知初始分布时的卡尔曼向前递推公式进行递推计算了。

30.9.2.3 初始化阶段方差阵简化公式

对\(t=1,\dots,d\)， 可以将\(P_{*,t}\)和\(P_{\infty,t}\)的递推写成一个统一的公式。 记 \[\begin{aligned} P_t^\dagger = \begin{bmatrix} P_{*,t} & P_{\infty,t} \end{bmatrix}, \quad L_t^\dagger = \begin{bmatrix} L_t^{(0)} & L_t^{(1)} \\ 0 & L_t^{(0)} \end{bmatrix}, \end{aligned}\] 取初值\(P_1^\dagger = \begin{bmatrix} P_{*} & P_{\infty} \end{bmatrix}\)， 对\(t=1,\dots,d\)有 \[\begin{aligned} P_{t+1}^\dagger = T_t P_t^\dagger [L_t^\dagger]^T + \begin{bmatrix} R_t Q_t R_t^T & 0 \end{bmatrix} . \end{aligned}\] 这与非发散先验时的公式形式相同。

若某个\(F_{\infty,t}=0\)， 则公式中 \[\begin{aligned} K_t^{(0)} =& T_t M_{*,t} F_{*,t}^{-1}, & L_t^{(0)} =& T_t - K_t^{(0)} Z_t, \quad L_t^{(1)} = 0 . \end{aligned}\]

30.9.3 精确初始化平滑算法

对\(t=n,n-1,\dots,d+1\)， 仍可以按原来的反向平滑算法进行平滑分布计算。

30.9.3.1 平滑的均值

对\(t=d,\dots,1\)， 考虑 \[\begin{aligned} \boldsymbol r_{t-1} =& Z_t^T F_t^{-1} \boldsymbol v_t + L_t^T \boldsymbol r_t . \end{aligned}\] 其中的\(F_t^{-1}\)和\(\boldsymbol v_t\), \(L_t\)都可以按\(\kappa^{-1}\)泰勒展开近似， 所以也有 \[\begin{aligned} \boldsymbol r_{t-1} =& \boldsymbol r_{t-1}^{(0)} + \kappa^{-1} \boldsymbol r_{t-1}^{(1)} + O(\kappa^{-2}) , \quad t=d,\dots,1 . \end{aligned}\] 用待定系数法可得， 当\(F_{\infty,t}\)正定时， \[\begin{aligned} \boldsymbol r_{t-1}^{(0)} =& [L_t^{(0)}]^T \boldsymbol r_{t}^{(0)}, \\ \boldsymbol r_{t-1}^{(1)} =& Z_t^T F_t^{(1)} \boldsymbol v_{t}^{(0)} + [L_t^{(0)}]^T \boldsymbol r_{t}^{(1)} + [L_t^{(1)}]^T \boldsymbol r_{t}^{(0)}, \quad t=d, \dots, 1 . \\ \boldsymbol r_{d}^{(0)} =& \boldsymbol r_{d}, \quad \boldsymbol r_{d}^{(1)} = \boldsymbol 0 . \end{aligned}\] 可以证明\(\text{Var}(\boldsymbol\alpha_t | \boldsymbol Y_n)\)有限， 从而由待定系数法可得 \[\begin{aligned} \hat{\boldsymbol\alpha}_t =E(\boldsymbol\alpha_t | \boldsymbol Y_n) = \boldsymbol a_t^{(0)} + P_{*,t} \boldsymbol r_{t-1}^{(0)} + P_{\infty,t} \boldsymbol r_{t-1}^{(1)}, \quad t=d,\dots,1 . \end{aligned}\] 可以合并写成 \[\begin{aligned} \boldsymbol r_{t-1}^\dagger =& \begin{pmatrix} \boldsymbol r_{t-1}^{(0)} \\ \boldsymbol r_{t-1}^{(1)} \end{pmatrix} = \begin{pmatrix} 0 \\ Z_t^T F_t^{(1)} \boldsymbol v_{t}^{(0)} \end{pmatrix} + [L_t^\dagger]^T \boldsymbol r_t^\dagger, \\ \hat{\boldsymbol\alpha}_t =& \boldsymbol a_t^{(0)} + P_t^\dagger \boldsymbol r_{t-1}^\dagger, \quad t=d,\dots, 1. \end{aligned}\] 其中初值\(\boldsymbol r_{d}^\dagger = \begin{pmatrix} \boldsymbol r_{d} \\ \boldsymbol 0 \end{pmatrix}\)， 这个递推形式与非发散先验时的公式形式相同。

这里略过\(P_{\infty,t} \neq 0\)但\(F_{\infty,t}=0\)的讨论， 这种情况比较少见。

30.9.3.2 平滑的方差阵

原来 \[\begin{aligned} V_t =& \text{Var}(\boldsymbol\alpha_t | \boldsymbol Y_n) = P_t - P_t N_{t-1} P_t, \\ N_n =& 0, \\ N_{t-1} =& Z_t^T F_t^{-1} Z_t + L_t^T N_t L_t . \end{aligned}\] 若\(F_{\infty,t}\)正定，则有 \[\begin{aligned} N_t =& N_t^{(0)} + \kappa^{-1} N_t^{(1)} + \kappa^{-2} N_t^{(2)} + O(\kappa^{-3}), \end{aligned}\] 用待定系数法可得递推公式为\(t=d,\dots,1\)时 \[\begin{aligned} N_d^{(0)} =& N_d, \quad N_d^{(1)} = N_d^{(2)} = 0, \\ N_{t-1}^{(0)} =& [L_t^{(0)}]^T N_t^{(0)} L_t^{(0)}, \\ N_{t-1}^{(1)} =& Z_t F_t^{(1)} Z_t + [L_t^{(0)}]^T N_t^{(1)} L_t^{(0)} \\ & + [L_t^{(1)}]^T N_t^{(0)} L_t^{(0)} + [L_t^{(0)}]^T N_t^{(0)} L_t^{(1)}, \\ N_{t-1}^{(2)} =& Z_t^T F_t^{(2)} Z_t + [L_t^{(0)}]^T N_t^{(2)} L_t^{(0)} \\ & + [L_t^{(0)}]^T N_t^{(1)} L_t^{(1)} + [L_t^{(1)}]^T N_t^{(1)} L_t^{(0)} \\ & + [L_t^{(0)}]^T N_t^{(0)} L_t^{(2)} + [L_t^{(2)}]^T N_t^{(0)} L_t^{(0)} \\ & + [L_t^{(1)}]^T N_t^{(0)} L_t^{(1)} . \end{aligned}\] 可以证明\(V_t\)都有限， 所以由待定系数法可得 \[\begin{aligned} V_t =& P_{*,t} - P_{*,t} N_{t-1}^{(0)} P_{*,t} \\ & - P_{*,t} N_{t-1}^{(1)} P_{\infty,t} - P_{\infty,t} N_{t-1}^{(1)} P_{*,t} \\ & - P_{\infty,t} N_{t-1}^{(2)} P_{\infty,t} . \end{aligned}\]

30.9.3.3 扰动项平滑

内容略过。

30.9.3.4 平滑抽样

内容略过。

30.9.4 精确初始化示例

30.9.4.1 结构时间序列模型

考虑局部趋势模型： \[\begin{aligned} y_t =& \mu_t + e_t, \\ \mu_{t+1} =& \mu_t + \nu_t + \xi_t, \\ \nu_{t+1} =& \nu_t + \zeta_t, \end{aligned}\] 写成状态空间模型，为 \[\begin{aligned} y_t =& (1\ 0) \begin{pmatrix} \mu_t \\ \nu_t \end{pmatrix} + e_t, \\ \begin{pmatrix} \mu_{t+1} \\ \nu_{t+1} \end{pmatrix} =& \begin{pmatrix} 1 & 1 \\ 0 & 1 \end{pmatrix} \begin{pmatrix} \mu_t \\ \nu_t \end{pmatrix} + \begin{pmatrix} \xi_t \\ \zeta_t \end{pmatrix} . \end{aligned}\] 各矩阵 \[\begin{aligned} Z_t =& (1, 0), \quad & T_t =& \begin{pmatrix} 1 & 1 \\ 0 & 1 \end{pmatrix} , \\ H_t =& \sigma_{\varepsilon}^2, \quad R_t = I_2, & Q_t =& \sigma_{\varepsilon}^2 \begin{pmatrix} q_\xi & 0 \\ 0 & q_\zeta \end{pmatrix} . \end{aligned}\]

\(t=1\)时\(\mu_1\)和\(\nu_1\)都是未知的， 作为\(\boldsymbol\delta\)，从而初值 \[\begin{aligned} \boldsymbol a_1 = \boldsymbol a_1^{(0)} = \boldsymbol 0, \quad P_{*,1} = 0, \quad P_{\infty,1} = I_2 . \end{aligned}\]

计算第一次更新： \[\begin{aligned} F_{\infty,1} =& Z_1 P_{\infty,1} Z_1^T = 1, & F_{*,1} =& Z_1 P_{*,1} Z_1^T + H_1 = \sigma_{\varepsilon}^2 , \\ M_{\infty,1} =& P_{\infty,1} Z_1^T = \begin{pmatrix} 1 \\ 0 \end{pmatrix}, & M_{*,1} =& P_{*,1} Z_1^T = \begin{pmatrix} 0 \\ 0 \end{pmatrix}, \\ F_1^{(0)} =& 0, \quad F_1^{(1)} = F_{\infty,1}^{-1} = 1, & F_1^{(2)} =& -F_{\infty,1}^{-1} F_{*,1} F_{\infty,1}^{-1} = -\sigma_{\varepsilon}^2, \\ K_1^{(0)} =& T_1 M_{\infty,1} F_1^{(1)} = \begin{pmatrix} 1 \\ 0 \end{pmatrix}, & K_1^{(1)} =& T_1 M_{*,1} F_1^{(1)} + T_1 M_{\infty,1} F_1^{(2)} \\ && =& -\sigma_{\varepsilon}^2 \begin{pmatrix} 1 \\ 0 \end{pmatrix}, \\ L_1^{(0)} =& T_1 - K_1^{(0)} Z_1 = \begin{pmatrix} 0 & 1 \\ 0 & 1 \end{pmatrix}, & L_1^{(1)} =& -K_1^{(1)} Z_1 = \sigma_{\varepsilon}^2 \begin{pmatrix} 1 & 0 \\ 0 & 0 \end{pmatrix}, \\ \boldsymbol v_1^{(0)} =& y_1 - Z_1 \boldsymbol a_1^{(0)} = y_1, \end{aligned}\] 进而 \[\begin{aligned} \boldsymbol a_2 =& \boldsymbol a_2^{(0)} = T_1 \boldsymbol a_1^{(0)} + K_1^{(0)} \boldsymbol v_1^{(0)} = \begin{pmatrix} y_1 \\ 0 \end{pmatrix}, \\ P_{*,2} =& T_1 P_{\infty,1} [L_1^{(1)}]^T + T_1 P_{*,1} [L_1^{(0)}]^T + R_t Q_t R_t \\ =& \sigma_{\varepsilon}^2 \begin{pmatrix} 1+q_{\xi} & 0 \\ 0 & q_{\zeta} \end{pmatrix}, \\ P_{\infty,2} =& T_1 P_{\infty,1} [L_1^{(0)}]^T = \begin{pmatrix} 1 & 1 \\ 1 & 1 \end{pmatrix} . \end{aligned}\]

继续下一轮更新可得 \[\begin{aligned} K_2^{(0)} =& \begin{pmatrix} 2 \\ 1 \end{pmatrix}, & K_2^{(1)} =& -\sigma_{\varepsilon}^2 \begin{pmatrix} 3+q_\xi \\ 2+q_\xi \end{pmatrix}, \\ L_2^{(0)} =& \begin{pmatrix} -1 & 1 \\ -1 & 1 \end{pmatrix}, & L_2^{(1)} =& \sigma_{\varepsilon}^2 \begin{pmatrix} 3+q_\xi & 0 \\ 2+q_\xi & 0 \end{pmatrix}, \end{aligned}\] 进而 \[\begin{aligned} \boldsymbol a_3 =& \begin{pmatrix} 2 y_2 - y_1 \\ y_2 - y_1 \end{pmatrix}, \\ P_{*,3} =& \sigma_{\varepsilon}^2 \begin{pmatrix} 5 + 2 q_{\xi} + q_{\zeta} & 3 + q_{\xi} + q_{\zeta} \\ 3 + q_{\xi} + q_{\zeta} & 2 + q_{\xi} + 2 q_{\zeta} \end{pmatrix}, \\ P_{\infty,3} =& \begin{pmatrix} 0 & 0 \\ 0 & 0 \end{pmatrix}, \end{aligned}\] 对\(t=3,\dots,n\)可以使用正常的卡尔曼向前滤波更新算法。

30.10 观测序列一元化

将\(p\)维的时间序列的状态空间模型设法改写成一元观测值的状态空间模型， 可以避免考虑精确发散先验中\(F_t\)不等于零也不满秩的情况， 同时可以提高计算效率， 也自然地容许不同时刻的观测值分量个数可变。

30.10.1 各分量条件独立的情形

设\(H_t\)是对角阵， 从而给定\(\boldsymbol\alpha_t\)时\(\boldsymbol y_t\)的各分量独立， 这样， 可以将各个分量看成是多个观测， 而不是一个观测的多个分量。 设\(\boldsymbol y_t = (y_{t,1}, \dots, y_{t,p_t})^T\)， 是\(p_t \times 1\)向量， 将所有观测值按如下次序排列成一元的观测值序列： \[ y_{1,1}, \dots, y_{1,p_1}, y_{2,1}, \dots, y_{n,p_n} . \] 仍设\(\boldsymbol\alpha_1 \sim \text{N}(\boldsymbol a_1, P_1)\)。 记\(\boldsymbol\varepsilon_t = (\varepsilon_{t,1}, \dots, \varepsilon_{t,p_t})^T\)， \[\begin{aligned} Z_t = \begin{pmatrix} Z_{t,1} \\ \vdots \\ Z_{t,p_t} \end{pmatrix}, \qquad H_t = \text{diag}(\sigma_{t,1}^2, \dots, \sigma_{t,p_t}^2), \end{aligned}\] 其中\(Z_{t,i}\)是\(1 \times m\)行向量。 则原来向量观测值的模型可以转换成如下的一元观测值的模型： \[\begin{aligned} y_{t,i} =& Z_{t,i} \boldsymbol\alpha_{t,i} + \varepsilon_{t,i}, \quad \varepsilon_{t,i} \sim \text{N}(0, \sigma_{t,i}^2), \quad i=1,\dots,p_t, \\ \boldsymbol\alpha_{t,i+1} =& \boldsymbol\alpha_{t,i}, \quad i=1,\dots, p_t-1, \\ \boldsymbol\alpha_{t+1,1} =& T_t \boldsymbol\alpha_{t,p_t} + R_t \boldsymbol\eta_t, \quad t=1,\dots,n . \end{aligned}\] 则\(\boldsymbol\alpha_{1,1} = \boldsymbol\alpha_1 \sim \text{N}(\boldsymbol a_1, P_1)\)， \(\boldsymbol\alpha_t = \boldsymbol\alpha_{t,1}\)， 令 \[\begin{aligned} \boldsymbol a_{t,1} =& E(\boldsymbol\alpha_{t,1} | \boldsymbol Y_{t-1}) = E(\boldsymbol\alpha_{t} | \boldsymbol Y_{t-1}) = \boldsymbol a_t, \\ P_{t,1} =& \text{Var}(\boldsymbol\alpha_{t,1} | \boldsymbol Y_{t-1}) = \text{Var}(\boldsymbol\alpha_{t} | \boldsymbol Y_{t-1}) = P_t, \\ \boldsymbol a_{t,i} =& E(\boldsymbol\alpha_{t,1} | \boldsymbol Y_{t-1}, y_{t,1}, \dots, y_{t,i-1}), \\ P_{t,i} =& \text{Var}(\boldsymbol\alpha_{t,1} | \boldsymbol Y_{t-1}, y_{t,1}, \dots, y_{t,i-1}), \quad i=2,\dots,p_t, \ t=1,\dots,n . \end{aligned}\]

向前滤波的计算公式分为在同一时刻不同分量的更新， 以及不同时刻的更新。

同一时刻的更新公式： \[\begin{aligned} v_{t,i} =& y_{t,i} - Z_{t,i} \boldsymbol a_{t,i}, \\ F_{t,i} =& Z_{t,i} P_{t,i} Z_{t,i}^T + \sigma_{t,i}^2, \\ K_{t,i} =& P_{t,i} Z_{t,i}^T F_{t,i}^{-1}, \\ \boldsymbol a_{t,i+1} =& \boldsymbol a_{t,i} + K_{t,i} v_{t,i}, \\ P_{t,i+1} =& P_{t,i} - K_{t,i} F_{t,i} K_{t,i}^T, \quad i=1,\dots,p_t, \ t=1,\dots, n . \end{aligned}\]

对每一个\(t\)， 在\(i=1,\dots,p_t\)的循环之前， \(\boldsymbol a_{t,1}\)保存了\(\boldsymbol a_t\)的值， \(P_{t,1}\)保存了\(P_t\)的值； 在\(i=1,\dots,p_t\)的循环之后， \(\boldsymbol a_{t,p_t+1}\)保存了\(\boldsymbol a_{t|t}\)的值， \(P_{t,p_t+1}\)保存了\(\Sigma_{t|t}\)的值， 即滤波分布\(p(\boldsymbol\alpha_t | \boldsymbol Y_t)\)。 \(i=1,\dots,p_t\)的循环起到了利用\(\boldsymbol y_t\)更新预报分布\(p(\boldsymbol\alpha_t | \boldsymbol Y_{t-1})\)的作用。

从\(t\)到\(t+1\)的更新公式为 \[\begin{aligned} \boldsymbol a_{t+1,1} =& T_t \boldsymbol a_{t, p_t+1}, \\ P_{t+1,1} =& T_t P_{t, p_t+1} T_t^T + R_t Q_t R_t^T, \ t=1,\dots,n . \end{aligned}\] 这一更新将滤波分布\(p(\boldsymbol\alpha_t | \boldsymbol Y_t)\)更新为预报分布\(p(\boldsymbol\alpha_{t+1} | \boldsymbol Y_t)\)， 这相当于\(\boldsymbol y_t\)缺失情况下的一步更新， 而\(\boldsymbol y_t\)的信息已经在前面关于\(i=1,\dots,p_t\)的循环利用过了。

注意\(F_{t,i}\)是标量， 在同时刻的更新公式中用到了\(F_{t,i}^{-1}\)， 如果\(F_{t,i}=0\)， 注意\(F_{t,i}\)是用\(\boldsymbol Y_{t-1}, y_{t,1}, \dots, y_{t,i-1}\)预测\(y_{t,i}\)的均方误差， \(F_{t,i}=0\)说明\(y_{t,i}\)是\(\boldsymbol Y_{t-1}, y_{t,1}, \dots, y_{t,i-1}\)的线性组合， 没有提供额外的信息， 所以\(F_{t,i}=0\)时递推公式只要变成 \[\begin{aligned} \boldsymbol a_{t,i+1} =& E(\boldsymbol\alpha_{t,i+1} | \boldsymbol Y_{t-1}, y_{t,1}, \dots, y_{t,i}) \\ =& E(\boldsymbol\alpha_{t,i+1} | \boldsymbol Y_{t-1}, y_{t,1}, \dots, y_{t,i-1}) \\ =& \boldsymbol a_{t,i}, \\ P_{t,i+1} =& P_{t,i} . \end{aligned}\]

对于反向递推平滑， 取\(\boldsymbol r_{n,p_n} = \boldsymbol 0\), \(N_{n,p_n} = 0\)， 同时刻反向递推 \[\begin{aligned} L_{t,i} =& I_m - K_{t,i} Z_{t,i}, \\ \boldsymbol r_{t,i-1} =& Z_{t,i}^T F_{t,i}^{-1} v_{t,i} + L_{t,i}^T \boldsymbol r_{t,i}, \\ N_{t,i-1} =& Z_{t,i}^T F_{t,i}^{-1} Z_{t,i} + L_{t,i}^T N_{t,i} L_{t,i}, \quad i=p_t, \dots, 1, \end{aligned}\] 从\(t\)到\(t-1\)递归为 \[\begin{aligned} \boldsymbol r_{t-1,p_{t-1}} =& T_{t-1}^T \boldsymbol r_{t,0}, \\ N_{t-1,p_{t-1}} =& T_{t-1}^T N_{t,0} T_{t-1}, \end{aligned}\] 这时 \[\begin{aligned} \boldsymbol r_{t,0} =& \boldsymbol r_{t-1}, & N_{t,0} =& N_{t-1}, \\ \boldsymbol a_t =& \boldsymbol a_{t,1}, & P_t =& P_{t,1}, \\ \hat{\boldsymbol\alpha}_t =& \boldsymbol a_t + P_t \boldsymbol r_{t-1}, & V_t =& P_t - P_t N_{t-1} P_t, \\ & t=n, \dots, 1 . \end{aligned}\]

30.10.2 各分量非条件独立的情形

若\(H_t\)不是对角阵， 上面的方法无法直接使用。 一种办法是将\(\boldsymbol\varepsilon_t\)整合进入状态向量\(\boldsymbol\alpha_t\)中， 这样观测方程误差变成了0， 符合方差阵对角阵要求。

另一种办法是对每个时刻\(t\)将\(H_t\)作LDL分解， 然后将\(\boldsymbol y_t\)和\(Z_t\)作变换， 使得新的观测变成方差阵对角形式。

30.11 参数的最大似然估计

状态空间模型中的参数， 主要是\(H_t\)和\(Q_t\)中的参数， 比如局部水平模型中的\(\sigma_e^2\), \(\sigma_\eta^2\)。 但是其它矩阵中也可以含有参数。 有时称这些参数为“超参数”。

因为观测值的联合分布是多元正态分布， 向前的卡尔曼滤波实际上是将很大的方差阵求逆问题， 转换成了对角阵求逆问题， 从而提高了运算效率， 确保了数值稳定性。

30.11.1 初始分布已知时的似然函数计算

设初始分布\(\boldsymbol\alpha_1 \sim \text{N}(\boldsymbol a_1, P_1)\), \(\boldsymbol a_1\), \(P_1\)已知。 似然函数为 \[ L(\boldsymbol Y_n) = p(\boldsymbol y_1) \prod_{t=2}^n p(\boldsymbol y_t | \boldsymbol Y_{t-1}) . \] 对数似然函数为 \[ \log L(\boldsymbol Y_n) = \sum_{t=1}^n \log p(\boldsymbol y_t | \boldsymbol Y_{t-1}) . \] 其中\(p(\boldsymbol y_1 | \boldsymbol Y_0)\)表示\(p(\boldsymbol y_1)\)。 注意这些条件分布都是多元正态分布， \(E(\boldsymbol y_t | \boldsymbol Y_{t-1}) = Z_t \boldsymbol a_t\)， \(\text{Var}(\boldsymbol y_t | \boldsymbol Y_{t-1}) = F_t\)， 所以 \[\begin{aligned} \log L(\boldsymbol Y_n) =& \sum_{t=1}^n \log p(\boldsymbol v_t) \\ =& -\frac{1}{2} np \log(2\pi) - \frac{1}{2} \sum_{t=1}^n \left\{ \log|F_t| + \boldsymbol v_t^T F_t^{-1} \boldsymbol v_t \right\} . \end{aligned}\] 若\(p>1\)， 这个公式需要所有\(F_t\)满秩。 如果\(p=1\)或者转化成一元观测值， 则\(F_t=0\)时只要忽略相应项。

30.11.2 发散初始化的似然函数计算

设对\(t=1,\dots,d\)，\(P_{\infty,t} \neq 0\)， \(t=d+1,\dots,n\)时\(P_{\infty,t} = 0\)。 如果沿用原来的似然函数， 会出现\(-\frac{1}{2} q \log\kappa\)这样的趋于负无穷的项， 所以似然函数应该改用发散似然函数 \[ \log L_d(\boldsymbol Y_n) = -\frac{1}{2} n p \log(2\pi) - \frac{1}{2} \sum_{t=1}^d w_t - \frac{1}{2} \sum_{t=d+1}^n \left\{ \log|F_t| + \boldsymbol v_t^T F_t^{-1} \boldsymbol v_t \right\} . \] 其中 \[ w_t = \begin{cases} \log|F_{\infty,t}|, & F_{\infty,t} > 0, \\ \log|F_{*,t}| + [\boldsymbol v_t^{(0)}]^T F_{*,t}^{-1} \boldsymbol v_t^{(0)}, & F_{\infty,t} = 0 . \end{cases} \]

30.11.3 一元化观测值的似然函数计算

一元化观测值以后的卡尔曼向前滤波算法产生\(v_{t,i}\)和\(F_{t,i}\)， 这时似然函数为

\[\begin{aligned} L(\boldsymbol Y_n) =& p(y_{1,1}) p(y_{1,2} | y_{1,1}) \dots p(y_{1,p_1} | y_{1,1}, \dots, y_{1,p_1-1}) \\ & p(p_{2,1} | y_{1,1}, \dots, y_{1,p_1}) \dots p(y_{n,n_p} | \boldsymbol Y_{n-1}, y_{n,1}, \dots, y_{n,n_p-1}) \\ =& p(v_{1,1}) p(v_{1,2}) \dots p(v_{1,p_1}) p(v_{2,1}) \dots p(v_{n,n_p}), \end{aligned}\] 于是 \[\begin{aligned} \log L(\boldsymbol Y_n) =& -\frac{1}{2} \sum_{t=1}^n \sum_{i=1}^{p_t} \iota_{t,i} \left[ \log(2\pi) + \log F_{t,i} + v_{t,i}^2 / F_{t,i} \right] . \end{aligned}\] 其中 \[\begin{aligned} \iota_{t,i} =& \begin{cases} 1, & F_{t,i} > 0, \\ 0, & F_{t,i} = 0 . \end{cases} \end{aligned}\]

如果使用发散先验， 对\(t=1,\dots,d\)， 计算\(v_{t,i}^{(0)}\), \(F_{\infty,t,i}\), \(F_{*,t,i}\)； 对\(t=d+1, \dots, n\)计算\(v_{t,i}\), \(F_{t,i}\)， 发散对数似然函数为 \[\begin{aligned} \log L(\boldsymbol Y_n) =& -\frac{1}{2} \log(2\pi) \sum_{t=1}^n \sum_{i=1}^{p_t} \iota_{t,i} -\frac{1}{2} \sum_{t=1}^d \sum_{i=1}^{p_t} w_{t,i} \\ & - \frac{1}{2} \sum_{t=d+1}^n \sum_{i=1}^{p_t} \iota_{t,i} \left[ \log F_{t,i} + v_{t,i}^2 / F_{t,i} \right] . \end{aligned}\] 其中 \[\begin{aligned} \iota_{t,i} =& \begin{cases} 1, & 1 \leq t \leq d \text{且} F_{*,t,i} > 0, \text{或} t>d \text{且} F_{t,i}>0 \\ 0, & \text{其它} . \end{cases} \\ w_{t,i} =& \begin{cases} \log|F_{\infty,t,i}|, & F_{\infty,t,i} > 0, \\ \log|F_{*,t,i}| + [\boldsymbol v_{t,i}^{(0)}]^2 / F_{*,t,i}, & F_{\infty,t} = 0 . \end{cases} \end{aligned}\]

30.12 滤波平滑的R程序

30.12.1 初始分布已知情形

设初始分布\(N(\boldsymbol a_1, P_1)\)已知， 模型矩阵\(Z_t, H_t, T_t, R_t, Q_t\)都不依赖于\(t\)。

30.12.1.1 似然函数计算

# 1.1 似然函数计算 ####
## 仅计算对数似然函数
kf_const_logL <- function(
    y,      # 输入观测值向量（长度n），或者$n \times p$矩阵
    Zmat,   # 观测$Z$矩阵，$p \times m$
    Tmat,   # 转移$T$矩阵，$m \times m$
    Hmat,   # 观测误差方差阵，$p \times p$，$p=1$时为标量
    Rmat = NULL,   # 状态方程误差变换矩阵，$m \times r$
    Qmat,   # 状态方程误差方差阵，$r \times r$
    a1,     # $t=1$时状态初始分布均值
    P1      # $t=1$时状态初始分布方差阵
    ){
  if(!is.matrix(y)){
    y <- cbind(y)
  }
  # 观测长度：
  n <- nrow(y)
  # 观测维数：
  p <- ncol(y)
  # 状态维数：
  m <- ncol(Zmat)
  # 状态方程误差维数：
  dim_err_sta <- nrow(Qmat)
  if(is.null(Rmat)) Rmat <- diag(m)
  
  Zt <- Zmat
  Tt <- Tmat
  Ht <- Hmat
  Rt <- Rmat
  Qt <- Qmat
  
  # 单个$v_t$向量
  vt <- numeric(p)
  # 单个$F_t$矩阵，$p \times p$
  Ft <- matrix(0, p, p)
  # 保存$F_t^{-1}$:
  Ftinv <- matrix(0, p, p)
  
  # 单个$a_t$向量
  at <- a1
  # 单个$P_t$矩阵
  Pt <- P1
  
  # 单个$K_t$矩阵，$m \times p$
  Kt <- matrix(0, m, p)
  # 单个$L_t$矩阵，$m \times m$
  Lt <- matrix(0, m, m)
  # t(Z_t)
  Zt_transp <- t(Zt)
  # $R_t Q_t R_t^T$
  RQRt <- Rt %*% Qt %*% t(Rt)
  
  # 进行卡尔曼滤波并累计计算对数似然函数值
  logL <- -0.5*n*p*log(2*pi)
  for(t in 1:n){
    vt[] <- y[t,] - Zt %*% at
    Ft[] <- Zt %*% Pt %*% Zt_transp + Ht
    if(p==1){ # $p=1$时不需要计算逆矩阵
      Ftinv[] <- 1/c(Ft)
      detFt <- c(Ft)
    } else {
      Ftinv[] <- solve(Ft)
      detFt <- det(Ft)
    }
    
    # 这里计算对数似然函数值
    if(p==1){
      logL <- logL - 0.5*(log(detFt) + c(vt)^2*Ftinv)
    } else{
      logL <- logL - 0.5*(
        log(detFt) + c(vt %*% Ftinv %*% vt)  )
    }
    
    # 向前一步预测
    Kt[] <- Tt %*% Pt %*% Zt_transp %*% Ftinv
    Lt[] <- Tt - Kt %*% Zt
    at <- Tt %*% at + Kt %*% vt
    Pt <- Tt %*% Pt %*% t(Lt) + RQRt
  }
  
  return(logL)
}

30.12.1.2 滤波平滑计算

## 1.2 向前滤波、向后平滑 ####
## 设初始分布$N(\boldsymbol a_1, P_1)$已知，
## 模型矩阵$Z_t, H_t, T_t, R_t, Q_t$都不依赖于$t$。
## 向前滤波得到$\boldsymbol a_t, P_t, \boldsymbol v_t, F_t, K_t, L_t$;
## 向后平滑得到$\boldsymbol r_{t-1}, N_t, \boldsymbol u_t, D_t$,
## 以及$E(\boldsymbol\alpha_t | \boldsymbol y_1, \dots, \boldsymbol y_n)$,
## $\text{Var}(\boldsymbol\alpha_t | \boldsymbol y_1, \dots, \boldsymbol y_n)$,
## 还有观测方程扰动项、系统方程扰动项的平滑分布。
kf_const_fwdbwd <- function(
    y,      # 输入观测值向量（长度n），或者$n \times p$矩阵
    Zmat,   # 观测$Z$矩阵，$p \times m$
    Tmat,   # 转移$T$矩阵，$m \times m$
    Hmat,   # 观测误差方差阵，$p \times p$，$p=1$时为标量
    Rmat = NULL,   # 状态方程误差变换矩阵，$m \times r$
    Qmat,   # 状态方程误差方差阵，$r \times r$
    a1,     # $t=1$时状态初始分布均值
    P1      # $t=1$时状态初始分布方差阵
){
  if(!is.matrix(y)){
    y <- cbind(y)
  }
  # 观测长度：
  n <- nrow(y)
  # 观测维数：
  p <- ncol(y)
  # 状态维数：
  m <- ncol(Zmat)
  # 状态方程误差维数：
  dim_err_sta <- nrow(Qmat)
  if(is.null(Rmat)) Rmat <- diag(m)
  
  Zt <- Zmat
  Tt <- Tmat
  Ht <- Hmat
  Rt <- Rmat
  Qt <- Qmat
  
  # 单个$v_t$向量
  vt <- numeric(p)
  # v[1:n], 保存为$p \times n$矩阵：
  v_arr <- matrix(0, p, n)
  
  # 单个$F_t$矩阵, $p \times p$
  Ft <- matrix(0, p, p)
  # 单个$F_t$逆矩阵, $p \times p$
  Ftinv <- matrix(0, p, p)
  # F[1:n], 保存为$p \times p \times n$的三维数组
  F_arr <- array(0, c(p, p, n))
  # F[1:n]逆矩阵, 保存为$p \times p \times n$的三维数组
  Finv_arr <- array(0, c(p, p, n))
  
  # 单个$a_t$向量
  at <- a1
  # a[1:n]，保存为$m \times n$矩阵：
  a_arr <- matrix(0, m, n)

  # 单个$P_t$矩阵
  Pt <- P1
  # P[1:n], 保存为$m \times m \times n$的三维数组
  P_arr <- array(0, c(m, m, n))

  # 状态滤波的均值向量1:n，保存为$m \times n$矩阵：
  a_filt_arr <- matrix(0, m, n)
  # 状态滤波的协方差阵1:n，保存为$m \times m \times n$三维数组
  P_filt_arr <- array(0, c(m, m, n))
  
  # 单个$K_t$矩阵，$m \times p$
  Kt <- matrix(0, m, p)
  # K[1:n]，保存为$m \times p \times n$的三维数组
  K_arr <- array(0, c(m, p, n))
  
  # 单个$L_t$矩阵，$m \times m$
  Lt <- matrix(0, m, m)
  # L[1:n], 保存为$m \times m \times n$的三维数组
  L_arr <- array(0, c(m, m, n))
  
  # t(Z_t)
  Zt_transp <- t(Zt)
  # $R_t Q_t R_t^T$
  RQRt <- Rt %*% Qt %*% t(Rt)
  
  # 进行向前的一步预报、滤波，并计算似然函数值
  logL <- -0.5*n*p*log(2*pi)
  for(t in 1:n){
    a_arr[,t] <- at
    P_arr[,,t] <- Pt
    
    vt[] <- y[t,] - Zt %*% at
    Ft[] <- Zt %*% Pt %*% Zt_transp + Ht
    v_arr[,t] <- vt
    F_arr[,,t] <- Ft
    if(p == 1){
      Ftinv[] <- 1 / c(Ft)
      detFt <- c(Ft)
    } else {
      Ftinv[] <- solve(Ft)
      detFt <- det(Ft)
    }
    Finv_arr[,,t] <- Ftinv
    
    # 这里计算对数似然函数值
    logL <- logL - 0.5*(
      log(detFt) + c(vt %*% Ftinv %*% vt)  )
    
    # 滤波
    a_filt_arr[,t] <- at + 
      Pt %*% Zt_transp %*% Ftinv %*% vt
    P_filt_arr[,,t] <- Pt - 
      Pt %*% Zt_transp %*% Ftinv %*% Zt %*% Pt
    
    # 向前一步预测
    Kt[] <- Tt %*% Pt %*% Zt_transp %*% Ftinv
    K_arr[,,t] <- Kt
    Lt[] <- Tt - Kt %*% Zt
    L_arr[,,t] <- Lt
    at <- Tt %*% at + Kt %*% vt
    Pt <- Tt %*% Pt %*% t(Lt) + RQRt
  }
  
  # 向后平滑所需存储结构：
  # 单个$r_t$，$m \times 1$
  rt <- numeric(m) # 保存$r_{t-1}$
  rtp1 <- rt # 保存$r_t$
  # r[0:(n-1)], 保存为$m \times n$矩阵
  # 第三下标表示从0开始的时刻
  r_arr <- matrix(0, m, n)
  # 单个$N_t$矩阵, $m \times m$
  Nt <- matrix(0, m, m)  # 保存$N_{t-1}$
  Ntp1 <- Nt # 保存$N_t$
  # N[0:(n-1)]
  # 第三下标表示从0开始的时刻
  N_arr <- array(0, c(m, m, n))
  
  # 状态平滑均值1:n
  a_sm_arr <- matrix(0, m, n)
  # 状态平滑方差阵1:n
  P_sm_arr <- array(0, c(m, m, n))
  
  # 观测方程扰动项平滑均值
  err_obs_mu_arr <- matrix(0, p, n)
  # 观测方程扰动项平滑方差阵
  err_obs_P_arr <- array(0, c(p, p, n))
  # 系统方程扰动项平滑均值
  err_sys_mu_arr <- matrix(0, dim_err_sta, n)
  # 系统方程扰动项平滑方差阵
  err_sys_P_arr <- array(
    0, c(dim_err_sta, dim_err_sta, n))
  
  # 向后平滑迭代
  rt[] <- 0
  Nt[] <- 0
  QtRtt <- Qt %*% t(Rt)
  for(t in n:1){
    rtp1[] <- rt 
    # $r_{t-1}$:
    rt[] <- Zt_transp %*% Finv_arr[,,t] %*% v_arr[,t] +
      t(L_arr[,,t]) %*% rt
    r_arr[,t] <- rt # 保存$r_{t-1}$
    
    Ntp1[] <- Nt 
    # $N_{t-1}$:
    Nt[] <- Zt_transp %*% Finv_arr[,,t] %*% Zt +
      t(L_arr[,,t]) %*% Nt %*% L_arr[,,t]
    N_arr[,,t] <- Nt # 保存$N_{t-1}$
    
    a_sm_arr[,t] <- a_arr[,t] + P_arr[,,t] %*% rt
    P_sm_arr[,,t] <- P_arr[,,t] - 
      P_arr[,,t] %*% Nt %*% P_arr[,,t] 
    
    err_obs_mu_arr[,t] <-
      Ht %*% (
        Finv_arr[,,t] %*% v_arr[,t] -
          t(K_arr[,,t]) %*% rtp1)
    err_obs_P_arr[,,t] <-
      Ht - Ht %*% (
        Finv_arr[,,t] + 
          t(K_arr[,,t]) %*% Ntp1 %*% K_arr[,,t]) %*% Ht
    
    err_sys_mu_arr[,t] <-
      QtRtt %*% rtp1
    err_sys_P_arr[,,t] <-
      Qt - QtRtt %*% Ntp1 %*% t(QtRtt)
  }
  
  res <- list(
    logLik = logL,
    a = a_arr, 
    P = P_arr,
    a_filt = a_filt_arr,
    P_filt = P_filt_arr,
    a_sm = a_sm_arr,
    P_sm = P_sm_arr,
    err_obs_mu = err_obs_mu_arr,
    err_obs_P = err_obs_P_arr,
    err_sys_mu = err_sys_mu_arr,
    err_sys_P = err_sys_P_arr
    )
  return(res)
}

30.12.2 各矩阵时变、初始分布已知、观测一元化情形

计算似然函数的程序。 需要将各个矩阵和观测值转换到一个列表中。

# 将输入转换成需要的统一格式，每个矩阵都是关于$t$的列表。
# 一元化要求$H_t$必须都是对角阵
# 输入y: 若$p=1$，输入观测值向量（长度n），
#   若$p>1$且固定，输入$n \times p$矩阵或数据框
#   若$p_t$随时间变化，输入为长度$n$的列表
kf_uni_input <- function(
    input = list(), # 用来保存已有转化结果
    y,      # 观测值
    Zmat,   # 观测$Z$矩阵，$p \times m$，可以是矩阵的列表
    Tmat,   # 转移$T$矩阵，$m \times m$，可以是矩阵的列表
    Hmat,   # 观测误差方差阵，$p \times p$，$p=1$时为标量，可以是矩阵的列表
    Rmat = NULL,   # 状态方程误差变换矩阵，$m \times r$，可以是矩阵的列表
    Qmat,   # 状态方程误差方差阵，$r \times r$，可以是矩阵的列表
    a1,     # $t=1$时状态初始分布均值
    P1      # $t=1$时状态初始分布方差阵
){
  
  # 将$Z_t$等矩阵转换为每个时间点一个矩阵的列表：
  make_mat_list <- function(Mat){
    # 以Mat保存$Z_t$为例
    if(is.matrix(Mat)){
      # Z_t 不随时间变化的情形
      Lis <- replicate(n, Mat, simplify=FALSE)
    } else if(is.array(Mat) && length(dim(Mat))==3) {
      # 保存为三维数组，Mat[,,t]为$Z_t$
      Lis <- apply(Mat, 3, rbind, simplify=FALSE)
    } else if(is.list(Mat)){
      # 已经是要求的格式
      Lis <- Mat
    } else {
      stop("Mat必须为矩阵、三维数组或者列表")
    }
    
    return(Lis)
  }
  
  if(!missing(y)){
    if(is.data.frame(y)) {
      # 适用于多个分量且维数固定情形
      y <- as.matrix(y)
    }
    
    if(is.matrix(y)) {
      # 适用于分量个数固定情形
      input$y <- apply(y, 1, c, simplify=FALSE)
    } else if(is.list(y)) {
      # 已经是列表形式
      input$y <- y
    } else {
      # 输入为长度n向量
      input$y <- as.list(y)
    }
    # 现在input$y是列表，y[[t]]保存了$t$时刻的各个分量为R向量形式
    input$n <- length(y)
    
    # 每个时刻的观测分量个数：
    input$pt <- vapply(y, length, 1L)
  }
  
  n <- input$n
  
  # $Z_t$矩阵：每个时刻保存为一个列表
  if(!missing(Zmat)){
    input$Z <- make_mat_list(Zmat)
    # 现在input$Z是列表，Z[[t]]为$Z_t$矩阵
  }
  
  # $T_t$矩阵：每个时刻保存为一个列表
  if(!missing(Tmat)){
    input$T <- make_mat_list(Tmat)
    # 现在input$T是列表，T[[t]]为$Z_t$矩阵
    input$m <- nrow(input$T[[1]])
  }
  
  # $H_t$矩阵：每个时刻保存为一个列表
  if(!missing(Hmat)){
    Ht_lis <- make_mat_list(Hmat)
    # 检查每个$H_t$必须都是对角阵
    if(!all(vapply(Ht_lis, Matrix::isDiagonal, TRUE))){
      stop("所有时刻的Hmat必须都是对角阵")
    }
    # H[[t]]保存$H_t$的对角元素
    input$H <- lapply(Ht_lis, diag)
  }
  
  # $R_t$，系统方程误差为$R_t \eta_t$
  if(is.null(input$R) || !missing(Rmat)){
    if(missing(Rmat)){
      Rmat <- diag(input$m)
    }
    input$R <- make_mat_list(Rmat)
    # 现在input$R是列表，R[[t]]为$R_t$
  } 
  
  # $Q_t$，系统扰动$\eta_t$的方差阵
  if(!missing(Qmat)){
    input$Q <- make_mat_list(Qmat)
    # 系统扰动$\eta$的维数，Durbin & Koopman(2012)中$r$
    input$dim_err_sta <- nrow(input$Q[[1]])
  }
  
  if(!missing(a1)){
    input$a1 <- a1
  }
  
  if(!missing(P1)){
    input$P1 <- P1
  }
  
  miss <- setdiff(
    c("y", "Z", "T", "H", "R", "Q", "a1", "P1"),
    names(input))
  if(length(miss) > 0){
    stop(paste("kf_uni_input: missing", miss))
  }
    
  return(input)
}

# 已知初始分布，所有输入已经列表化，
# 用一元化方法向前递推计算对数似然函数值
# x: 其中包含y, Z, T, H, R, Q, a1, P1等元素，
#    每个元素是关于时间$t$的列表。
kf_uni_logL <- function(x){
  n <- x$n
  m <- x$m
  pt <- x$pt
  
  # 判断$F_{t,i}=0$的标准
  Fti_tol <- 1E-8
  
  ati <- x$a1
  Pti <- x$P1
  Zti <- x$Z[[1]][1,]
  Kti <- Zti[]
  logL <- 0
  num_innov <- 0 # 累计$F_{t,i}>0$个数
  for(t in 1:n){
    # 当前ati的值是$\boldsymbol a_t$,
    # Pti的值是$P_t$
    for(i in 1:pt[t]){
      # 同一时间点的更新
      Zti[] <- x$Z[[t]][i,]
      vti <- x$y[[t]][i] - sum(Zti * ati)
      Fti <- c(Zti %*% Pti %*% Zti) + x$H[[t]][i]
      if(abs(Fti) > Fti_tol){
        num_innov <- num_innov + 1
        logL <- logL - 0.5 * (
          log(Fti) + vti^2 / Fti )
        
        Kti[] <- (Pti %*% Zti) / Fti
        ati[] <- ati + Kti * vti
        Pti[] <- Pti - outer(Kti, Kti)*Fti
      } # else $y_{t,i}$不包含新的信息，不用更新
    }
    # 当前ati和和Pti的值代表滤波分布
    # $p(\boldsymbol\alpha_t | \boldsymbol Y_t)$的参数
    
    # 从t到t+1的更新
    ati[] <- x$T[[t]] %*% ati
    Pti[] <- x$T[[t]] %*% Pti %*% t(x$T[[t]]) + 
      x$R[[t]] %*% x$Q[[t]] %*% t(x$R[[t]])
    # 当前ati和和Pti的值代表预报分布
    # $p(\boldsymbol\alpha_{t+1} | \boldsymbol Y_t)$的参数
  }
  logL <- logL - 0.5*log(2*pi)*num_innov
  
  return(logL)
}

30.12.3 各矩阵时变、初始分布发散、观测一元化情形

计算似然函数的程序。 需要将各个矩阵和观测值转换到一个统一的列表中保存。

# 将输入转换成需要的统一格式
# 要求$H_t$必须都是对角阵
# 输入:
#  y: 若$p=1$，输入观测值向量（长度n），
#     若$p>1$且固定，输入$n \times p$矩阵或数据框
#     若$p_t$随时间变化，输入为长度$n$的列表
#  P1inf: 输入时表示发散先验部分$P_{\infty}$，不输入表示没有发散先验
#  P1: 没有P1inf时就是$\boldsymbol\alpha_1$方差阵，否则是$P_{\star}$
kf_uni_input <- function(
    input = list(), # 用来保存已有转化结果
    y,      # 观测值
    Zmat,   # 观测$Z$矩阵，$p \times m$
    Tmat,   # 转移$T$矩阵，$m \times m$
    Hmat,   # 观测误差方差阵，$p \times p$，$p=1$时为标量
    Rmat = NULL,   # 状态方程误差变换矩阵，$m \times r$
    Qmat,   # 状态方程误差方差阵，$r \times r$
    a1,     # $t=1$时状态初始分布均值
    P1,     # $t=1$时状态初始分布方差阵，非发散部分
    P1inf   # 初始分布方差阵发散部分
){
  
  # 将$Z_t$等矩阵转换为列表的函数
  make_mat_list <- function(Mat){
    # 以Mat保存$Z_t$为例
    if(is.matrix(Mat)){
      # Z_t 不随时间变化的情形
      Lis <- replicate(n, Mat, simplify=FALSE)
    } else if(is.array(Mat) && length(dim(Mat))==3) {
      # 保存为三维数组，Mat[,,t]为$Z_t$
      Lis <- apply(Mat, 3, rbind, simplify=FALSE)
    } else if(is.list(Mat)){
      # 已经是要求的格式
      Lis <- Mat
    } else {
      stop("Mat必须为矩阵、三维数组或者列表")
    }
    
    return(Lis)
  }
  
  if(!missing(y)){
    if(is.data.frame(y)) {
      # 适用于多个分量且维数固定情形
      input$y <- as.matrix(y)
    }
    
    if(is.matrix(y)) {
      # 适用于分量个数固定情形
      input$y <- apply(y, 1, c, simplify=FALSE)
    } else if(is.list(y)) {
      # 已经是列表形式
      input$y <- y
    } else {
      # 输入为长度n向量
      input$y <- as.list(y)
    }
    # 现在input$y是列表，y[[t]]保存了$t$时刻的各个分量为R向量形式
    input$n <- length(y)
    
    # 每个时刻的观测分量个数：
    input$pt <- vapply(y, length, 1L)
  }
  
  n <- input$n
  
  # $Z_t$矩阵：每个时刻保存为一个列表
  if(!missing(Zmat)){
    input$Z <- make_mat_list(Zmat)
    # 现在input$Z是列表，Z[[t]]为$Z_t$矩阵
  }
  
  # $T_t$矩阵：每个时刻保存为一个列表
  if(!missing(Tmat)){
    input$T <- make_mat_list(Tmat)
    # 现在input$T是列表，T[[t]]为$Z_t$矩阵
    input$m <- nrow(input$T[[1]])
  }
  
  # $H_t$矩阵：每个时刻保存为一个列表
  if(!missing(Hmat)){
    Ht_lis <- make_mat_list(Hmat)
    # 检查每个$H_t$必须都是对角阵
    if(!all(vapply(Ht_lis, Matrix::isDiagonal, TRUE))){
      stop("所有时刻的Hmat必须都是对角阵")
    }
    # H[[t]]保存$H_t$的对角元素
    input$H <- lapply(Ht_lis, diag)
  }
  
  # $R_t$，系统方程误差为$R_t \eta_t$
  if(is.null(input$R) || !missing(Rmat)){
    if(missing(Rmat)){
      Rmat <- diag(input$m)
    }
    input$R <- make_mat_list(Rmat)
    # 现在input$R是列表，R[[t]]为$R_t$
  } 
  
  # $Q_t$，系统扰动$\eta_t$的方差阵
  if(!missing(Qmat)){
    input$Q <- make_mat_list(Qmat)
    # 系统扰动$\eta$的维数，Durbin & Koopman(2012)中$r$
    input$dim_err_sta <- nrow(input$Q[[1]])
  }
  
  if(!missing(a1)){
    input$a1 <- a1
  }
  
  if(!missing(P1)){
    input$P1 <- P1
  }
  
  if(!missing(P1inf)){
    input$P1inf <- P1inf
  }
  
  miss <- setdiff(
    c("y", "Z", "T", "H", "R", "Q", "a1", "P1"),
    names(input))
  if(length(miss) > 0){
    stop(paste("kf_uni_input: missing", miss))
  }
  
  return(input)
}

# 已知允许发散先验精确初始化，所有输入已经列表化，
# 用一元化方法向前递推计算发散先验对数似然函数值
# 注意$H_t$矩阵是对角阵
kf_uni_dif_logL <- function(x){
  # 判断$F_{t,i}=0$和$P_{\infty,t,i}=0$的标准
  Fti_tol <- 1E-8
  
  n <- x$n
  m <- x$m
  pt <- x$pt
  
  is_diffuse <- is.null(x$P1inf)
  if(is_diffuse){
    x$P1inf <- matrix(0, m, m)
  }
  
  logL <- 0
  num_innov <- 0 # 累计$F_{t,i}>0$个数
  
  # ati: 发散时保存$a^{(0)}_{t,i}$，以后保存$a_{t,i}$
  ati <- x$a1 
  # Pti: 发散时保存$P_{*,t,i}$，以后保存$P_{t,i}$
  Pti <- x$P1
  # Ptii: 发散时保存$P_{\infty,t,i}$，以后保持为0矩阵
  Ptii <- x$P1inf
  # Ftii: 保存$F_{\infty,t,i}$标量
  # Zti: $Z_t$第$i$行
  Zti <- x$Z[[1]][1,]
  # K0ti: $K_{t,i}^{(0)}$
  K0ti <- Zti[]
  # K1ti: $K_{t,i}^{(1)}$
  K1ti <- Zti[]
  # Mti: $M_{*,t,i}$
  Mti <- Zti[]
  # Mtii: $M_{\infty,t,i}$
  Mtii <- Zti[]
  # L0ti: $L^{(0)}_{t,i}$
  L0ti <- diag(m)
  # L1ti: $L^{(1)}_{t,i}$
  L1ti <- diag(m)
  Tt <- x$T[[1]]
  
  d <- 0 # 记住$P_{\infty,t}$不等于零的$t$的个数
  for(t in 1:n){ # 发散阶段
    if(max(abs(diag(Ptii))) < Fti_tol){
      break
    } else {
      d <- d + 1
    }
    
    # 当前$ati$保存了$\boldsymbol a_t$,
    # $Pti$保存了$P_{*,t}$, $Ptii$保存了$P_{\infty,t}$
    # 即$\boldsymbol\alpha_t$在$\boldsymbol Y_{t-1}$下的条件分布

    for(i in 1:pt[t]){
      # 下面是按照有无穷方差部分向前递推，
      # 同一时间点$t$的多个观测值
      Zti[] <- x$Z[[t]][i,]
      Fti <- c(Zti %*% Pti %*% t(Zti)) + x$H[[t]][i]
      Ftii <- c(Zti %*% Ptii %*% t(Zti))
      if(Fti > Fti_tol) {
        num_innov <- num_innov + 1
      }
      
      if(Ftii > Fti_tol){ # $F_{\infty,t,i} > 0$?
        Mti[] <- Pti %*% Zti
        Mtii[] <- Ptii %*% Zti
        F0ti <- 0
        F1ti <- 1/Ftii
        F2ti <- -Fti / Ftii^2
        K0ti[] <- Mtii * F1ti
        K1ti[] <- Mti * F1ti + Mtii * F2ti
        L0ti[] <- diag(m) - outer(K0ti, Zti)
        L1ti[] <- -outer(K1ti, Zti)
        v0ti <- x$y[[t]][i] - sum(Zti * ati)
        
        ati[] <- ati + K0ti*v0ti
        Pti[] <- Ptii %*% t(L1ti) + Pti %*% t(L0ti)
        Ptii[] <- Ptii %*% t(L0ti)
        
        # 似然
        if(Fti > Fti_tol){ # $F_{*,t,i} > 0$时才有贡献
          logL <- logL - 0.5*log(Ftii)
        }
        
      } else { # $F_{\infty,t,i} = 0$
        Mti[] <- Pti %*% Zti
        K0ti[] <- Mti / Fti
        L0ti[] <- diag(m) - outer(K0ti, Zti)
        v0ti <- x$y[[t]][i] - sum(Zti * ati)
        
        ati[] <- ati + K0ti * v0ti
        Pti[] <- Pti %*% t(L0ti)
        Ptii[] <- Ptii %*% t(L0ti)
        
        # 似然
        if(Fti > Fti_tol){ # $F_{*,t,i} > 0$?
          logL <- logL - 0.5*(log(Fti) + v0ti^2/Fti)
        }
      }
      
    } # for i
    
    # 当前$ati$和$Pti$, $Ptii$保存了$\boldsymbol\alpha_t$
    # 在$\boldsymbol Y_{t}$下的条件分布,
    # 即滤波分布

    # 从$t$到$t+1$步，时间$(t,p_t+1)$到$(t+1,1)$
    # 相当于观测缺失，没有新息，
    # $K_{t,p_t+1}=0$, $L_{t,p_t+1}=T_t$
    # 只有状态转移与扰动
    # 均值只需要状态转移，没有新息可以修正
    ati[] <- Tt %*% ati 
    Pti[] <- Tt %*% Pti %*% t(Tt) + 
      x$R[[t]] %*% x$Q[[t]] %*% t(x$R[[t]])
    Ptii[] <- Tt %*% Ptii %*% t(Tt)
    
    # 当前$ati$保存了$\boldsymbol a_{t+1}$,
    # $Pti$保存了$P_{*,t+1}$, $Ptii$保存了$P_{\infty,t}$，
    # 即$\boldsymbol\alpha_{t+1}$在$\boldsymbol Y_t$下的条件分布
  } # for t, 到t=d为止
  
  Kti <- numeric(m)
  for(t in (d+1):n){
    
    # 当前$ati$保存了$\boldsymbol a_t$,
    # $Pti$保存了$P_{*,t}$, $Ptii$保存了$P_{\infty,t}$
    # 即$\boldsymbol\alpha_t$在$\boldsymbol Y_{t-1}$下的条件分布
    
    for(i in 1:pt[t]){
      # 同一时间点的更新
      Zti[] <- x$Z[[t]][i,]
      vti <- x$y[[t]][i] - sum(Zti * ati)
      Fti <- c(Zti %*% Pti %*% Zti) + x$H[[t]][i]
      if(abs(Fti) > Fti_tol){
        num_innov <- num_innov + 1
        logL <- logL - 0.5 * (
          log(Fti) + vti^2 / Fti )
        
        Kti[] <- (Pti %*% Zti) / Fti
        ati[] <- ati + Kti * vti
        Pti[] <- Pti - outer(Kti, Kti)*Fti
      } # else $y_{t,i}$不包含新的信息，不用更新
    }
    
    # 当前$ati$和$Pti$保存了$\boldsymbol\alpha_t$
    # 在$\boldsymbol Y_{t}$下的条件分布,
    # 即滤波分布
    
    # 从t到t+1的更新, 看成从$(t, p_t+1)$到$(t+1,1)$的时间变化
    # 没有新息，$K_{t,p_t+1}=0$, $v_{t,p_t+1}=0$, $L_{t,p_t+1}=T_t$
    Tt[] <- x$T[[t]]
    ati[] <- Tt %*% ati
    Pti[] <- Tt %*% Pti %*% t(Tt) + 
      x$R[[t]] %*% x$Q[[t]] %*% t(x$R[[t]])
    
    # 当前$ati$保存了$\boldsymbol a_{t+1}$,
    # $Pti$保存了$P_{*,t+1}$, $Ptii$为零矩阵，
    # 即$\boldsymbol\alpha_{t+1}$在$\boldsymbol Y_t$下的条件分布
  }
  logL <- logL - 0.5*log(2*pi)*num_innov
  
  return(logL)
}

B 参考文献

Durbin, James, and Siem Jan Koopman. 2012. Time Series Analysis by State Space Methods. 2nd ed. Oxford University Press.