[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}
となる. 係数がゼロになるかならないか, という極端な選択になるので, 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 事前分布を使うならば原理上は MCMC で計算が可能である.

モデル平均化法

spike-and-slab によってベイズ的な方法で変数選択を行っているため, 係数パラメータは事後分布から得ることになる. ここで, パラメータを点推定するのではなく, 事後分布から複数回抽出するため, モデルが複数通りでき, よって予測値も複数通り得られ, ここから予測値の事後分布を推定することができる. この方法は Madigan & Raftery (1994, Model Selection and Accounting for Model Uncertainty in Graphical Models Using Occam’s Window) によって, 点推定で予測するよりも, 理論上, 優れていることが示されている*3.

bsts の利点

このように, ベイズ的な方法で時系列モデリングをするため, モデルの設定が柔軟になるという利点がある: MCMC は複雑な分布に対応できるため, 非ガウシアン状態空間モデルにも対応できる. また, 共軛事前分布を用いると, 事後分布が陽に得られるためギブスサンプラーが使えるから, 計算も比較的高速になる.

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() 関数を使うだけである*4.

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

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

実験

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

参考文献


  • George, E. I., & Mcculloch, R. E. (1997). Approaches for Bayesian Variable Selection. Statistica Sinica, 7, 339–373.
  • 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:が, 私はまだ精読していない……

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

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