渋谷駅前で働くデータサイエンティストのブログ

元祖「六本木で働くデータサイエンティスト」です / 道玄坂→銀座→東京→六本木→渋谷駅前

Rでベイジアン動的線形モデルを学ぶ(4)：季節要素のあるローカルレベル・モデル

R DLM 時系列分析 統計学

色々と興味が発散していて違う話題ばかりしてますが、これもまだ全然終わってないので粛々と進めようと思います。ということで今回は季節調整のお話。Commandeur本の進行に合わせて、季節調整ありただしトレンドなしというモデルでいきます。もちろんテキストはこちらの2冊と、R実践例をまとめたPDF book。


状態空間時系列分析入門

状態空間時系列分析入門

  • 作者: J.J.F.コマンダー,S.J.クープマン,Jacques J.F. Commandeur,Sime Jan Koopman,和合肇
  • 出版社/メーカー: シーエーピー出版
  • 発売日: 2008/09
  • メディア: 単行本
  • 購入: 2人 クリック: 4回
  • この商品を含むブログを見る

Rによるベイジアン動的線形モデル (統計ライブラリー)

Rによるベイジアン動的線形モデル (統計ライブラリー)


前回同様、詰まったらPDF bookを参照しつつRコードの分からないところを補足しながら、という感じでだらだらやっていきましょう。


季節調整


単なる季節調整の話題なら以前stl関数でLOESSベースでやる方法を紹介してますが、ここではきちんとモデリングして季節調整する方法を紹介します。式としてはCommandeur本p.34に従い以下のように立てます。ちなみにこれは四半期データの場合です。


y_t = \mu_t + \gamma_t + \epsilon_t, \epsilon_t \sim NID(0,\sigma^2_\epsilon)

\mu_{t+1} = \mu_t + \xi_t, \xi_t \sim NID(0,\sigma^2_\xi)

\gamma_{1,t+1} = - \gamma_{1,t} - \gamma_{2,t} - \gamma_{3,t} + \omega_t, \omega_t \sim NID(0,\sigma^2_\omega)

\gamma_{2,t+1} = \gamma_{1,t}

\gamma_{3,t+1} = \gamma_{2,t}


ここで\gamma_t = \gamma_{1,t}は季節要素を表します。攪乱項\omega_tは季節要素が時間ごとにランダムに変動することを表します。そして初期値 \mu_1, \gamma_{1,1}, \gamma_{2,1}, \gamma_{3,1}は固定して、未知係数として最尤推定する対象になります。\gamma同士の式は確定的関係なので、攪乱項は入らないという解釈になっています。


Commandeur本p.34にもありますが、この式の形から分かるように季節周期がsの場合(s-1)本の状態方程式が必要とされます。なので、周期が長ければ長くなるほど状態方程式の本数も増えるという面倒なことになります。。。{dlm}なら簡単ですが、例えばWinBUGS / JAGS / Stanで同じようなモデリングをやる時には、相当に面倒くさいことになりそうです。


ちなみに、{dlm}でやる場合はdlmModSeas関数で簡単に季節要素を追加できます。なので、やり方としてはdlmModPolyでローカルレベル・モデルやローカル線形トレンドモデルを定義し、そこに単なる線形和としてdlmModSeas関数で季節要素を加える、という手順を取ることになります。


英国ドライバーのデータでやってみる


このシリーズ記事ご覧になってきた人なら、明らかに例のUK driverのデータに年ごとの季節要素（つまり周期12）があることに気付いているんじゃないかと思うわけで（笑）。ということでやってみます。Commandeur本では段階を踏んで色々やってますが、面倒なのでいきなり4.2節「確率的レベルと確率的季節要素」だけやります。4.1節と4.3節はスキップするので悪しからず。


とりあえずUK driverのデータをd1という名前でインポートしてあるとします（やり方はこのシリーズの過去記事を参照のこと）。後はPDF bookにもあるように以下の通りやるだけです。

> library(dlm)
> build1<-function(param){
+     model<-dlmModPoly(order=1)+dlmModSeas(frequency=12) # ローカルレベル・モデルと季節要素の線形和
+     V(model)<-exp(param[1]) # 元のdVでこれはモデル全体での値
+     diag(W(model))[1:2]<-exp(param[2:3]) # 同じく元のdWでこれもモデル全体での値
+     return(model) # 普通に返値としてモデルを返す
+ }
> fit1<-dlmMLE(d1,rep(0,3),build1) # 最尤推定する
> model1<-build1(fit1$par) # 推定したパラメータを与えてモデルを確定させる
> V(model1)
            [,1]
