非農業部門雇用者数を状態空間モデルで予測する

概要

・状態空間モデルの練習も兼ねて非農業部門雇用者数(の変化)をモデル化し、1期先予測を行なう。
・使うモデルはローカルレベルモデル

ローカルレベルモデルとは

学んだばかりのものをドヤ顔で紹介するわけですが、ローカルレベルモデルとは以下のように定式化されます。
$$
y_t = \mu_t + v_t \\
\mu_t = \mu_{t-1} + w_t
$$
ここで\( v_t, w_t \)は正規分布に従い
$$
v_t \sim N \left( 0, V \right) \\
w_t \sim N \left( 0, W \right)
$$
と表されます。解釈としては、\( \mu_t \)が\( t \)時点での観測不可能な状態、これがランダムウォークに従い、また実際に観測される値(今回は非農業部門雇用者数(NFP)の変化が該当する)は、その状態から観測誤差( \( v_t \) )を伴って現れるということですね。

このローカルレベルモデルは、状態空間モデルの中の動的線形モデル(DLM: Dynamic Linear Models)というクラスに分類されます。状態がダイナミック(定数でない)ということですね。Rではdlmという素晴らしいライブラリがありますので、こちらを使用させていただきます。

データ取得

例の如く、セントルイス連銀先生からデータを取得します。

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

diffを取り、plotすると以下のようになります(xtsオブジェクトに変換しています)。
dataが取得したデータを格納している変数です。

140301_1

NFP自体は1939年1月からの系列ですので、diffを取ると、1939年2月からの系列となります。
リーマン・ショック期や、第二次世界大戦中は大きく減少しています。

パラメータ推定

上述したDLMには正規分布の分散という未知パラメータ\( V, W \)がありますので、まずはこれを推定します。
参考文献(後述)のP.150あたりを参考にすると

dlmMLE関数で、初期値(0.1)、パラメータの下限(未知パラメータが分散なので負の値は取らない)を与えた上で最適化してくれます。

『DLMの尤度関数には、多数の局所最大が存在する可能性がある』(P.147)
『特に、要素convergenceは、常に確認する必要がある。この値が非零であれば、最小値への収束が達成されていない知らせとなる。』(P.148)

とありますので、先ほどの結果は、最小値は達成されているが、局所最大の可能性があるということですね。実際初期値を変えて推定し直すと、異なったパラメータが推定されます。が、そこまで大きく違いもないようなので、今回はこのまま使います。

カルマンフィルタ

先ほどの推定されたパラメータをもとにカルマンフィルタを動作させます。カルマンフィルタにより各時点での状態\( \mu_t \)を推定できます。

実際に推定結果をプロットすると以下

140301_2

黒線が実際の観測値、赤線がカルマンフィルタにより推定された状態\( \mu_t \)です。観測誤差分を除いた比較的滑らかな線が引かれています。これが観測不可能な状態の推定結果ということですね。

各時点の状態を推定することにより、1期先の観測値も予測することが出来ます。

1期先予測

2014年1月時点での値が152.277なので、2月のNFPの予測も152.277となりますね(正規分布の期待値は0なので)。
うん、何か当たらない気がする笑

どの程度予測力があるのか

何やら小難しいモデルを使っていますが、そもそもこのローカルレベルモデルは予測力あるんですかというのは当然の疑問。業界でよく使ってる3ヶ月平均とかでも良いんじゃねーのと思いますわな。
ということで予測力の比較を行います。

そこで使われるのがTheilのUというもので詳しくは参考文献P.100を参照。TheilのU自体は以下のようにして計算されます。
$$
U = \sqrt{ \frac{ \sum_{t=2}^n \left( y_t – f_t \right)^2 }{ \sum_{t=2}^n \left( y_t – y_{t-1} \right)^2 } }
$$

変化のないモデル(今期の値を来期の値の予測値とするモデル)でのMSEと、比較したいモデルのMSEの比をとっています。Uの値が小さいほど、当該モデルのMSEが小さいことを意味しますので、予測力も高いということですね。

まずはDLMの場合のTheilのU

カルマンフィルタがデータに適応する猶予を与えるため、最初の10サンプルは除外した上で計算しています。で、いくつかというと

ということで、変化のないモデルに対して、67%のMSEで予測を行なっていることが分かります。予測誤差が33%減少しているイメージですかね。

移動平均の計算自体は手前味噌ですが過去作った関数を流用します。

[R]テクニカル指標の計算
http://nekopuni.holy.jp/?p=974

3ヶ月移動平均と、6ヶ月移動平均も計算しておきます。

DLMのときと同じサンプルに対する予測力を比較するため、こちらも最初の10サンプルを除外しています。
で、いくつになるかというと

ということで、3ヶ月移動平均を用いた予測、6ヶ月移動平均を用いた予測とも、DLMよりは劣る予測力ということが分かります。とはいっても3ヶ月移動平均なんかは3%ほどの予測力の違いですので、コストパフォーマンスを考えたら移動平均でも良いんじゃないかという気がしないでもない結果に・・汗

観測値、DLM、移動平均のプロット

最後に、どんな感じでそれぞれのモデルが計算しているかを見ておきます。

140301_3

黒線が観測値、赤線がDLM、緑線が3ヶ月移動平均です。やはり2013年12月の74,000が効いていますね笑

参考文献

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