9.7 fableによるARIMAモデル化

どうやってARIMA()はモデルを選択しているのか?

fableパッケージのARIMA()関数は、単位根検定とAICc最小化、そして最尤法を組み合わせてARIMAモデルを得るHyndman-Khandakarアルゴリズム (Hyndman & Khandakar, 2008) を使っています。ARIMA()の引数を指定することで、アルゴリズムの多くの変形が使えます。ここで述べるのは、デフォルトでの動作です。

自動ARIMAモデル化のためのHyndman-Khandakarアルゴリズム
  1. 差分の次数\(0 \le d\le 2\)をKPSS検定を繰り返して決める。
  1. 次に、\(d\)回差分を取ったデータでAICcを最小化して、\(p\)\(q\)の値を選ぶ。アルゴリズムは\(p\)\(q\)の全ての可能な組み合わせを考えるのではなく、ステップワイズ探索を使ってモデル空間を横断する。
  1. 最初に適合させる4つのモデルは、
    • ARIMA\((0,d,0)\)
    • ARIMA\((2,d,2)\)
    • ARIMA\((1,d,0)\)
    • ARIMA\((0,d,1)\)
    \(d=2\)でなければ、定数を含める。\(d \le 1\)なら、もう1つ追加で適合させるモデルは、
    • 定数なしのARIMA\((0,d,0)\)
  1. (a)のステップで適合させた中で(最小のAICc値を持つ)最良モデルを「当座のモデル」とする。
  1. 当座のモデルの変形を以下のようにして考える。
    • 当座のモデルから、\(\pm1\)だけ\(p\)\(q\)の両方か一方を変えてみる。
    • 当座のモデルから、\(c\)を除いたり、含めたりしてみる。
    ここまでで最良と考えられるモデル(現在のモデルかこれらの変形の1つ)が新しい当座のモデルになる。
  1. (c)のステップをAICcがそれ以上小さくならなくなるまで繰り返す。
Hyndman-Khandakarのステップワイズ探索プロセスの例示

図 9.11: Hyndman-Khandakarのステップワイズ探索プロセスの例示

9.11は、Hyndman-KhandakarアルゴリズムがどうやってARMA次数の空間を横断しているか、ある事例で示したものです。グリッドは、左上のARMA(\(0,0\))から始まって、AR次数は縦軸に沿って、MA次数は横軸に沿って、それぞれ増えていく格好で、ARMA(\(p,q\))の組み合わせをカバーしています。

オレンジ色のセルは、アルゴリズムが最初に検討したモデルの集合を示しています。この例では、これらのモデルの中では、ARMA(2,2)モデルが最小のAICc値でした。これを「当座のモデル」と言い、黒色の丸で示しています。アルゴリズムは次に、青色の矢印で示した近隣のモデルを探索します。もしより良いモデルが見つかれば、それが新しい「当座のモデル」になります。この例では、新しい「当座のモデル」はARMA(3,3)モデルです。アルゴリズムはこのような格好で、より良いモデルが見つからなくなるまで続けます。この例では、返ってきたモデルはARMA(4,2)でした。

デフォルトのプロセスでは、より良いモデルを見つけると、近隣モデル全ての探索を完了せずに、直ちに新しい「当座のモデル」へ切り替えます。greedy=FALSEと指定すれば、完全な近隣探索を行います。

デフォルトのプロセスはまた、探索を高速化するため、いくつかの近似値を使っています。引数でapproximation=FALSEと指定すると、これらの近似値を避けることができます。これらの近似値のため、あるいは、ステップワイズ手順を使っているため、AICcが最小のモデルを見つけ損ねることがあり得ます。引数でstepwise=FALSEと指定すると、モデルのずっと大きな集合を探索します。引数の完全な解説は、ヘルプを参照してください。

モデル化の手順

(非季節性)時系列データにARIMAモデルを適合させる際は、以下の手順が有益な一般的アプローチを提供してくれます。

  1. データをプロットして、異常値を特定する。
  2. 必要なら、(Box-Cox変換を使って)データを変換し、分散を安定化させる。
  3. データが非定常なら、定常になるまでデータの一つ目差分を取る。
  4. ACF/PACFを検討する: ARIMA(\(p,d,0\))モデルとARIMA(\(0,d,q\))モデルのどちらが適切か?
  5. 選んだモデルを試す、そして、AICcを使ってより良いモデルを探索する。
  6. 選んだモデルからの残差を、ACFプロットやportmanteau検定でチェックする。残差がホワイトノイズらしくなければ、モデルを修正して試す。
  7. 残差がホワイトノイズらしく見えるようになったら、予測を計算する。

