Rでベイジアン動的線形モデルを学ぶ(3)：ローカル線形トレンドモデル - 渋谷駅前で働くデータサイエンティストのブログ

相変わらずグダグダな上に挙句の果てに既にRでやっちゃった例をまとめたPDF bookまであると判明してモチベーションだだ下がりなんですが、備忘録も兼ねてめげずに続けます。もちろんテキストは相変わらずこちらの2冊。

状態空間時系列分析入門

作者: J.J.F.コマンダー,S.J.クープマン,Jacques J.F. Commandeur,Sime Jan Koopman,和合肇
出版社/メーカー: シーエーピー出版
発売日: 2008/09
メディア: 単行本
購入: 2人クリック: 4回
この商品を含むブログを見る

Rによるベイジアン動的線形モデル (統計ライブラリー)

作者: G.ペトリス,S.ペトローネ,P.カンパニョーリ,和合肇,萩原淳一郎
出版社/メーカー: 朝倉書店
発売日: 2013/05/08
メディア: 単行本
この商品を含むブログを見る

後は冒頭に挙げたPDF bookを参照しながら、パラメータ設定のところあたりの分からないところを補足しながらだらだらやっていこうかと思います。

ローカル線形トレンドモデル

モデル式自体は実は既出で、以下の通りです。

$Y_t = \mu_t + \nu_t, \nu_t \sim \cal{N} (0,V)$ （観測方程式）

$\mu_{t} = \mu_{t-1} + \beta_{t-1} + \omega_{t,1}, \omega_{t,1} \sim \cal{N} (0,\sigma_\mu^2)$ （状態方程式）

$\beta_{t} = \beta_{t-1} + \omega_{t,2}, \omega_{t,2} \sim \cal{N} (0,\sigma_\beta^2)$ （トレンド方程式）

これを行列表現に直したものが下記。

$Y_t = F_t \theta_t + \nu_t, \nu_t \sim \cal{N} (0,V_t)$ （行列化した観測方程式）

$\theta_{t+1} = G_t \theta_t + \omega_t, \omega_t \sim \cal{N} (0,W_t)$ （行列化した状態方程式）

ただし

$\theta_t = \begin{bmatrix} \mu_t \\ \beta_t \end{bmatrix}, G = \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix}, W = \begin{bmatrix} \sigma_\mu^2 & 0 \\ 0 & \sigma_\beta^2 \end{bmatrix}, F = \begin{bmatrix} 1 & 0 \end{bmatrix}$

これだけ見てもよく分からないのでローカルレベル・モデルとの違いを端的に書くと、トレンド方程式＆状態方程式のところに線形トレンド項 $\beta_t$ が追加されている点が違います。この項によって中短期のトレンドを表現しようというわけです。

ちなみに誤差項を見れば分かりますが、このモデルでは観測方程式に1つ、状態方程式に1つ、トレンド方程式に1つ、誤差項のパラメータが存在するので合計3つの未知パラメータがあります。なので、最尤法でパラメータ推定する際は未知パラメータを最大3つまで推定することになります。

トレンドの意味合いについては、既に以前沖本本シリーズ記事で単位根過程＆トレンド過程を取り上げた際に色々コメントしたかと思いますので、そちらを読んでくださいということで。

ここでカルマンフィルタとカルマンスムーザ（平滑化）の話をば

そう言えば、前回前々回と何気なくカルマンフィルタとdlmFilter関数の話をしていて何の説明もつけていなかったので、自分の備忘録として書いておきます。Petris本2.7節には

状態ベクトルを推定するためには、条件付き密度 $\pi(\phi_s | y_{1:t})$ を計算する。ここで、フィルタリング(s = tの場合)、状態予測(s > t)、そして平滑化(s < t)の問題を区別する。

と書いてあるので、ここで一旦それぞれの説明を引用してそれぞれの概念がどういうものかを確認しておきます。

フィルタリング

フィルタリングの問題では、データは時間順に到着すると想定されている。これは多くの応用問題にあてはまる。例えば、移動中の目標を追跡する問題でもそうだし、金利の期間構造を毎日推定しなければならず、翌日の市場で新しいデータが観測された際に現在の推定値を更新するような金融分野での応用でもそうである。

これらの場合においては、時点t（[いま]）までの観測値に基づいて状態ベクトルの現在の値を推定し、時点t + 1で新しいデータが利用可能になると推定値と予測値を更新するような手続きが欲しい。

