9.9 季節性ARIMAモデル

ここまでは、非季節性データと非季節性ARIMAモデルだけを見てきました。しかし、ARIMAモデルは、幅広い季節性データのモデル化もできます。

季節性ARIMAモデルは、ここまで見てきたARIMAモデルに季節項を追加して形成します。以下のように書きます。

ARIMA \(\underbrace{(p, d, q)}\) \(\underbrace{(P, D, Q)_{m}}\)
\({\uparrow}\) \({\uparrow}\)
モデルの モデルの
非季節パート 季節パート

ただし、\(m =\)季節周期(例えば、一年当たりの観測値の数)です。モデル中の季節パートには大文字を、非季節パートには小文字を使うことにします。

モデルの季節パートは、モデルの非季節パートと似た項から成りますが、季節周期の後方シフトが関わっています。例えば、ARIMA(1,1,1)(1,1,1)\(_{4}\)モデル(定数なし)は、四半期データ(\(m=4\))のためのモデルですが、以下のように書けます。 \[ (1 - \phi_{1}B)~(1 - \Phi_{1}B^{4}) (1 - B) (1 - B^{4})y_{t} = (1 + \theta_{1}B)~ (1 + \Theta_{1}B^{4})\varepsilon_{t} \]

追加した季節項を、非季節項に掛けるだけです。

ACF/PACF

ARモデル、もしくは、MAモデルであれば、その季節パートは、PACFとACFの季節ラグから分かります。例えば、ARIMA(0,0,0)(0,0,1)\(_{12}\)モデルでは、以下のようになります。

  • ACFのラグ12にスパイクがあり、その他に有意なスパイクはなく、
  • PACFの季節ラグ(つまり、ラグ12、24、36、…)は、指数関数的に減衰している。

同様に、ARIMA(0,0,0)(1,0,0)\(_{12}\)モデルでは、以下のようになります。

  • ACFの季節ラグは、指数関数的に減衰しており、
  • PACFのラグ12だけに、有意なスパイクがある。

季節性ARIMAモデルで季節パートの適切な次数を考える際は、季節ラグに注意を限定すべきです。

モデル化の手順は非季節性データに対するものとほぼ同じですが、モデルの非季節パートだけでなく、ARとMAの季節項も選ぶ必要がある点が違っています。プロセスは事例で示した方が良さそうです。

事例: 米国余暇・娯楽産業の月次就業者数

9.18に示した、米国余暇・娯楽産業の月次就業者数の2001年1月から2019年9月までのデータを使って、季節性ARIMAモデルを述べていきます。

leisure <- us_employment %>%
  filter(Title == "Leisure and Hospitality",
         year(Month) > 2000) %>%
  mutate(Employed = Employed/1000) %>%
  select(Month, Employed)
autoplot(leisure, Employed) +
  labs(title = "米国: 余暇・娯楽産業",
       y="就業者数 (百万人)")
米国余暇・娯楽産業の月次就業者数、2001-2019年

図 9.18: 米国余暇・娯楽産業の月次就業者数、2001-2019年

データには強い季節性と非線形なトレンドがあり、明らかに非定常なので、まず季節差分を取ります。季節差分データを図9.19に示します。

leisure %>%
  gg_tsdisplay(difference(Employed, 12),
               plot_type='partial', lag=36) +
  labs(title="季節差分", y="")
米国余暇・娯楽産業の月次就業者数の季節差分

図 9.19: 米国余暇・娯楽産業の月次就業者数の季節差分

まだ明らかに非定常なので、さらに一つ目差分を取ったものが、図9.20です。

leisure %>%
  gg_tsdisplay(difference(Employed, 12) %>% difference(),
               plot_type='partial', lag=36) +
  labs(title = "二階差分", y="")
米国余暇・娯楽産業の月次就業者数の二階差分

図 9.20: 米国余暇・娯楽産業の月次就業者数の二階差分

