9.9 Modelli ARIMA stagionali

Finora, abbiamo limitato la nostra attenzione a serie di dati non stagionali e a modelli ARIMA non stagionali. Tuttavia, i modelli ARIMA sono anche in grado di modellare componenti stagionali presenti nei dati.

Un modello ARIMA stagionale è formato includendo termini stagionali aggiuntivi nei modelli ARIMA che abbiamo visto finora. Può scriversi come segue:

ARIMA \(\underbrace{(p, d, q)}\) \(\underbrace{(P, D, Q)_{m}}\)
\({\uparrow}\) \({\uparrow}\)
Non-seasonal part Seasonal part
of the model of the model

dove \(m =\) è il periodo stagionale (cioè, il numero di osservazioni per anno). Nel seguito, useremo la notazione maiuscola per indicare le componenti stagionali del modello, e la notazione minuscola per le componenti non stagionali.

La parte stagionale del modello consiste in termini che sono simili alle componenti non stagionali del modello, ma coinvolgono ritardi del periodo stagionale. Per esempio, un modello ARIMA(1,1,1)(1,1,1)\(_{4}\) (senza costante) è un modello per dati trimestrali (\(m=4\)), e può essere scritto come \[ (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}. \]

I temini stagionali aggiuntivi sono semplicemente moltiplicati per i termini non stagionali.

ACF/PACF

La componente stagionale di un modello AR o MA si può vedere ai ritardi stagionali della PACF e ACF. Per esempio, un modello ARIMA(0,0,0)(0,0,1)\(_{12}\) avrà:

  • un picco al ritardo 12 della ACF e nessun altro picco significativo;
  • un decrescimento esponenziale ai ritardi stagionali della PACF (cioè, ai ritardi 12, 24, 36, …).

In modo simile, un modello ARIMA(0,0,0)(1,0,0)\(_{12}\) mostrerà:

  • un decrescimento esponenziale ai ritardi stagionali della ACF;
  • un singolo picco significativo al ritardo 12 della PACF.

Per identificare gli ordini stagionali appropriati per un modello ARIMA stagionale, è bene limitare l’attenzione ai ritardi stagionali.

La procedura di modellazione è quasi la stessa dei dati non stagionali, se non che dobbiamo selezionare le componenti AR e MA stagionali oltre alle componenti non stagionali del modello. Il processo è meglio illustrato tramite alcuni esempi.

Esempio: Occupazione mensile negli Stati Uniti per il tempo libero e l’ospitalità

Descriveremo, ora, come identificare un modello ARIMA stagionale utilizzando i dati mensili dell’occupazione statunitense per i posti di lavoro nel settore del tempo libero e dell’ospitalità da gennaio 2001 a settembre 2019, mostrati nella figura 9.18.

leisure <- us_employment %>%
  filter(Title == "Leisure and Hospitality",
         year(Month) > 2000) %>%
  mutate(Employed = Employed/1000) %>%
  select(Month, Employed)
autoplot(leisure, Employed) +
  labs(title = "Occupazione mensile USA: tempo libero e ospitalità",
       y="Numero di persone (milioni)")
Occupazione mensile USA per tempo libero e ospitalità, 2001-2019.

Figura 9.18: Occupazione mensile USA per tempo libero e ospitalità, 2001-2019.

I dati sono chiaramente non stazionari, con una forte componente stagionale e un trend non lineare, pertanto dobbiamo utilizzare una differenza stagionale. La differenza stagionale della serie è mostrata in figura 9.19.

leisure %>%
  gg_tsdisplay(difference(Employed, 12),
               plot_type='partial', lag=36) +
  labs(title="Differenza stagionale", y="")
Occupazione mensile USA per tempo libero e ospitalità: differenza stagionale.

Figura 9.19: Occupazione mensile USA per tempo libero e ospitalità: differenza stagionale.

Appare evidente che ci siano altre componenti non stazionarie, pertanto sulla serie (già differenziata stagionalmente) viene applicata una differenza prima. La serie ulteriormente differenziata è riportata in figura 9.20.

leisure %>%
  gg_tsdisplay(difference(Employed, 12) %>% difference(),
               plot_type='partial', lag=36) +
  labs(title = "Doppia differenziazione", y="")
Occupazione mensile USA per tempo libero e ospitalità: doppia differenziazione.

Figura 9.20: Occupazione mensile USA per tempo libero e ospitalità: doppia differenziazione.

In nostro obiettivo, ora, è quello di trovare un modello ARIMA appropriato basandosi sulle funzioni ACFe PACF riportate in figura 9.20. In picco significativo al ritardo 2 nella ACF suggerisce una componente MA(2) non stagionale. Il picco significativo al ritardo 12 nella ACF suggerisce una componente MA(1) stagionale. Conseguentemente, iniziamo con un modello ARIMA(0,1,2)(0,1,1)\(_{12}\), che indica la presenza di una differenza prima, una differenza stagionale, una componente MA(2) non stagionale e una MA(1) stagionale. Se avessimo iniziato, opsservando la PACF, avremmo selezionato un modello ARIMA(2,1,0)(0,1,1)\(_{12}\) — usando la PACF per identificare la componente non stagionale e l’ACF per identificare la componente stagionale. Abbiamo anche consideranto un modello identificato in modo automatico. Ponendo stepwise=FALSE e approximation=FALSE, richiediamo del lavoro extra ad R per trovare un buon modello. Ovviamente questo richiede più tempo, ma avendo una sola serie, ciò non costituisce un grande problema.

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.

