Rパッケージ紹介ばかりが続いていて恐縮ですが。。。最近になってこんなものがFacebookからリリースされていたのを知りました。
これはこれで使いやすそうだなと思ったんですが、実はGoogleからも同様のMCMCサンプリングベースの時系列分析向けCRANパッケージ{bsts}がしばらく前から出ていたりします。見た感じ日本ではほとんど知られていないように見受けられるのですが、どんなものなんでしょうか？
ということで、一応他社フレームワークの紹介をするよりはこちらのフレームワークの紹介を先にする方が筋かな*1と思ったもので（笑）、ここでは遅ればせながら{bsts}パッケージの紹介をすることといたします。
そもそも{bsts}とは
名前が示す通り、Bayesian Structural Time Series（ベイズ構造時系列）の時系列分析を行うパッケージです。愚かにもVignetteがないので詳しい説明自体がofficialには殆どないんですが、要は時系列データのベイジアンモデリングを回せるパッケージだと思えば合ってます。このパッケージの特筆すべき点は、パッケージ自体でMCMCサンプリングを行える点です*2。
作者のSteveからは「Stanなど汎用MCMCサンプラーを回すほど手間をかけるような話ではない時系列『予測』に使うべきなのが{bsts}」とのコメントを貰っていて、実際回帰なども気軽に突っ込めるのでよほど複雑な時系列データでない限りは{bsts}を使うと便利なんじゃないかと思ってます。
{ggplot2}を使い慣れている人には見慣れた感じかもしれませんが、モデルへの各成分の追加はstate.specification変数に後から後からAdd...という関数で付け足していく形で表現できます。回帰については最後のmodel関数のところで指定する方式になっています。また、xts形式データであれば日付カラムに基づいて週末・祝日効果もモデルに組み込むなどの関数も用意されていてなかなか使い勝手が良いです。
個々のAdd...関数が何を表現しているかについては、このブログでも以前取り上げた動的線形モデル＋状態空間モデルの記事で紹介した各用語がそのまま当てはまるので、こちらを参考にしていただいても良いかもしれません。
基本的には{forecast}パッケージへのアンチテーゼ？として作られたものゆえ、回帰そのものよりも予測にかなり特化した仕様になっています。学習データへの当てはまりを見たりするのにはあまり適していないようなのでご注意ください。
パッケージのサンプルデータセットで試してみる
同根のサンプルデータセットとしてGoogleの株価時系列データが入っているので、これを使って気軽に試すことができます。
> library(bsts) > data(goog) > train <- as.vector(goog)[1:1372] > test <- as.vector(goog)[1372:1428] > ss <- AddLocalLevel(list(), train) > ss <- AddSeasonal(ss, train, nseasons = 7) > model <- bsts(train, ss, niter=1000) # 略 > pred <- predict(model, horizon=56, burn=100) > plot(pred, ylim=c(250,800), xlim=c(1,1428)) > par(new=T) > plot(test, ylim=c(250,800), xlim=c(-1372,56), type='l', lwd=3, col='red', axes=F) > legend('topleft', legend=c('Predicted','Actual'), col=c('blue','red'), lty=1, lwd=3, cex=1.5)
株価予想なのにあまり合ってない気もするんですが、95% confidence intervalには収まっているので良しとしろってことなんでしょうかね（笑）。基本的にpredict関数が返すのはbsts.predictionというクラスで原則として予測のみ行う点に注意が必要です（つまり回帰の当てはめを見るには一手間要る）が、こんな感じで使えるということはお分かりいただけだかと思います。
実際に真のモデルが分かっているデータセットで試してみる
今度は回帰を含むモデリングをやってみましょう。ここでは以前Stanで時系列回帰をやったのとは少し異なるデータセットを使います。とは言っても、以下のまとめ記事の7番目で使ったものと同種のものです。
このデータセットは線形回帰部分とトレンドと周期7の季節調整を含むモデルから生成したものです*3。よってそのようにAddLocalLinearTrend及びAddSeasonal関数で設定しておきます。
> d <- read.csv('https://raw.githubusercontent.com/ozt-ca/tjo.hatenablog.samples/master/r_samples/public_lib/jp/sample_dlm_season_trend.csv', sep='\t') > d_train <- d[1:80,] > d_test <- d[81:100,] > library(bsts) > ss <- AddLocalLinearTrend(list(),d_train$cv) > ss <- AddSeasonal(ss, d_train$cv, nseasons = 7) > model <- bsts(cv~., state.specification = ss, niter=1000, data=d_train) # 略 > pred <- predict(model, burn=100, newdata=d_test[,-1]) > plot(pred, ylim=c(0, 1500), xlim=c(0,100), xlab='', ylab='') > par(new=T) > plot(d_test$cv, ylim=c(0,1500), xlim=c(-80,20), axes=F, type='l', col='red', lwd=3, lty=1, xlab='', ylab='') > legend('topleft', legend=c('Predicted','Actual'), col=c('blue','red'), lty=1, lwd=3, cex=1.5)
テストデータはちゃんと95% confidence intervalにも収まっていますし、MCMCサンプリングで得られた事後分布の「濃い」範囲にも収まっています。予測精度としてはなかなか悪くなさそうです。こんな感じで使えますよー、という宣伝でした（笑）。
