Osmundsen, K.K., Kleppe, T.S., Ogland A. (2019) MCMC for Markov-switching models—Gibbs sampling vs. marginalized likelihood. Communications in Statistics – Simulation and Computation.
仕事の都合で、マルコフ・スイッチングVARモデルについてあれこれ考えているのだけれど、組みたいモデルを推定できるソフトがなかなか見当たらない。RではかつてMSBVARというパッケージがあったのだけれど、開発が止まったようで、CRANからもremoveされている。
これは自力で書けっていうこと？ 到底そんな学力はない。Stanでどうにかする？ でも離散潜在変数ってめっちゃ難しいんじゃなかったっけ… と思い悩みながらwebの海をさまよっていたら、Stanでの実装を示している人がいた。いわく、詳細はこの論文を読め。はい、読みます。
1. イントロダクション
[…マルコフ・スイッチング・モデルについての説明…]
Sims & Xha (2006) はマルコフ・スイッチング・ベクトル自己回帰(MS-VAR)モデルを使って財政政策の構造的ブレークについて検討した。またBloom(2009)は経済活動の不確実性ショックをモデル化した。経済学ではこれらがきっかけでMS-VARモデルが使われるようになり、[…いくつか紹介]。
さて、MS-VARモデルの推定方法には、一般的に使われるGibbsサンプリングと、周辺尤度のMCMC推定とがある。隠れたレジーム変数をサンプリングするか、周辺化して尤度関数から消去するかの違いである。後者は状態空間モデルでよく使われている擬似周辺尤度みたいなものである。MHやHMCのようなより一般的なMCMC手続きを使うことができる。
本論文ではどっちの効率が良いかを調べる。
2. 方法論
2.1 マルコフ・スイッチング・ベクトル自己回帰モデル
話を簡単にするために、以下ではずっとラグ1で考えますね。観察を\(\mathbf{Y}_t \in \mathbb{R}^d, t \in (1,2,3,…n)\)とする。時点2以降について以下とする。$$ \mathbf{Y}_t = \phi_{(S_t)} \mathbf{Y}_{t-1} + \mathbf{\mu}_{(S_t)} + \epsilon_t$$ $$ \epsilon_t \sim N(0, \Sigma_{(S_t)}) $$ \(S_t \in (1,2,\ldots,m)\)は潜在状態である。\(\phi_{(S)}\)は自己回帰係数行列、\(\mathbf{\mu}_{(S)}\)は平均ベクトル。
状態は時間等質な一次マルコフ連鎖に従うものとする。遷移行列を\(\mathbf{Q}\)とする。遷移行列というのは\(m \times m\)の正方行列で、要素は\(p_{ij} = p(S_t = j | S_{t-1} = i)\)ね。
\(\{S_t\}\)は既約で非周期的であり、よって定常分布\(\delta\)を持つとする。最初の観察は既知としモデル化しない。
[ああ、なんてクリアな定式化だろうか。すべての論文がこんな感じであってほしい。というか、この前に読んだVLSTARってやつ、あれわかりにくかったよな。いま気が付いた]
尤度関数について考えよう。パラメータをまとめて $$ \Theta = (\mathbf{Q}, \phi_{(1)}, \ldots, \phi_{(m)}, \mu_{(1)}, \ldots, \mu_{(m)}, \Sigma_{(1)}, \ldots, \Sigma_{(m)}) $$ とする。ここから\(\mathbf{Q}\)を抜いた奴を\(\Theta_{-Q}\)とする。各状態の下での観察の確率を対角にいれた\(m \times m\)対角行列 $$ \mathbf{P}(\mathbf{Y_t}) = diag \left( \{p(\mathbf{Y}_t | \mathbf{Y}_{t-1}, \Theta_{-Q}, S_t= j)\}_{j=1}^{m} \right) $$ をつくっておく。
さて、尤度関数はどうなるかというと… 要素がすべて1の行ベクトルを\(\mathbf{1}\)として、$$ l(\Theta) = p(\mathbf{Y}_{2:n}|\Theta, \mathbf{Y_1}) = \delta \mathbf{P}(\mathbf{Y}_2) \mathbf{Q} \mathbf{P}(\mathbf{Y}_2) \mathbf{Q} \mathbf{P}(\mathbf{Y}_3) \cdots \mathbf{Q}(\mathbf{Y}_n) \mathbf{1}^\top $$ [あー、そうか！ マルコフ連鎖が時間等質で定常分布を持つから、尤度関数では\(S_t\)を消してめっちゃ簡単に書けるのか。なるほどー。言ってる意味がやっとわかったよ]
2.2 \(\Theta\)のサンプリング
事前分布を\(p(\Theta)\)として、事後分布は $$ p(\Theta | Y_{1:n}) \propto l(\Theta) p(\theta) = p(\Theta) \times \int p(Y_{2:n}|S_{2:n}, \Theta, \mathbf{Y}_1) p(S_{2:n}|\Theta) d S_{2:n} $$ と書ける。どうやって推定するか。ふたつ方法がある。
- \(S_{2:n}\)について周辺化する。アイデア自体は新しくない。本論文ではHMCでやる。
- Gibbsサンプリングで\(\Theta\)と\(S_{2:n}\)を同時にサンプリングする。いまよく使われている方法。
2.2.1 周辺尤度に基づくMCMC
以下のアルゴリズムを用いる。[擬似コードが書いてあるんだけど、日本語でメモする]
- 時点2について、状態 \(i\) ごとに以下を求める。$$ \alpha_{2,i} = \log \delta_i + \log p(\mathbf{Y}_t | \mathbf{Y}_{t-1}, \Theta_{-Q}, S_t = 1) $$ [ああなるほど、時点2だけは対数尤度の求め方が違うのか]
- 時点2について、状態 \(i\) ごとに以下を求める。$$ \alpha_{3,i} = \log \left( \sum_{j=1}^m \exp(\alpha_{2,i} + \log \mathbf{Q}_{j,i} + \log p(\mathbf{Y}_t | \mathbf{Y}_{t-1}, \Theta_{-Q}, S_t = i)) \right) $$
- 時点3, 4, …について繰り返す。
- 以下を求める。$$ l(\Theta) = \log p(\mathbf{Y}_{2:n} | \Theta, \mathbf{Y}_1) = \log \sum_{i=1}^m \exp(\alpha_{n,i}) $$ [あれれ？ これ、状態で条件づけた尤度の和を最後の時点だけについて求めて対数にしてない？ 変だなあ… これでいいの？ ]
なお、この方法では状態についての推論は行わないわけだが、あとでViterbiアルゴリズムを使って推定できる。[よくわからんけど、ま、パラメータの事後分布と\(\mathbf{Y}_t\)を使って時点2から順に推定していくとか、たぶんそんな感じだろうな]
2.2.2 Gibbsサンプリング
[まるごとスキップ。要は比較対象の話だからとばしていいや、と自分に言い訳し… 難しい話苦手なんです]
2.2.3 実装
周辺化アプローチはStanでやる。Viterbiアルゴリズムはgenerate quantitiesブロックで書ける。
GibbsサンプラーはJAGSでやる。他に、Rcppパッケージを使って自分で実装したのでそれも試してみる。
2.4 ラベルスイッチング
教師なしのマルコフ・スイッチングモデルでは、ラベル・スイッチングによって事後分布が多峰になるのがふつうである。解釈のためにはなんらかの方法で状態を識別する必要がある。ふつうはなにかしらのパラメータについて順序制約をつける。特にベイズ推定の場合、制約が事後分布のジオメトリを無視しているとラベルが一意にならなくなるので、制約を注意深くかける必要がある。いっぽう、MCMCのなかではなにもせず、あとからどうにかすりゃいいじゃんという説もある。本研究もこの説に従い、Rのlabel.switchingパッケージを使ってあとからなんとかする。もっともStanの場合、うまい制約をつけるとサンプリングが速くなるかもしれない。
3. 数値比較
[すいません… ここがこの論文の本題なんだけど、私、全然関心なくて… まるごとスキップ]
4. 天然ガスと原油価格の同時ダイナミクス
[数値例。というか、こうやってわざわざ示す例が2変量のVARなんだ… やっぱし、計算は結構大変なんだろうな…]
5. 考察
全体として、周辺化アプローチのほうが頑健な結果を示し効率も良いことがわかった。JAGSは簡単だがめっちゃ効率が悪く、自分で実装したのは効率が超よかったけど書くのが大変だった。
云々。
—————–
というわけで、命じられるままに論文に目を通したうえで、Stanコードに戻ったんだけど… ううう、超めんどくさい。やっぱし俺、社会人に向いてない。働くのに向いてない。猫とか石とかに生まれればよかった。