状态空间模型一般包括一个量测方程 (mesurement equation) 和一个转移方程 (transition equation)，后者描述了状态如何与观测向量之间相互作用。
一个学习文档 (存放的是谷歌云盘)； 以下内容几乎来源于郑挺国老师博士论文。
Kalman Filter：
状态空间模型的一般结构
\(x_t\) $ , r \times 1$ 它被称为系统状态或状态向量。给定 \(x_{t−1}\) 的实现，\(x_t\) 的条件均值等于 \(\Phi_{t}(x_{t−1})\)。函数 \(\Phi_{t}(·)\) 的下标 \(t\) 表示此函数可能依赖于时间。
新生向量 \(v_t\) 是一个独立的 \(r \times 1\) 维白噪声过程(White Noise)。
\[\mathbf{x}_{t} = \Phi_{t} \left(\mathbf{x}_{t-1}\right)+\mathbf{v}_{t} \tag{1} \label{eq1}
\]
量测方程将一个 \(N \times 1\) 维观测向量表示为同期状态 \(x_t\) 和误差项 \(u_t\) 的一个 (可能时间相依的) 函数，即
\[\mathbf{y}_{t}=\Psi_{t}\left(\mathbf{x}_{t}\right)+\mathbf{u}_{t} \tag{2} \label{eq2}
\]
这里 \(u_t\) 也是一个独立的 \(N \times 1\) 维白噪声过程。函数 \(\Psi_{t}(·)\) 度量了观测向量与转移向量之间的相关性。
对以上转移方程和量测方程，外生或前定变量可能引入到函数 \(\Phi_t (·)\) 和 \(\Psi_{t}(·)\) 中，而新生向量 \(v_t\) 和 \(u_t\) 可能为高斯分布也可能为非高斯分布，
为以下讨论方便，这里分别令 \(\mathcal{X}_{t}\) 和 $\mathcal{Y}_{t} $ 表示 \(\left(\mathbf{x}_{1}^{\prime}, \mathbf{x}_{2}^{\prime}, \ldots, \mathbf{x}_{t}^{\prime}\right)^{\prime}\)和 \(\left(\mathbf{y}_{1}^{\prime}, \mathbf{y}_{2}^{\prime}, \ldots, \mathbf{y}_{t}^{\prime}\right)^{\prime}\)，\(\theta\) 表示该状态空间模型中含有的待估参数向量。
确切滤波
公式 \eqref{eq1} 和公式 \eqref{eq2} 直接计算是非常难的，在具有 T 个观测值下，联合分布（Joint Distribution）：
\[p\left(\mathcal{Y}_{T}, \mathcal{X}_{T} ; \theta\right)=\prod_{t=1} p\left(\mathbf{y}_{t} \mid \mathcal{Y}_{t-1}, \mathcal{X}_{t} ; \theta\right) p\left(\mathbf{x}_{t} \mid \mathcal{Y}_{t-1}, \mathcal{X}_{t-1} ; \theta\right)
\]
似然函数（Likelihood fuction）
\[p\left(\mathcal{Y}_{T} ; \theta\right)=\int \prod_{t=1}^{T} p\left(\mathbf{y}_{t} \mid \mathcal{Y}_{t-1}, \mathcal{X}_{t} ; \theta\right) p\left(\mathbf{x}_{t} \mid \mathcal{Y}_{t-1}, \mathcal{X}_{t-1} ; \theta\right) \prod_{t=1}^{T} d \mathbf{x}_{t}
\]
这是一个积分维数等于样本个数与未观测变量 \(x_t\) 对应维数乘积的高维积分，因此这在计算上实际是不可行的。
令 \(p(x_{t−1}|Y_{t−1})\) 为 \(t = 1, 2, \cdots , T\) 次迭代输入的滤子密度，对第 1 次迭代可以考虑 \(x_0\) 的无条件分布 \(p(x_0)\)。
算法 1 (一般确切滤波).
- 在 t = 1，初始化输入密度 \(p(x_0|Y_0)\) ;
- 计算 \(x_{t}\) 和 \(x_{t−1}\) 的联合条件分布，它可分解为输入密度和转移密度的乘积：
\[p\left(\mathbf{x}_{t}, \mathbf{x}_{t-1} \mid \mathcal{Y}_{t-1}\right)=p\left(\mathbf{x}_{t} \mid \mathbf{x}_{t-1}\right) p\left(\mathbf{x}_{t-1} \mid \mathcal{Y}_{t-1}\right)
\]
- 由边际化，计算 \(x_t\)的预测密度 (prediction density) ：
\[p\left(\mathbf{x}_{t} \mid \mathcal{Y}_{t-1}\right)=\int p\left(\mathbf{x}_{t}, \mathbf{x}_{t-1} \mid \mathcal{Y}_{t-1}\right) d \mathbf{x}_{t-1}
\]
- 计算 \(y_t\) 和 \(x_t\) 的联合密度，它可以分解为量测密度（measurement density）和预测密度的乘积，
\[p\left(\mathbf{y}_{t}, \mathbf{x}_{t} \mid \mathcal{Y}_{t-1}\right)=p\left(\mathbf{y}_{t} \mid \mathbf{x}_{t}\right) p\left(\mathbf{x}_{t} \mid \mathcal{Y}_{t-1}\right)
\]
\[p\left(\mathbf{y}_{t} \mid \mathcal{Y}_{t-1}\right)=\int p\left(\mathbf{y}_{t}, \mathbf{x}_{t} \mid \mathcal{Y}_{t-1}\right) d \mathbf{x}_{t}
\]
这个式子特别有用，它可以帮助我们轻松地构建模型的似然函数。
- 由条件概率公式，计算输出的滤子密度 (filtering density) ：
\[p\left(\mathbf{x}_{t} \mid \mathcal{Y}_{t}\right)=\frac{p\left(\mathbf{y}_{t}, \mathbf{x}_{t} \mid \mathcal{Y}_{t-1}\right)}{p\left(\mathbf{y}_{t} \mid \mathcal{Y}_{t-1}\right)}
\]
- 如果 $ t < T $，设 $ t = t + 1 $，并返回第 (ii) 步，否则结束。
只有在非常特殊的情形中，我们才可能从此一般滤波算法中获得解析递归算法，例如高斯、线性假设下的状态空间模型可以采用卡尔曼滤波，
在高斯线性情形中，初始输入滤子密度 $p (x_0 | Y_0) $、量测密度和转移密度假设为高斯，在每步算法，高斯性质被保存下来，于是所有输出也都是高斯的。
线性高斯模型及其估计
假设一般性状态空间模型 公式 \eqref{eq1} 和公式 \eqref{eq2} 中函数 \(\Phi_t (·)\) 和 \(\Psi_t (·)\)都是线性的，且误差向量 \(u_t\) 和 \(v_t\) 均服从于高斯正态分布。满足这些假设的模型就称为线性高斯状态空间
模型。
模型结构
转移方程：
\[\mathbf{x}_{t}=\Gamma \mathbf{x}_{t-1}+\alpha+\mathbf{v}_{t}
\]
量测方程：
\[\mathbf{y}_{t}=H \mathbf{x}_{t}+\mu+\mathbf{u}_{t}
\]
where
\[\mathbf{u}_{t} \sim i i d \mathcal{N}(0, Q) \quad \mathbf{v}_{t} \sim i i d \mathcal{N}(0, R)
\]
初始条件设为
\[\mathbf{x}_{0} \sim \mathcal{N}\left(\bar{a}_{0}, \bar{P}_{0}\right)
\]
初始条件一般采用 \(x_t\) 平稳条件的无条件均值和无条件方差，也即是稳态条件下求得的均值和方差。
对所有的t
\[E\left(\mathbf{v}_{t} \mathbf{x}_{0}^{\prime}\right)=0 \quad \quad E\left(\mathbf{u}_{t} \mathbf{x}_{0}^{\prime}\right)=0
\]
上述模型的一些量 \(\alpha , \mu , \bar{a}_0 , \Gamma , H, Q , R和 \bar{P}_0\)均假设为相应维数的参数向量和矩阵，因此该线性高斯状态空间模型又是时齐次的。
卡尔曼滤波
对上述模型，它满足转移密度 \(p (x_t | x_{t−1})\) 和量测密度 \(p (y_t | x_t)\) 都是正态的。可以证明预测密度和滤子密度也是正态的:
\[ \begin{aligned}
\mathbf{x}_{t} \mid \mathcal{Y}_{t-1} & \sim \mathcal{N}\left(a_{t \mid t-1}, \Sigma_{t \mid t-1}\right) \\
\mathbf{x}_{t} \mid \mathcal{Y}_{t} & \sim \mathcal{N}\left(a_{t \mid t}, \Sigma_{t \mid t}\right) \\
\mathbf{y}_{t} \mid \mathcal{Y}_{t-1} & \sim \mathcal{N}\left(\mathbf{y}_{t \mid t-1}, F_{t \mid t-1}\right)
\end{aligned} \]
因此，我们必须找到条件均值 \(a_{t | t −1}\) , \(a_{t | t}\) , \(y_{t | t −1}\)和条件方差协方差矩阵 \(\sum_{t \mid t-1}, \Sigma_{t \mid t}, F_{t \mid t-1}\),这些量可以通过运行卡尔曼滤波递归地获得。
算法 2 (卡尔曼滤波).
- 初始化( \(t =1\) ) ：设 \(a_{0 | 0} =\bar{a}_0\) ，\(\Omega_{0 | 0} =\bar{P}_0\) ;
- 预测( 从 t −1 到 t ) ：给定 $a_{t −1 | t −1} $ 和 \(\Omega_{t −1 | t −1}\) ，计算
\[ \begin{aligned}
a_{t \mid t-1} &=\Gamma a_{t-1 \mid t-1}+\alpha \\
\Sigma_{t \mid t-1} &=\Gamma \Sigma_{t-1 \mid t-1} \Gamma^{\prime}+R \\
\mathbf{y}_{t \mid t-1} &=H a_{t \mid t-1}+\mu \\
F_{t \mid t-1} &=H \Sigma_{t \mid t-1} H^{\prime}+Q
\end{aligned} \]
\[ \begin{aligned}
K_{t} &=\Sigma_{t \mid t-1} H^{\prime} F_{t \mid t-1}^{-1} \\
a_{t \mid t} &=a_{t \mid t-1}+K_{t}\left(\mathbf{y}_{t}-\mathbf{y}_{t \mid t-1}\right) \\
\Sigma_{t \mid t} &=\Sigma_{t \mid t-1}-K_{t} H \Sigma_{t \mid t-1}
\end{aligned} \]
- 如果 \(t < T\) ，设 \(t =t +1\)，并返回第 （ii） 步，否则结束。
\(K_t\)为卡尔曼增益，它决定了分配给包含在预测误差中关于 \(x_t\) 新信息的权重。关于卡尔曼滤波的推导更为详细的讨论，参见 Hamilton(1994) 。
平滑滤波
平滑滤波主要用于基于所有可获取信息，即 \(\mathcal{Y}_T\) ，直接对 \(x_t\) 进行推断，并提供了平滑密度\(p(x_t|\mathcal{Y}_T)\) 。
令 \(p ( x_{t +1} | \mathcal{Y}_T )\) 为迭代输入密度，\(t =T−1, T−2, . . . , 1\)。我们可以基于信息集 \(\mathcal{Y}_t\)，将\(x_{t +1}\) , \(x_t\) 的联合密度分解为转移密度和滤子密度的乘积：
\[p\left(\mathbf{x}_{t+1}, \mathbf{x}_{t} \mid \mathcal{Y}_{t}\right)=p\left(\mathbf{x}_{t+1} \mid \mathbf{x}_{t}\right) p\left(\mathbf{x}_{t} \mid \mathcal{Y}_{t}\right)
\]
基于得自滤波的预测密度，我们获得以下条件密度：
\[p\left(\mathbf{x}_{t} \mid \mathbf{x}_{t+1}, \mathcal{Y}_{t}\right)=\frac{p\left(\mathbf{x}_{t+1}, \mathbf{x}_{t} \mid \mathcal{Y}_{t}\right)}{p\left(\mathbf{x}_{t+1} \mid \mathcal{Y}_{t}\right)}
\]
其中 \(p ( x_{t +1} | \mathcal{Y}_t )\) 为 t +1 时刻的预测密度，\(p ( x_{t} | \mathcal{Y}_t)\) 为 t 时刻的滤子密度，两者可由确切滤波迭代过程生成，而转移密度 \(p ( x_{t +1} | x_{t} )\) 在本文混合模型中被假设为已知正态密度。
基于信息集 \(\mathcal{Y}_T\) ，\(x_{t +1}\) , \(x_t\) 的联合密度由条件密度 \(p (x_t | x_{t +1} , Y_T)\) 和算法输入密度 \(p ( x_{t +1} | Y_{T})\)的乘积给出。
信息集 \(x_{t +1}\) , \(Y_T\) 包含在信息集 \(x_{t +1}\) , \(Y_t\) , \(\mathbf{u}_{t +1}^{T}\) , \(\mathbf{v}^{T}{t +2}\) 中，其中 \(\mathbf{u}_{t+1}^{T}=\left(\mathbf{u}_{t+1}, \ldots, \mathbf{u}_{T}\right)^{\prime}\)，\(\mathbf{v}_{t+1}^{T}=\left(\mathbf{v}_{t+1}, \ldots, \mathbf{v}_{T}\right)^{\prime}\)
假如 \(\mathbf{u}_{t +1}^{T}\) 和 \(\mathbf{v}^{T}{t +2}\) 独立于 \(x_{t+1}\) , \(x_t\) ,\(Y_t\) ，我们就可以推断\(p\left(\mathbf{x}_{t} \mid \mathbf{x}_{t+1}, \mathcal{Y}_{T}\right) = p\left(\mathbf{x}_{t} \mid \mathbf{x}_{t+1}, \mathcal{Y}_{t}\right)\)。
\[p\left(\mathbf{x}_{t+1}, \mathbf{x}_{t} \mid \mathcal{Y}_{T}\right)=p\left(\mathbf{x}_{t} \mid \mathbf{x}_{t+1}, \mathcal{Y}_{T}\right) p\left(\mathbf{x}_{t+1} \mid \mathcal{Y}_{T}\right)=p\left(\mathbf{x}_{t} \mid \mathbf{x}_{t+1}, \mathcal{Y}_{t}\right) p\left(\mathbf{x}_{t+1} \mid \mathcal{Y}_{T}\right)
\]
最后，由边际化我们获得 \(x_t\) 的平滑密度，即输出密度为：
\[p\left(\mathbf{x}_{t} \mid \mathcal{Y}_{T}\right)=\int p\left(\mathbf{x}_{t+1}, \mathbf{x}_{t} \mid \mathcal{Y}_{T}\right) d \mathbf{x}_{t+1}=\int p\left(\mathbf{x}_{t} \mid \mathbf{x}_{t+1}, \mathcal{Y}_{t}\right) p\left(\mathbf{x}_{t+1} \mid \mathcal{Y}_{T}\right) d \mathbf{x}_{t+1}
\]
注意只有在线性高斯情形中，可能获得解析逆向递归。以下介绍线性高斯情形下的卡尔曼平滑估计。
平滑估计值\(\left\{a_{t \mid T}, t=1, \ldots, T\right\}\),为我们提供了关于 \(x_t\) 更为精确的推断。为获得平滑估计值，我们必须首先运行卡尔曼滤波，并保存相应预测估计值 \(\{a_{t | t-1} \}\) 和滤子估计值 \(\{a_{t | t} \}\)，以及它们的均方误差 ( MSE ) 矩阵 \(\{\Omega_{t | t-1}\}\) 和 \(\{\Omega_{t | t}\}\)。
\[ \begin{aligned}
a_{t \mid T} &=a_{t \mid t}+\Sigma_{t \mid t} T^{\prime} \Sigma_{t+1 \mid t}^{-1}\left(a_{t+1 \mid T}-a_{t+1 \mid t}\right) \\
\Sigma_{t \mid T} &=\Sigma_{t \mid t}+\Sigma_{t \mid t} T^{\prime} \Sigma_{t+1 \mid t}^{-1}\left(\Sigma_{t+1 \mid T}-\Sigma_{t+1 \mid t}\right)\left(\Sigma_{t+1 \mid t}^{-1}\right)^{\prime} T \Sigma_{t \mid t}^{\prime}
\end{aligned} \]
这里 \(a_{T| T}\)和 \(\Omega_{T| T}\)为平滑算法的初始值，是卡尔曼滤波在最后一次迭代 ( 即在 T时刻) 获得的滤子估计值。
在线性高斯模型中，平滑估计值 \(a_{t | T}\)和 \(Σ_{t | T}\) 分别为 \(x_t\) 基于 \(Y_T\)的均值和方差，且
\[\mathbf{x}_{t} \mid \mathcal{Y}_{T} \sim \mathcal{N}\left(a_{t \mid T}, \Sigma_{t \mid T}\right)
\]
极大似然估计
在线性高斯状态空间模型中，参数 \(\theta\) 的极大似然 (ML) 估计特别简单，因为卡尔曼滤波可以用来形成极大似然估计中可能使用到的一些统计量。
在正态假设下，\(y_t\) 基于 \(\mathcal{Y_{t−1}}\) 的分布是一个均值为 \(y_{t|t−1}\)、方差协方差矩阵为 \(F_{t|t−1}\) 的 N 维正态分布。因此，\(y_{t}\) 的条件密度可表示为:
\[p\left(\mathbf{y}_{t} \mid Y^{t-1} ; \theta\right)=\frac{1}{(2 \pi)^{N / 2}\left|F_{t \mid t-1}\right|^{1 / 2}} \exp \left(-\frac{1}{2} \eta_{t \mid t-1}^{\prime} F_{t \mid t-1}^{-1} \eta_{t \mid t-1}\right)
\]
where \(\eta_{t \mid t-1} = y_{t} - y_{t|t-1}\)
\[ \begin{aligned}
\ell(\theta) &=\log \mathcal{L}(\theta)=\sum_{t=1}^{T} \log p\left(\mathbf{y}_{t} \mid \mathcal{Y}_{t-1} ; \theta\right) \\
&=-\frac{N T}{2} \log 2 \pi-\frac{1}{2} \sum_{t=1}^{T} \log \left|F_{t \mid t-1}\right|-\frac{1}{2} \sum_{t=1}^{T} \eta_{t \mid t-1}^{\prime} F_{t \mid t-1}^{-1} \eta_{t \mid t-1}
\end{aligned} \tag{3} \label{eq3}\]
实际上，\(\eta_{t \mid t-1}, F_{t \mid t-1}\)正好可由卡尔曼滤波简单获得。
算法 3 (线性高斯状态空间模型的极大似然估计)
- 选取参数 \(\theta\) 的一个初始值，即 \(\theta_0\) ；
- 运行卡尔曼滤波，并保存序列\(\left\{\eta_{t \mid t-1}\left(\theta_{0}\right)\right\}\) 和 \(\left\{F_{t \mid t-1}\left(\theta_{0}\right)\right\}\)
- 利用序列\(\left\{\eta_{t \mid t-1}\left(\theta_{0}\right)\right\}\) 和 \(\left\{F_{t \mid t-1}\left(\theta_{0}\right)\right\}\)计算 \eqref
- 用优化程序重复上述步骤，直至 \eqref{eq3} 达到最大化。
通过算法3，可以获得极大似然的估计值:
\[\hat{\theta}=\underset{\theta \in \Theta}{\arg \max } \ell(\theta)
\]
梯度 (gradients) 可以通过采用有限差分值计算。但是，通过利用解析梯度可以显著地稳定优化过程。因子向量的第 \(i\) 元素为：
\[\frac{\partial \ell(\theta)}{\partial \theta_{i}}=-\frac{1}{2} \sum_{t=1}^{T}\left[\operatorname{tr}\left(\left(F_{t \mid t-1}^{-1} \frac{\partial F_{t \mid t-1}}{\partial \theta_{i}}\right)\left(I-F_{t \mid t-1}^{-1} \eta_{t} \eta_{t}^{\prime}\right)\right)-\frac{\partial \eta_{t}^{\prime}}{\partial \theta_{i}} F_{t \mid t-1}^{-1} \eta_{t}\right]
\]
如何计算\(\frac{\partial \eta_{t}^{\prime}}{\partial \theta_{i}}\) 和 \(\frac{\partial F_{t}^{\prime}}{\partial \theta_{i}}\)
在一定条件下，极大似然估计 \(\theta\) 是渐近正态分布的:
\[\sqrt{T}\left(\hat{\theta}-\theta^{0}\right) \stackrel{d}{\rightarrow} \mathcal{N}\left(0, \mathcal{I}_{H T}^{-1}\right)\quad \quad where \ \ \mathcal{I}_{H T}=-\frac{1}{T} E\left(\left.\frac{\partial^{2} \ell(\theta)}{\partial \theta \partial \theta^{\prime}}\right|_{\theta=\theta^{0}}\right)
\]
对此渐近分布的有效正则条件包含以下条件：真实参数向量 \(\theta^0\) 位于参数空间内部，转移方程平稳，且参数是可识别的,
信息矩阵可一致地由以下方程估计，即
\[\hat{\mathcal{I}}_{H T}=-\left.\frac{1}{T} \frac{\partial^{2} \ell(\theta)}{\partial \theta \partial \theta^{\prime}}\right|_{\theta=\hat{\theta}}=-\left.\frac{1}{T} \sum_{t=1}^{T} \frac{\partial^{2} \log p\left(\mathbf{y}_{t} \mid \mathcal{Y}_{t-1} ; \theta\right)}{\partial \theta \partial \theta^{\prime}}\right|_{\theta=\hat{\theta}}
\]
经验研究中给出的 \(\hat{\theta}\) 估计方差协方差矩阵为 :
\[\widehat{\operatorname{Var}(\hat{\theta})}=\frac{1}{T}\left(\hat{\mathcal{I}}_{H T}\right)^{-1}
\]
如果在线性模型中，状态新生和度量误差不是高斯分布，那么人们仍然能通过 (错误地) 假设正态性、借助卡尔曼滤波计算对数似然值及关于 \(\theta\) 最大化来获得模型参数估计值。这种方法
称为伪极大似然估计。在某些条件下，它将仍然产生一致性估计量，这种估计量是渐近正态分布的。