9.20のACFとPACFに基づいて、適切なARIMAモデルを見つけるのが狙いです。ACFのラグ2での有意なスパイクは、非季節MA(2)成分を示唆しています。ACFのラグ12での有意なスパイクは、季節MA(1)成分を示唆しています。ですから、1回の一つ目差分、1回の季節差分、非季節MA(2)、季節MA(1)を指すARIMA(0,1,2)(0,1,1)\(_{12}\)モデルをまず選びます。PACFから始めていたら、つまり、PACFを使ってモデルの非季節パートを選び、ACFを使ってモデルの季節パートを選ぶと、ARIMA(2,1,0)(0,1,1)\(_{12}\)モデルを選ぶことになります。自動選択モデルも含めましょう。 stepwise=FALSEapproximation=FALSEを設定すると、Rにもう一仕事して良いモデルを見つけさせることになります。そうするとずっと時間がかかりますが、モデル化する系列が1つだけなら、追加でかかる時間は問題になりません。

fit <- leisure %>%
  model(
    arima012011 = ARIMA(Employed ~ pdq(0,1,2) + PDQ(0,1,1)),
    arima210011 = ARIMA(Employed ~ pdq(2,1,0) + PDQ(0,1,1)),
    auto = ARIMA(Employed, stepwise = FALSE, approx = FALSE)
  )
fit %>% pivot_longer(everything(), names_to = "Model name",
                     values_to = "Orders")
#> # A mable: 3 x 2
#> # Key:     Model name [3]
#>   `Model name`                    Orders
#>   <chr>                          <model>
#> 1 arima012011  <ARIMA(0,1,2)(0,1,1)[12]>
#> 2 arima210011  <ARIMA(2,1,0)(0,1,1)[12]>
#> 3 auto         <ARIMA(2,1,0)(1,1,1)[12]>
glance(fit) %>% arrange(AICc) %>% select(.model:BIC)
#> # A tibble: 3 × 6
#>   .model       sigma2 log_lik   AIC  AICc   BIC
#>   <chr>         <dbl>   <dbl> <dbl> <dbl> <dbl>
#> 1 auto        0.00142    395. -780. -780. -763.
#> 2 arima210011 0.00145    392. -776. -776. -763.
#> 3 arima012011 0.00146    391. -775. -775. -761.

ARIMA()関数は、unitroot_nsdiffs()を使って(行う季節差分の回数)\(D\)を、 unitroot_ndiffs()を使って(行う一つ目差分の回数)\(d\)を、それらが指定されていない場合、決めます。その他のパラメータ(\(p,q,P, Q\))は全て、非季節性ARIMAモデルと同様、AICcの最小化により決めます。

これら3つの当てはめたモデルのAICc値は、自動選択モデルがやや良いですが、似たようなものになりました。私たちの二番目の「当て推量」であるARIMA(2,1,0)(0,1,1)\(_{12}\)は、結果的には自動選択モデルのARIMA(2,1,0)(1,1,1)\(_{12}\)にかなり近かったわけです。

最良モデルの残差を図9.21に示します。

fit %>% select(auto) %>% gg_tsresiduals(lag=36)
当てはめたモデルARIMA(2,1,0)(1,1,1)\(_{12}\)からの残差

図 9.21: 当てはめたモデルARIMA(2,1,0)(1,1,1)\(_{12}\)からの残差

ACFで有意なスパイクは36のうち小さいのが1つ(ラグ11)だけで、ホワイトノイズと矛盾しません。確認のため、注意深くモデル内のパラメータ数に一致するよう自由度を設定して、Ljung-Box検定を使います。

augment(fit) %>%
  filter(.model == "auto") %>%
  features(.innov, ljung_box, lag=24, dof=4)
#> # A tibble: 1 × 3
#>   .model lb_stat lb_pvalue
#>   <chr>    <dbl>     <dbl>
#> 1 auto      16.6     0.680

p値は大きく、残差がホワイトノイズに似ていることが確認できました。

こうして、必要なチェックに合格した予測に使える季節性ARIMAモデルがやっと得られました。モデルから次の3年間を予測して図9.22に示します。予測は季節パターンをとても上手く取り込んでおり、増加トレンドは最近のパターンを延長しています。予測のトレンドは二階差分がもたらしたものです。

forecast(fit, h=36) %>%
  filter(.model=='auto') %>%
  autoplot(leisure) +
  labs(title = "米国: 余暇・娯楽産業",
       y="就業者数 (百万人)", level = "区間予測")
米国余暇・娯楽産業の月次就業者数のARIMA(2,1,0)(1,1,1)\(_{12}\)モデルを使った予測。80%と95%の区間予測も示す