フィルタリング問題を解決するためには、条件付き密度 $\pi(\theta_t | y_{1:t})$ を計算する。DLMでは、新しいデータが利用可能になった際に、カルマン・フィルタが状態ベクトルの現在の推定を更新する式を提供し、フィルタリング密度は $\pi(\theta_t | y_{1:t})$ から $\pi(\theta_{t+1} | y_{1:t+1})$ に移る。

基本的には「現在を含むこれまでの」データから「現在の」状態を推定する手続きなので、基本的には前向きのパスに基づいて現在の（外からは見えない状態の）真の値が欲しい時に使うもの、ということのようです。Petris本の2.7.1 - 2.7.2節に理論的な説明が載っています。

平滑化

一方、平滑化問題もしくは回顧的分析は、データ $y_1, \cdots, y_t$ が与えられた下で時点 $1, \cdots, t$ における状態系列を推定することからなる。多くの応用で望まれるのは、一定期間時系列を観測してから、過去に遡ってその観測値の基礎になるシステムの振る舞いを検討することである。

例えば、経済学の研究では、ある国における一定年数分の消費や国内総生産の時系列があり、過去に遡ってシステムの社会経済的な振る舞いを理解することに関心が持たれるかもしれない。

平滑化問題を解決するためには、 $y_{1:t}$ が与えられた下で $\theta_{1:t}$ の条件付き分布を計算する。フィルタリングの場合と同様に、平滑化は逐次的なアルゴリズムとして実装できる。

フィルタリングとは異なり、後ろ向きのパスに基づいて過去の全データ履歴を遡りそれらを全て勘案することでより正確な全期間の真の値を得ようとするのが平滑化、ということでいいんでしょうか。実際、その後の記述を見ると「全ての状態の履歴を後ろ向きに推定する」とあります。Petris本2.7.4節に理論的説明が載っています。

状態予測

実際問題として、時系列分析では予測が主な課題となることも多い。ここでは、状態の推定は、将来の観測値を予測するための単なる一ステップとなる。

すなわち、データ $y_{1:t}$ に基づいて次の観測値 $Y_{t+1}$ を予測する一期先予測では、状態ベクトルの次の値 $\theta_{t+1}$ を最初に推定し、この推定値に基づいて $Y_{t+1}$ の予測値を計算する。

状態の一期先予測密度は $\pi(\theta_{t+1} | y_{1:t})$ であり、これは $\theta_t$ のフィルタリング密度に基づいている。状態の一期先予測密度から、観測値の一期先予測密度 $\pi(y_{t+1} | y_{1:t})$ を得る。

これまでと違うのは、状態ではなく出力（観測値）を求めるという操作を伴うところですね。基本的にはPetris本2.7.2節で述べられているように、カルマンフィルタの副産物として得られるものなので見かけ上やってることはフィルタリングと変わらないとも言えます。

ただし、Commandeur本8.4節p.90でも述べられているように、状態と観測値とを分離していることによって予測値を得る仕組みは若干直観的なイメージとは異なる気がするので、気になる人は確認しておいた方が良いかも。

以上の3つの推定についてはCommandeur本では8.4節でコンパクトにまとめて説明されているので、読み比べてみるのも良いかもしれません。

Commandeur本のサンプルデータでやってみる

では、実際にやってみましょう。今回はCommandeur本の第3章からです。ただし使うデータは相変わらずUK driverとFinland fatalityとかいう2つの物騒なデータです（笑）。

ローカル線形トレンドモデルは一応「二次」のモデルということになっているので、dlmModPoly関数でモデルを立てる際にはorder = 2（ただしデフォルト値でもある）と与える必要があります。そこを踏まえて、まずはUK driverの方に当てはめてみましょう。

ところで、今回はCommandeur本のあんちょこがあるのでこれにやり方も合わせることにします。

UK driverのデータの場合

これはCommandeur本3.2節の「確率的レベル＋確率的傾き」のケースであることに留意。

# 既にUK driverのlogとしてd1、Finland fatalityのlogとしてd2があるものとします
> level0<-d1[1]
> slope0<-mean(diff(d1))
# ベクトルm0を固定するために値を決めておく
# なおm0のレベルはd1生データの最初の値、トレンドは全体の傾きの平均に固定
> build1<-function(parm){dlmModPoly(order=2,dV=exp(parm[1]),dW=exp(parm[2:3]),m0=c(level0,slope0),C0=2*diag(2))}
# 3つの未知パラメータは全て最尤法で推定する、ただしC0（m0の誤差分散の二乗）は2に固定
> fit1<-dlmMLE(d1,rep(0,3),build1)
> fit1$par
[1]  -6.157164  -4.412262 -24.910858
> model1<-build1(fit1$par)
> V(model1)
            [,1]
[1,] 0.002118253
# Commandeur本p.24の値と一致

