渋谷駅前で働くデータサイエンティストのブログ

元祖「六本木で働くデータサイエンティスト」です / 道玄坂→銀座→東京→六本木→渋谷駅前

Rで計量時系列分析:AR, MA, ARMA, ARIMAモデル, 予測

前回の記事では計量時系列分析とは何ぞや?みたいなところをやりましたので、今回はいろはのイともいえるARIMAまわりから始めていこうと思います。


ということで改めて、使用テキストはいつものこちらです。


経済・ファイナンスデータの計量時系列分析 (統計ライブラリー)

経済・ファイナンスデータの計量時系列分析 (統計ライブラリー)


以下タイトルにのっとってRで各モデルの挙動を見ながらやっていきます。


必要なRパッケージ


{forecast}をインストールして展開して下さい。Rそのものの初心者向け説明はここでは全面的に割愛するので、適宜何かしらの初心者向け説明をご参照あれ。


今回のモデルで目指すもの


前回の記事では、要は「自己相関が大事よー」という話を何度もしました。ということは、時系列モデリングで最も大事なこともやはり「自己相関をモデリングする」ことだとも言えます。


例えば、ここでは1次の自己相関を持つ時系列y_tについて考えてみます(沖本本p.19)。ここでy_ty_{t-1}とが相関を持つようなモデルを構築しようと思ったら、y_ty_{t-1}のそれぞれのモデルに共通の成分を含めるのが手っ取り早いはずです。つまり、


 \left \{ \begin{array}{l} ~ y_t = a + b \\ y_{t-1} = b + c \end{array} \right.


のようにモデリングすれば、共通の成分bを通じてy_ty_{t-1}とが自己相関を持つことが期待されます。またもう一つの方法としては、y_tのモデルにy_{t-1}を含めてしまうというものが考えられます。つまり、


 y_t = a y_{t-1} + b


のようにモデリングするということです。これならy_ty_{t-1}とが自己相関を持つのは明らかですね。これら2つのうち、前者はMA(Moving Average: 移動平均)の考え方、後者はAR(Autoregressive: 自己回帰)の考え方を表しています。これこそが今回のお題です。


MA(Moving Average:移動平均)過程


 \left \{ \begin{array} ~ y_t = a + b \\ y_{t-1} = b + c \end{array} \right.


のように、y_ty_{t-1}とが共通成分を持っているのがMA過程です。で、前回の記事を踏まえて考えると、多分それはホワイトノイズあたりが適切だなーと思われるわけです。そこで実際にそう仮定して1次MA過程をモデリングすると、


 y_t = \mu + \epsilon_t + \theta_1 \epsilon_{t-1}, \\ \epsilon_t \sim W.N.(\sigma^2)

※参考
 y_{t-1} = \mu + \epsilon_{t-1} + \theta_1 \epsilon_{t-2}


ホワイトノイズ\epsilon_t \sim W.N.(\sigma^2)が共通成分となって、これをベースとして自己相関が決まるという仕組みです。


ちなみにこの過程から生じる時系列にははっきりとした特徴があり、\theta_1 > 0だとy_ty_{t-1}が同じ方向に動く傾向があるためグラフは滑らかになり、\theta_1 < 0だとy_ty_{t-1}が逆の方向に動く傾向があるためグラフはギザギザになる、ということが言えます。


ところでRでこれを表現しようと思ったら、一番手っ取り早いのがarima.sim(){base}関数でMA次数&係数を指定してランダム生成することです。

> x1_ma<-arima.sim(n=200,list(ma=1)) #MA次数は1、係数も1
> plot(x1_ma)
> x2_ma<-arima.sim(n=200,list(ma=-1)) #MA次数は1、係数は-1
> plot(x2_ma)


f:id:TJO:20130711152740p:plain

f:id:TJO:20130711152751p:plain


見ての通り、\theta_1 > 0(上)だとやや滑らかな、\theta_1 < 0(下)だとややギザギザなグラフになることが分かりますね。今回は1次MA過程のみ定式化しましたが、全く同じ理屈でn次まで拡張できます。


なお、MA過程の厳密な性質を知る上でMA特性方程式について知っておくと、後で役立つかもしれません。


 1 + \theta_1 z + \theta_2 z^2 + \cdots + \theta_p z^p = 0


この方程式の全ての解の絶対値が1より大きい時、MA過程は「反転可能」と言われます。


AR(Autoregressive:自己回帰)過程


 y_t = a y_{t-1} + b


一方で、上の式のように直接y_ty_{t-1}との間に何かしらの関係性を持たせてやることで、つまり過程が自身の過去に回帰された形で、自己相関を表現しようとするのがAR過程です。1次AR過程を定式化すると、


 y_t = c + \phi_1 y_{t-1} + \epsilon_t,\\ \epsilon_t \sim W.N.(\sigma^2)


と表されます。この式の形から何となく分かった人もいるかと思いますが、 |\phi_1| < 1であればこのAR過程は定常になりますが、そうでなければどんどん値が増大して発散します(explosive:爆発的と称する)。ちなみにMA過程と同じく\phi_1 > 0だとy_ty_{t-1}が同じ方向に動く傾向があるためグラフは滑らかになり、\phi_1 < 0だとy_ty_{t-1}が逆の方向に動く傾向があるためグラフはギザギザになるということが言えます。


この辺の振る舞いをRで見てみましょう。

> x1_ar<-arima.sim(n=200,list(ar=0.5)) # AR次数は1、係数は0.5
> plot(x1_ar)
> x2_ar<-arima.sim(n=200,list(ar=-0.5)) # AR次数は1、係数は-0.5
> plot(x2_ar)
> x3_ar<-arima.sim(n=200,list(ar=1)) # AR次数は1、係数を1にしてみると…
Error in arima.sim(n = 200, list(ar = 1)) : 
  'ar' part of model is not stationary
# AR係数が1以上だとエラーが出る


f:id:TJO:20130711161323p:plain

f:id:TJO:20130711161332p:plain


MA過程と同じようになっていることが見て取れます。なお、MA過程のケースと同じく、簡単にn次AR過程まで拡張していくことが可能です。


AR過程にもMA過程と同じく、AR特性方程式があります。


 1 - \phi_1 z - \phi_2 z^2 - \cdots - \phi_p z^p = 0


この方程式の全ての解の絶対値が1より小さい大きい時、AR過程は定常過程となります。


ARMA(Autoregressive and Moving Average:自己回帰移動平均)過程


MA過程もAR過程も、特に競合するパートが数式の中にないので、容易に2つを合わせて定式化することができます。これがARMA過程です。(p,q)次ARMA過程は


y_t = c + \phi_1 y_{t-1} + \cdots + \phi_p y_{t-p} + \epsilon_t + \theta_1 \epsilon_{t-1} + \cdots + \theta_q \epsilon_{t-q},\\ \epsilon_t \sim W.N.(\sigma^2)


と表せます。なお、ARMA過程の反転可能性ですが、沖本本p.38から引用するとMA過程の反転可能性については

MA過程がAR(∞)過程に書き直せるとき、MA過程は反転可能といわれる。


とあります。実際に問題になるのはARMA過程全体の反転可能性なんですが、沖本本pp.39-40にもある通りARMA過程は本来定常AR過程と定常MA過程の足し算なので、全体としては通常なら定常過程になります。なので、単純にMA過程の部分の反転可能性が示せれば、ARMA過程全体も反転可能だということになります。


ARMA過程のRでの表し方はちょっと一旦置いといて、サクッと次にいっちゃいましょう。


ARIMA(Autoregressive, Integrated and Moving Average:自己回帰和分移動平均)過程


詳しい説明はまた単位根過程の話題になった時に書きますが、世の中には「直前の値に定常過程に従う何かの値を足すことで現在の値になる」時系列って結構多いんですよね。代表例は株価などのファイナンスデータです。そういう、1階階差が自分自身の足し算として表されるような過程を和分(Integrated)過程と呼びます。


これは、実はARMA過程にそのまま付け足すことができます。ただし、数式としては表現するのが面倒*1なので、数式は割愛。定義については沖本本p.106(ちょっと飛びます)にある通り、

d階差分をとった系列が定常かつ反転可能なARMA(p,q)過程に従う過程は次数(p,d,q)の自己回帰和分移動平均過程もしくはARIMA(p,d,q)と呼ばれる。


基本的にはARMA過程とほぼ同じで、しかも基本的にはd = 1と相場が決まっています*2。なので、普通はARIMA(p,1,q)過程だけ知っていれば十分です。


Rではarima.sim(){stats}関数で、関数名の通りARIMA次数(p,d,q)と係数を与えることでランダムに発生させることができます。

> x1_arima<-arima.sim(n=200,list(order=c(2,0,1),ar=c(0.2,-0.1),ma=0.1))
# 次数orderでARIMA(2,0,1)を指定、その後さらに個々の係数を指定
> plot(x1_arima)
> x2_arima<-arima.sim(n=200,list(order=c(1,1,2),ar=0.3,ma=c(0.3,-0.4)))
# 次数orderでARIMA(1,1,2)を指定、その後さらに個々の係数を指定
> plot(x2_arima)


f:id:TJO:20130711165609p:plain

f:id:TJO:20130711165616p:plain


一目見て、ARIMA(2,0,1)過程のグラフは定常で、ARIMA(1,1,2)過程のグラフは非定常であることが分かるんじゃないかと思います。


ARIMAモデルの推定


沖本本pp.40-53にある通り、最小二乗法(OLS)もしくは最尤法で推定するのが基本です。係数&次数の選択はAICでやるのが普通です。この辺はRだとauto.arima(){forecast}関数で大変お手軽にできます(最尤法を使用)。探索的方法とAICによるモデル選択によって最適なARIMA過程の係数と次数(p,d,q)を推定してくれます。


ただし、あくまでも探索的にAICを計算していくだけなので割と真の値以外のところに収まってしまうことが少なくないです。事前に何かしらの手掛かりを持っていて、特に次数(p,d,q)の探索範囲を絞れる場合は、自分で絞ってしまった方が確実です。

> auto.arima(x1_arima,max.p=2,max.q=2,stepwise=T,trace=T)
# AR次数の上限を2、MA次数の上限を2に制限
# さらに逐次探索ではなく階段式の高速(端折り)探索を選択
# 個々の探索におけるAICを表示するようにした

 ARIMA(2,0,2) with non-zero mean : 578.1241
 ARIMA(0,0,0) with non-zero mean : 592.2657
 ARIMA(1,0,0) with non-zero mean : 575.1168
 ARIMA(0,0,1) with non-zero mean : 574.037
 ARIMA(1,0,1) with non-zero mean : 576.7794
 ARIMA(0,0,2) with non-zero mean : 576.091
 ARIMA(1,0,2) with non-zero mean : 576.4368
 ARIMA(0,0,1) with zero mean     : 572.0696
 ARIMA(1,0,1) with zero mean     : 574.7596
 ARIMA(0,0,0) with zero mean     : 590.3755
 ARIMA(0,0,2) with zero mean     : 574.1006
 ARIMA(1,0,2) with zero mean     : 574.3327

 Best model: ARIMA(0,0,1) with zero mean # 設定通りにならなかった(笑)

Series: x1_arima 
ARIMA(0,0,1) with zero mean     

Coefficients:
         ma1
      0.3222
s.e.  0.0659

sigma^2 estimated as 1.001:  log likelihood=-283.99
AIC=571.98   AICc=572.04   BIC=578.58

> auto.arima(x2_arima,d=1,stepwise=T,trace=T)
# 和分次数を1(つまり単位根過程)に設定

 ARIMA(2,1,2) with drift         : 569.3479
 ARIMA(0,1,0) with drift         : 628.4033
 ARIMA(1,1,0) with drift         : 596.3629
 ARIMA(0,1,1) with drift         : 570.4191
 ARIMA(1,1,2) with drift         : 568.7291
 ARIMA(1,1,1) with drift         : 568.2689
 ARIMA(1,1,1)                    : 567.6839
 ARIMA(0,1,1)                    : 569.0287
 ARIMA(2,1,1)                    : 568.7358
 ARIMA(1,1,0)                    : 595.1957
 ARIMA(1,1,2)                    : 568.0267
 ARIMA(0,1,0)                    : 628.2451
 ARIMA(2,1,2)                    : 568.1838

 Best model: ARIMA(1,1,1) # これも設定通りじゃないし(笑)

Series: x2_arima 
ARIMA(1,1,1)                    

Coefficients:
          ar1     ma1
      -0.1649  0.7782
s.e.   0.1055  0.0721

sigma^2 estimated as 0.9675:  log likelihood=-280.83
AIC=567.66   AICc=567.78   BIC=577.55


なお、stepwise=Fにすることである程度網羅的に推定してくれます*3。ところで、見ての通り次数推定は結構外れてます(笑)。これは割とauto.arimaではよくあることなので、できればちゃんと一つ一つ手作業でも検証できると良いところ。そこで、arima(){forecast}関数*4AIC(){stats}関数を使って

> x1_arima.model.estimated<-Arima(x1_arima,order=c(0,0,1))
# auto.arimaで推定した結果でモデリング
> x1_arima.model.original<-Arima(x1_arima,order=c(2,0,1))
# 実際の設定通りにモデリング
> AIC(x1_arima.model.estimated)
[1] 573.8909
> AIC(x1_arima.model.original)
[1] 578.6755

> x2_arima.model.estimated<-Arima(x2_arima,order=c(1,1,1))
# auto.arimaで推定した結果でモデリング
> x2_arima.model.original<-Arima(x2_arima,order=c(1,1,2))
# 実際の設定通りにモデリング
> AIC(x2_arima.model.estimated)
[1] 567.6576
> AIC(x2_arima.model.original)
[1] 569.364


というように、細かく突き詰めて調べることもできます。今回のケースでは割とAICが拮抗しているので、誤判定されても仕方ないという感じですが、場合によってはこれでより適切なモデルが分かることもある*5ので、その都度しっかりチェックする習慣にしても良いでしょう。


ARIMAモデルからの予測


基本的には既にモデル式があるので、それを使って1期ずつ増やしていって予測すれば良いわけなんですが。当たり前ですが、予測値もばらつきます。時系列モデリングの世界ではMSE(Mean Squared Error:平均二乗誤差)をばらつきの目安として用いることになっています。これらのモデル式は、h期先の値を予測するものとして


純粋予測値
 y_{t+h} = (1 + \phi_1 + \phi_1^2 + \cdots + \phi_1^{h-1})c + \phi_1^h y_t
 + \epsilon_{t+h} + \phi_1 \epsilon_{t+h-1} + \phi_1^2 \epsilon_{t+h-2} + \cdots + \phi_1^{h-1} \epsilon_{t+1}
 = c {\displaystyle \sum_{k=1}^{h}}\phi_1^{k-1} + \phi_1^h y_t + {\displaystyle \sum_{k=0}^{h-1} \phi_1^k \epsilon_{t+h-k}}


最適予測値
 \hat{y}_{t+h|t} = c {\displaystyle \sum_{k=1}^{h}} \phi_1^{k-1} + \phi_1^h y_t = \frac{(1-\phi_1^h)c}{1-\phi_1} + \phi_1^h y_t


MSE
 MSE(\hat{y}_{t+h|t}) = E \left( {\displaystyle \sum_{k=0}^{h-1}} \phi_1^{k} \epsilon_{t+h-k} \right)^2 = \sigma^2 {\displaystyle\sum_{k=0}^{h-1}} \phi_1^{2k} = \frac{(1-\phi_1^{2h}) \sigma^2}{1-\phi_1^2}


95%信頼区間
 \left( \hat{y}_{t+h|t} - 1.96 \sqrt{MSE(\hat{y}_{t+h|t})}, \hat{y}_{t+h|t} + 1.96 \sqrt{MSE(\hat{y}_{t+h|t})} \right)


ということで、|\phi|<1であることを考えるとhが大きくなればなるほどMSEも大きくなる=末広がりになるということが分かるかと思います。ちなみに信頼区間の計算は現実には計算負荷の問題から困難とされ、通常はカルマンフィルタで推定するか、モンテカルロ・シミュレーションで疑似標本を発生させることで近似値を求めるようです。


で、Rではforecast(){forecast}関数で予測と信頼区間推定をそのまんまかけることができます。

> x1_forecast<-Arima(x1_arima,order=c(2,0,1))
> x2_forecast<-Arima(x2_arima,order=c(1,1,2))
> plot(forecast(x1_forecast,level=c(50,95),h=30))
> plot(forecast(x2_forecast,level=c(50,95),h=30))
# 信頼区間を50%, 95%で設定して表示(前者の方が狭い)
# 向こう30区間の予測値と信頼区間を表示


f:id:TJO:20130711175754p:plain

f:id:TJO:20130711175803p:plain


もちろん、auto.arima(){forecast}関数が返したモデルに対して予測を計算することもできます。

> x1_model<-auto.arima(x1_arima,max.p=2,max.q=2,stepwise=T,trace=T)

 ARIMA(2,0,2) with non-zero mean : 578.1241
 ARIMA(0,0,0) with non-zero mean : 592.2657
 ARIMA(1,0,0) with non-zero mean : 575.1168
 ARIMA(0,0,1) with non-zero mean : 574.037
 ARIMA(1,0,1) with non-zero mean : 576.7794
 ARIMA(0,0,2) with non-zero mean : 576.091
 ARIMA(1,0,2) with non-zero mean : 576.4368
 ARIMA(0,0,1) with zero mean     : 572.0696
 ARIMA(1,0,1) with zero mean     : 574.7596
 ARIMA(0,0,0) with zero mean     : 590.3755
 ARIMA(0,0,2) with zero mean     : 574.1006
 ARIMA(1,0,2) with zero mean     : 574.3327

 Best model: ARIMA(0,0,1) with zero mean     

> summary(x1_model)
Series: x1_arima 
ARIMA(0,0,1) with zero mean     

Coefficients:
         ma1
      0.3222
s.e.  0.0659

sigma^2 estimated as 1.001:  log likelihood=-283.99
AIC=571.98   AICc=572.04   BIC=578.58

Training set error measures:
                     ME     RMSE       MAE      MPE     MAPE      MASE
Training set 0.02058361 1.000744 0.8001445 43366.89 48853.11 0.8174317
                    ACF1
Training set 0.005273937
> x2_model<-auto.arima(x2_arima,max.p=2,d=1,max.q=2,stepwise=T,trace=T)

 ARIMA(2,1,2) with drift         : 569.3479
 ARIMA(0,1,0) with drift         : 628.4033
 ARIMA(1,1,0) with drift         : 596.3629
 ARIMA(0,1,1) with drift         : 570.4191
 ARIMA(1,1,2) with drift         : 568.7291
 ARIMA(1,1,1) with drift         : 568.2689
 ARIMA(1,1,1)                    : 567.6839
 ARIMA(0,1,1)                    : 569.0287
 ARIMA(2,1,1)                    : 568.7358
 ARIMA(1,1,0)                    : 595.1957
 ARIMA(1,1,2)                    : 568.0267
 ARIMA(0,1,0)                    : 628.2451
 ARIMA(2,1,2)                    : 568.1838

 Best model: ARIMA(1,1,1)                    

> summary(x2_model)
Series: x2_arima 
ARIMA(1,1,1)                    

Coefficients:
          ar1     ma1
      -0.1649  0.7782
s.e.   0.1055  0.0721

sigma^2 estimated as 0.9675:  log likelihood=-280.83
AIC=567.66   AICc=567.78   BIC=577.55

Training set error measures:
                      ME      RMSE       MAE       MPE     MAPE
Training set -0.07201374 0.9811826 0.7681955 -989.5171 2507.598
                  MASE        ACF1
Training set 0.8381687 -0.00697138

> plot(forecast(x1_model,level=c(50,95),h=30))
> plot(forecast(x2_model,level=c(50,95),h=30))


f:id:TJO:20130711180305p:plain

f:id:TJO:20130711180332p:plain


上の2枚と下の2枚とを見比べると、実は予測範囲としてはあまり変わらないことが分かります。ARIMA係数&次数推定で本当に大事なのは、線形トレンドの有無と和分次数が正確に推定できるかどうかだという話題になるんですが、これはまた改めて。

*1:本気で表現したらスーパー長ったらしい式が出来上がります

*2:というよりd > 1になるケースはあまり知られてないはず

*3:その代わり時間はかかる

*4:arima(){stats}関数のwrapper

*5:見ての通りauto.arimaのAICAIC(){stats}の結果と食い違うことがある