図 9.22: 米国余暇・娯楽産業の月次就業者数のARIMA(2,1,0)(1,1,1)\(_{12}\)モデルを使った予測。80%と95%の区間予測も示す

事例: オーストラリアでのコルチコステロイド薬売上高

二番目の事例では、オーストラリアでのコルチコステロイド薬の月次売上高の予測を試みます。解剖学的治療のための医薬品分類スキームでは、H02薬として知られているものです。

h02 <- PBS %>%
  filter(ATC2 == "H02") %>%
  summarise(Cost = sum(Cost)/1e6)
h02 %>%
  mutate(log(Cost)) %>%
  pivot_longer(-Month) %>%
  mutate(name = recode(name, "Cost" = "売上高", "log(Cost)" = "対数変換後")) %>% 
  ggplot(aes(x = Month, y = value)) +
  geom_line() +
  facet_grid(name ~ ., scales = "free_y") +
  labs(y="", title="コルチコステロイド薬の売上高 (H02)")
オーストラリアでのコルチコステロイド薬の売上高 (百万オーストラリア・ドル/月)。下のパネルは対数変換後のデータ

図 9.23: オーストラリアでのコルチコステロイド薬の売上高 (百万オーストラリア・ドル/月)。下のパネルは対数変換後のデータ

1991年7月から2008年6月までのデータをプロットしたのが、図9.23です。水準につれて分散が少し増加しているので、分散を安定化させるため対数を取ります。

データには強い季節性があり、明らかに非定常なので、季節差分を取ります。季節差分を取ったデータを図9.24に示します。この段階では、もう一段差分を取るべきか否か明らかではありません。取らないことにしますが、選択は明白ではありません。

最後の2つ、3つの観測値は、それ以前のデータとは違う(より変動が大きい)ように見えます。売上の報告が遅れた場合、データが遡って修正されることがあるせいかもしれません。

h02 %>% gg_tsdisplay(difference(log(Cost), 12),
                     plot_type='partial', lag_max = 24)
オーストラリアでのコルチコステロイド薬の売上高の季節差分 (百万オーストラリア・ドル/月)

図 9.24: オーストラリアでのコルチコステロイド薬の売上高の季節差分 (百万オーストラリア・ドル/月)

季節差分データのプロットでは、PACFのラグ12と24にスパイクがありますが、ACFの季節ラグにはスパイクはありません。これは、季節AR(2)項を示唆しています。非季節性ラグでは、PACFに3つ有意なスパイクがあり、AR(3)項の可能性を示唆しています。ACFのパターンはどの単純なモデルも示唆していません。

ということで、この初期の分析からは、このデータで可能性のあるモデルはARIMA(3,0,0)(2,1,0)\(_{12}\)になります。このモデルとそのいくつかの変形を適合させ、AICc値を計算したものが表9.2です。

表 9.2: H02月次売上高データに適用したさまざまなARIMAモデルのAICc値
モデル AICc
ARIMA(3,0,1)(0,1,2)\(_{12}\) -485.5
ARIMA(3,0,1)(1,1,1)\(_{12}\) -484.2
ARIMA(3,0,1)(0,1,1)\(_{12}\) -483.7
ARIMA(3,0,1)(2,1,0)\(_{12}\) -476.3
ARIMA(3,0,0)(2,1,0)\(_{12}\) -475.1
ARIMA(3,0,2)(2,1,0)\(_{12}\) -474.9
ARIMA(3,0,1)(1,1,0)\(_{12}\) -463.4

これらのモデルのうち最良(つまり、AICc値が最小)なのはARIMA(3,0,1)(0,1,2)\(_{12}\)モデルです。このモデルからのイノベーション残差を図9.25に示します。

fit <- h02 %>%
  model(ARIMA(log(Cost) ~ 0 + pdq(3,0,1) + PDQ(0,1,2)))
fit %>% gg_tsresiduals(lag_max=36)
H02の月次売上高データに適用したARIMA(3,0,1)(0,1,2)$_{12}$モデルからのイノベーション残差

図 9.25: H02の月次売上高データに適用したARIMA(3,0,1)(0,1,2)\(_{12}\)モデルからのイノベーション残差

augment(fit) %>%
  features(.innov, ljung_box, lag = 36, dof = 6)
