こんにちは。なるせです。
以前のエントリーで、状態空間モデルのパラメータ推定方法のひとつとしてマルコフ連鎖モンテカルロ(MCMC)を実践しました。

今回は別の方法として、カルマンフィルターによる推定をやってみます。
準備
想定するモデル
前回と同様、次の観測方程式\((1)\)と状態方程式\((2)(3)(4)\)をもつ状態空間モデルを扱います。
y_{t} &= \beta_{0,t}+\beta_{1,t}x_{1,t}+\beta_{2,t}x_{2,t}+\epsilon_{t} \tag{1} \\
\beta_{0,t} &= \beta_{0,t-1}+\epsilon_{0,t} \tag{2} \\
\beta_{1,t} &= \beta_{1,t-1}+\epsilon_{1,t} \tag{3} \\
\beta_{2,t} &= \beta_{2,t-1}+\epsilon_{2,t} \tag{4}
\end{align}
データ
# パラメタの設定 # サンプルサイズ N <- 200 # 観測誤差の標準偏差 observationErrorSD <- 20 # 切片の過程誤差の標準偏差 processErrorSD <- 10 # 係数1,2の過程誤差の標準偏差 coefErrorSD1 <- 0.5 coefErrorSD2 <- 0.2 # 各種シミュレーションデータの生成 # 説明変数をシミュレーションで作る set.seed(1) x1 <- rnorm(n=N, mean=10, sd=10) set.seed(2) x2 <- rnorm(n=N, mean=20, sd=5) # 係数1,2のシミュレーション set.seed(12) beta1 <- cumsum(rnorm(n=N, mean=0, sd=coefErrorSD1)) + 10 set.seed(13) beta2 <- cumsum(rnorm(n=N, mean=0, sd=coefErrorSD2)) + 5 plot(beta1, main="時間によるbeta1の変化", xlab="day") plot(beta2, main="時間によるbeta2の変化", xlab="day") # beta0のシミュレーション set.seed(3) beta0 <- cumsum(rnorm(n=N, mean=0, sd=processErrorSD)) plot(beta0, main="時間によるbeta0の変化", xlab="day") # responseのシミュレーション set.seed(4) response <- beta0 + beta1*x1 + beta2*x2 + rnorm(n=N, mean=0, sd=observationErrorSD) par(mfrow=c(3,2)) plot(x1, main="x1のシミュレーション結果", xlab="day",type="l") plot(x2, main="x2のシミュレーション結果", xlab="day",type="l") plot(beta0, main="beta0のシミュレーション結果", xlab="day",type="l") plot(beta1, main="beta1のシミュレーション結果", xlab="day",type="l") plot(beta2, main="beta2のシミュレーション結果", xlab="day",type="l") plot(response, main="responseのシミュレーション結果", xlab="day",type="l") par(mfrow=c(1,1))
パッケージ
カルマンフィルターによる回帰係数推定をやってみます。
KFASというパッケージを使用します。
変化する回帰係数の推定ということで、隼本から第5部-10章の広告効果の推定を参考にしました。
一部、説明変数が2つということで推定結果の取り出しに躓いたため、KFASパッケージのドキュメントを参照しました。
https://cran.r-project.org/web/packages/KFAS/KFAS.pdf
まず、必要なパッケージを使えるようにしておきます。
# KFASによる2時変係数モデルのパラメータ推定 -------------------- # 必要なパッケージ library(KFAS) library(ggplot2) library(gridExtra)
gridExtraパッケージはggplotによるグラフを並べて表示するために使います。
ggplotに対してはpar(mfrow())は使えないんですね。
モデルの構造と推定
モデルの構造を決める
まずはモデルの構造を決めます。
SSModelの中に、状態空間モデルの構造を記述します。
# Step1:モデルの構造を決める build_reg <- SSModel( H = NA, response ~ SSMtrend(degree = 1, Q = NA) + SSMregression( ~ x1 + x2 , Q = diag(rep(NA, 2))) )
Hは観測方程式の誤差分散、Qは状態方程式の誤差分散で、どちらも推定する対象なのでNAとします。
切片はSSMtrendで記述し、degree =1とすることでローカルレベルモデル扱いになります。degree = 2だとトレンド付きになります。
回帰部分はSSMregressionで記述します。説明変数が2つなので、行列の対角成分にNAを2つ並べます。
パラメータを推定する
次にパラメータの推定です。
4つの未知数(観測誤差の分散、切片と説明変数2つの分散)の初期値を1としてfitSSM関数にて推定します。
# Step2 パラメータ推定 fit_reg <- fitSSM(build_reg, inits = c(1, 1, 1, 1))
今回の場合、初期値は回帰係数1、回帰係数2、切片、観測誤差の順に並んでいます1。
この1行だけでfit_regの中に推定結果が格納されます。
結果の確認
それでは、推定された回帰係数を信頼区間付きで見てみましょう。
まずは切片から。
切片を取り出すにはpredict()関数でstatus ="level"を指定します。
# 変化する回帰係数の図示 ------------------------------------ # 信頼区間付きの回帰係数0 ----------------------------------- interval_beta0 <- predict(fit_reg$model, states = "level", interval = "confidence", level = 0.95) # データの整形 beta0_df <- cbind( data.frame(time = 1:200, beta0 = beta0), as.data.frame(interval_beta0) ) # グラフを格納 g1 <- ggplot(data = beta0_df, aes(x = time, y = beta0)) + labs(title="カルマンフィルターによるbeta0の推定値") + geom_point(alpha = 0.6, size = 0.9) + geom_line(aes(y = fit), size = 1.2) + geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.3)
後のためにグラフは変数に格納しました。
# 図示 g1
点がシミュレーションで生成した切片の推移、実線が平滑化された推定値、灰色の範囲が推定値の99%信頼区間です。
概ね正しく推定されているようです。
続いて回帰係数を見てみます。
隼本では広告効果の回帰係数をstatus = "regression"とすることで取り出していますが、今回のモデルで同様に指定すると2つの回帰係数を同時に取り出してしまうため別の指定方法が必要です。
ドキュメントを見ると次のような記載があります。
Which states are used in computing the predictions. Either a numeric vector containing the indices of the corresponding states, or a character vector defining the types of the corresponding states. Possible choices are "all", "level", "slope", "trend", "regression", "arima", "custom", "cycle" or "seasonal", where "trend" extracts all states relating to trend. These can be combined. Default is "all".
対応するインデックスで指定することができると書かれているので、回帰係数1はstatus = 1として取り出します。
また、格納される値は回帰係数と説明変数の積となっているため、説明変数で除する必要があります。
その際、説明変数の符号によって上限と下限が入れ替わってしまうため調整するコードを追加しています。
# 信頼区間付きの回帰係数1 ---------------------------------------------- interval_beta1 <- predict(fit_reg$model, states = 1, interval = "confidence", level = 0.95) # 説明変数で除して係数のみにする interval_beta1_adj <- data.frame( "fit" = t(apply(interval_beta1/x1,1,sort))[,2], "lwr" = t(apply(interval_beta1/x1,1,sort))[,1], "upr" = t(apply(interval_beta1/x1,1,sort))[,3] ) # データの整形 beta1_df <- cbind( data.frame(time = 1:200, beta1 = beta1), as.data.frame(interval_beta1_adj) ) # グラフを格納 g2 <- ggplot(data = beta1_df, aes(x = time, y = beta1)) + labs(title="カルマンフィルターによるbeta1の推定値") + geom_point(alpha = 0.6, size = 0.9) + geom_line(aes(y = fit), size = 1.2) + geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.3)
同様に回帰係数2はstatus = 2として取り出せます。
# 信頼区間付きの回帰係数2 ----------------------------------------------- interval_beta2 <- predict(fit_reg$model, states = 2, interval = "confidence", level = 0.95) interval_beta2_adj <- data.frame( "fit" = t(apply(interval_beta2/x2,1,sort))[,2], "lwr" = t(apply(interval_beta2/x2,1,sort))[,1], "upr" = t(apply(interval_beta2/x2,1,sort))[,3] ) # データの整形 beta2_df <- cbind( data.frame(time = 1:200, beta2 = beta2), as.data.frame(interval_beta2_adj) ) # グラフを格納 g3 <- ggplot(data = beta2_df, aes(x = time, y = beta2)) + labs(title="カルマンフィルターによるbeta2の推定値") + geom_point(alpha = 0.6, size = 0.9) + geom_line(aes(y = fit), size = 1.2) + geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.3)
さて、2つの推定値を見てみましょう。
gridExtra::grid.arrange(g2, g3, nrow = 2)
上手く推定できているようです。
MCMCによる推定との比較
MCMCによる推定値と比べるとどうでしょうか。
先のエントリーで実行したものをg4 g5 g6に格納し、並べて表示します。
gridExtra::grid.arrange(g1, g4, g2, g5, g3, g6, ncol = 2)
ほぼ同等の結果が得られました。実際の回帰係数からの乖離の仕方も同じように見えます。
2変数の動的線形モデルの場合、どちらの方法でも同じような結果が得られるため、計算の速いカルマンフィルターを使用したほうがいいのかもしれません。
しかし、全く異なる手法でこれだけ同じ結果が得られるって、統計学すごい。
- fit_reg$model$Qやfit_reg$model$Hで中身を調べると確認できます。↩
コメント