狀態空間模型一般包括一個量測方程 (mesurement equation) 和一個轉移方程 (transition equation)，后者描述了狀態如何與觀測向量之間相互作用。
學習文檔：https://drive.google.com/file/d/1R9WGousjJlFXz-dBQLZlRhZNqJvHU3wR/view
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{eq3}
- 用優化程序重復上述步驟，直至 \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\) 最大化來獲得模型參數估計值。這種方法
稱為偽極大似然估計。在某些條件下，它將仍然產生一致性估計量，這種估計量是漸近正態分布的。