Hyndman-Khandakarアルゴリズムが面倒を見るのは、3–5のステップだけです。ですから、それを使っても、他のステップは依然として自分で面倒を見る必要があります。

9.12は、このプロセスを要約したものです。

ARIMAモデルを使った予測のための一般的プロセス

図 9.12: ARIMAモデルを使った予測のための一般的プロセス

事例: 中央アフリカ共和国の輸出

9.13に示す、中央アフリカ共和国の輸出にこの手順を適用してみましょう。

global_economy %>%
  filter(Code == "CAF") %>%
  autoplot(Exports) +
  labs(title="中央アフリカ共和国の輸出",
       y="GDP比 (%)")
中央アフリカ共和国の輸出、GDP比

図 9.13: 中央アフリカ共和国の輸出、GDP比

  1. 時間プロットは、全体として減少傾向なので、いくらか非定常性なようです。1994年の改善は、軍事政権を打倒した新政権が初期の成功を収めたためです。その後、政情不安が経済を一層落ち込ませてしまいましたが。

  2. 分散が変化している証拠はないので、Box-Cox変換はしません。

  3. 非定常性に対処するため、データの一つ目差分を取ります。差分データを図9.14に示します。

    global_economy %>%
      filter(Code == "CAF") %>%
      gg_tsdisplay(difference(Exports), plot_type='partial')
    中央アフリカ共和国の輸出の差分データでの時間プロット、ACFプロット、PACFプロット

    図 9.14: 中央アフリカ共和国の輸出の差分データでの時間プロット、ACFプロット、PACFプロット

    今や、定常に見えます。

  4. 9.14に示したPACFは、AR(2)モデルを示唆しています。ですから、最初のモデル候補はARIMA(2,1,0)になります。ACFはMA(3)モデルを示唆しています。ですから、次の候補はARIMA(0,1,3)になります。

  5. ARIMA(2,1,0)モデルとARIMA(0,1,3)モデルの両方、および、2つの自動化されたモデル選択(一つはデフォルトのステップワイズを使うもの、もう一つは働き者でより大きなモデル空間を探索するもの)を適合させます。

    caf_fit <- global_economy %>%
      filter(Code == "CAF") %>%
      model(arima210 = ARIMA(Exports ~ pdq(2,1,0)),
            arima013 = ARIMA(Exports ~ pdq(0,1,3)),
            stepwise = ARIMA(Exports),
            search = ARIMA(Exports, stepwise=FALSE))
    
    caf_fit %>% pivot_longer(!Country, names_to = "Model name",
                             values_to = "Orders")
    #> # A mable: 4 x 3
    #> # Key:     Country, Model name [4]
    #>   Country                  `Model name`         Orders
    #>   <fct>                    <chr>               <model>
    #> 1 Central African Republic arima210     <ARIMA(2,1,0)>
    #> 2 Central African Republic arima013     <ARIMA(0,1,3)>
    #> 3 Central African Republic stepwise     <ARIMA(2,1,2)>
    #> 4 Central African Republic search       <ARIMA(3,1,0)>
    glance(caf_fit) %>% arrange(AICc) %>% select(.model:BIC)
    #> # A tibble: 4 × 6
    #>   .model   sigma2 log_lik   AIC  AICc   BIC
    #>   <chr>     <dbl>   <dbl> <dbl> <dbl> <dbl>
    #> 1 search     6.52   -133.  274.  275.  282.
    #> 2 arima210   6.71   -134.  275.  275.  281.
    #> 3 arima013   6.54   -133.  274.  275.  282.
    #> 4 stepwise   6.42   -132.  274.  275.  284.

    4つのモデルのAICc値はほぼ同じになりました。当てはめたモデルのうち、完全探索のARIMA(3,1,0)のAICc値が最小で、ARIMA(2,1,0)とARIMA(0,1,3)がたいして変わらない値で続きます。後者の2つは、ACFプロットとPACFプロットから候補としたものです。自動化ステップワイズによる選択はARIMA(2,1,2)モデルを特定しましたが、4つのモデルの中でAICc値は最大でした。

  6. ARIMA(3,1,0)モデルからの残差のACFプロットでは、全ての自己相関は臨界値内にあり、残差の動きはホワイトノイズのようだと言えそうです。

    caf_fit %>%
      select(search) %>%
      gg_tsresiduals()
    ARIMA(3,1,0)モデルからの残差のプロット

    図 9.15: ARIMA(3,1,0)モデルからの残差のプロット

    portmanteau検定は大きなp値を返しており、これまた、残差はホワイトノイズだと示唆しています。

    augment(caf_fit) %>%
      filter(.model=='search') %>%
      features(.innov, ljung_box, lag = 10, dof = 3)
    #> # A tibble: 1 × 4
    #>   Country                  .model lb_stat lb_pvalue
    #>   <fct>                    <chr>    <dbl>     <dbl>
    #> 1 Central African Republic search    5.75     0.569
  7. 9.16に、選択したモデルからの予測を示します。

    caf_fit %>%
      forecast(h=5) %>%
      filter(.model=='search') %>%
      autoplot(global_economy) +
      labs(title="中央アフリカ共和国の輸出",
       y="GDP比 (%)", level = "区間予測")
    中央アフリカ共和国の輸出の予測

    図 9.16: 中央アフリカ共和国の輸出の予測

    予測の平均は、ランダムウォーク(ARIMA(0,1,0)と同値です)から得るであろう予測に、とても似ていますね。この例では、AR項とMA項を追加しても点予測にたいした違いは出ませんでしたが、区間予測はランダムウォーク・モデルよりもずっと狭くなっています。

