読者です 読者をやめる 読者になる 読者になる

六本木で働くデータサイエンティストのブログ

元祖「銀座で働くデータサイエンティスト」です / 道玄坂→銀座→東京→六本木

Rでベイジアン動的線形モデルを学ぶ(5):説明変数のあるローカル・レベル・モデル

R DLM 時系列分析

何か月1回しか書かなくなりつつあるこのシリーズですが、中には@さんのようにツッコミ倒すのを楽しみにして下さっている方もおられるようなので、久しぶりに更新してみます。


もちろん参考文献は以下の2冊 + PDF book。お題はCommandeur本の第5章「説明変数のあるローカル・レベル・モデル」です。


状態空間時系列分析入門

状態空間時系列分析入門

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

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


ちなみにずっと{dlm}パッケージでやってますが、以前の記事でも紹介したようにもちろんStanなどMCMC系ソフトウェアでモデリング&パラメータ推定をかけることは可能です。{dlm}パッケージの挙動が気に入らない!という人はStanやBUGSでやるというのも良いかと思います。


説明変数


仰々しく説明変数とか言ってますが、要はこれまで内生変数だけ扱ってきましたがこいつに外生変数を加えて、回帰モデルにしてやろうというお話です。Commandeur本の第5章冒頭49ページには説明変数を加えた観測方程式として


y_t = \mu_t + \displaystyle \sum^{k}_{j=1} \beta_{jt} x_{jt} + \epsilon_t (5.1)


という式になると説明されています。ここでx_{jt}は連続な予測変数で、\beta_{jt}j = 1, \cdots, kに対して未知の回帰ウェイト、あるいは回帰係数*1\beta_t = \beta_{1t}の1つの予測変数に対して、モデルはt = 1, \cdots, nに対して


y_t = \mu_t + \beta_t x_t + \epsilon_t, \epsilon_t \sim NID(0, \sigma^2_{\epsilon})

\mu_{t+1} = \mu_t + \xi_t, \xi_t \sim NID(0, \sigma^2_{\xi})

\beta_{t+1} = \beta_t + \tau_t, \tau_t \sim NID(0, \sigma^2{\tau}) (5.2)


となります。ちなみに重回帰としてk個の説明変数をモデル化するにはさらにk本の状態方程式が個々の説明変数に対して1セット必要になります。なお上記3本目の撹乱項\tau_tは回帰関係を安定化させるために0に固定するのが通常だとか。


基本的には上記モデルの3本の式の撹乱項のパラメータ\sigma^2を固定したりフリーにするなりして、最尤法でモデル推定を行う流れになります。


確定的レベルと説明変数


要は\xi_tをゼロに固定するモデルです。ここで(5.2)式のように回帰モデルを立てた場合、Commandeur本pp.50-51でも指摘されているようにごく普通の古典的な線形回帰分析に相当します。


となると、PDF bookにあるように実はこれは{dlm}パッケージを用いるまでもなく、普通にlm関数で線形回帰すればおしまいなのでした。以下そのコード。ただしいつものUK driver KSIをd1として読み込み、対数変換済みtsオブジェクトにしておきます。

# まずUKdriverKSIデータを読み込む
> d1<-read.table("OxCodeIntroStateSpaceBook/Chapter_5/UKdriversKSI.txt",skip=1)
> colnames(d1)<-"logKSI"
> d1<-ts(log(d1),start=c(1969),frequency=12)
# 単純な時間変化で線形回帰したいのでただの1, 2, ...という数列を作る
> t<-1:nrow(d1)
# まずその線形回帰を計算する
> summary(lm(d1~t))

Call:
lm(formula = d1 ~ t)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.33649 -0.10850 -0.01256  0.10368  0.40749 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  7.5458427  0.0219747 343.387  < 2e-16 ***
t           -0.0014480  0.0001975  -7.333 6.31e-12 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1517 on 190 degrees of freedom
Multiple R-squared:  0.2206,	Adjusted R-squared:  0.2165 
F-statistic: 53.77 on 1 and 190 DF,  p-value: 6.313e-12
# 切片&偏回帰係数ともCommandeur本p.51の通りになった


