前回の記事では計量時系列分析とは何ぞや?みたいなところをやりましたので、今回はいろはのイともいえるARIMAまわりから始めていこうと思います。
ということで改めて、使用テキストはいつものこちらです。
経済・ファイナンスデータの計量時系列分析 (統計ライブラリー)
- 作者: 沖本竜義
- 出版社/メーカー: 朝倉書店
- 発売日: 2010/02/01
- メディア: 単行本
- 購入: 4人 クリック: 101回
- この商品を含むブログ (6件) を見る
以下タイトルにのっとってRで各モデルの挙動を見ながらやっていきます。
必要なRパッケージ
{forecast}をインストールして展開して下さい。Rそのものの初心者向け説明はここでは全面的に割愛するので、適宜何かしらの初心者向け説明をご参照あれ。
今回のモデルで目指すもの
前回の記事では、要は「自己相関が大事よー」という話を何度もしました。ということは、時系列モデリングで最も大事なこともやはり「自己相関をモデリングする」ことだとも言えます。
例えば、ここでは1次の自己相関を持つ時系列について考えてみます(沖本本p.19)。ここでととが相関を持つようなモデルを構築しようと思ったら、とのそれぞれのモデルに共通の成分を含めるのが手っ取り早いはずです。つまり、
のようにモデリングすれば、共通の成分bを通じてととが自己相関を持つことが期待されます。またもう一つの方法としては、のモデルにを含めてしまうというものが考えられます。つまり、
のようにモデリングするということです。これならととが自己相関を持つのは明らかですね。これら2つのうち、前者はMA(Moving Average: 移動平均)の考え方、後者はAR(Autoregressive: 自己回帰)の考え方を表しています。これこそが今回のお題です。
MA(Moving Average:移動平均)過程
のように、ととが共通成分を持っているのがMA過程です。で、前回の記事を踏まえて考えると、多分それはホワイトノイズあたりが適切だなーと思われるわけです。そこで実際にそう仮定して1次MA過程をモデリングすると、
※参考
ホワイトノイズが共通成分となって、これをベースとして自己相関が決まるという仕組みです。
ちなみにこの過程から生じる時系列にははっきりとした特徴があり、だととが同じ方向に動く傾向があるためグラフは滑らかになり、だととが逆の方向に動く傾向があるためグラフはギザギザになる、ということが言えます。
ところで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)
見ての通り、(上)だとやや滑らかな、(下)だとややギザギザなグラフになることが分かりますね。今回は1次MA過程のみ定式化しましたが、全く同じ理屈でn次まで拡張できます。
なお、MA過程の厳密な性質を知る上でMA特性方程式について知っておくと、後で役立つかもしれません。
この方程式の全ての解の絶対値が1より大きい時、MA過程は「反転可能」と言われます。
AR(Autoregressive:自己回帰)過程
一方で、上の式のように直接ととの間に何かしらの関係性を持たせてやることで、つまり過程が自身の過去に回帰された形で、自己相関を表現しようとするのがAR過程です。1次AR過程を定式化すると、
と表されます。この式の形から何となく分かった人もいるかと思いますが、であればこのAR過程は定常になりますが、そうでなければどんどん値が増大して発散します(explosive:爆発的と称する)。ちなみにMA過程と同じくだととが同じ方向に動く傾向があるためグラフは滑らかになり、だととが逆の方向に動く傾向があるためグラフはギザギザになるということが言えます。
この辺の振る舞いを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以上だとエラーが出る
MA過程と同じようになっていることが見て取れます。なお、MA過程のケースと同じく、簡単にn次AR過程まで拡張していくことが可能です。
AR過程にもMA過程と同じく、AR特性方程式があります。
この方程式の全ての解の絶対値が1より小さい大きい時、AR過程は定常過程となります。
ARMA(Autoregressive and Moving Average:自己回帰移動平均)過程
MA過程もAR過程も、特に競合するパートが数式の中にないので、容易に2つを合わせて定式化することができます。これがARMA過程です。(p,q)次ARMA過程は
と表せます。なお、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)
一目見て、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}関数*4とAIC(){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期先の値を予測するものとして
純粋予測値
最適予測値
MSE
95%信頼区間
ということで、であることを考えると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区間の予測値と信頼区間を表示
もちろん、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))
上の2枚と下の2枚とを見比べると、実は予測範囲としてはあまり変わらないことが分かります。ARIMA係数&次数推定で本当に大事なのは、線形トレンドの有無と和分次数が正確に推定できるかどうかだという話題になるんですが、これはまた改めて。