#> # A tibble: 1 × 3
#>   .model                                             lb_stat lb_pvalue
#>   <chr>                                                <dbl>     <dbl>
#> 1 ARIMA(log(Cost) ~ 0 + pdq(3, 0, 1) + PDQ(0, 1, 2))    50.7    0.0104

ACFに2つ、3つ有意なスパイクがあり、このモデルはLjung-Box検定を合格しません。それでもこのモデルを予測に使うことはできますが、残差に相関があるため、区間予測は不正確かもしれません。

次に、自動化ARIMAアルゴリズムを使って試してみます。全ての引数をデフォルト値のままにしてARIMA()を実行すると、ARIMA(2,1,0)(0,1,1)\(_{12}\)モデルになりました。stepwise=FALSEapproximation=FALSEとしてARIMA()を実行すると、ARIMA(2,1,3)(0,1,1)\(_{12}\)モデルになりました。しかし、どちらのモデルも36ラグでのLjung-Box検定を合格しません。全ての検定を合格するモデルを見つけるのが不可能なことがあるのです。

テストセットでの評価

これまでに当てはめたモデルのいくつかを、データの最後の2年分のテストセットを使って、比較してみましょう。では、1991年7月から2006年6月までのデータを使ってモデルを適合させ、2006年7月から2008年6月までの売上高を予測します。結果は、表9.3に要約しました。

表 9.3: H02月次売上高データに適用したさまざまなARIMAモデルの2006年7月から2008年6月までのテストセットでのRMSE値
.model RMSE
ARIMA(3,0,1)(1,1,1)\(_{12}\) 0.0619
ARIMA(3,0,1)(0,1,2)\(_{12}\) 0.0621
ARIMA(2,1,1)(0,1,1)\(_{12}\) 0.0622
ARIMA(2,1,2)(0,1,1)\(_{12}\) 0.0623
ARIMA(2,1,4)(0,1,1)\(_{12}\) 0.0627
ARIMA(2,1,3)(0,1,1)\(_{12}\) 0.0628
ARIMA(3,0,1)(0,1,1)\(_{12}\) 0.0630
ARIMA(3,0,2)(0,1,1)\(_{12}\) 0.0630
ARIMA(2,1,0)(0,1,1)\(_{12}\) 0.0630
ARIMA(3,0,1)(0,1,3)\(_{12}\) 0.0630
ARIMA(3,0,3)(0,1,1)\(_{12}\) 0.0631
ARIMA(3,0,2)(2,1,0)\(_{12}\) 0.0651
ARIMA(3,0,1)(2,1,0)\(_{12}\) 0.0653
ARIMA(2,1,0)(1,1,0)\(_{12}\) 0.0666
ARIMA(3,0,1)(1,1,0)\(_{12}\) 0.0666
ARIMA(3,0,0)(2,1,0)\(_{12}\) 0.0668

テストセットでのRMSE値に基づくと、手動で選んだモデルが最良モデルに近いですが、ARIMA()で自動で選んだモデルもそれほど離れていません。

AICc値を使ったモデル比較が可能なのは、全てのモデルの差分次数が同じである場合だけです。他方、テストセットを使った比較は、予測がどうやって生成されたかと関係なく、いつも有効です。ですから、上記の表には、季節差分だけのモデルと一つ目差分と季節差分の両方を行うモデルとを含めることができます。その前のAICc値を含む表では、季節差分は取ったが一つ目差分は取っていないモデルだけを比較していました。

ここで検討したモデルに、全ての残差検定を合格したものはありません。実際には、たとえ全ての検定を合格できなくても、見つけた最良モデルを使うことになるのが普通です。

(テストセットでのRMSE値が二番目に小さく、季節差分だけのモデルの中ではAICc値が最小だった)ARIMA(3,0,1)(0,1,2)\(_{12}\)モデルからの予測を図9.26に示します。

h02 %>%
  model(ARIMA(log(Cost) ~ 0 + pdq(3,0,1) + PDQ(0,1,2))) %>%
  forecast() %>%
  autoplot(h02) +
  labs(y="オーストラリア・ドル (百万)",
       title="コルチコステロイド薬 (H02) の売上高", level = "区間予測")
H02月次売上高データに適用したARIMA(3,0,1)(0,1,2)$_{12}$モデルからの予測

図 9.26: H02月次売上高データに適用したARIMA(3,0,1)(0,1,2)\(_{12}\)モデルからの予測