ill-identified diary

所属組織の見解などとは一切関係ありません

[R] bsts パッケージの使い方

概要

  1. ベイジアン構造時系列モデリングを行う bsts パッケージは最近リリースされたばかりである.
  2. このパッケージはまだ情報が少ないため, ここで包括的に説明することを試みる

某有名IT企業の某有名データサイエンティストの方が bsts パッケージの入門記事
tjo.hatenablog.com
を書いていたが, より詳細で広範な説明がほしいところであった. 奇しくも bsts の開発者もこの方と同じ企業に所属しているので, 何か特別な思い入れがあるのかとも思ったが, 本人はなかなかより詳細な記事を書く気配がないため, しびれを切らして自分が書くことにした (あてつけではない).

この記事の構成は次の通り. まず, 第2セクションではベイジアン構造時系列モデルとは何か, 理論的な解説をする. 実際にはこのセクションはほぼ S. Scott & Varian (2013, Bayesian Variable Selection for Nowcasting Economic Time Series), S. L. Scott & Varian (2014, Predicting the present with Bayesian structural time series) の要約である. 第3セクションでは, 実際に bsts パッケージを操作できるように構文を解説する. 第4セクションでは, 具体例に基づいて実演する.

bsts の基本理論

bsts 開発のモチベーションは, 時系列予測の中でも, 随時, 短期間の将来の状態を予測する “nowcasting” の実現にある. 従来的な時系列予測と異なるのは, 時系列の多さと, それとは相対的にデータの期間が短い*1ことである. 膨大なデータから予測モデル構築に有効な変数を全て人の手で取捨選択するのは困難である*2. このような状況で, 予測モデルを構築するにあたって重要な機能は

  1. 時系列モデル特有の要素 (定常トレンドや周期トレンドの要素)

  2. 変数選択

  3. モデル平均化法 (model averaging)

の3つである.

時系列モデル特有の要素

(1) は, 状態空間モデルをカルマンフィルタで計算すれば, 容易である. 状態空間モデルを用いれば, ローカルレベルトレンドや季節性トレンドと回帰モデル(係数が変動する dynamic regression を含む) をミックスしたモデルもカルマンフィルタで計算が可能になる. これは以前にいくつかの記事で言及したので, これ以上の説明は省略する.
ill-identified.hatenablog.com
ill-identified.hatenablog.com

変数選択

状態空間モデル,


\begin{align}
\boldsymbol{y}_t &= \mathbf{Z}_t\boldsymbol{\alpha}_t + \boldsymbol{\varepsilon}_t,\\
\boldsymbol{\alpha}_{t+1} &= \mathbf{T}_t\boldsymbol{\alpha}_t + \mathbf{R}\boldsymbol{v}_t
\end{align}

の上段の式が観測方程式, 下段が遷移方程式 (状態方程式) である. 観測方程式の係数行列 \mathbf{Z}_t は, 時間によって変化させることもできる. つまり, 毎期 t ごとに値が変わってもモデルの計算は可能である. よって, \mathbf{Z}_t の成分に, 毎期で値の変わる説明変数を与えることで, 回帰モデルとして扱うこともできる. この場合, 対応する状態変数 \boldsymbol{\alpha}_t の成分が回帰係数に対応する. 状態変数も値を変化させることができるが, この変数選択の場面では, 全期間を通して値が一定の場合のみを考える.

しかし, 先に書いたように, 今は説明変数の候補が非常に多い状況を考えている. 無意味な説明変数を増やすことは計算の困難さや過剰適合 (オーバーフィット) にもつながるので, 必要最小限の変数だけを選択しなければならない. 回帰モデルの変数選択といえば, Lasso のような罰則付き回帰, あるいは主成分回帰 (PCR) など, 様々な方法があるが, bstsGeorge & Mcculloch (1997, Approaches for Bayesian Variable Selection) によって提案された spike-and-slab 事前分布を使うというベイズ的な方法を採用している. 説明変数の係数 \beta_{i} について, \begin{align}
\gamma_{i}:= & 1(\beta_{i}\neq0)\end{align} という指示変数を用意する. \gamma_{i} は係数 \beta_{i} がゼロでないなら1, ゼロならゼロとなる. \beta_{i} がゼロでない説明変数だけが, 実際に使われることになるから, まずは\gamma_{i} の事前分布を考える. \gamma_{i}=1 となる確率を \pi_{i} とおくと, \boldsymbol{\gamma}:=\begin{bmatrix}\gamma_{1} & \cdots & \gamma_{K}\end{bmatrix} の事前分布は

