DLMを用いた動的回帰

概要

・回帰分析における回帰係数が動的に変動する動的回帰というものをやってみる
・というかほとんどネタ消化モード
・データはいつものNFPとADP。ADPのほうが先に公表されるので、ADPを使って同月のNFPが予測できたら嬉しい

データの準備

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となります。
ココらへんはいつもの作業

140420_1

モデル式

動的回帰ですので回帰係数が動的に変動します。
イメージとしては以下の様な感じですかねえ

$$
{\rm NFP}_t = \beta_{1,t} + \beta_{2,t} \cdot {\rm ADP}_t + v_t \\
v_t \sim N \left( 0, \sigma_v^2 \right)
$$

また、回帰係数はランダムウォークに従うという仮定でして、

$$
\beta_{1, t} = \beta_{1, t-1} + w_{1, t} \\
\beta_{2, t} = \beta_{2, t-1} + w_{2, t} \\
w_{i, t} \sim N \left(0, \sigma_{w,i}^2 \right)
$$

まあこんな感じと。
一般的には回帰係数がランダムウォークでない形になるかと思うのですが、Rのパッケージdlmがランダムウォークにしか対応していないようなので、今回はこの形でご勘弁ということで。

未知パラメータ推定

今回は回帰式における\( v \)の標準偏差\( \sigma_v \)、回帰係数ランダムウォークの標準偏差 \( \sigma_{w, 1}, \sigma_{w, 2} \)が未知となりますのでこちらを最尤推定

動的回帰において、dlmモデルを構築する関数はdlmModRegとなりまして、第一引数には説明変数のマトリックスを代入します。
そして、各未知パラメータは非負制約がありますのでexpをかませています。
ドキュメントによると重回帰もいけるようなので誰かお願いします。

あとはdlmMLEで最尤推定します。
dlmMLEの第一引数には非説明変数であるNFPのデータを代入します。

で、推定された数字をもとにdlmを構築し、パラメータを確認すると以上のようになります。
分散が7177.945ですんで、標準偏差が84.7と。元の数字が雇用者数の伸びで100とか200とかの値であることを考えるとだいぶブレますな…(― ―)(― ―)(― ―)

カルマンスムーザで回帰係数が過去どのように推移していたかを見る

平滑化ですな、

140420_2

プロットすると以上のようになりまして、上が\(\beta_{1, t} \)、下が\( \beta_{2, t} \)の推移となります。
まず\( \beta_{2, t} \)から先に見ると、これ推移しているように見えますが縦軸みるとスケールがかなり細かいことになっていましてこれはほとんど0.957から動いていませんな。ということで、ADPの96%掛けがNFPに大体なりますと。

で、次に切片項の\(\beta_{1, t} \)を見ると、2000年前半は値が10を超えていまして、NFPのほうがADP*0.96よりも上方にでる傾向があったということですかねえ。まあ回帰式の標準偏差が84.7だったんでアレですが。それと直近は\(\beta_{1, t} \)の値も0に近くなっていて、NFPとADP*0.96の差が無くなってきていることを示唆しているということなんですかね。まあ回帰式の標準偏差(略

1期先回帰

で、お馴染みの1期先回帰でありますが、まだADPが発表されていないので予測はできませんが回帰係数がランダムウォークなので回帰係数の期待値自体は求めることができまして、

まずデータの方が

となっておりまして、回帰係数のほうが

となりますので、
$$
{\rm NFP}_{Apr-2014} = 1.22 + 0.96 \cdot {\rm ADP}_{Apr-2014}
$$

という感じですかな。まあ回帰式の(略

参考文献

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

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

概要

・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

異なる祝日に従う複数の時系列データをまとめる+ggplot2を用いたプロット

概要

・複数の時系列データをまとめる+プロットするためのほぼ自分用メモ
・例えばNYSEに上場している株式と、東証に上場している株式では米国祝日と日本祝日の日にそれぞれ欠損値がある。
・ggplot2の練習も兼ねて複数の時系列をまとめてプロットする

データの準備

いちいち株価を取得するのは少し面倒なので自前でテストデータを準備

最後の2行はデフォルトのplotを使っていまして、実際のプロット結果としては以下の様な感じ。

140402_1

このままだと日付と関連付けられていないため、一旦日付を紐付けxtsオブジェクトへと変換

どんなデータとなっているかというと

このような感じで、stockAは4月6日が祝日だっためデータが欠損していますが、stockBは平日だったためデータが存在。
4月12日は逆パターンといった感じ。

複数の時系列データをまとめる

といっても1行で済むのですが、

merge関数で複数のxtsオブジェクトをまとめます。
この関数はデフォルトでjoin=”outer”という設定がされておりまして、SQLをイメージするとわかりやすいかもしれないですが、外部結合ってやつですかね。
mergeだけ行った場合は以下の様な感じに。
片方しか存在しない日付に関しては、相方にNAが格納されています。

さらにその外側でna.locf関数をかませることで、NAの部分を直近の値で上書きするようにしてくれます。
エクセルでいうVlookupの第四引数をTRUEにした場合でしょうか

こうすると何が嬉しいかというと、リターン系列を簡単に計算することができる。

ret1はadditive return(という言葉があるかは知らんが汗)の計算でして、要は今日の株価から昨日の株価を引いたものが今日のリターンとするもの。
ret2はいわゆる普通のリターンの計算でして、今日の株価/昨日の株価-1という計算ですな。
arithmeticという引数をTrue(デフォルト)にすると、単純なdiffをとってくれる一方で、Falseにすると、比率で計算してくれるというわけですな。

祝日に関してはリターン0が格納されているので問題無いですね。

プロットしてみる

で、このデータをまとめてプロットしてみます。
普通(?)のやり方としては

140402_2

一方で、ggplot2を使う場合は

140402_3

facets=NULLを指定することで1つの画面に重ね書きしてくれます。

参考

Package ‘xts’
http://cran.r-project.org/web/packages/xts/xts.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