9.5 非季節性ARIMAモデル

差分化に自己回帰モデルや移動平均モデルを組み合わると、非季節性ARIMAモデルになります。ARIMAは AutoRegressive Integrated Moving Average の略語です。(この文脈では、“integration”は差分化の逆の意味です。)完全なモデルは以下のように書けます。 \[\begin{equation} y'_{t} = c + \phi_{1}y'_{t-1} + \cdots + \phi_{p}y'_{t-p} + \theta_{1}\varepsilon_{t-1} + \cdots + \theta_{q}\varepsilon_{t-q} + \varepsilon_{t} \tag{9.1} \end{equation}\] ただし、\(y'_{t}\)は差分系列(1回を超えて差分化されている場合もあります)。右辺の「予測変数」には、\(y_t\)のラグと誤差のラグの両方があります。これをARIMA(\(p, d, q\))モデルと言います。ただし、

\(p=\) 自己回帰パートの次数
\(d=\) 一つ目差分化の次数
\(q=\) 移動平均パートの次数

です。自己回帰モデルと移動平均モデルで使ったのと同じ定常性制約と反転可能性制約が、ARIMAモデルにも当てはまります。

これまで議論してきたモデルの多くは、表9.1に示すように、ARIMAモデルの特殊ケースに当たります。

表 9.1: ARIMAモデルの特殊ケース
ホワイトノイズ 定数なしのARIMA(0,0,0)
ランダムウォーク 定数なしのARIMA(0,1,0)
ドリフト付きランダムウォーク 定数ありのARIMA(0,1,0)
自己回帰 ARIMA(\(p\),0,0)
移動平均 ARIMA(0,0,\(q\))

このように部品を組み合わせてより複雑なモデルを作り始めると、後方シフト表記の方がずっと簡単になります。例えば、(9.1)式は、後方シフト表記では以下のように書けます。 \[\begin{equation} \tag{9.2} \begin{array}{c c c c} (1-\phi_1B - \cdots - \phi_p B^p) & (1-B)^d y_{t} &= &c + (1 + \theta_1 B + \cdots + \theta_q B^q)\varepsilon_t\\ {\uparrow} & {\uparrow} & &{\uparrow}\\ \text{AR($p$)} & \text{$d$次の差分化} & & \text{MA($q$)}\\ \end{array} \end{equation}\]

\(p\)\(d\)\(q\)の適切な値を選択するのは困難かもしれません。しかし、fableパッケージのARIMA()関数を使うと、自動的に選択してくれます。9.7節で、この関数がどうやって選択しているのか、自分でこれらの値を選択するにはどうしたら良いか、を学びます。

事例: エジプトの輸出

9.7は、エジプトの輸出をGDP比で1960年から2017年まで表したものです。

global_economy %>%
  filter(Code == "EGY") %>%
  autoplot(Exports) +
  labs(y = "GDP (%)", title = "エジプトの輸出")
エジプトの年次輸出、1960年以降のGDP比

図 9.7: エジプトの年次輸出、1960年以降のGDP比

以下のRコードを実行すると、非季節性ARIMAモデルが1つ自動的に選ばれます。

fit <- global_economy %>%
  filter(Code == "EGY") %>%
  model(ARIMA(Exports))
report(fit)
#> Series: Exports 
#> Model: ARIMA(2,0,1) w/ mean 
#> 
#> Coefficients:
#>          ar1      ar2      ma1  constant
#>       1.6764  -0.8034  -0.6896    2.5623
#> s.e.  0.1111   0.0928   0.1492    0.1161
#> 
#> sigma^2 estimated as 8.046:  log likelihood=-141.6
#> AIC=293.1   AICc=294.3   BIC=303.4

選ばれたのは、ARIMA(2,0,1)モデルでした。 \[ y_t = 2.56 + 1.68 y_{t-1} -0.80 y_{t-2} -0.69 \varepsilon_{t-1} + \varepsilon_{t} \] ただし、\(\varepsilon_t\)はホワイトノイズで、その標準偏差は\(2.837 = \sqrt{8.046}\)です。モデルからの予測を図9.8に示します。過去数十年明らかに存在するエジプト経済の循環を拾い上げていますね。

fit %>% forecast(h=10) %>%
  autoplot(global_economy) +
  labs(y = "GDP (%)", title = "エジプトの輸出", level = "区間予測")
エジプト輸出の予測

図 9.8: エジプト輸出の予測

ARIMAモデルの理解

ARIMA()関数は便利ですが、何にせよ自動化されたものは少々危険です。ですから、モデル選択の際は自動化手順に頼るとしても、モデルの動きの幾分かは理解しておく価値があります。

定数\(c\)は、これらモデルから得られる長期予測に大きな影響があります。

  • \(c=0\)\(d=0\)だと、長期予測はゼロになる。
  • \(c=0\)\(d=1\)だと、長期予測はゼロでない定数になる。
  • \(c=0\)\(d=2\)だと、長期予測は直線になる。
  • \(c\ne0\)\(d=0\)だと、長期予測はデータの平均になる。
  • \(c\ne0\)\(d=1\)だと、長期予測は直線になる。
  • \(c\ne0\)\(d=2\)だと、長期予測は2次曲線トレンドになる。(これはお薦めできないので、fableは許しません。)

\(d\)の値は、区間予測にも影響します。\(d\)の値が大きいほど、より急速に区間予測の幅が広がります。\(d=0\)だと、長期予測の標準偏差は過去のデータの標準偏差になるので、区間予測は本質的に同一になります。