\begin{align}
p(\boldsymbol{\gamma};\boldsymbol{\pi})\sim & \prod_{i}\pi_{i}^{\gamma_{i}}(1-\pi_{i})^{1-\gamma_{i}}\end{align}
となる.  \pi_i を変更することで, 全て変数のうちだいたい何割が採用されるかを簡単に調整できる. 係数がゼロになるかならないか, という極端な選択になるので, spike-and-slab の spike (突起) 部分である. 上記の式はベルヌーイ分布を使っているが, MCMC を使えば, ベルヌーイ分布以外の分布であっても計算できる. \boldsymbol{\gamma} を所与とすると, あとは非ゼロの係数 \beta_{i} だけを考えれば良い. これを \boldsymbol{\beta}_{\gamma} と表すと, (条件つき) 事前分布を
\begin{align}
\boldsymbol{\beta}_{\gamma}\mid\boldsymbol{\gamma},\sigma^{2}\sim & \mathcal{N}(\boldsymbol{b}_{\gamma},\sigma^{2}\boldsymbol{\Omega}_{\gamma}),\\
\sigma^{2}\sim & \Gamma^{-1}(\frac{\mathit{df}}{2},\frac{\mathit{ss}}{2})\end{align}

と設定する. ただし, \Gamma^{-1} は逆ガンマ分布を表す. この事前分布は平坦な弱情報分布になるので, spike-and-slab の slab (平板) 部分である. 以上を合成した,

\begin{align}
p(\boldsymbol{\beta},\boldsymbol{\gamma},\sigma^{2})= & p(\boldsymbol{\beta}_{\gamma}\mid\boldsymbol{\gamma},\sigma^{2})p(\sigma^{2}\mid\boldsymbol{\gamma})p(\boldsymbol{\gamma})\end{align}
が spike-and-slab 分布である. \boldsymbol{\beta}_\gamma, \boldsymbol{\Omega}_\gamma \gamma で添字されているように, \gamma の値に依存して決まる. \gamma=0 ならば, 平均ゼロ, 分散もゼロかゼロに近い値となり, 係数がゼロと推定されやすくなる*3. これが係数の事前分布に用いられる. このような事前分布を採用することで, 説明変数のうち, 適合しないものは係数がゼロになりやすくなり, 適合する説明変数の係数はよくある裾の厚い事前分布になる.

(2018/8/9: spike-slab 事前分布の意義についてよくわからない説明していたので修正)

ではなぜ, よく用いられる L1, L2 正則化を使わず, この事前分布を採用したかだ. 機械学習の L1, L2 正則化は, それぞれベイズ統計ではラプラス分布と正規分布を事前分布に設定するのと同じで, 特に L1 正則化 (=Lasso) は変数選択についてよく似た性質を持つし, 計算ももっと簡単なはずだ. なぜこのような複雑な事前分布を考えたのか. 本文では L1 正則化に比べてどこが優れているかを明示していない*4ので, 私なりの考えを述べる. まず1つは, 既に書いたようにハイパーパラメータ  \pi_i の調整でどの程度の変数が排除されるのか直感的にわかりやすいという点があるだろう. ここでは変数が観測数よりも多いような状況も想定しているので, スケールパラメータをどの程度増減させれば変数がどれだけ排除されるかすぐにわからない L1 正則化よりも使い勝手がよい, ということなのだと思う. 一方で, 一旦採用されてしまえば, L2 正則化と同様に,「緩やかな」縮小推定がなされ, 過剰に係数の大きさが抑制されることがなくなる. しかし, このような事前分布では, どの変数が選ばれるかが確率に依存しやすく, モデルの予測値にも大きなばらつきが出ると思われる. その点については, 次に紹介するモデル平均化法で対処している.

モデル平均化法

(2018/8/10: モデル平均化の意味を最近まで分かっていなかったことを告白しつつ, 説明を刷新)

bsts ではベイジアンモデル平均化 (BMA) というテクニックが使われている. これについて本文で参照されているページは現在は存在しなかったので, 関連が深いと思われる Hoteling et al. (1999, Bayesian Model Averaging: A Tutorial) を参考にする*5. BMA は, 要するに機械学習でいうアンサンブル学習と同じアイディアである. バイアスは小さいがバリアンスの大きい複数のモデルの予測値を平均し, バリアンスを低減することで予測精度を上げるというテクニックは, ランダムフォレストやバギング, スタッキングなどのアンサンブル学習のテクニックと同じである.