Rでの定数の理解

非季節性ARIMAモデルは以下のように書けます。 \[\begin{equation} \tag{9.3} (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, \end{equation}\] あるいは、同値として、 \[\begin{equation} \tag{9.4} (1-\phi_1B - \cdots - \phi_p B^p)(1-B)^d (y_t - \mu t^d/d!) = (1 + \theta_1 B + \cdots + \theta_q B^q)\varepsilon_t \end{equation}\] ただし、\(c = \mu(1-\phi_1 - \cdots - \phi_p )\)\(\mu\)\((1-B)^d y_t\)の平均です。fableパッケージは(9.3)式のパラメータ化を使っていますが、他のほとんどのRパッケージは(9.4)式を使っています。

このように、非季節性ARIMAモデルに定数を含めることは、予測に\(d\)次多項式のトレンドを持たせることと同値です。(定数を省くと、予測は\(d-1\)次多項式のトレンドを含むことになります。)\(d=0\)の場合は、\(\mu\)\(y_t\)の平均になる特殊ケースになります。

デフォルトでは、ARIMA()関数は定数を含めるべきか否か自動的に決めます。\(d=0\)\(d=1\)の場合、AICc値が改善するなら定数を含めます。\(d>1\)の場合、常に定数を省きます。2次やそれより高次のトレンドは、予測の際は特に危険だからです。

定数を含めるか否か指定するには、(lm()における切片と同様に)モデル式に01を書き足します。例えば、定数を含むARIMAモデルを自動選択するためには、ARIMA(y ~ 1 + ...)とします。同様に、ARIMA(y ~ 0 + ...)とすれば、定数を省けます。

特性根のプロット

(本節はより高度な内容なので、飛ばして構いません。)

(9.3)式は以下のように書き直せます。 \[\phi(B) (1-B)^d y_t = c + \theta(B) \varepsilon_t\] ただし、\(\phi(B)= (1-\phi_1B - \cdots - \phi_p B^p)\)\(B\)\(p\)次多項式で、\(\theta(B) = (1 + \theta_1 B + \cdots + \theta_q B^q)\)\(B\)\(q\)次多項式です。

モデルの定常性制約は\(\phi(B)\)\(p\)個の複素根が単位円の外側にあることで、反転可能性制約は\(\theta(B)\)\(q\)個の複素根が単位円の外側にあることです。ですから、複素単位円との関係で根をプロットすれば、モデルが反転可能性や定常性に近いのか見ることができます。

根の替わりに、根の逆数をプロットすれば、制約は全て単位円の内側にあることになるので、より簡単になります。Rでは簡単にできます。中央アフリカ共和国の輸出へ当てはめたモデルARIMA(3,1,0)では、図9.17が得られます。

gg_arma(caf_fit %>% select(Country, search))
中央アフリカ共和国の輸出へ当てはめたモデルARIMA(3,1,0)の特性根の逆数

図 9.17: 中央アフリカ共和国の輸出へ当てはめたモデルARIMA(3,1,0)の特性根の逆数

プロットの3つのオレンジ色の点は、それぞれ多項式\(\phi(B)\)の根に対応しています。3つとも単位円の内側にあります。fableは当てはめたモデルを定常かつ反転可能なものにしているので、予想通りです。単位円に近い根が1つでもあると数値的に不安定かもしれないので、そのモデルは予測に向いていないことになります。

ARIMA()関数は単位円の外側に根の逆数があるモデルを決して返しません。また、ARIMA()関数が自動選択するモデルが単位円に近い根を含むこともありません。従って、ARIMA()が返すモデルよりもAICc値が良いモデルを見つけ出せることが時々ありますが、そうしたモデルには問題がある可能性があります。

参考文献

Hyndman, R. J., & Khandakar, Y. (2008). Automatic time series forecasting: The forecast package for R. Journal of Statistical Software, 27(1), 1–22. [DOI]