La funzione ARIMA() usa unitroot_nsdiffs() per determinare il valore di \(D\) (il numero di differenze stagionali da utilizzare), e unitroot_ndiffs() per determinare il valore di \(d\) (il numero di differenze non stagionali), quando tali valori non sono specificati. La selezione degli altri ordini del modello (\(p,q,P\) e \(Q\)) è determinata minimizzando il valore del criterio di informazione automatico AICc, come abbiamo visto per i modelli non stagionali.

I tre modelli identificati hanno valori molto simili del criterio AICc, con il modello identificato in modo automatico leggermente migliore. La nostra seconda “scelta”, il modello ARIMA(2,1,0)(0,1,1)\(_{12}\), è risultato essere molto vicino al modello selezionato automaticamente con ARIMA(2,1,0)(1,1,1)\(_{12}\).

I grafici relativi ai residui per il modello migliore sono riportati nella figura 9.21.

fit %>% select(auto) %>% gg_tsresiduals(lag=36)
Residui del modello ARIMA(2,1,0)(1,1,1)\(_{12}\).

Figura 9.21: Residui del modello ARIMA(2,1,0)(1,1,1)\(_{12}\).

Sebbene ci sia un piccolo, ma significativo, picco (al ritardo 11) sui 36 considerati l’ipotesi che i residui siano white noise può considerarsi verificata. Per esserne certi, usiamo il test di Ljung-Box, facendo attenzione a impostare i gradi di libertà in modo che corrispondano al numero di parametri nel modello.

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

I grandi valori dei p-value confermano che i residui si comportano come white noise.

Abbiamo quindi identificato un modello ARIMA stagionale che supera i controlli richiesti ed è pronto per le previsioni. Le previsioni del modello per i successivi tre anni sono riportate nella figura 9.22. Grazie alla doppia differenziazione, le previsioni seguono molto bene il trend e l’andamento stagionale.

forecast(fit, h=36) %>%
  filter(.model=='auto') %>%
  autoplot(leisure) +
  labs(title = "Occupazione USA: tempo libero e ospitalità",
       y="Numero di persone (milioni)")
Previsioni per l’occupazione mensile negli Stati Uniti per il tempo libero e l’ospitalità usando il modello ARIMA(2,1,0)(1,1,1)\(_{12}\). Sono riportati gli intervalli di previsione all’80% e al 95%.

Figura 9.22: Previsioni per l’occupazione mensile negli Stati Uniti per il tempo libero e l’ospitalità usando il modello ARIMA(2,1,0)(1,1,1)\(_{12}\). Sono riportati gli intervalli di previsione all’80% e al 95%.

Esempio: vendita di corticosteroidi in Australia

Come secondo esempio, cerchiamo di prevedere le vendite mensili di farmaci corticosteroidi in Australia. Questi farmaci sono anche conosciuti come farmaci H02 (secondo lo schema di classificazione Anatomico Terapeutico Chimico).

h02 <- PBS %>%
  filter(ATC2 == "H02") %>%
  summarise(Cost = sum(Cost)/1e6)
h02 %>%
  mutate(log(Cost)) %>%
  pivot_longer(-Month) %>%
  ggplot(aes(x = Month, y = value)) +
  geom_line() +
  facet_grid(name ~ ., scales = "free_y") +
  labs(y="", title="Vendite di corticosteroidi (H02)")
Vendite di corticosteroidi in Australia (in milioni di script per mese). Il grafico inferiore riporta il logaritmo dei dati..

Figura 9.23: Vendite di corticosteroidi in Australia (in milioni di script per mese). Il grafico inferiore riporta il logaritmo dei dati..

I dati, dal luglio 1991 al giugno 2008, sono riportati nella figura 9.23. Poiché è evidente una leggera crescita della variabilità, prendiamo il logartitmo dei dati per stabilizzare la varianza.

La serie contiene una forte componente stagionale ed è, ovviamente, non stazionaria, pertanto è necessario usare una differenza stagionale. La serie differenziata stagionalmente è riportata nella figura 9.24. Non è chiaro, a questo punto, se dovrebbe essere utilizzata anche una differenza prima o no. Noi decidiamo di non usarla, ma la scelta non è evidente.

Le ultime osservazioni sembrano essere diverse (più variabili) dai dati precedenti. Questo può essere dovuto al fatto che i dati sono talvolta revisionati quando le vendite precedenti sono riportate in ritardo.

h02 %>% gg_tsdisplay(difference(log(Cost), 12),
                     plot_type='partial', lag_max = 24)
Serie stagionalmente differenziata delle vendite di corticosteriodi in Australia (milioni di script per mese).