今回のモデルは予測分布をモンテカルロサンプリングで得られるので, 1つの y_{t+1} に対して複数の予測値を与えられる. しかし BMA はそういう話ではなく, モデル自体を複数作成することを考える. ベイズモデルは, モデルの事後確率*6を定義できるため, これをウエイトにして各モデルの予測値の加重平均を取ることで, 単純平均よりも当てはまりを良くしている. これが BMA の基本アイディアである.

ただし実際には, あまりに性能の悪いモデルが含まれていると平均化は上手くいかない. Madigan & Raftery (1994) では, 相対的に確率の低い値で, かつ一定の水準以下のモデルを排除する, 複雑度のより大きいモデルを排除する, といった情報量規準に似たルールでモデルを取捨選することが提案されている. さらに別の方法として, メトロポリス-ヘイスティングスアルゴリズムを使ってモデルを複数サンプリングすることで, モデル平均を直接計算する MC3 アルゴリズムというものも提案されている*7.

bsts の利点

ベイズ統計モデルとして定式化すると, 事後分布が複雑であっても計算できるマルコフ連鎖モンテカルロ法 (MCMC) で対応できる. MCMC は計算が遅いと言われているが, 一部のパラメータに共軛事前分布を用いて事後分布が陽に得られるなら, ギブスサンプラーが使えるから, 計算も比較的高速になる. spike-and-slab 分布は一見複雑に見えるが, 実は一部をギブスサンプラーでサンプリングできる.

bsts の構文と使用法

状態空間モデルとして記述されたモデルを扱うため, 状態空間モデルの特徴を記録するオブジェクトが必要である. bsts では list オブジェクトが用いられ, 記述を追加するには Add で始まる各種関数を用いる. 例えば, ローカルレベルモデルならば,

ss <- AddLocalLevel(state.specification=list(), data$y)

となる. これで ss という list オブジェクトに, ローカルレベルモデルの記述がなされた. AddLocalLevel() 関数の第1引数は, 状態空間モデルを指定するオブジェクトであるが, ここでは1から作ったため, 空の list を与えている. 既存のリストオブジェクトにローカルレベル要素を追加したい場合は, ここに追加先のオブジェクトを指定する. 第2引数 data$y は観測された時系列データである. この ss に4期間周期の季節性トレンドを追加するならば,

ss <- AddLocalLevel(state.specification=list(), data$y)
ss <- AddSeasonal(ss, data$y, nseasons=4)

というふうに記述する. 観測データは毎回与える必要がある. その他にも, 休日効果を追加する関数や, trigonometric trend を追加する関数もある. 変化点を与えるものは今のところないようだ.

記述が終われば, 次に予測を行う. 先に書いたように, bsts パッケージは完全にベイズ的な方法でなされるため, dlmkfas と違いパラメータ推定と予測が bsts() という1つの関数で行われる. bsts() では, 観測方程式の分布族や MCMC の試行回数, (必要ならば) 回帰する説明変数の指定などを指定できる. なお, ここで指定した説明変数には, 自動で spike-and-slab による変数選択が適用される.

bsts() 関数の返した結果の格納された bsts オブジェクトを, plot() 関数に与えることで結果を様々なグラフで表示できる. たとえば, "coefficients" または "coef" オプションを追加した場合は, spike-and-slab によって選別された説明変数が, モデルに含まれる確率がプロットされる (棒グラフに影がついているものは係数がマイナスであることを意味する.).

f:id:ill-identified:20170907235529p:plain

さらに, "components" または "comp" ならば時系列をローカルトレンド, 季節トレンド, といった部品に分解したものがプロットされる.

f:id:ill-identified:20170907235412p:plain

予測を求めたいなら, これもシンプルで, predict() 関数と plot() 関数を使うだけである*8.

f:id:ill-identified:20170907235511p:plain

さらに, 複数のモデルを用意した場合, CompareBstsModels() 関数で複数のモデルの誤差を比較できる.

実験