\(d=0\)\(c\ne0\)だった図9.8に、この動きが見られます。区間予測は予測期間の最後の方ではほとんど同じ幅になっており、最後の点予測はデータの平均に近くなっています。

\(p\)の値は、データが循環を示す場合に重要です。循環的な予測を得るためには、\(p\ge2\)、かつ、パラメータがいくつかの追加条件を満たす必要があります。AR(2)モデルでは\(\phi_1^2+4\phi_2<0\)なら(エジプトの輸出モデルではそうなっています)、循環的動きが生じます。その場合、循環の平均周期は以下のようになります。16 \[ \frac{2\pi}{\text{arc cos}(-\phi_1(1-\phi_2)/(4\phi_2))} \]

ACFプロットとPACFプロット

データに適切な\(p\)\(q\)の値を決めるのに、時間プロットだけからは判断できないのが普通です。しかし、ACFプロット、および、密接に関連するPACFプロットも使えば、判断できる場合もあります。

ACFプロットは、異なる\(k\)値での\(y_t\)\(y_{t-k}\)の関係を測った自己回帰を表すことを思い起こしましょう。さて、もし\(y_t\)\(y_{t-1}\)に相関があれば、\(y_{t-1}\)\(y_{t-2}\)にも相関があるに違いありません。しかし、\(y_t\)\(y_{t-2}\)に相関があったとしても、\(y_{t-2}\)\(y_t\)の予測に使える情報が含まれているとは限りません。単に両方とも\(y_{t-1}\)につながっているために相関があるだけかもしれないからです。

この問題を克服するために、偏自己相関を使うことができます。これは、ラグ\(1, 2, 3, \dots, k - 1\)の影響を除去した後の\(y_{t}\)\(y_{t-k}\)の関係を測ります。ですから、ラグ1の偏自己相関はラグ1の自己相関と同じです。間に除去すべきものがないからです。偏自己相関は、自己回帰モデルの最後の項の係数としてそれぞれ推計できます。具体的には、ラグ\(k\)の偏自己相関\(\alpha_k\)は、AR(\(k\))モデルの\(\phi_k\)の推計値に等しいです。実際には、\(\alpha_k\)を計算するのに、これらの自己回帰モデルを一つ一つ全て適合させるよりも効率的なアルゴリズムがあり、同じ結果が得られます。

9.99.10は、図9.7で示したエジプトの輸出データのACFプロットとPACFプロットです。偏自己相関の臨界値は、通常の自己相関と同じ\(\pm 1.96/\sqrt{T}\)で、図9.10にあるようにプロットに描くのが普通です。

global_economy %>%
  filter(Code == "EGY") %>%
  ACF(Exports) %>%
  autoplot()
エジプト輸出のACF

図 9.9: エジプト輸出のACF

global_economy %>%
  filter(Code == "EGY") %>%
  PACF(Exports) %>%
  autoplot()
エジプト輸出のPACF

図 9.10: エジプト輸出のPACF

時間プロット、ACFプロット、PACFプロットを1つのコマンドで作成する便利なやり方は、plot_type = "partial"を指定してgg_tsdisplay()関数を使うことです。

データがARIMA(\(p\),\(d\),0)かARIMA(0,\(d\),\(q\))のモデルから来ている場合には、\(p\)\(q\)の値を決めるのに、ACFプロットとPACFプロットが助けになり得ます。\(p\)\(q\)の両方が正の場合には、両プロットは\(p\)\(q\)の適切な値を見つける助けになりません。

差分データのACFプロットとPACFプロットが以下のパターンを示すなら、データはARIMA(\(p\),\(d\),0)モデルに従っているのかもしれません。

  • ACFが指数関数的な減衰、もしくは、サインカーブを描いており、
  • PACFに有意なスパイクがラグ\(p\)にあるが、ラグ\(p\)超にはない。

差分データのACFプロットとPACFプロットが以下のパターンを示すなら、データはARIMA(0,\(d\),\(q\))モデルに従っているのかもしれません。

  • PACFが指数関数的な減衰、もしくは、サインカーブを描いており、
  • ACFに有意なスパイクがラグ\(q\)にあるが、ラグ\(q\)超にはない。

9.9のACFでは減衰するサインカーブのパターンがあり、図9.10のPACFではラグ4に最後の有意なスパイクがあります。これは、ARIMA(4,0,0)モデルからのデータに期待するものと一致します。

fit2 <- global_economy %>%
  filter(Code == "EGY") %>%
  model(ARIMA(Exports ~ pdq(4,0,0)))
report(fit2)
#> Series: Exports 
#> Model: ARIMA(4,0,0) w/ mean 
#> 
#> Coefficients:
#>          ar1      ar2     ar3      ar4  constant
#>       0.9861  -0.1715  0.1807  -0.3283    6.6922
#> s.e.  0.1247   0.1865  0.1865   0.1273    0.3562
#> 
#> sigma^2 estimated as 7.885:  log likelihood=-140.5
#> AIC=293.1   AICc=294.7   BIC=305.4

このモデルは、ARIMA()が特定したARIMA(2,0,1)モデルより、ほんの少し悪いだけです。(AICc値が294.70と、294.29と比べて少し大きくなっています。)

ARIMA()が探索できるpdq()の値を、特定の値の範囲内に指定することもできます。例えば、\(p\in\{1,2,3\}\)\(q\in\{0,1,2\}\)\(d=1\)のうちから最良なARIMAモデルを見つけるには、ARIMA(y ~ pdq(p=1:3, d=1, q=0:2))と指定できます。


  1. arc cosはコサインの逆関数です。あなたの電卓にもあるはずです。acosとかcos\(^{-1}\)のラベルになっているかもしれません。↩︎