MCMCサンプリング(非農業部門雇用者数編)

概要

・前回(http://nekopuni.holy.jp/?p=983)、非農業部門雇用者数(NFP)に対してローカルレベルモデルを適用し、そのパラメータは最尤推定により推定した。
・今回は、巷で話題沸騰のMCMCを用いて、パラメータの不確実性を考慮した推定を行なう。

最尤推定の問題点

最尤法でパラメータ \( \psi \)を推定できるからいいジャマイカというのが尤もな感覚なのですが、
『フィルタリング漸化式や平滑化漸化式にMLE\( \hat{ \psi } \)を使うのは一般的な慣習ではあるが、\( \psi \)の不確実性を適切に考慮しづらいという欠点がある。』(P.150)

最尤推定のアイディアとしては、現在取得できるデータに最もフィットするパラメータを選択するというものですから、真のパラメータとは外れたところで推定されている可能性もあり、そのパラメータのバラ付き具合も含めて推定してあげたほうが紳士的なんじゃないですかということですかね(ベイズ推定)

データ取得

前回同様(というかいつもですな)、NFPのデータはセントルイス連銀先生から。

All Employees: Total nonfarm (PAYEMS) – FRED – St. Louis Fed
http://research.stlouisfed.org/fred2/series/PAYEMS/

diffを取ったものを変数nfpに格納しておきます。

パラメータ(分散)の分布を仮定

ローカルレベルモデルも再掲

$$
y_t = \mu_t + v_t \\
\mu_t = \mu_{t-1} + w_t
$$
ここで\( v_t, w_t \)は正規分布に従い
$$
v_t \sim \left( 0, V \right) \\
w_t \sim \left( 0, W \right)
$$

今回推定する未知パラメータは\( V, W\)の2つです。
ここでパラメータ(分散)がガンマ分布に従うという仮定を置きます。
正確には、分散の逆数(以下、精度と言う)がガンマ分布に従うというものでして、

$$
V_t^{-1} = \phi_y \sim G \left( \alpha_y, \beta_y \right) \\
W_t^{-1} = \phi_\mu \sim G \left( \alpha_\mu, \beta_\mu \right)
$$

とします。

実際にMCMC推定を行なう

dlmパッケージを用いることで簡単にMCMCサンプリングを行なうことが出来ます。
中身でやっていることは、ギブスサンプリングでして、

\( \pi \left( \theta_{0:T} | y_{1:T}, \psi \right) \)から\( \theta_{0:T} \)
をサンプリングし、ここでサンプリングされた\( \theta_{0:T} \)を用いて
\( \pi \left( \psi | y_{1:T}, \theta_{0:T} \right) \)から \( \psi \)
をサンプリング

というのを何回も行なうというもの(らしい)です。片方を固定して、もう片方をサンプリングってイメージですかね。
で、どうやってコードで書くかというと

となります。
引数の意味としては
・a.y: \( \phi_y \)の事前平均
・b.y: \( \phi_y \)の事前分散
・a.theta: \( \phi_\mu \)の事前平均
・b.theta: \( \phi_\mu \)の事前分散
・n.sample: サンプリングする数
・thin: 何回おきにサンプルを保存するか
・save.states: 状態\( \mu \)を保存するか
となるようです。
事前平均と事前分布でなく、ガンマ分布の事前パラメータを指定することも可能で、その場合はshape.y, rate.y, shape.theta, rate.thetaとなるようです。ガンマ分布の事前平均・分散のほうが直感的に入れやすいということなんですかね。
thin = 1とすることで1個おきのサンプルのみ最終的に保存されます。すなわち、実際には44000回サンプリングは行なわれているということですね。

で、この関数は実行すると無茶苦茶時間がかかります。(Mac book proで1時間半ほど…)

また、このdlmGibbsDIG関数ですが、公式のドキュメントにもあるように、
『The d-inverse-gamma model is a consistent univariate DLM with unknown observation variance, diagonal system variance with unknown diagonal entries.』
となっておりますので、多変量DLMや、\( W \)が対角行列以外を想定した場合には対応していないということですかね。

結果

サンプルの生データをプロットしてみますと

140321_1

140321_2

といった具合でサンプルされたようです。
まあ一口にパラメータと言ってもこれだけバラつきますよというお話ですね。

ここで、最初の2000サンプルをバーンイン期間として除外し、パラメータの移動エルゴード平均を計算すると以下。
バーンイン期間を除いた上で、500個目から移動エルゴード平均を計算しています。次第にある程度の水準で収束してそうな気配がしておりますな。
(このバーンイン期間をどの程度取ればよいのかが謎ですが。)

140321_3

140321_4

また、バーンイン期間を除く、MCMCサンプルの平均をとると

となりまして、\( V \)の推定値が22042.9、\( W \)の推定値が4117.3となりました。
カッコ内の数字はモンテカルロ標準偏差の推定値。

最尤推定との比較

ちなみに最尤推定による結果は

となりました(前回記事投稿時から、NFPのデータが更新されているので推定しなおしております)。
まあこんなもんかなあという感じ

この手法による欠点

・dlmGibbsDIG関数を一般化できない
前述したように、多変量DLMには使えないですし、分散行列も対角行列でなければ対応してない。
(単変量DLMで分散行列に非対角要素が存在するケースを考えなければならない場合はそんなに無いかもしれないですが…)

・オンライン推定できない
新たな観測値が観測された場合、再度MCMCを行う必要がある。
非常に計算コストがかかるので、カルマンフィルタのように逐次的に計算できる方法のほうが良い
この問題点に関しては逐次モンテカルロ法により解決できる(らしい)。
ということで当面の課題としては、多変量の推定ができるようになるということとオンライン推定ができるようになることということで(― ―)(― ―)(― ―)

参考文献

Amazon.co.jp: Rによるベイジアン動的線形モデル (統計ライブラリー): G.ペトリス, S.ペトローネ, P.カンパニョーリ, 和合 肇, 萩原 淳一郎: 本
http://www.amazon.co.jp/dp/4254127960

Giovanni Petris: dlm: an R package for Bayesian analysis of Dynamic Linear Models
http://cran.r-project.org/web/packages/dlm/vignettes/dlm.pdf

Package ‘dlm’
http://cran.r-project.org/web/packages/dlm/dlm.pdf