R Time Series Analysis 時系列解析(10) ARIMAモデル2 検証

前回に引き続き、ARIMAモデル。今回はARIMA(p,d,q)を利用することを前提として、そのパラメータの決定をする。Rではarima()関数がパラメータを最尤推定してくれる。

データの下準備は以下。株価のデータを使ってみる。

stock <- read.csv("http://k-db.com/site/jikeiretsu.aspx?c=7201-T&hyouji=&download=csv",skip=1,header=TRUE)

stockday <- matrix(stock[,1],ncol=1)
stockdaynum <- matrix(seq(1,nrow(stock),by=1),ncol=1)
stockindex <- cbind(stockday,stockdaynum)

delete <- c(2,3,4,6,7)
stock <- stock[,-delete]
stock <- ts(stock,start=1,freq=1)
stock <- stock[order(stock[,1]),]

stocksub <- matrix(rep(0,(nrow(stock)-1)*2),ncol=2)
for (i in 1:(nrow(stock)-1)) {
  stocksub[i,2] <- stock[i+1,2]-stock[i,2]
  stocksub[i,1] <- i
}
  
stocksub <- ts(stocksub,start=1,freq=1)

par(mfrow=c(2,1),ps=15)
ts.plot(stock[,2],type="l",col=1,lty=1,
        main="日産自動車(7201)の株価(終値)120626-130701",
        ylab="株価(円)")
kekka <- acf(stock[,2],main="自己相関係数")

f:id:nakhirot:20130702020341p:plain

1回階差をとった後に"定常"に近いとみなし、階差d=1としてarima関数を使う。ARモデルの次数pとMAモデルの次数qはとりあえず2としておく。

(arima212 <- arima(stock[,2],order=c(2,1,2),transform.pars=FALSE))

以下のような結果を返してくれる。

Call:
arima(x = stock[, 2], order = c(2, 1, 2), transform.pars = FALSE)

Coefficients:
          ar1      ar2     ma1     ma2
      -1.4276  -0.7219  1.4359  0.8017
s.e.   0.2715   0.2361  0.2396  0.1835

sigma^2 estimated as 398.3:  log likelihood = -1098.82,  aic = 2207.63

上記から、以下のように表すことが出来るとわかる。

${\triangle^{2}y_{t} = -1.4276\triangle^{2}y_{t-1} -0.7219\triangle^{2}y_{t-p} + a_{t} + 1.4359a_{t-1} + 0.8017a_{t-q},\:\:\:\:a_t \sim N(0, 398.3)}$

s.eは各係数の標準誤差。これで、時系列を満たすARIMAモデルを推定出来た。あとは、推定結果の評価が必要。

パラメータの検定

ARIMAモデル

${\triangle^{d}y_{t} = m + \phi_{1}\triangle^{d}y_{t-1} + \ldots + \phi_{p}\triangle^{d}y_{t-p} + a_{t} + \theta_{1} a_{t-1} + \ldots + \theta _{q}a_{t-q},\:\:\:\:a_t \sim N(0, \sigma^2)}$

のパラメータ

${\beta \equiv (\phi_{1}, \ldots , \phi_{p}, \theta_{1}, \ldots ,\theta _{q} )}$

は(近似的に)多変量正規分布に従うことが知られている。帰無仮説としてこれらのパラメータが0であるとして、t検定を行って有意性を判断する。

b <- arima212$coef　#係数
V <- arima212$var.coef　#標準誤差
t <- numeric(4)

for (j in 1:4) {
  t[j] <- b[j]/sqrt(V[j,j])
}
names(t) <- c("t_ar1","t_ar2","t_ma1","t_ma2")
{(t<0)&(pnorm(t)<0.05)} | {(t>=0)&pnorm(t)>0.95)}
#文字化け対策のため一番外側の"("を"{"に変えています

t_ar1 t_ar2 t_ma1 t_ma2 
TRUE TRUE TRUE TRUE

となって、係数が0であるという帰無仮説は棄却される。

残差

モデルと実績値の残差がホワイトノイズatとみなせるかどうか。つまりは、残差がランダムなatと見なせれば良いし、そうでない場合は異なる2時点間に関係を持つ(モデルに含まれない)他の項が存在する事になる。これをチェックする代表的な検定にLjung-Box検定というものがあるらしい。RではBox.test関数を使う。

p.value <- numeric(20)

for(i in 1:20) {
  result <- Box.test(arima212$resid,lag=i,type="Ljung-Box")
  p.value[i] <- result$p.value
}

plot(p.value,main="Ljung-Box test p value",xlab="lag",ylim=c(0,1))
abline(h=0.05)

帰無仮説は残差の系列に自己相関が無いこと。この元でp.valueが常に0.05を上回る結果となったため、残差系列は(20次までは)ホワイトノイズと見なせる。

f:id:nakhirot:20130703024152p:plain

定常性

ARモデルの時にも少し触れたが定常性は係数の条件で表現が出来る。一般に、ARIMAモデルの定常性の条件は、「特性方程式

${\lambda^{p} - \hat{\phi_{1}}\lambda^{p-1} - \ldots - \hat{\phi_{p}} = 0}$

の解の絶対値が1より小さいこと」となるので、これをpolyrootという関数を使ってチェックする。

test <- abs(polyroot(c(-b[2],-b[1],1)))
(test < 1)

結果は、

TRUE TRUE

となるので、定常性の条件も満たされるとわかる。

以上3診断をパス出来るかどうか、複数のモデルが3診断をパスした場合は、arima()の結果で返ってくるAICの値を参考にして決めるのがオーソドックスらしい。1