Rの状態空間モデルdlmパッケージの使い方

Share Button

素晴らしい諸先輩方の解説を見てぼくもなんとかここまで使えるようになりました。
基本的な使い方は下の参考のところに挙げた諸先輩方の解説ページや文献など読んでもらうとして,このページはあまりウェブで見つからなかったことを中心とした忘備録です。
 

モデリング

別段新しいことはありませんが,重回帰というか説明変数を複数投入してのdlmModRegの使用例が日本語であまり見当たらなかったので,まあなんかウェブへの知識の蓄積(笑)に貢献できればみたいな。
そこだけちゃんと書きますが,ちゃんと書くといってもただ説明変数のところを変数が並んだMatrixに変えるだけです。
2011年1月からの月次データdatみたいなものを仮定してですね。(元データを示さずにサンプルプログラムを公開するのってなんか不思議な感じですが,いい感じのが見つかったらまたそれでやります)
 
#準備
start=c(2011,1)
freq=12

#従属変数
dep<-ts(dat[,1],start=start,frequency=freq)

#説明変数
var1<-ts(dat[,2],start=start,frequency=freq)
var2<-ts(dat[,3],start=start,frequency=freq)
var3<-ts(dat[,4],start=start,frequency=freq)
vars<-cbind(var1,var2,var3)
 
dlmModRegでの説明変数のところにマトリックスのvarsをぶちこんで,あとは何列目から何列目まで使うのかを指定してあげればよろしいわけです。
いちいち列の数とか指定するのめんどいので勝手に列の数は読んでもらいます。
build <- function(u) {
  reg <- dlmModReg(vars[,(1:length(vars[1,]))],dV=exp(u[1]),dW=exp(u[2:(length(vars[1,])+2)]))
  return(reg)
}
outMLE<-dlmMLE(window(dep,start=start),parm=rep(0,(length(vars[1,])+2)),build,method="Nelder-Mead")
outMLE$convergence
mod<-build(outMLE$par)
outF<-dlmFilter(dep,mod)
outS<-dlmSmooth(outF)
 
できあがりですね。
 

パラメータはどう格納されている?

このモデルだと,outF$mやoutS$sにはIntercept, var1〜var3の順で回帰係数が出力されます。つまり,
intercept <-outS$s[,1]
beta1     <-outS$s[,2]
beta2     <-outS$s[,3]
beta3     <-outS$s[,4]
したがってモデリング結果の出力は
est<-window((intercept + beta1*var1 + beta2*var2 + beta3*var3),start=start)
ですね。
諸先輩方の解説にもあるように頭にゴミがつきますので,dropFirstなりで落としてあげてください。
ぼくはいつもwindowで切り取ってます。
 
元プロットと重ねてみると,
leg1<-rgb(0.8,0.5,0.1,alpha=0.6)
leg2<-rgb(0.70, 0.13, 0.13, alpha=0.8)
plot(dep,type = "l", ylab="", col = leg1, lwd=5)
lines(est,type="l", ylab="",lty=1, lwd=2, col=leg2)
legend("bottomleft",legend=c("Original","Model"),lty=c(1,1),lwd=c(5,2),col=c(leg1,leg2),bg = "transparent")
とまあこんな感じですね。
 

切片や特定の回帰係数のみ一定にしたい場合

まあ必ずしもすべての回帰係数が時変である必要性はそんなないってことも多いんじゃないですかね。むしろコンスタントなものも入れたいとか。そりゃそうだ。
モデル設定部分のdWのとこをいじるだけですね。初めのモデルでは下のような感じですが,
dW=exp(u[2:(length(vars[1,])+2)]))
これuを2から順番にvarsの要素数分だけ並べているわけですが,中身は
dW=c(exp(u[2]), #切片の分散
     exp(u[3]), #var1回帰係数分散
     exp(u[4]), #var2回帰係数分散
     exp(u[5])  #var3回帰係数分散
     )
って感じで並んでるので,var1だけを時変にして他のパラメータをコンスタントにする場合には
dW=c(0,    #切片の分散
     exp(u[2]), #var1回帰係数分散
     0,    #var2回帰係数分散
     0     #var3回帰係数分散
    )
まあこうなるわけですね。
outS$sのパラメータ見比べてみるとたぶん一定になってるんじゃないですかね。プロットすると初期に極小さく変動している場合がありますが。
 
これで自在にモデルが組めるじゃないですかやったー!
 

dlmForecastのこと

んでまあ当然予測がやりたくなるわけですが,このモデルでdlmForecast走らせると何が起きるかというと
Error: dlmForecast only works with constant models
はい使えなーい。
回帰係数が時変のモデルで使えなかったら意味なーい。
 
広い広い世界に目を向けるとdlmFilter()のカルマンフィルタでうまくやれるらしいんですが,今のところぼくは全くうまくいきませんので解説お待ちしております。
 

Lapack

あとこれ使ってるとちょいちょいLapack Routine Errorで止まるんですけど,その辺全然解明できてないです。モデリングでうまく収束しなかったとかそもそもモデルに無理があるとかそんな感じのようには見えるんですけど。
 
参考:
DLM カテゴリーの記事一覧 – 銀座で働くデータサイエンティストのブログ
dlmの使い方 | 時系列分析 | Logics of Blue
 ・Rによるベイジアン動的線形モデル (統計ライブラリー)  ・状態空間時系列分析入門  
広告:
ニート院生に愛の手(職)を差し伸べ世界からまたひとりの無職を救おうキャンペーンも展開しています。
 
 


Share Button

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です

This site uses Akismet to reduce spam. Learn how your comment data is processed.