Rのカルマンフィルタライブラリ、KFASの結果解釈備忘録

この本を読んで勉強していますがすぐに忘れてしまうので備忘録として．

線形ガウス状態空間モデルの状態方程式，観測方程式は，

$x_t = T_t x_{t-1} + R_t \xi_t, \qquad \xi_t \sim N(0, Q_t) \\ y_t = Z_t x_t + \epsilon_t, \qquad \epsilon_t \sim N(0, H_t)$

ローカルレベルモデル

ローカルレベルモデルの状態方程式，観測方程式は，
$\mu_t = \mu_{t-1} + w_t, \qquad w_t \sim N(0, \sigma^2_w) \\ y_t = \mu_t + \nu_t, \qquad \nu_t \sim N(0, \sigma^2_\nu)$

なので一般形との対応は，
$T_t = 1, \quad R_t = 1, \quad Z_t = 1 \\ x_t = \mu_t, \quad Q_t = \sigma^2_w, \quad H_t = \sigma^2_\nu$

ローカル線形トレンドモデル

ローカル線形トレンドモデルの状態方程式，観測方程式は，
$\delta_t = \delta_{t-1} + \zeta_t, \qquad \zeta_t \sim N(0, \sigma^2_\zeta) \\ \mu_t = \mu_{t-1} + \delta_{t-1} + w_t, \qquad w_t \sim N(0, \sigma^2_w) \\ y_t = \mu_t + \nu_t, \qquad \nu_t \sim N(0, \sigma^2_\nu)$

なので一般形との対応は，
$T = \begin{pmatrix} 1 & 1 \\ 0 & 1 \end{pmatrix}, \quad R_t = \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix}, \quad Z_t = \begin{pmatrix} 1 & 0 \end{pmatrix} \\ x_t = \begin{pmatrix} \mu_t \\ \delta_t \end{pmatrix}, \quad Q_t = \begin{pmatrix} \sigma^2_w & 0 \\ 0 & \sigma^2_\zeta \end{pmatrix}, \quad H_t = \sigma^2_\nu$

教科書第5部9章のサンプルを追試してみます．

mod <- SSModel(sales_train ~ SSMtrend(2, Q=c(list(NA), list(NA))), H=NA)
fit <- fitSSM(mod, rep(1, 3))
kfs <- KFS(fit$model, filtering = c("state", "mean"), smoothing = c("state", "mean"))

観測誤差の分散は約24となり信の値25に近い推定値が得られています．

> fit$model$H
, , 1

         [,1]
[1,] 23.63745

同様に，過程誤差の分散は，水準の変動が0.88（真の値1），トレンドの変動が0.00046（真の値0）となり，どちらも真の値に近い推定値が得られています．

> fit$model$Q
, , 1

          [,1]         [,2]
[1,] 0.8833319 0.0000000000
[2,] 0.0000000 0.0004648067

その他の要素も上記定式化の通りであることを確認できます．

> fit$model$T
, , 1

      level slope
level     1     1
slope     0     1

> fit$model$R
, , 1

      [,1] [,2]
level    1    0
slope    0    1

> fit$model$Z
, , 1

     level slope
[1,]     1     0

推定された平滑化状態はkfs$alphahatに格納されており，levelが水準（ $\mu_t$ ），slopleがトレンド（ $\delta_t$ ）です．

> head(kfs$alphahat, n=3)
        level     slope
[1,] 7.592715 0.1979556
[2,] 7.573414 0.1980699
[3,] 7.817400 0.1981601

ここで，教科書p.255にも記載がありますが，ここで推定したトレンドは状態の推定値であって観測値ではないことに注意が必要です．
上で定式化したように，観測値はこの累積和になっています．

kfs$aをみると最後（t=401）の時点での水準値とトレンド値を得られます．

> tail(kfs$a, n=2)
          level       slope
[400,] 47.12591 -0.09044255
[401,] 44.60478 -0.13952919

ここから，t=402の推定値は，
$y_{402} = \mu_{402} = \mu_{401} + \delta_{401} = 44.60478 - 0.13952919 = 44.46525$
と求まります．同様に，
$Y_{403} = 44.46525 - 0.13952919 = 44.32572$
となります．これらはpredictの結果と一致します．

> predict(fit$model, interval="prediction", level=0.9, n.ahead=3)
Time Series:
Start = 401 
End = 403 
Frequency = 1 
         fit      lwr      upr
401 44.60478 35.70135 53.50821
402 44.46525 35.39030 53.54020
403 44.32572 35.07563 53.57581

予測区間については，KFASのソースコードから，
$\hat{\mu} \pm qnorm((1 - level) / 2) * \sqrt{V_\mu + H}$
と求めているのは分かるのですが， $V_\mu$ が何なのかまだ理解できていません...

追記

こちらの本によると $V_\mu$ は平滑化状態分散というそうです．

ダミー変数を用いた周期成分

例えば週次変動を考慮する場合，6つのダミー変数を用意します．このとき一般式は次のようになります．
$T_t = \begin{pmatrix} 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & -1 & -1 & -1 & -1 & -1 & -1 \\ 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \end{pmatrix}$
$R_t = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix} \\ Z_t = \begin{pmatrix} 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \end{pmatrix} \\ x_t = \begin{pmatrix} \mu_t \\ \delta_t \\ \gamma_1 \\ \gamma_2 \\ \gamma_3 \\ \gamma_4 \\ \gamma_5 \\ \gamma_6 \end{pmatrix} \\ Q_t = \begin{pmatrix} \sigma^2_w & 0 & 0 \\ 0 & \sigma^2_\zeta & 0 \\ 0 & 0 & \sigma^2_\eta \end{pmatrix} \\ H_t = \sigma^2_\nu$