計量経済学メモ:ボックス=ジェンキンスのモデル選択法、R

・本稿の内容
ボックスジェンキンスの時系列モデル選択法に関するメモです。最初にモデル選択法のプロセスを簡単に確認したあとに、Enders[2014](新谷・藪訳[2019])の練習問題とRを用いて、実際にモデル選択を実施してみます。

・本文

Ⅰ:ボックス=ジェンキンスのモデル選択法の流れ
Ⅱ:Rでボックス=ジェンキンスのモデル選択法を実践

Ⅰ:ボックス=ジェンキンスのモデル選択法の流れ

以下の3段階でモデルを選択する。

1.童貞同定
対象の時系列データを可視化し、自己相関(以下、 $ACF$ )と偏自己相関(以下、 $PACF$ )を確認する。理論上の $ACF$ 、 $PACF$ と比較し、推定するモデルの候補を決める。

2.推定
候補に挙がったモデルを推定し、データの当てはまりや、係数の推定値を分析する。定常で当てはまりが良く、倹約的なモデルを選択する。

3.診断
推定したモデルの残差がホワイトノイズであるどうかを調べる。ホワイトノイズであれば、推定したモデルは適切なモデルとみなすことができる。*1

Ⅱ:Rでボックス=ジェンキンスのモデル選択法を実践

新谷・藪訳[2019]の第2章の練習問題問11を例にボックス=ジェンキンスのモデル選択法を実践していく。*2

以下で使用するデータはここにアップされているSIM2.XLSの系列 $y3$ である。
$y3$ は

$y_t = 0.7y_{t-1}-0.49y_{t-2}+\varepsilon_{t}$

から発生させたものである。

・問11の(a)

以下のコードを実行し、y3系列をグラフ化する。

#使用するライブラリ
library("tidyverse")
library("rugarch")

#データ読み込み
data <- read.xls("C:/Users/81803/Desktop/SIM2.xls")

#y3系列を抽出
Y3 <- data[,"y3"]

#tidyverseのggplotでグラフを作成
data%>%
  ggplot()+
  geom_line(aes(x = X,y= y3))+
  labs(title = "y3系列")+
  theme(
    plot.title    = element_text(hjust = 0.5),   
  )

図1; $y3$ 系列のグラフ

図1をみると、平均と分散がおおむね安定しているように見える。*3
次に以下のコードを実行し、 $ACF$ と $PACF$ を求める。

#ACF,PACF
acf.Y3 <- acf(Y3,lag.max=20,tck=.02,xlab="",ylab="",main="")
acf.Y3
pacf.Y3 <- pacf(Y3,lag.max=20,tck=.02,xlab="",ylab="",main="")
pacf.Y3

図2:ACF

図3:PACF

$AR(2)$ の理論上の $ACF$ は0に減衰し*4、理論上の $PACF$ は3次以上で0となる。図2を見ると、 $ACF$ は振動しながら0に収束しているように見える。図3を見ると3次以上でほぼ0であり、 $AR(2)$ モデルが候補のモデルとして考えられる。

・問11の(b)

ここでは、(a)の結果、 $AR(1)$ モデルを推定するべきだと判断したとして、 $AR(1)$ モデルを推定し、そのモデルを診断してみる。
以下のコードを実行し、 $AR(1)$ モデルを推定する。

#AR(1)モデルを推定
spec.ar1 <- arfimaspec(mean.model = list(armaOrder=c(1,0),include.mean=FALSE))
fit.ar1 <- arfimafit(spec=spec.ar1,data=Y3)
fit.ar1

実行結果の一部を以下に抜粋する。

図4:推定結果

実行結果より、推定されたモデルは

$y_t = \underset{(5.29)}{0.467}y_{t-1}+e_{t}$

である。※()は $t$ 値。図4には記載されていないが、 $AIC$ (赤池情報量基準)の値は0.695である。

次に残差の診断を行うために、以下のコードを実行する。

