はじめに
ARIMAモデルの復習として、以前のARモデルをRで試してみるの記事同様に記録を残します。
前回同様、基本的にはRと時系列(2)に沿って実施しているだけですが、
Nileデータについても追加でやっています。
また、パラメータ推定については実装はせずにforecastパッケージのauto.arimaを利用しました。
ARMAモデルとARIMAモデルについて
ARはAutoRegressiveのことで自己回帰のことでした。
これに誤差の移動平均（Moving Average）$\sum_{j=1}^{q}{b_je_{t-j}}$を加えたものが、ARMAモデルです。
数式にすると次です。
y_t = \sum_{i=1}^{p}{a_iy_{i-1}} + e_t + \sum_{j=1}^{q}{b_je_{t-j}}
さらにy_tのd階の差分$\Delta^dy_d$のARMAモデルを自己回帰和分移動平均
(ARIMA:AutoRegressive Integrated Moving Average)モデルと言います。
ARIMA($p$,$d$,$q$)で表し、式では次のようになります。
y_t - y_{t-d} = \sum_{i=1}^{p}{a_iy_{i-1}} + e_t + \sum_{j=1}^{q}{b_je_{t-j}}
上式を$y_t$で整理すると、$y_t$の項に$y_{t-d}$に足される形なので、和分となっていると思います。
なお、差分で表現するのはトレンド除去の対応など時系列データの整形でよくあることなので、
差分でデータ整形しARMAモデル適用=ARIMAモデルという理解で多分良いのかなと思います。
ARIMAモデルのパラメータ推定方法
ARモデルの時はユールウォーカー法で必要なパラメータを予測しましたが、
今回のARIMA($p$,$d$,$q$)は、各パラメータのすべての組み合わせの中から
AICなどの情報量基準を使い、最も小さい組み合わせを選ぶようです。
参考としているRと時系列(2)では、別途実装しパラメータ推定を行っていますが、
forecastパッケージのauto.arimaを使えば、次のように計算を簡単に実行できました。
> library(forecast)
> auto.arima(lh, ic="aic", stepwise=T, trace=T)
ARIMA(2,0,2) with non-zero mean : 66.42642
ARIMA(0,0,0) with non-zero mean : 82.09291
ARIMA(1,0,0) with non-zero mean : 64.75832
ARIMA(0,0,1) with non-zero mean : 68.10389
ARIMA(0,0,0) with zero mean : 224.6837
ARIMA(2,0,0) with non-zero mean : 64.50375
ARIMA(3,0,0) with non-zero mean : 64.18482
ARIMA(4,0,0) with non-zero mean : 65.84091
ARIMA(3,0,1) with non-zero mean : 64.47047
ARIMA(2,0,1) with non-zero mean : 65.20321
ARIMA(4,0,1) with non-zero mean : 66.41819
ARIMA(3,0,0) with zero mean : 80.68744
Best model: ARIMA(3,0,0) with non-zero mean
Series: lh
ARIMA(3,0,0) with non-zero mean
Coefficients:
ar1 ar2 ar3 mean
0.6448 -0.0634 -0.2198 2.3931
s.e. 0.1394 0.1668 0.1421 0.0963
sigma^2 estimated as 0.1949: log likelihood=-27.09
AIC=64.18 AICc=65.61 BIC=73.54
同じようにNileデータもやってみました。
> auto.arima(Nile, ic="aic", stepwise=T, trace=T)
ARIMA(2,1,2) with drift : Inf
ARIMA(0,1,0) with drift : 1298.645
ARIMA(1,1,0) with drift : 1283.346
ARIMA(0,1,1) with drift : 1270.309
ARIMA(0,1,0) : 1296.697
ARIMA(1,1,1) with drift : 1267.637
ARIMA(2,1,1) with drift : 1269.132
ARIMA(1,1,2) with drift : 1269.134
ARIMA(0,1,2) with drift : 1268.544
ARIMA(2,1,0) with drift : 1279.282
ARIMA(1,1,1) : 1267.255
ARIMA(0,1,1) : 1269.091
ARIMA(1,1,0) : 1281.48
ARIMA(2,1,1) : 1268.896
ARIMA(1,1,2) : 1268.923
ARIMA(0,1,2) : 1267.957
ARIMA(2,1,0) : 1277.484
ARIMA(2,1,2) : Inf
Best model: ARIMA(1,1,1)
Series: Nile
ARIMA(1,1,1)
Coefficients:
ar1 ma1
0.2544 -0.8741
s.e. 0.1194 0.0605
sigma^2 estimated as 20177: log likelihood=-630.63
AIC=1267.25 AICc=1267.51 BIC=1275.04
ただし、Rで計量時系列分析：AR, MA, ARMA, ARIMAモデル, 予測によると、
このパラメータ推定ですが探索的に行なっていることから、
真のパラメータ以外のところをベストモデルと判定してしまうことがあるようです。
そのため、ある程度あたりがついている場合、探索範囲を狭めるのが良いようです。
推定されたパラメータを元に予測してみる
せっかくauto.arimaでパラメータ推定したので予測をしたいと思います。
推定したパラメータでarimaモデルを作ってから、predict関数で予測します。前回と同じです。
ということで、まずはlhデータを。
> lh.arima<-arima(lh,order=c(3,0,0))
> lh.pr<-predict(lh.arima,n.ahead=10)
> lh.SE1<-lh.pr$pred+2*lh.pr$se
> lh.SE2<-lh.pr$pred-2*lh.pr$se
> ts.plot(lh,lh.pr$pred,lh.SE1,lh.SE2,gpars=list(lt=c(1,2,3,3),col=c(1,2,4,4)))
次に、Nileデータを。
> Nile.arima<-arima(Nile,order=c(1,1,1))
> Nile.pr<-predict(Nile.arima,n.ahead=10)
> Nile.SE1<-Nile.pr$pred+2*Nile.pr$se
> Nile.SE2<-Nile.pr$pred-2*Nile.pr$se
> ts.plot(Nile,Nile.pr$pred,Nile.SE1,Nile.SE2,gpars=list(lt=c(1,2,3,3),col=c(1,2,4,4)))
ということで、簡単ですがARIMAモデルのパラメータ推定と予測を行ってみました。
なお、ARIMAモデルの診断についてはtsdiagを使うそうですが、
これについては後日この記事に追記するか、別の記事として改めて取り上げたいと思います。