FKFパッケージを使ってカルマンフィルタ（Kalman Filter）を実行する

カルマンフィルタ自体は以下の状態空間モデルを前提としている。

$\alpha_{t + 1}= d_{t}+T_{t}\alpha_ t + H_{t}\eta_t$ ...状態方程式
$y_{t}= c_{t}+Z_{t}\alpha_ t + G_{t}\eps_t$ ...観測方程式
（ $\eta_{t}$ , $\eps_{t}$ は標準正規分布に従う乱数）

観測可能（取得可能）なデータは $y_t$ で、それを使って状態方程式（状態変数 $\alpha_t$ ）の推定や $d_t$ , $c_t$ 等の各パラメーターを推定したり、状態変数の予測をしたりする。粒子フィルタの線型版。

なにはともあれRでカルマンフィルタを使うためのパッケージをインストールしておく。
（他にもsspirというのもあるけど、ぱっと見の使いやすさからこれ使ってる）

install.packages("FKF")

簡単な例として以下のようなモデルを考える。

$\alpha_{t + 1}= \alpha_ t + 2.0\eta_t$
$y_{t} = \alpha_ t + \eps_t$

このモデルに合わせて、データを作成。 $\alpha_0$ は0とした。 $P_0$ は $\alpha_0$ の推計精度を分散で指定する。ここではモデルパラメーター $H_t$ , $G_t$ が不明として、それを推定することを考える。（実際にはそれぞれ2と1）

library(FKF)
## alpha[t+1] = alpha[t] + eta[t], eta[t] ~ N(0, 2)
## y[t]       = alpha[t] + eps[t], eps[t] ~ N(0, 1)
alpha <- cumsum(2.0 * rnorm(100)) 
y     <- alpha + rnorm(100) 
dt <- matrix(0)
ct <- matrix(0)
Zt <- matrix(1)
Tt <- matrix(1)
a0 <- 0
P0 <- matrix(0.01)

$H_t$ , $G_t$ が不明なので、それを推定するためにカルマンフィルタから計算される尤度を最大にするような最適化を行う。最適化に使う目的関数を以下のように書いておく。

# objective function
objective <- function(par, ...){
  -fkf(HHt = matrix(par[1]), GGt = matrix(par[2]), ...)$logLik
}

HHt,GGtと書いている箇所がそれぞれ $H_t$ , $G_t$ の2乗、つまり標準偏差ではなく分散に対応する点に注意。

最適化はRに入ってるoptim関数を使う。HHt,GGtの初期値はそれぞれ適当に置く、今回は答えを知ってるのでそれにしておいた。ytには「1行×データ数分の列」の行列形式でデータを入れないといけないのでrbindしている。

## Estimate parameters:
fit.fkf <- optim(c(HHt = 4, GGt = 1), fn = objective,
  yt = rbind(y), a0 = a0, P0 = P0, dt = dt, ct = ct, Zt = Zt, Tt = Tt)

最適化した結果のパラメーターは$parで取得できる。

>fit.fkf$par
      HHt       GGt 
4.2356580 0.8431528

これの平方根を取ったものが、 $H_t$ , $G_t$ に対応する。

次に推計したパラメーターを使って、カルマンフィルタを実施した結果を取得する。

## Filter data with estimated parameters:
fkf.obj <- fkf(a0, P0, dt, ct, Tt, Zt, HHt = matrix(fit.fkf$par[1]),
  GGt = matrix(fit.fkf$par[2]), yt = rbind(y))

このfkf関数の返却値はfkfクラスのオブジェクト（S3）になっていて、それぞれの要素は
(n：データ数,m：状態数の次元,d：観測データの次元)

att：フィルターされた状態ベクトル $E(\alpha_t|y_t)$ 。m×n行列
at：予測状態ベクトル $E(\alpha_t|y_{t-1})$ 。m×(n+1)行列
Ptt：フィルターされた状態ベクトルの分散 $var(\alpha_t|y_t)$ 。m×m×n行列
Pt：予測状態ベクトルの分散 $var(\alpha_t|y_{t-1})$ 。m×m×(n+1)行列
vt：予測誤差( $y_t - c_t - Z_t E(\alpha_t|y_{t-1})$ )。d×n行列
Ft：予測誤差の分散。d×d×n行列
Kt：カルマンゲイン。m×d×n行列
logLik：尤度値
status LAPACKのdpotri と　dpotrf変数の状態が含まれている. (0; 0) だとうまくいったということ
sys.time：経過時間

という意味合い。今回の場合だと

>str(fkf.obj)
List of 10
 $ att     : num [1, 1:100] 2.75 2.41 2.55 2.42 1.63 ...
 $ at      : num [1, 1:101] 1.1 2.75 2.41 2.55 2.42 ...
 $ Ptt     : num [1, 1, 1:100] 1.63 1.22 1.19 1.19 1.19 ...
 $ Pt      : num [1, 1, 1:101] 100 4.66 4.26 4.23 4.22 ...
 $ vt      : num [1, 1:100] 1.678 -0.457 0.183 -0.175 -1.097 ...
 $ Ft      : num [1, 1, 1:100] 101.66 6.32 5.91 5.88 5.88 ...
 $ Kt      : num [1, 1, 1:100] 0.984 0.738 0.72 0.718 0.718 ...
 $ logLik  : num -229
 $ status  : int [1:2] 0 0
 $ sys.time:Class 'proc_time'  Named num [1:5] 0 0 0 NA NA
  .. ..- attr(*, "names")= chr [1:5] "user.self" "sys.self" "elapsed" "user.child" ...
 - attr(*, "class")= chr "fkf"

というようになっている。

推定した状態変数と元データを重ねてプロットしてみると

## Plot data together with filtered value:
plot(y)
lines(fkf.obj$att[1, ],col = "blue")
legend("bottomright", c("original data", "estimation"), col = c("black", "blue"), lty = 1)