前の記事でもリンクさせていただいているが、サイト 「状態空間時系列分析入門」をRで再現する では以下のテキストを {dlm}, {KFAS} で再現されており非常にありがたい。これらのライブラリの使い方については リンク先を読めば困らない感じだ。

自分も勉強のために似たことやりたい、、でも同じことやるのもなあ、、と考えた結果 同テキストの内容 {rstan} を使ってやってみた。

補足 Stan には状態空間表現用の関数 gaussian_dlm_obs_log ( 利用例 ) があるのだが、自分は使ったことがない。7章までのモデルは全て漸化式で表現されているため、それらを Stan のモデルとして記述した。

github リポジトリ

各スクリプト / モデルは GitHub にあげた。現在 7 章までのモデルと kick する R スクリプト, markdown ファイルを置いている。

github.com

モデルの一覧

各章のモデルについて、その実行結果を RPubs にあげている。各モデルの概要とリンクを以下に記載する。自分のモデリング能力不足のため いくつか課題があると考えており、特に怪しいものには * をつけた。

• 第1章 はじめに
  • 線形回帰
• 第2章 ローカル・レベル・モデル
  • 確定的レベル
  • 確率的レベル
  • ローカル・レベル・モデルとノルウェイの事故
• 第3章 ローカル線形トレンド・モデル
  • 確率的レベルと確率的傾き *
  • 確率的レベルと確定的傾き
  • ローカル線形トレンド・モデルとフィンランドの事故
• 第4章 季節要素のあるローカル・レベル・モデル
  • 確定的レベルと確定的季節要素
  • 確率的レベルと確率的季節要素
  • ローカル・レベルと季節モデルと英国インフレーション
• 第5章 説明変数のあるローカル・レベル・モデル
  • 確定的レベルと(確定的)説明変数
  • 確率的レベルと(確定的)説明変数 *
• 第6章 干渉変数のあるローカル・レベル・モデル
  • 確定的レベルと(確定的)干渉変数
  • 確率的レベルと(確定的)干渉変数
• 第7章 英国シートベルト法とインフレーション・モデル
  • 確定的レベルと確定的季節要素 *
  • 確率的レベルと確率的季節要素 *
  • 英国インフレーション・モデル

かんたんな説明

必要パッケージ

まず、実行には以下のCRAN非公開パッケージが必要である。後者はどうでもよいが {pforeach} はいろいろ捗るのでインストールしていない方はただちに入れるべきだ。

• hoxo-m/pforeach: Stan 並列化に利用。
• sinhrks/ggfortify: 結果のプロットに利用。

データの読み込み

データは著者のサポートサイト、An Introduction to State Space Time Series Analysis から入手できる。各データは適宜 ts インスタンス化しておく。

Stan の実行

Stan の呼び出しは全て以下のように並列化して行う。

model_file <- 'モデルを記述したファイルのパス'

# モデルのコンパイルを実行
stan_fit <- stan(file = model_file, chains = 0)

# Stan を並列実行する。実行結果に sflist2stanfit を適用し、stanfit インスタンスを得る
fit <- pforeach(i = 1:4, .final = sflist2stanfit)({
  stan(fit = stan_fit, data = standata, chains = 1, seed = i)
})

また、モデルが収束したかどうかは 全ての Rhat が 1.1 未満となっているかどうかで確認する。そのままでは収束しない / しにくいモデルには、初期値 / 事前分布 / 制約を入れて収束させた。

is.converged <- function(stanfit) {
  summarized <- summary(stanfit)  
  all(summarized$summary[, 'Rhat'] < 1.1)
}

stopifnot(is.converged(fit))

また、テキストに最尤推定の結果が記載されている場合は、以下のようにして近い値になっているかを確認した。値がずれる場合は stopifnot を外し、差異を出力する。

# is.almost.fitted の定義は各リンク先を参照
stopifnot(is.almost.fitted(mu[[1]], 7.4157))

結果のプロット

テキストのグラフを再現する形で行う。使っているパッケージの説明はこちら。

• RPubs - Plotting Time Series with ggplot2 and ggfortify
• RPubs - Plotting State Space Time Series with ggplot2 and ggfortify

課題

以下の2つについてはもう少し Stan に習熟できたら直したい。

複数の確率状態をもつモデルをうまく収束させたい

このようなモデルは 収束しにくく、また 最尤推定した状態推定値 ( テキスト、{dlm}、{KFAS}の結果 ) とずれやすい。特に 現状 制約を入れて収束させているモデルについて、事前分布 / 初期値 / 尤度関数の範囲でなんとかしたい。

テキストの尤度関数を再現したい

テキストに記載の式 (概要は以下記事参照) に揃えようとすると、自分の記述では収束しにくくなってしまう。

sinhrks.hatenablog.com