ところで、Commandeur本p.51と付録AにあるようにこのUK driver KSIデータは石油価格との間に密接な関係があります。これはある意味当たり前で、ガソリン価格が上がればみんなガス代ケチって車に乗らなくなる→交通事故死傷者も減るし、逆なら死傷者は増えるわけです。これも当然ただの線形回帰で関係性を推定することができます。

# 石油価格データを読み込む
> x<-read.table("OxCodeIntroStateSpaceBook/Chapter_5/logUKpetrolprice.txt",skip=1)
> x<-ts(x,start=c(1969),frequency=12)

> summary(lm(d1~x))

Call:
lm(formula = d1 ~ x)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.37612 -0.09896 -0.01594  0.09077  0.36594 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  5.87873    0.20889  28.142  < 2e-16 ***
x           -0.67166    0.09173  -7.322 6.74e-12 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1517 on 190 degrees of freedom
Multiple R-squared:  0.2201,	Adjusted R-squared:  0.216 
F-statistic: 53.61 on 1 and 190 DF,  p-value: 6.742e-12
# 切片&偏回帰係数ともCommandeur本p.51の通りになった

> matplot(cbind(d1,predict(lm(d1~x),newdata=x)),type='l',lty=1,lwd=2)
# 図5.1をプロットしてみる

f:id:TJO:20150130160036p:plain


これ自体はどちらかというと動的線形モデルとかどうでもいいような、ただの線形回帰です。なので、確定的レベルではなく確率的レベルに問題設定を変更して解き直さないと醍醐味が分からない感があります。


確率的レベルと説明変数


今度は撹乱項をフリーにして、{dlm}パッケージの諸関数を使って最尤推定することを考えます。

# dlmModeReg関数で単回帰モデルを線形に加える
> build1<-function(params){
+     model<-dlmModPoly(order=1)+dlmModReg(x,addInt = FALSE)
+     V(model)<-exp(params[1])
+     diag(W(model))[1]<-exp(params[2])
+     model
+ }
> fit1<-dlmMLE(d1,rep(0,2),build1) # まず推定
> model1<-build1(fit1$par) # 最尤推定量を突っ込み直してモデル構築
> V(model1)
            [,1]
[1,] 0.002347963
> diag(W(model1))[1]
[1] 0.01166743
# この辺のパラメータはCommandeur本p.55とほぼ一致
> filt1<-dlmFilter(d1,model1) # カルマンフィルタ
> smooth1<-dlmSmooth(filt1) # カルマンスムージング
> smooth1$s[1,1]
[1] 6.82036 # mu1でCommandeur本p.55と一致
> smooth1$s[1,2]
[1] -0.2610497 # beta1でCommandeur本p.55と一致
> matplot(cbind(d1,dropFirst(smooth1$s[,1])),type='l',lty=1,lwd=2)
# じゃあプロットしてみると。。。

f:id:TJO:20150130162531p:plain


あれー。。。PDF bookと全く同じ不一致になってるorz これ、どう見ても単回帰モデルの側に切片が入ってないからだと思うんですが、ちょっと調べた範囲ではこれを{dlm}パッケージ及びdlmModReg関数の範疇で解決する方法が見つからず。。。うーむ、どうしよう(汗)。

追記(2016年12月27日)

という指摘をいただいたので、最後のmatplotのところだけ直しました。

> matplot(cbind(d1,dropFirst(smooth1$s[,1]+smooth1$s[,2]*x)),type='l',lty=1,lwd=2,ylab='')

f:id:TJO:20161227142843p:plain

はい、おかげさまで期待した通りになりました。有難うございます。


ちなみに第4章で季節調整をやっているのに今回やっていないのを見れば分かる通り、当たり前ですが第5章のモデルは季節調整がないので不適切です。この辺は次章及び次々章以降で取り上げますよ、ということで。

*1:一応重回帰に発展可能なように立てられている