R Time Series Analysis 時系列解析(9) ARIMAモデル1 紹介

ARモデルとMAモデルを組み合わせたものがARMAモデル。p次のARモデルAR(p)とq次のMAモデルMA(q)を組み合わせてARMA(p,q)と書く。以下のようになる。

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

パラメータを増やすとモデルの自由度が増える一方、パラメータに起因する誤差も増加するが、Box-JenkinsによるとARMAモデルのp,qは1～3程度で実用上問題がないらしい。

以上、ARモデル、MAモデル、ARMAモデルを勉強してきたが、現実のデータはなかなか当てはまらず、そうみなせるように変換を施すようにしなければならない。その発想で生まれたのがARIMA（Autoregressive integrated moving average process、自己回帰和分平均過程)モデルと理解した。

なお、「変換」のうち頻繁に使われるのは差分をとる操作。

${\triangle^{1}y_{t} \equiv y_{t} - y_{t-1}$

2回差分をとることを

${\triangle^{2}y_{t}=\triangle^{1}y_{t}-\triangle^{1}y_{t-1}$

と書き、一般化してd階について

${\triangle^{d}y_{t}=\triangle^{d-1}y_{t}-\triangle^{d-1}y_{t-1}$

と書くと、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)}$

d回差分をとることで定常過程となる場合にその確率を「d次の和分過程」と呼び、d回差分をとったものがARMA(p,q)過程となるものがARIMA(p,d,q)過程とのこと。

Rにはarima.sim( )関数というものがあり、ARIMA過程に従った標本データを作成することができる。

#パラメータの設定
d <- 1 #階差
no <- 200 #標本数
sd <- 1.5 #ホワイトノイズの標準偏差
phi <- c(0.5,0.3) #ARモデルの係数
theta <- c(0.3,0.2) #MAモデルの係数

#arima.sim()による標本の作成
p <- length(phi)
q <- length(theta)
mdl <- list(order=c(p,d,q),ar=phi,ma=theta)
y <- arima.sim(n=no,model=mdl,sd=sd

#グラフタイトル
title1 <- paste("ARIMA Model(",p,",",d,",",q,")")
title2 <- "ar ="
for(j in 1:p) {title2 <- paste(title2,phi[j]," ")}
title3 <- "ma ="
for(j in 1:q) {title3 <- paste(title3,theta[j]," ")}
title2 <- paste(title2,title3)
title <- c(title1,title2)

plot(y,type="l",main=title)

ARIMAモデルによる標本系列が出来た。

f:id:nakhirot:20130630190253p:plain

いろいろパラメータを変えてやってみる。

f:id:nakhirot:20130630191622p:plain

気づいたこと。（仮説）

階差dを増加させるとモデルは"滑らか"になる　⇒　あべこべな（dを増加させた方が振動が激しいグラフになる）気がするが後で考えよう。
MAモデルの係数(の絶対値) > ARモデルの係数(の絶対値)とすると、モデルのグラフはギザギザ(振動)が大きく表れる。

最後に、d＝1の場合のモデルが1階階差をとると、定常過程になることを確認する。

par(mfrow=c(2,1))
Dy <- rep(0,200)
for (i in 1:(length(y)-1)) {
  Dy[i] <- y[i+1]-y[i]
}
plot(Dy,main=paste(title1,"(Difference sequence)"),type="l")
acf(Dy)

確かに、定常過程過程（平均一定、分散一定、自己相関係数一定）となっているような様子が見受けられる。

f:id:nakhirot:20130630195103p:plain