R Time Series Analysis 時系列解析(14) SARIMAモデル

SARIMAモデルとは

ARIMAモデルの変化形であるSARIMAモデルの確認を忘れていた。SARIMAモデルとは、通常のARIMA(p,d,q)と季節階差に関するARIMA(k,l,m)を合わせたモデル。季節階差とは、例えば「前年同期との差」のこと。1周期がsである場合に、季節階差を

$\triangle_{s} y_{t} \equiv y_{t} - y_{t-s} = (1 - B^s) y_{t}$

と書くことが出来る。l回季節階差をとった場合は $\triangle_{s}^{l} y_{t}$ と書く。 $\triangle_{s}^{l} y_{t}$ がAR過程に従うのであれば、

$\Phi_{s} (B) \equiv 1 - \phi_{s1} B^s - \phi_{s2} B^2s - \cdots - \phi_{sk} B^ks$

として、

$\Phi_{s} (B) \triangle_{s}^{l} y_{t} = a_{t}$ (1)

と書ける。MA過程部分を

$\Theta_{s} (B) \equiv 1 + \theta_{s1} B^s + \theta_{s2} B^2s + \cdots + \theta_{sk} B^ks$

として(1)に追加すれば、

$\Phi_{s} (B) \triangle_{s}^{l} y_{t} = \Theta_{s} (B) a_{t}$ (2)

と書ける。(2)は $\triangle_{s}^{l} y_{t}$ がARIMA(k,l,m)に従うモデル。このモデルは $t, t-s, t-2s, \cdots$ 間の関係のみを含んでいるため、これに通常の $y_{t}$ に関するARIMA(p,d,q)を重ねる。

$\Phi(B) \Phi_{s} (B) \triangle^{d} \triangle_{s}^{l} y_{t} = \Theta (B) \Theta_{s} (B) a_{t}$ (3)

もちろん、

$\Phi(B) = 1 - \phi_{1} B - \phi_{2} B^2 - \cdots - \phi_{p} B^p$

$\Theta(B) = 1 + \theta_{1} B - \theta_{2} B^2 + \cdots + \theta_{q} B^q$

$B y_{t} = y_{t-1}, B a_{t} = a_{t-1}$

という意味。…なので、(3)は季節階差がARIMA(k,l,m)に従う系列が、ARIMA(p,d,q)に従うモデルを意味する。これをSARIMA(seasonal ARIMA)と呼ぶ。

RでSARIMAモデルを使ってみる

Rではarima()関数の引数seasonalに季節階差のARIMAモデルのk,l,mおよび周期を与えてやることでSARIMAモデルを使うことが出来る。季節性が無いと話にならないので、野菜のデータを使うことにした。東京都中央卸売市場の月別入荷量の推移らしい。ここではダウンロードしたデータのうち、玉ねぎのみを残し、縦軸がYear, 横軸がMonthのデータ（平成元年～23年）を用いる。

まず、下準備とグラフ図示、自己相関係数のチェック。

yasai <- read.csv("yasai.csv",header=TRUE,row.names=1)
yasai_dat <- rep(0,nrow(yasai)*ncol(yasai))
k <- 1

for (i in 1:nrow(yasai)){
  for (j in 1:ncol(yasai)) {
   yasai_dat[k] <- yasai[i,j]
   k <- k + 1
  }
}

yasai_ts <- ts(yasai_dat,start=c(1989,1),frequency=12)
par(mfrow=c(2,1))
plot(yasai_ts)
acf(yasai_ts)
par(mfrow=c(1,1))

面白い程に波打っている。例年、入荷量が最大になる月はほぼ同じで、前年の結果との相関が強いことがわかる。当たり前といえば当たり前。

f:id:nakhirot:20130716025743p:plain

とりあえず検定うんぬんの話は置いておき、SARIMAを当てはめる。おそらく一回階差をとると線形トレンドが除去されて定常過程になりそうなのでARIMA(1,1,1)とする(p,q=1は慣例)。季節階差についても（この場合は前年同期との差）線形下降トレンドを描いているのでARIMA(1,1,1)とする(これもp,q=1は慣例)。

kisetsu <- list(order=c(1,1,1),period=12)
(sarima111 <- arima(yasai_ts,order=c(1,1,1),seasonal=kisetsu,transform.pars=FALSE)) 
#引数transform.parsはデフォルトでTRUEで、AR部分の係数については、定常性の条件を満たす範囲で尤度を最大化する。

結果が出てきた。

Call:
arima(x = yasai_ts, order = c(1, 1, 1), seasonal = kisetsu, transform.pars = FALSE)

Coefficients:
         ar1      ma1    sar1     sma1
      0.2956  -0.9174  0.0592  -0.7878
s.e.  0.0750   0.0400  0.0788   0.0478

sigma^2 estimated as 850404:  log likelihood = -2174.83,  aic = 4359.67

…というわけで、結果は

$\Phi(B) \Phi_{s} (B) \triangle^{d} \triangle_{s}^{l} y_{t} = \Theta (B) \Theta_{s} (B) a_{t}$

$\Phi(B) = (1 - 0.2956B)$

$\Phi_{s} (B) = (1 - 0.0592B)$

$\Theta(B) = (1 - 0.9174B)$

$\Theta_{s} (B) = (1 - 0.7878B)$

$\triangle_{s}^{12} \triangle y_{t} = (1-B^{12})(1-B) y_{t}$

となることがわかった。ついでに、モデルの残差について調べる。

par(mfrow=c(2,1))
resi <- sarima111$resid
yasai_tshat <- yasai_ts - resi
plot(yasai_ts,type="l",main="東京卸売市場たまねぎ入荷量(トン)")
lines(yasai_tshat,lty=2,col=2)

plot(resi,type="l",col=4,main="残差系列")
abline(h=mean(resi))
par(mfrow=c(1,1))

結果は下記の通り。残差はARIMAのノイズ項によって発生する。下記グラフの赤線がモデルによる推定値となる。下のグラフはノイズ項による影響。

f:id:nakhirot:20130716025425p:plain