前回KFASでうまくいかなかったが、時変係数をランダムウォークとするのであれば、
最終時点の係数を予測に使えばいいわけで（期待値は不変だから）、
わざわざKFASを使わなくてもいいか、ということでdlmでやってみる。
100%こちらを参考に実践
説明変数を2変数に増やしたパターンでもちゃんと推定できるか確認。
パッケージさえインストールしていれば以下コピペでそのままグラフまで描ける。
↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓
# パラメタの設定
# サンプルサイズ
N <- 500
# 観測誤差の標準偏差
observationErrorSD <- 20
# 切片の過程誤差の標準偏差
processErrorSD <- 10
# 係数1の過程誤差の標準偏差
coefErrorSD1 <- 0.5
# 係数2の過程誤差の標準偏差
coefErrorSD2 <- 0.2
# 各種シミュレーションデータの生成
# 説明変数をシミュレーションで作る
set.seed(1)
explanatory1 <- rnorm(n=N, mean=10, sd=10)
set.seed(2)
explanatory2 <- rnorm(n=N, mean=20, sd=5)
# 各種シミュレーションデータの作成
# slopeのシミュレーション
set.seed(12)
slope1 <- cumsum(rnorm(n=N, mean=0, sd=coefErrorSD1)) + 10
set.seed(13)
slope2 <- cumsum(rnorm(n=N, mean=0, sd=coefErrorSD2)) + 5
plot(slope1, main="時間によるslope1の変化", xlab="day")
plot(slope2, main="時間によるslope2の変化", xlab="day")
# 各種シミュレーションデータの作成
# interceptのシミュレーション
set.seed(3)
intercept <- cumsum(rnorm(n=N, mean=0, sd=processErrorSD))
plot(intercept, main="時間によるinterceptの変化", xlab="day")
# responseのシミュレーション
set.seed(4)
response <- intercept + explanatory1*slope1 + explanatory2*slope2 + rnorm(n=N, mean=0, sd=observationErrorSD)
plot(response, main="responseのシミュレーション結果", xlab="day")
# 解析の準備
# ==================================================
# dlmパッケージによる、状態空間モデルの推定
# ==================================================
# パッケージの読み込み
library(dlm)
# モデルの推定
# Step1 モデルの型を決める
buildDlmReg <- function(theta){
dlmModReg(
X=cbind(explanatory1,explanatory2),
dV=exp(theta[1]),
dW=c(exp(theta[2]), exp(theta[3]), exp(theta[4]))
)
}
# モデルの推定
# Step2 パラメタ推定
fitDlmReg <- dlmMLE(
response,
parm=c(2, 1, 1, 1),
buildDlmReg,
method = "SANN"
)
# 結果の確認
# 推定されたパラメタの確認
fitDlmReg
sqrt(exp(fitDlmReg$par))
# モデルの推定
# 推定された分散を使って、モデルを組みなおす
modDlmReg <- buildDlmReg(fitDlmReg$par)
# モデルの推定
# Step3
# カルマンフィルター
filterDlmReg <- dlmFilter(response, modDlmReg)
# Step4
# スムージング
smoothDlmReg <- dlmSmooth(filterDlmReg)
# dlmの出力
modDlmReg
# ==================================================
# 推定結果の図示と確認
# ==================================================
# 推定されたresponseの状態
estimatedLevel <-
dropFirst(smoothDlmReg$s)[,1] +
explanatory1*dropFirst(smoothDlmReg$s)[,2] +
explanatory2*dropFirst(smoothDlmReg$s)[,3]
# 図示
par(mfrow=c(2,2))
# 元データ
plot(response, col=1, main="response", pch=21, xlab="day")
# 推定された状態
lines(estimatedLevel, col=4)
# 観測誤差の入っていない、正しい値
lines(intercept + explanatory1*slope1 + explanatory2*slope2, col=2)
# 凡例
legend("topleft", legend=c("スムージングの結果","正しい値"), lty=1, col=c(4,2))
# 推定された切片
plot(dropFirst(filterDlmReg$m)[,1], type="l", col=2, xlab="day", ylab="intercept")
lines(dropFirst(smoothDlmReg$s)[,1], type="l", col=4)
# 正しい切片
lines(intercept, col=1)
# 凡例
legend("topleft", legend=c("フィルタリングの結果","スムージングの結果","正しい値"), lty=1, col=c(2,4,1))
# 推定された傾き1
plot(dropFirst(filterDlmReg$m)[,2], type="l", ylim=c(0, 15), col=2, xlab="day", ylab="slope1")
lines(dropFirst(smoothDlmReg$s)[,2], type="l", ylim=c(0, 15), col=4)
# 正しい傾き
lines(slope1, col=1)
# 凡例
legend("topleft", legend=c("フィルタリングの結果","スムージングの結果","正しい値"), lty=1, col=c(2,4,1))
# 推定された傾き2
plot(dropFirst(filterDlmReg$m)[,3], type="l", ylim=c(0, 15), col=2, xlab="day", ylab="slope2")
lines(dropFirst(smoothDlmReg$s)[,3], type="l", ylim=c(0, 15), col=4)
# 正しい傾き
lines(slope2, col=1)
# 凡例
legend("topleft", legend=c("フィルタリングの結果","スムージングの結果","正しい値"), lty=1, col=c(2,4,1))
# 図を戻す
par(mfrow=c(1,1))
↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑
Rのコードはここまで
うまくいった。
次は季節性を加えたモデルをやってみる。