10.2 fableでのARIMA誤差を持つ回帰

ARIMA()関数は、式に外生の予測変数がある場合、ARIMA誤差を持つ回帰モデルを適合させます。9.5節で紹介したように、pdq()特殊関数でARIMAモデルの次数を指定します。差分化を指定すると、モデル推計の前に回帰モデルの変数全てを差分化します。例えば、以下のコマンド

ARIMA(y ~ x + pdq(1,1,0))

は、モデル\(y_t' = \beta_1 x'_t + \eta'_t\)を適合します。ただし、\(\eta'_t = \phi_1 \eta'_{t-1} + \varepsilon_t\)はAR(1)誤差です。これは、以下のモデルと同値です。 \[ y_t = \beta_0 + \beta_1 x_t + \eta_t \] ただし、\(\eta_t\)はARIMA(1,1,0)誤差です。差分のモデルでは定数項が消えていますね。差分のモデルに定数を含めるには、モデル式に1を書き足します。

ARIMA()関数に、誤差の最良ARIMAモデルを自動選択させることもできます。pdq()特殊関数を指定しないと、自動選択します。最小2乗法を使って回帰モデルを推計し、残差にKPSS検定を行って、差分化が必要か否か決めています。差分化が必要なら、変数全てを差分化し、最尤法を使ってモデルを再推計します。たとえ差分化した変数で推計しても、最終モデルは元の変数で表します。

最終モデルではAICcが計算されるので、この値を使って最良の予測変数の組み合わせを決めることができます。つまり、考えるべき予測変数のサブセット全てに対してこの手順を繰り返して、AICc値が最小なモデルを選択するわけです。

事例: 米国の個人消費と可処分所得

10.1は、1970年から2019年第2四半期までの個人消費と可処分所得の四半期伸び率です。可処分所得の伸び率に基づいて個人消費の伸び率を予測したいとします。可処分所得の伸び率は、個人消費の伸び率に必ずしも瞬時に結び付くわけではありません。(例えば、失職後、新しい環境に合わせて支出を減らすのに2、3カ月かかるかもしれません。)しかし、ここではその複雑さは無視して、可処分所得の平均的伸び率の個人消費の平均的伸び率への瞬時の影響を測ることを試みます。

us_change %>%
  pivot_longer(c(Consumption, Income),
               names_to = "var", values_to = "value") %>%
  mutate(var = recode(var, "Consumption" = "個人消費", "Income" = "可処分所得")) %>% 
  ggplot(aes(x = Quarter, y = value)) +
  geom_line() +
  facet_grid(vars(var), scales = "free_y") +
  labs(title = "米国の個人消費と可処分所得",
       y = "四半期伸び率 (%)")
米国の個人消費と可処分所得の四半期伸び率(%)、1970 Q1から2019 Q2まで

図 10.1: 米国の個人消費と可処分所得の四半期伸び率(%)、1970 Q1から2019 Q2まで

fit <- us_change %>%
  model(ARIMA(Consumption ~ Income))
report(fit)
#> Series: Consumption 
#> Model: LM w/ ARIMA(1,0,2) errors 
#> 
#> Coefficients:
#>          ar1      ma1     ma2  Income  intercept
#>       0.7070  -0.6172  0.2066  0.1976     0.5949
#> s.e.  0.1068   0.1218  0.0741  0.0462     0.0850
#> 
#> sigma^2 estimated as 0.3113:  log likelihood=-163
#> AIC=338.1   AICc=338.5   BIC=357.8

(個人消費と可処分所得はそのままではなく、伸び率(%)で見ているので)データは既に明らかに定常です。ですから、差分化の必要はありません。当てはめたモデルは、以下の通りです。 \[\begin{align*} y_t &= 0.595 + 0.198 x_t + \eta_t \\ \eta_t &= 0.707 \eta_{t-1} + \varepsilon_t -0.617 \varepsilon_{t-1} + 0.207 \varepsilon_{t-2}\\ \varepsilon_t &\sim \text{NID}(0,0.311) \end{align*}\]

residuals()関数を使うと、\(\eta_t\)\(\varepsilon_t\)両系列の推計値を回収できます。

bind_rows(
    "回帰残差" =
        as_tibble(residuals(fit, type = "regression")),
    "ARIMA残差" =
        as_tibble(residuals(fit, type = "innovation")),
    .id = "type"
  ) %>%
  mutate(
    type = factor(type, levels=c(
      "回帰残差", "ARIMA残差"))
  ) %>%
  ggplot(aes(x = Quarter, y = .resid)) +
  geom_line() +
  facet_grid(vars(type))
当てはめたモデルからの回帰残差($\eta_t$)とARIMA残差($\varepsilon_t$)

図 10.2: 当てはめたモデルからの回帰残差(\(\eta_t\))とARIMA残差(\(\varepsilon_t\))

ホワイトノイズに似ているべきなのは、ARIMA誤差の推計値(イノベーション残差)です。

fit %>% gg_tsresiduals()
イノベーション残差(つまり、ARIMA誤差の推計値)は、ホワイトノイズと有意に違わない

図 10.3: イノベーション残差(つまり、ARIMA誤差の推計値)は、ホワイトノイズと有意に違わない

augment(fit) %>%
  features(.innov, ljung_box, dof = 5, lag = 8)
#> # A tibble: 1 × 3
#>   .model                      lb_stat lb_pvalue
#>   <chr>                         <dbl>     <dbl>
#> 1 ARIMA(Consumption ~ Income)    5.21     0.157