#リュン=ボックス検定
res.ar1 <- fit.ar1@fit$residuals
Box.test(res.ar1,lag=1,type="Ljung-Box")
Box.test(res.ar1,lag=2,type="Ljung-Box")
Box.test(res.ar1,lag=3,type="Ljung-Box")
#以下、適当なラグ次数まで同様に実施。

実行結果は以下である。

図5:リュン=ボックス検定の結果

図5を見るとすべてのラグ次数でQ統計量の値(X-squared)が大きい(つまり、 $p$ 値が小さい。)ため、「残差に自己相関が無い」という帰無仮説を棄却できる。よって、系列相関が存在する可能性があるため、 $AR(1)$ モデルは不適切なモデルであると判断できる。*5

・問11の(c)

問11の(b)で実施したことを $ARMA(1,1)$ で同様に実施する。本稿では省略するが、 $ARMA(1,1)$ も $AR(1)$ と同様に不適切なモデルであることが確認できる。

・問11の(d)
ここでは、(a)の結果、 $AR(2)$ モデルを推定するべきだと判断したとして、 $AR(2)$ モデルを推定し、そのモデルを診断してみる。
以下のコードを実行し、 $AR(2)$ モデルを推定する。

#AR(2)モデルを推定
spec.ar2 <- arfimaspec(mean.model = list(armaOrder=c(2,0),include.mean=FALSE))
fit.ar2 <- arfimafit(spec=spec.ar2,data=Y3)
fit.ar2

実行結果の一部を以下に抜粋する。

図6:推定結果

実行結果より、推定されたモデルは

$y_t = \underset{(7.88)}{0.692}y_{t-1}\underset{(-5.47)}{-0.48}y_{t-2}+e_{t}$

である。※()は $t$ 値。図6には記載されていないが、 $AIC$ (赤池情報量基準)の値は0.454である。
推定された係数を見ると、真のデータ生成過程である

$y_t = 0.7y_{t-1}-0.49y_{t-2}+\varepsilon_{t}$

の係数とほぼ同じ値になっている。さらに、AICの値が問11(b)で推定した $AR(1)$ モデルよりも小さくなっており、データのあてはまりも改善していると言える。*6

次に残差の診断を行うために、以下のコードを実行する。

#リュン=ボックス検定
res.ar2 <- fit.ar2@fit$residuals
Box.test(res.ar2,lag=1,type="Ljung-Box")
Box.test(res.ar2,lag=2,type="Ljung-Box")
Box.test(res.ar2,lag=3,type="Ljung-Box")
#以下、適当なラグ次数まで同様に実施。

実行結果は以下である。

図7:リュン=ボックス検定の結果

図7を見るとすべてのラグ次数でQ統計量の値(X-squared)が小さい(つまり、 $p$ 値が大きい。)ため、「残差に自己相関が無い」という帰無仮説を棄却できない。よって、系列相関が存在するとは言えず、 $AR(2)$ モデルは適切なモデルであると判断できる。

・参考文献
ウォルター・エンダース著,新谷元嗣・薮友良訳 (2019)『実証のための計量時系列分析』有斐閣

・参考サイト
ウォルター・エンダース著,新谷元嗣・薮友良訳 (2019)のサポートページwww.fbc.keio.ac.jp
David Gabauer氏のホームページsites.google.com
パッケージ「rugarch」のドキュメント
https://cran.r-project.org/web/packages/rugarch/rugarch.pdf

*1:推定したモデルの残差に自己相関[系列相関]が残っているなら、その自己相関を説明できるより良いモデルが存在するという意味です。

*2:問題文の全文掲載は怒られそうなので控えます・・・

*3:本来は異常値や構造変化が生じている可能性、変数が非定常過程である可能性などをしかるべき方法でチェックするべきですが、本稿では省略します。

*4: $y_t = 0.7y_{t-1}-0.49y_{t-2}+\varepsilon_{t}$ の場合は振動しながら0に収束します。詳しくは新谷・藪訳[2019]のP68～71を参照してください。