Figura 9.24: Serie stagionalmente differenziata delle vendite di corticosteriodi in Australia (milioni di script per mese).

Nei grafici relativi alla serie stagionalmente differenziata, sono evidenti due picchi significativi nella PACF ai ritardi 12 e 24, ma nella ACF non ci sono picchi significativi ai ritardi stagionali. Questo suggerisce una componente AR(2) stagionale. Ai ritardi non stagionali, troviamo tre picchi significativi nella PACF, che suggeriscono una possibile componente AR(3). La funzione ACF non suggerisce alcun semplice modello.

Conseguentemente, questa analisi iniziale suggerisce che un possibile modello per i nostri dati sia un ARIMA(3,0,0)(2,1,0)\(_{12}\). Adattiamo, quindi, questo modello ai dati assieme ad altri modelli che variano poco rispetto a questo, e ne calcoliamo i valori del criterio AICc che sono riportati nella Tabella 9.2.

Tabella 9.2: Valori del criterio AICc per i diversi modelli ARIMA adattati alla serie delle vendite mensili H02.
Model 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

Di questi modelli, il migliore (cioè quello col minor valore del criterio AICc) è il modello ARIMA(3,0,1)(0,1,2)\(_{12}\). I residui di tale modello sono riportati nella figura 9.25.

fit <- h02 %>%
  model(ARIMA(log(Cost) ~ 0 + pdq(3,0,1) + PDQ(0,1,2)))
fit %>% gg_tsresiduals(lag_max=36)
Residui del modello ARIMA(3,0,1)(0,1,2)$_{12}$ adattato alla serie delle vendite mensili di corticosteroidi.

Figura 9.25: Residui del modello ARIMA(3,0,1)(0,1,2)\(_{12}\) adattato alla serie delle vendite mensili di corticosteroidi.

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

Ci sono alcuni picchi significativi nella ACF e i residui non superano il test di Ljung-Box. Il modello può ancora essere usato per prevedere, ma gli intervalli di previsione potrebbero non essere accurati a causa dei residui correlati.

Proviamo ad usare l’algoritmo ARIMA automatico. L’esecuzione di ARIMA() lasciando tutti gli argomenti ai loro valori di default ha selezionato un modello ARIMA(2,1,0)(0,1,1)\(_{12}\). Eseguendo ARIMA() con stepwise=FALSE e approximation=FALSE si ottiene un modello ARIMA(2,1,3)(0,1,1)\(_{12}\). Tuttavia, entrambi i modelli non superano il test di Ljung-Box per 36 ritardi. A volte non è proprio possibile trovare un modello che superi tutti i test.

Valutazione su un insieme di prova:

Confrontiamo ora alcuni dei modelli adattati finora, usando un sottoinsieme di prova set costituito dai dati degli ultimi due anni. In pratica, adattiamo i modelli considerati usando solo i dati da luglio 1991 a giugno 2006, e prevediamo, quindi,le vendite degli script per il periodo luglio 2006 – giugno 2008. I risultati sono riassunti nella tabella 9.3.

Tabella 9.3: Valori del RMSE per diversi modelli ARIMA models applicati alla serie mensile delle vendite H02 nel periodo Luglio 2006 – Giugno 2008 .
.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

Sull’insieme di dati di prova, i modelli scelti manualmente hanno un RMSE molto prossimo al valore del miglior modello, ma anche i modelli scelti automaticamente con ARIMA() non presentano valori molto distanti.

Quando i modelli vengono confrontati usando i valori del criterio AICc, è importante che abbiano tutti lo stesso ordine di differenziazione. Tuttavia, quando confrontiamo più modelli usando un sottoinsieme di dati di prova, when comparing models using a test set, qualsiasi sia il modo in cui sono state calcolate le previsioni — i confronti sono sempre validi. Di conseguenza, nella tabella qui sopra, possiamo includere alcuni modelli con la sola differenza stagionale e altri modelli che contengono sia la differenza prima sia quella stagionale, mentre nella tabella precedente contenente i valori AICc, abbiamo confrontato solo i modelli con la differenza stagionale, ma senza la differenza prima.

Nessuno dei modelli qui considerati supera tutti i test dei residui. In pratica, useremo il miglior modello che possiamo trovare, anche se non supera tutti i test.

Le previsioni dal modello ARIMA(3,0,1)(0,1,2)\(_{12}\) (che ha il secondo valore più basso del RMSE nell’insieme di prova, e il miglior valore del criterio AICc fra i modelli con solo la differenza stagionale) sono mostrate nella figura 9.26.

h02 %>%
  model(ARIMA(log(Cost) ~ 0 + pdq(3,0,1) + PDQ(0,1,2))) %>%
  forecast() %>%
  autoplot(h02) +
  labs(y=" $AU (millions)",
       title="Vendite di corticosteroidi (H02)")
Previsioni dal modello ARIMA(3,0,1)(0,1,2)$_{12}$ applicato alla serie mensile delle vendite di farmaci H02.

Figura 9.26: Previsioni dal modello ARIMA(3,0,1)(0,1,2)\(_{12}\) applicato alla serie mensile delle vendite di farmaci H02.