現行バージョンでは, 周期トレンドはカレンダー効果と固定周期のトレンドしか指定できないため, 日次単位のデータや四半期単位のデータは有効だが, 時間単位でありながら数カ月あるような長期データなどを扱うのは難しい. また, 前回の動的市場反応モデルも, どの変数を使うべきかが理論モデルからある程度決まっているので, bsts の恩恵があまりない. この要件をみたす実験用データを探すのは面倒なので, 結局, 公式の解説で使われているのと同じデータを用いることにした. このデータには米国の毎週の失業者数*9が記録されている. 失業者数のような政府統計はふつう, まず不正確な速報が発表されてから, しばらくして確報が発表される. そこで, より早く最新の失業者数を知るために, 失業者数の予測モデルを作るモチベーションが発生する. 失業者数の発表は木曜日で, グーグルのログは最新でも2日前のものしか利用できない.

公式では, 3種類のモデルを設定し比較しているが, 元記事そのまんまでは面白くないので, 適当に2つほどモデル候補を追加してみる.

  1. ローカルレベル + 季節トレンド(52週間(=1年)で1周期のトレンド)

  2. モデル (1) + 説明変数

  3. モデル (2) の係数の事前分布を緩和 (係数がゼロになりにくい)

  4. (2) に2階の自己回帰項 AR(2) を追加する (係数には spike-and-slab が適用される).

  5. (4) の回帰を dynamic regression にする (spike-and-slab は適用できない).

プログラムは以下のとおりである.


gist4bb6fdc0714122bec7b12d2f1422e7c8


モデル (5) は, 係数が時変パラメータであるため, "coef" が使えない. かわりに "dynamic" ("dyn") オプションをつけて平滑化推定量をプロットする.

f:id:ill-identified:20170908201640p:plain

CompareBstsModels() を使った各モデルの比較結果は以下のようになる. 2つの時系列グラフが出力され, そのうち, 上段は各モデルの1期先予測値の誤差の絶対値の累積値で, 下段は観測値のスケールを調整した時系列である. このグラフにより, 時系列のどのあたりで予測精度が悪くなっているか, その原因は原系列の構造変化ではないか, といったことを端的に知ることができる.

f:id:ill-identified:20170908000348p:plain

モデル (4) は元記事で提案された (2), (3) の精度に近いが, モデル (5) は, 最初は誤差が小さいものの, 2009年 (おそらくリーマン・ショック) を境に誤差が積み重なっている. おそらく過剰適合しているのだろう.

蛇足: 状態空間モデルつながりで, KFAS についてもそのうち解説記事を書こうかと思う. 同じく状態空間モデルで有名な dlm は既に のように充実した解説が存在する.

追記: KFASパッケージの解説を書いた.
ill-identified.hatenablog.com

参考文献


  • George, E. I., & Mcculloch, R. E. (1997). Approaches for Bayesian Variable Selection. Statistica Sinica, 7, 339–373.
  • Hoteling, J. A., Madigan, D., Raftery, A. E., & Volinsky, C. T. (1999). Bayesian Model Averaging: A Tutotial. Statistical Science, 14(4), 382–401.
  • Madigan, D., & Raftery, A. E. (1994). Model Selection and Accounting for Model Uncertainty in Graphical Models Using Occam’s Window. Journal of the American Statistical Association, 89(428), 1535. doi:https://doi.org/10.2307/2291017
  • Scott, S. L., & Varian, H. R. (2014). Predicting the present with Bayesian structural time series. International Journal of Mathematical Modelling and Numerical Optimisation, 5(1/2), 4. doi:https://doi.org/10.1504/IJMMNO.2014.059942
  • Scott, S., & Varian, H. (2013). Bayesian Variable Selection for Nowcasting Economic Time Series. NBER Working Paper Series, 1–22. Retrieved from http://www.nber.org/papers/w19567.pdf

*1:Google の検索ログは2004年からしかない.

*2:これを著者らは 'fat regression' と呼んでいる.

*3:といっても, この説明からすると, \beta_\gamma はどっちにしろゼロになる?

*4:George & Mcculloch (1997) でも言及されていない

*5:なお, ぐぐると『パターン認識機械学習』の勉強会スライドがいくつもヒットするので, こちらでもBMAが解説されているようだ.

*6:ベイジアン情報量規準 (BIC) はこの確率の対数を近似している

*7:アイディアが直球すぎて笑ってしまった.

*8:ただし, 回帰モデルを含む場合は, 予測のために将来の説明変数データも必要になる

*9:正確には initial claims number; 失業者による初回の失業給付申請件数. また, データセットに収録されている数値は未加工のカウントデータではなく標準化したものである