この記事はカルマンフィルタの勉強 ＃備忘録_1の続きとなります．
ローカルレベルモデル
はじめに
単変量時系列の観測値を$y_1,\cdots,y_n$とする．観測値は真の水準$\alpha_t$に正規ホワイトノイズ$\epsilon_t$が乗ることで
y_t=\alpha_t+\epsilon_t,\quad \epsilon\sim N(0,\sigma_\epsilon^2),\quad t=1,\cdots,n\tag{1}
のように生成されていると仮定する．$N(0,\sigma_\epsilon^2)$は平均ゼロ，分散$\sigma_\epsilon^2$の正規分布を表す．
真の水準$\alpha_t$は観測されない潜在変数であり状態と呼ばれる．観測誤差$\epsilon_t$は観測値撹乱項と呼ばれる．
状態$\alpha_t$は時点間の変動量を表す正規ホワイトノイズ$\eta_t$により
\alpha_{t+1}=\alpha_{t}+\eta_{t},\quad \eta_t\sim N(0,\sigma_\eta^2),\quad t=1,\cdots,n-1\tag{2}
と更新されていくと仮定する．状態の変動量$\eta_t$は状態撹乱項と呼ばれる．初期時点の状態$\alpha_1$の分布は
\alpha_t\sim N(a_1,P_1)\tag{3}
の正規分布に従うと仮定する．
ここまで式$(1),(2),(3)$により分布が与えられた$\epsilon_1,\cdots,\epsilon_n,\quad \eta_1,\cdots,\eta_{n-1},\quad \alpha_1$について互いが独立であることを仮定する．
これにより，撹乱項は過去の状態，観測値に対して独立となる．よって，状態$\alpha_t$，観測値$y_t$を含む確率変数の任意の線形結合の組は多変量正規分布に従う．
式$(1),(2),(3)$により定義された時系列モデルはローカルレベルモデルと呼ばれる．また，式$(1)$は観測方程式，式$(2)$は状態方程式と呼ばれる．
状態の推定と観測値の予測
状態空間モデルでは，与えられた時系列の観測値$y$に基づいて状態$\alpha$を推定し，将来の観測値を予測，欠測値の補間することが主な目的となる．ここでは，撹乱分散項，初期状態が従う正規分布の平均，分散が既知であると仮定する．
状態空間モデルの推定問題は大きく3つある．
1つ目は，途中時点$t$までの観測値$Y=\{y_1,\cdots,y_t\}$に基づいて当該時点の状態$\alpha_t$を推移するフィルタリング．
2つ目は，与えられた全ての観測値$y=\{y_1,\cdots,y_n\}$を用いて各時点の状態$\alpha_1,\cdots,\alpha_n$を推定する平滑化．また，欠測値の補間も平滑化の処理の中で行われる．
3つ目は，フィルタリングの最終時点$t=n$からさらに先の時点へと進めた，将来の状態，観測値に対する予測．
これらの処理は計算量が膨大になるため，計算時間的に実行不可能になることがある．
線形状態空間モデルのフィルタリングのための計算効率の高いアルゴリズムとしてカルマンフィルタがある．このアルゴリズムは逐次計算に基づくものである．
カルマンフィルタ
カルマンフィルタは与えられた初期状態$\alpha_1$の平均$a_1$，分散$P_1$から出発して，各時点$t$に対する状態$\alpha_t$の1期先予測$a_t=E(\alpha_t\mid Y_{t-1})$とその予測誤差分散$P_t=Var(\alpha_t\mid Y_{t-1})$，状態$\alpha_t$のフィルタ化推定量$a_{t\mid t}=E(\alpha_t\mid Y_{t})$とその推定誤差分散$P_{t\mid t}=Var(\alpha_t\mid Y_{t})$を交互に求めていく逐次計算アルゴリズムである．
ローカルレベルモデルでは，状態の1期先予測$a_t$がそのまま観測値の1期先予測$E(y_t\mid Y_{t-1})=E(\alpha_t+\epsilon_t\mid Y_{t-1})=a_t$となる．ここで，観測値$y_t$の1期先予測誤差$v_t=y_t-a_t$を1期先予測誤差と呼ぶ．さらに，$v_t$の条件付き分散を$F_t=Var(v_t\mid Y_{t-1})$と表し1期先予測誤差分散と呼ぶ．
既知である$a_1=E(\alpha_1)$，$P_1=Var(\alpha_1)$から$a_{1\mid1}=E(\alpha_1\mid y_1)$，$P_{1\mid1}=Var(\alpha_1\mid y_1)$を求めると$v_1$について
\begin{eqnarray}
E(v_1)&=&E(y_1-a_1)=E(\alpha_1+\epsilon_1-a_1)=a_1-a_1=0\\
F_1=Var(v_1)&=&Var(\alpha_1+\epsilon_1-a_1)=Var(\alpha_1)+Var(\epsilon)=P_1+\sigma_\epsilon^2\\
Cov(\alpha_1,v_1)&=&Cov(\alpha_1,\alpha_1+\epsilon_1-a_1)\\
&=&Var(\alpha_1)+Cov(\alpha_1,\epsilon_1)+Cov(\alpha_1,a_1)=P_1
\end{eqnarray}
($Cov(\alpha_1,\epsilon_1)=0$は$\alpha_t,\epsilon_t$は互いに独立，$Cov(\alpha_1,a_1)=0$は$a_1$が定数だから)
が成り立つ．
ここで，$v_1=y_1-a_1$は観測値$y_1$から定数$a_1$を引いただけのものなので，$E(\cdot\mid y_1)=E(\cdot\mid v_1)$が成り立つので多変量(2変量)正規分布の条件付き平均，条件付き分散と同様に求めることができ，
\begin{eqnarray}
a_{1\mid1}&=&E(\alpha_1\mid v_1)=E(\alpha_1)+Cov(\alpha_1,v_1)Var(v_1)^{-1}v_1\\
&=&a_1+P_1F_1^{-1}v_1=a_1+K_1v_1\\
P_{1\mid1}&=&Var(\alpha_1\mid v_1)=Var(\alpha_1)-Cov(\alpha_1,v_1)^2Var(v_1)^{-1}\\
&=&P_1-P_1^2F_1^{-1}=P_1(1-K_1)=P_1L_1
\end{eqnarray}
となる．ただし，$K_1=P_1/F_1, L_1=1-K_1$とおいた．
$a_{1\mid1},P_{1\mid1}$から$a_2,P_2$を求めるのは簡単である．(ローカルレベルモデルでは状態の1期先予測がそのまま観測値の1期先予測になるから)
\begin{eqnarray}
a_2&=&E(\alpha_2\mid y_1)=E(\alpha_1+\eta_1\mid y_1)=a_{1\mid1}\\
P_2&=&Var(\alpha_2\mid y_1)=Var(\alpha_1+\eta_1\mid y_1)=P_{1\mid1}+\sigma_\eta^2
\end{eqnarray}
これらを一般の時点$t=2,\cdots,n$について導出する．$Y_{t-1}=\{y_1,\cdots,y_{t-1}\}$が与えられた，$Y_{t-1}=\{y_1,\cdots,y_n\}$が与えられた下で，$v_t$について
\begin{eqnarray}
E(v_t\mid Y_{t-1})&=&E(\alpha_t+\epsilon_t-a_t\mid Y_{t-1})=a_t-a_t=0\\
F_t=Var(v_t\mid Y_{t-1})&=&Var(\alpha_t+\epsilon_t-a_t\mid Y_{t-1})\\
&=&Var(\alpha_t\mid Y_{t-1})+Var(\epsilon_t\mid Y_{t-1})=P_t+\sigma_\epsilon^2\\
Cov(\alpha_t,v_t\mid Y_{t-1})&=&Cov(\alpha_t+\epsilon_t-a_t\mid Y_{t-1})\\
&=&Var(\alpha_t\mid Y_{t-1})+Cov(\alpha_t,\epsilon_t\mid Y_{t-1})=P_t
\end{eqnarray}
が成り立つ．ここで，観測値$Y_{t-1}$が与えられた下で$a_t=E(\alpha_t\mid Y_{t-1})$は定数となるため，$y_t$と$v_t=y_t-a_t$は1対1に対応し，$E(\cdot\mid Y_t)=E(\cdot\mid v_t,Y_{t-1})$がいえる．よって，
\begin{eqnarray}
a_{t\mid t}&=&E(\alpha_t\mid v_t,Y_{t-1})\\
&=&E(\alpha_t\mid Y_{t-1})+Cov(\alpha_t,v_t\mid Y_{t-1})Var(v_t\mid Y_{t-1})^{-1}v_t\\
&=&a_t+P_tF_t^{-1}v_t=a_t+K_tv_t\\
P_{t\mid t}&=&Var(\alpha_t\mid v_t,Y_{t-1})\\
&=&Var(\alpha_t\mid Y_{t-1})-Cov(\alpha_t,v_t\mid Y_{t-1})^2Var(v_t\mid Y_{t-1})^{-1}\\
&=&P_t-P_t^2F_t^{-1}=P_t(1-K_t)=P_tL_t
\end{eqnarray}
が得られる．ただし，$K_t=P_t/F_t,L_t=1-K_t$とおいた．
これは$Y_{t-1}$が与えられた下での1期先の状態$\alpha_t$の推定値$a_t$が，新たな観測$v_t$が与えられたことにより$K_tv_t$だけ補正され，誤差分散$P_t$が$L_t=1-K_t$倍に縮まったと解釈できる．
$a_{t\mid t},P_{t\mid t}$が与えられた下で，$a_{t+1},P_{t+1}$は
\begin{eqnarray}
a_{t+1}&=&E(\alpha_{t+1}\mid y_t)=E(\alpha_t+\eta_t\mid Y_t)=a_{t\mid t}\\
P_{t+1}&=&Var(\alpha_{t+1}\mid y_t)=Var(\alpha_t+\eta_t\mid y_t)=P_{t\mid t}+\sigma_\eta^2
\end{eqnarray}
となる．
これらをまとめると，カルマンフィルタにおけるフィルタ化推定量と1期先予測，1期先予測誤差の逐次計算アルゴリズムは$t=1,\cdots,n$に対して
\begin{eqnarray}
v_t&=&y_t-a_t,\qquad F_t=P_t+\sigma_\epsilon^2,\\
a_{t\mid t}&=&a_t+K_tv_t,\quad P_{t\mid t}=P_tL_t,\\
a_{t+1}&=&a_{t\mid t},\qquad\quad P_{t+1}=P_{t\mid t}
\end{eqnarray}
ただし，$K_t=P_t/F_t,L_t=1-K_t$であり，$K_t$はカルマンゲインと呼ばれる．
また，上記の式を再帰的に計算すると，
\begin{eqnarray}
a_{t+1}&=&a_{t\mid t}=a_t+K_tv_t\\
&=&a_{t-1\mid t-1}+K_tv_t\\
&=&a_{t-1}+K_tv_t+K_{t-1}v_{t-1}\\
&=&\cdots=a_1+\sum_{j=1}^tK_jv_j
\end{eqnarray}
となる．