多変量動的線形モデルにおけるパラメータ最尤推定

概要

・Rのdlmパッケージを用いて多変量(今回は2変数)の動的線形モデルのパラメータ最尤推定を行う
・対象となるデータはNFP及びADP
・1期予測もしてみる(3月末までに入手可能な情報をもとに計算しています)

データの取得

恒例のセントルイス連銀先生

NFPのデータはこちら
All Employees: Total nonfarm – FRED – St. Louis Fed
https://research.stlouisfed.org/fred2/series/PAYEMS/

ADPのデータはこちら
Total Nonfarm Private Payroll Employment – FRED – St. Louis Fed
https://research.stlouisfed.org/fred2/series/nppttl/

diffをとって、両系列が入手可能な2001年5月から一応データの確認。
黒線がNFP, 赤線がADPとなります。

140405_1

モデルの概要

今回は最終的に和分ランダムウォークを2変量へと拡張します。
まずは、2変量のDLMを構築するわけですが、といっても概形は1変量とあまり変わりません。

$$
\left(
\begin{array}{c}
y_{1,t} \\
y_{2,t} \\
\end{array}
\right)
=
\left(
\begin{array}{cccc}
1 & 0 & 0 & 0 \\
0 & 1 & 0 & 0 \\
\end{array}
\right)
\left(
\begin{array}{c}
\mu_{1, t} \\
\mu_{2, t} \\
\beta_{1, t} \\
\beta_{2, t} \\
\end{array}
\right)
+
\left(
\begin{array}{c}
v_{1, t} \\
v_{2, t} \\
\end{array}
\right)
$$

$$
\left(
\begin{array}{c}
\mu_{1, t} \\
\mu_{2, t} \\
\beta_{1, t} \\
\beta_{2, t} \\
\end{array}
\right)
=
\left(
\begin{array}{cccc}
1 & 0 & 1 & 0\\
0 & 1 & 0 & 1\\
0 & 0 & 1 & 0\\
0 & 0 & 0 & 1\\
\end{array}
\right)
\left(
\begin{array}{c}
\mu_{1, t-1} \\
\mu_{2, t-1} \\
\beta_{1, t-1} \\
\beta_{2, t-1} \\
\end{array}
\right)
+
\left(
\begin{array}{c}
w_{1, t} \\
w_{2, t} \\
w_{3, t} \\
w_{4, t} \\
\end{array}
\right)
$$

ここで\( v_t, w_t \)は正規分布に従い
$$
v_t \sim \left( 0, V \right) \\
w_t \sim \left( 0, W \right)
$$
というような感じです。ここで\( V, W \)は分散共分散行列となります。

もう少し書き下すと
$$
V=
\left(
\begin{array}{cc}
\sigma_{1, y}^2 & \rho_{1,2;y} \sigma_{1, y} \sigma_{2, y} \\
\rho_{1,2;y} \sigma_{1, y} \sigma_{2, y} & \sigma_{2, y}^2 \\
\end{array}
\right) \\
W=
\left(
\begin{array}{cc}
W_\mu & 0 \\
0 & W_\beta \\
\end{array}
\right)
$$
となります。
また、
$$
W_\mu=
\left(
\begin{array}{cc}
\sigma_{1, \mu}^2 & \rho_{1,2;\mu} \sigma_{1, \mu} \sigma_{2, \mu} \\
\rho_{1,2;\mu} \sigma_{1, \mu} \sigma_{2, \mu} & \sigma_{2, \mu}^2 \\
\end{array}
\right) \\
W_\beta=
\left(
\begin{array}{cc}
\sigma_{1, \beta}^2 & \rho_{1,2;\beta} \sigma_{1, \beta} \sigma_{2, \beta} \\
\rho_{1,2;\beta} \sigma_{1, \beta} \sigma_{2, \beta} & \sigma_{2, \beta}^2 \\
\end{array}
\right)
$$
となります。

で\( W_\mu, W_\beta \)の中身なんですが、全部の成分がある場合は線形成長モデル、\( W_\mu \)がゼロ行列の場合には和分ランダムウォークとなります。
線形成長モデルの場合には、観測行列のパラメータも含めて9個の未知パラメータがあるわけですが、パラメータは少ないほうが良いということで、今回は和分ランダムウォークということで未知パラメータは6個となります。
また、\( \mu \)はローカルレベルと呼ばれましたが、\( \beta \)はローカル成長率と呼ばれるようです。