> W(model1)
           [,1]        [,2]
[1,] 0.01212771 0.00000e+00
[2,] 0.00000000 1.51828e-11
# Commandeur本p.24の値と2つとも一致

> filt1<-dlmFilter(d1,model1)
# カルマンフィルタでフィルタリング
> smooth1<-dlmSmooth(filt1)
# 平滑化
> res1<-c(residuals(filt1,sd=F))
# 残差も出しておく
> matplot(cbind(d1,dropFirst(smooth1$s[,1])),type='l',lwd=c(2,3),lty=c(1,2),col=c('red','blue'),ylab='')
# 生データと平滑化後系列をプロットする
> plot(res1,type='l',lty=3,lwd=2,ylab=')
# 残差をプロットする
> plot(dropFirst(smooth1$s[,2]),type='l',lty=3,lwd=3,ylab='')
# トレンドをプロットする

f:id:TJO:20141015174939p:plain

f:id:TJO:20141015174949p:plain

f:id:TJO:20141015174956p:plain

Commandeur本図3.1, 3.2, 3.3とほぼ同じプロットになっていることを確認してください。これは状態のばらつきはゼロと仮定し、その代わり観測とトレンドにはばらつきがあると仮定するモデルですね。得られる $\mathbf{\theta}_t$ は2次元なので、必然的にカルマンフィルタ＆平滑化の返り値も2次元になります。

なお、残差のチェックをした結果は以下の通り。

> shapiro.test(res1)

	Shapiro-Wilk normality test

data:  res1
W = 0.9666, p-value = 0.0001563
# ダメ

> Box.test(res1,lag=15,type="Ljung")

	Box-Ljung test

data:  res1
X-squared = 101.8532, df = 15, p-value = 5.773e-15
# ダメ

どちらもパスせず、ということでUK driverのデータに対してはこれでもまだ当てはまりは良くないという結論になりました。では、Commandeur本3.3節の「確率的レベル＋確定的傾き」の場合ではどうでしょうか？

> build1m<-function(parm){dlmModPoly(order=2,dV=exp(parm[1]),dW=c(exp(parm[2]),0))}
# dWの1番目を推定（確率的レベル）、2番目をゼロで固定（確定的傾き）
> fit1m<-dlmMLE(d1,rep(0,2),build1m)
> model1m<-build1m(fit1m$par)
> V(model1m)
            [,1]
[1,] 0.002118081
> W(model1m)
           [,1] [,2]
[1,] 0.01212834    0
[2,] 0.00000000    0

プロットは省略しますが、Commandeur本p.27記載のパラメータと一致することが分かります。

Finland fatalityのデータの場合

今度は3.4節のローカル線形トレンドモデル、即ち「確定的レベル＋確率的傾き」のケースです。

> build2<-function(parm){dlmModPoly(order=2,dV=exp(parm[1]),dW=c(0,exp(parm[2])))}
# dWの1番目をゼロに固定（確定的レベル）、2番目を最尤法で推定（確率的傾き）
> fit2<-dlmMLE(d2,rep(0,2),build2)
> model2<-build2(fit2$par)
> V(model2)
            [,1]
[1,] 0.003200851
# Commandeur本p.29の値と一致
> W(model2)
     [,1]        [,2]
[1,]    0 0.000000000
[2,]    0 0.001533121
# Commandeur本p.29の値とほぼ一致
> filt2<-dlmFilter(d2,model2)
# カルマンフィルタ
> smooth2<-dlmSmooth(filt2)
# 平滑化
> par(mfrow=c(2,1))
> matplot(cbind(d2,dropFirst(smooth2$s[,1])),type='l',lwd=c(2,3),lty=c(1,2),col=c('red','blue'),ylab='')
> plot(dropFirst(smooth2$s[,2]),type='l',lty=3,lwd=3,ylab='')
> segments(0,0,35,0,lty=3)
# 図3.5のようなものをプロット
> par(mfrow=c(1,1))
> plot(res2,type='l',lty=3,lwd=2,ylab='')
# 図3.6のようなものをプロット

f:id:TJO:20141015190654p:plain

f:id:TJO:20141015190706p:plain

図3.5, 3.6を再現した結果です。でも残差だけ合ってない気がする。。。ちなみに残差のチェック結果は以下の通り。

> shapiro.test(res2)

	Shapiro-Wilk normality test

data:  res2
W = 0.9714, p-value = 0.5005
# 良いっぽい

> Box.test(res2,lag=15,type="Ljung")

	Box-Ljung test

data:  res2
X-squared = 10.045, df = 15, p-value = 0.8169
# 良いっぽい