10.2 Regressione con errori di tipo ARMA usando fable

La funzione ARIMA() adatterà un modello di regressione con errori ARIMA se i regressori esogeni sono inclusi nella formula. Come introdotto nel paragrafo 9.5, la funzione speciale pdq() specifica l’ordine del modello per il termine di errore ARIMA. Se si specifica la differenziazione, questa è applicata a tutte le variabili nel modello di regressione prima che il modello sia stimato. Per esempio, il comando

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

stima il modello \(y_t' = \beta_1 x'_t + \eta'_t\), dove \(\eta'_t = \phi_1 \eta'_{t-1} + \varepsilon_t\) è un errore AR(1). Questo è equivalente al modello \[ y_t = \beta_0 + \beta_1 x_t + \eta_t, \] dove \(\eta_t\) è un errore ARIMA(1,1,0). Si noti che il termine costante scompare a causa della differenziazione. Per includere una costante nel modello differenziato, si deve aggiungere 1 alla formula del modello.

La funzione ARIMA() può anche essere usata per selezionare il miglior modello ARIMA per gli errori. Questo viene fatto non specificando la funzione speciale pdq(). La necessità o meno di differenziare viene determinata applicando un test KPSS ai residui del modello di regressione stimato utilizzando i minimi quadrati ordinari. Se è richiesta la differenziazione, allora tutte le variabili vengono differenziate ed il modello viene nuovamente stimato usando la stima di massima verosimiglianza. Il modello finale sarà espresso in termini di variabili originarie, anche se è stato stimato utilizzando variabili differenziate.

L’indice AICc è calcolato per il modello finale, e questo valore può essere usato per determinare i migliori predittori. Questo significa che la procedura dovrebbe essere ripetuta per tutti i sottoinsiemi di predittori da considerare, e sarà poi selezionato il modello con il più basso valore AICc.

Esempio: consumi individuali e reddito negli Stati Uniti

La figura 10.1 mostra le variazioni trimestrali della spesa per consumi personali e del reddito personale disponibile dal primo trimestre (Q1) del 1970 al secondo trimestre (Q2) del 2019. Si è interessati a prevedere i cambiamenti nella spesa a partire dai cambiamenti nel reddito. Un cambiamento nel reddito non si traduce necessariamente in un cambiamento istantaneo nel consumo (ad esempio, dopo la perdita del lavoro, ad una persona potrebbero essere necessari alcuni mesi per ridurre le spese per tenere conto delle nuove circostanze). Tuttavia, in questo esempio si ignorerà questa complessità e si cercherà di misurare l’effetto istantaneo del cambiamento medio del reddito sul cambiamento medio delle spese di consumo.

us_change %>%
  pivot_longer(c(Consumption, Income),
               names_to = "var", values_to = "value") %>%
  ggplot(aes(x = Quarter, y = value)) +
  geom_line() +
  facet_grid(vars(var), scales = "free_y") +
  labs(title = "Consumi e reddito personale negli Stati Uniti",
       x = "Trimestre", y = "Variazioni % trimestrali")
Variazioni percentuali nelle spese personali per consumo e nel reddito personale disponibile negli Stati Uniti, dal primo trimestre (Q1) del 1970 al secondo trimestre (Q2) del 2019.

Figura 10.1: Variazioni percentuali nelle spese personali per consumo e nel reddito personale disponibile negli Stati Uniti, dal primo trimestre (Q1) del 1970 al secondo trimestre (Q2) del 2019.

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

I dati sono chiaramente già stazionari (tenuto conto che si stanno considerando le variazioni percentuali piuttosto che la spesa e il reddito grezzi), quindi non c’è bisogno di alcuna differenziazione. Il modello stimato è \[\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*}\]

È possibile recuperare le stime di entrambe le serie \(\eta_t\) e \(\varepsilon_t\) usando la funzione residuals().

bind_rows(
    `Regression residuals` =
        as_tibble(residuals(fit, type = "regression")),
    `ARIMA residuals` =
        as_tibble(residuals(fit, type = "innovation")),
    .id = "type"
  ) %>%
  mutate(
    type = factor(type, levels=c(
      "Regression residuals", "ARIMA residuals"))
  ) %>%
  ggplot(aes(x = Quarter, y = .resid)) +
  geom_line() +
  labs(x = "Trimestre", y = "Residui") +
  facet_grid(vars(type))
Residui della regressione ($\eta_t$) -- pannello in alto, e residui del modello ARIMA ($\varepsilon_t$) -- pannello in basso, sul modello stimato.

Figura 10.2: Residui della regressione (\(\eta_t\)) – pannello in alto, e residui del modello ARIMA (\(\varepsilon_t\)) – pannello in basso, sul modello stimato.

In questo caso sono gli errori stimati ARIMA (i residui dell’innovazione) che dovrebbero assomigliare a una serie di rumore bianco.

fit %>% gg_tsresiduals()
I residui di innovazione (i.e., gli errori ARIMA stimati) non sono significativamente differenti da un white noise.

Figura 10.3: I residui di innovazione (i.e., gli errori ARIMA stimati) non sono significativamente differenti da un white noise.

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