最尤推定の準備

まず最尤推定用に、DLMモデルを構築する関数を定義します。

引数paramsは長さ6のベクトルで、1番めから順番に、
1. \( y_1 \)の観測誤差の分散
2. \( y_2 \)の観測誤差の分散
3. \( y_1, y_2 \)の観測誤差の相関係数
4. \( \theta_1 \)の遷移誤差の分散
5. \( \theta_2 \)の遷移誤差の分散
6. \( \theta_1, \theta_2 \)の遷移誤差の相関係数
を表します。ただし、分散は非負制約、相関係数は[-1, 1]の制約があります。分散に関してはエクスポネンシャル、相関係数に関してはtanhをかませることでその制約を満たすようにしています。
分散共分散行列のうち、非対角要素に関しては、共分散を直接推定するのではなく、相関係数を推定することで間接的に共分散を推定しています。

このへんの推定方法に関しては、参考文献2のChapter4.R、参考文献3のdlmMLEの辺りなんかも参考になります。

ちなみにtanhの形としては以下。

140105_2

『始めに(一変量)線形成長モデルを作成し、次にクロネッカー積を用いて行列\( F \)と\( G \)を再定義している。このアプローチであれば、\( F \)と\( G \)の個別の要素を手入力するより打ち間違いを減らしやすい』(P.132より)

2変量のローカルレベルモデルであればそこまで打ち間違いは発生しないと思いますが、3変数以上の場合や、状態のオーダーが2以上になった場合には分散共分散行列がカオスになりますのでクロネッカー積を使ったほうが間違えにくいということですね。

最尤推定

前節で作った関数をもとに、実際にパラメータ推定を行ないます。

第一引数に元データ、第二引数にパラメータの初期値、第三引数にさきほど定義したbuildModel関数を用います。
2行目で収束しているかどうか、3行目で実際に推定されたパラメータを表示させています。結果を見ますと、

ということでconvergeceは0となってりますのでとりあえずの収束はしているようです。
ただし、ここで表示しているのはパラメータの仮の姿ということで、実際に欲しい分散はエクスポネンシャルなどをかませたあとの数字となります。

実際に推定されたパラメータを用いてモデルを構築します。
その後、実際の観測誤差、遷移誤差の分散共分散行列を表示させると以下。

カルマンフィルタ + 1期先予測

上記で推定されたパラメータをもとにカルマンフィルタを動作させます。

model.filtにカルマンフィルタ結果、model.filt.stateに推定されたNFP, ADPの状態を保存しておきます。
dropFirstしている理由は、$mに事前分布も格納されているためです。

\( F \)の形にありますように、1期先のNFP, ADPの期待値は、今期の状態\( \mu_1, \mu_2 \)と等しくなります。

で、観測値と状態\( \mu_1, \mu_2, \beta_1, \beta_2 \)をそれぞれプロットしてみますと以下のようになります。
それぞれ、黒線が実測値、赤線が\( \mu \)、濃いグレーの面グラフが\( \beta \)となります。
(面グラフを重ねてかくためだけにggplot2を使ってます…)

140105_3

140405_4

で、直近の状態を出力しますと

となりまして、3月のNFP期待値が139k, ADP期待値が115kとなります。
直近の\( \beta \)なんかを見ても、ローカル成長率はマイナスの値となっておりますし、観測値の成長率としては以前マイナス方向を示唆しているというわけですな(― ―)

ま、3月の実績値はADPが191kで、NFPが192kなんですけどね!!!

再送-UPDATE 1-米ADP民間雇用者数、3月は19.1万人増 寒波の影響薄らぐ兆候示す | マネーニュース | 最新経済ニュース | Reuters
http://jp.reuters.com/article/marketsNews/idJPL4N0MU33Q20140402

UPDATE 2-米雇用統計、3月の雇用者数19.2万人増 寒波からの回復鮮明 | マネーニュース | 最新経済ニュース | Reuters
http://jp.reuters.com/article/marketsNews/idJPL4N0MW31S20140404

参考文献

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

Website for “Dynamic Linear Models with R”
http://definetti.uark.edu/~gpetris/dlm/#chapter4

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

R Financial Time Series Plotting
http://timelyportfolio.github.io/rCharts_time_series/history.html

R: Convenience Functions for Plotting zoo Objects with ggplot2
http://hosho.ees.hokudai.ac.jp/~kubo/Rdoc/library/zoo/html/ggplot2.zoo.html