[1,] 0.003514183
# 統合されたdV
> W(model1)
              [,1]        [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
 [1,] 0.0009455539 0.00000e+00    0    0    0    0    0    0    0     0     0     0
 [2,] 0.0000000000 9.14866e-10    0    0    0    0    0    0    0     0     0     0
 [3,] 0.0000000000 0.00000e+00    0    0    0    0    0    0    0     0     0     0
 [4,] 0.0000000000 0.00000e+00    0    0    0    0    0    0    0     0     0     0
 [5,] 0.0000000000 0.00000e+00    0    0    0    0    0    0    0     0     0     0
 [6,] 0.0000000000 0.00000e+00    0    0    0    0    0    0    0     0     0     0
 [7,] 0.0000000000 0.00000e+00    0    0    0    0    0    0    0     0     0     0
 [8,] 0.0000000000 0.00000e+00    0    0    0    0    0    0    0     0     0     0
 [9,] 0.0000000000 0.00000e+00    0    0    0    0    0    0    0     0     0     0
[10,] 0.0000000000 0.00000e+00    0    0    0    0    0    0    0     0     0     0
[11,] 0.0000000000 0.00000e+00    0    0    0    0    0    0    0     0     0     0
[12,] 0.0000000000 0.00000e+00    0    0    0    0    0    0    0     0     0     0
# 季節要素なのでちゃんと12×12になっている
> diag(W(model1))[1:2]
[1] 9.455539e-04 9.148660e-10
# 統合されたdW、でもCommandeur本と微妙に違う気が。。。
> filt1<-dlmFilter(d1,model1)
# カルマンフィルタ
> smooth1<-dlmSmooth(filt1)
# カルマンスムージング
> sm1<-dropFirst(smooth1$s)
# 先頭のゴミを落とす
> mu1<-c(sm1[,1])
# 状態値
> nu1<-c(sm1[,2])
# 季節要素
> res1<-c(residuals(filt1,sd=F))
# モデル残差
> matplot(cbind(d1,mu1),type='l',lwd=c(1,3),lty=1,col=c('black','red'))
> plot(nu1,type='l')
> plot(res1,type='l')
# 各種プロット

f:id:TJO:20141120221741p:plain

f:id:TJO:20141120221829p:plain

f:id:TJO:20141120221838p:plain


大体図4.6, 4.7, 4.9っぽい図が。。。あれ、4.9っぽい図のつもりの3番目だけ変だな。何でだろう（汗）。ともあれ、これでとりあえず季節調整済みの状態値を得ることができました。


これで、UKドライバーのデータのうち季節要素による上下動まではほぼ説明できるようになりました。後は、最後の方に見られる「カックン」と下がっている謎の変動を説明するだけです。これまた後ろの章でのお楽しみ、ということで。


英国インフレーションのデータでやってみる


第4章では英国インフレーションというデータが初めて出てきます。これも後ろの章では結構面白いテーマの対象となっていくんですが、ひとまずここでは上の英国ドライバーのデータ同様にやってみましょう。原典の公式サイトから落としてきたデータの中の、Chapter4のフォルダに入ってます。一応d2という名前でインポートしておきましょう。


ちょっとプロットしてみれば分かりますが、これは周期4の季節データです。

f:id:TJO:20141120230223p:plain

ということで、以下のようにやってみます。

> build2<-function(param){
+     model<-dlmModPoly(order=1)+dlmModSeas(frequency=4)
+     V(model)<-exp(param[1])
+     diag(W(model))[1:2]<-exp(param[2:3])
+     return(model)
+ }
> fit2<-dlmMLE(d2,rep(0,3),build2)
> filt2<-dlmFilter(d2,model2)
> smooth2<-dlmSmooth(filt2)
> sm2<-dropFirst(smooth2$s)
> mu2<-c(sm2[,1])
> nu2<-c(sm2[,2])
> res2<-c(residuals(filt2,sd=F))
> matplot(cbind(d2,mu2),type='l',lwd=c(1,3),lty=1,col=c('black','red'))
> plot(nu2,type='l')
> plot(res2,type='l')

f:id:TJO:20141120230323p:plain

f:id:TJO:20141120230332p:plain

f:id:TJO:20141120230339p:plain


こんな感じで図4.10の3つのプロットと大体同じ結果になりました。2つ目のプロットに見えるように、季節要素の推定があまり良くないことが窺えます。ということで、ここをやはり後ろの章で改善していく感じになります。


最後に


そう言えばこの{dlm}での推定結果をStanのコードで再現しようとすると、ものによっては何故かうまくいかないことがあるんですよね。。。何でだろう。特に「確定的○○」の時にうまくいってない気が。これは@先生に伺わなければ（他力本願）。


追記


我らがStanユーザーのドン@さんがStanでのやり方をブログにupして下さいましたー。


ということで、これでStandでもやれそうです。@さん、有難うございました！