9.5 Modelli ARIMA non stagionali

Se combiniamo la differenziazione con modelli autoregressivi e/o a media mobile, otteniamo il modello ARIMA non stagionale. ARIMA è un acronimo per AutoRegressive Integrated Moving Average (in questo contesto, “integrazione” è l’operazione inversa della differenziazione). Il modello con tutte le componenti (autoregressiva, integrata e a media mobile) può essere scritto come \[\begin{equation} y'_{t} = c + \phi_{1}y'_{t-1} + \cdots + \phi_{p}y'_{t-p} + \theta_{1}\varepsilon_{t-1} + \cdots + \theta_{q}\varepsilon_{t-q} + \varepsilon_{t}, \tag{9.1} \end{equation}\] dove \(y'_{t}\) è la serie differenziata (che può essere differenziata anche più di una volta). I “predittori” sulla parte destra includono sia valori ritardati di \(y_t\) sia errori ritardati. Denotiamo questo modello con ARIMA(\(p, d, q\)), dove

\(p=\) ordine della componente autoregressiva;
\(d=\) numero di differenziazioni;
\(q=\) ordine della componente a media mobile.

Le stesse condizioni di stazionarietà e invertibilità che si applicano a modelli puramente autoregressivi e puramente a media mobile si applicano anche ad un modello ARIMA.

Molti dei modelli che abbiamo già discusso sono casi particolari del modello ARIMA, come mostrato nella Tabella 9.1.

Tabella 9.1: Casi speciali di modelli ARIMA.
White noise ARIMA(0,0,0) senza costante
Random walk ARIMA(0,1,0) senza costante
Random walk con drift ARIMA(0,1,0) con costante
Autoregressivo ARIMA(\(p\),0,0)
Media mobile ARIMA(0,0,\(q\))

Nel momento in cui combiniamo diverse componenti per ottenere modelli più complicati, è molto più conveniente utilizzare la notazione con l’operatore ritardo. Per esempio, l’equazione (9.1) può essere scritta usando l’operatore ritardo come \[\begin{equation} \tag{9.2} \begin{array}{c c c c} (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\\ {\uparrow} & {\uparrow} & &{\uparrow}\\ \text{AR($p$)} & \text{$d$ differences} & & \text{MA($q$)}\\ \end{array} \end{equation}\]

La selezione di valori appropriati per gli ordini \(p\), \(d\) e \(q\) può essere difficile. Tuttavia, la funzione ARIMA() del pacchetto fable lo farà per noi in modo automatico. Nella Sezione 9.7, impareremo come lavora questa funzione, insieme ad alcuni metodi per scegliere da soli questi valori.

Esempio: Esportazioni egiziane

La figura 9.7 mostra le esportazioni annuali egiziane come una percentuale del GDP nel periodo 1960-2017.

global_economy %>%
  filter(Code == "EGY") %>%
  autoplot(Exports) +
  labs(y = "% del GDP", title = "Esportazioni egiziane")
Esportazioni annuali egiziane come percentuale del GDP, anni 1960-2017.

Figura 9.7: Esportazioni annuali egiziane come percentuale del GDP, anni 1960-2017.

Il seguente codice R identifica, in modo automatico, un modello ARIMA non stagionale.

fit <- global_economy %>%
  filter(Code == "EGY") %>%
  model(ARIMA(Exports))
report(fit)
#> Series: Exports 
#> Model: ARIMA(2,0,1) w/ mean 
#> 
#> Coefficients:
#>          ar1      ar2      ma1  constant
#>       1.6764  -0.8034  -0.6896    2.5623
#> s.e.  0.1111   0.0928   0.1492    0.1161
#> 
#> sigma^2 estimated as 8.046:  log likelihood=-141.6
#> AIC=293.1   AICc=294.3   BIC=303.4

Tale modello è un modello ARIMA(2,0,1): \[ y_t = 2.56 + 1.68 y_{t-1} -0.80 y_{t-2} -0.69 \varepsilon_{t-1} + \varepsilon_{t}, \] dove \(\varepsilon_t\) è white noise con una deviazione standard pari a \(2.837 = \sqrt{8.046}\). Le previsioni ottenute utilizzando tale modello sono mostrate nella figura 9.8. Si osservi come negli ultimi decenni la serie presenti dei cicli molto evidenti.

fit %>% forecast(h=10) %>%
  autoplot(global_economy) +
  labs(y = "% of GDP", title = "Esportazioni egiziane")
Previsioni per le esportazioni egiziane.

Figura 9.8: Previsioni per le esportazioni egiziane.

Comprendere i modelli ARIMA

La funzione ARIMA() è utile, ma, come tutte le procedure automatiche, può essere un po’ fuorviante, e vale la pena cercare di capire qualcosa sul comportamento dei modelli, anche quando ci si affida a una procedura automatica che sceglie il modello per noi

La costante \(c\) ha un effetto molto importante sulle previsioni di lungo periodo ottenute da questo tipo di modelli.

  • Se \(c=0\) e \(d=0\), le previsioni di lungo periodo tenderanno a zero.
  • Se \(c=0\) e \(d=1\), le previsioni di lungo periodo tenderanno a una costante diversa da zero.
  • Se \(c=0\) e \(d=2\), le previsioni di lungo periodo seguiranno una retta.
  • Se \(c\ne0\) e \(d=0\), le previsioni di lungo periodo tenderanno alla media dei dati.
  • Se \(c\ne0\) e \(d=1\), le previsioni di lungo periodo seguiranno una retta.
  • Se \(c\ne0\) e \(d=2\), le previsioni di lungo periodo seguiranno un trend quadratico. (Tuttavia, ciò non è raccomandato e, infatti, fable non lo permette.)

Il valore di \(d\) influisce anche sugli intervalli di previsione — più alto è il valore di \(d\), più rapidamente gli intervalli di previsione aumentano di dimensione. Per \(d=0\), la deviazione standard della previsione a lungo termine sarà pari alla deviazione standard dei dati storici, quindi gli intervalli di previsione avranno, essenzialmente, tutti la stessa ampiezza, indipendentemente dall’orizzonte di previsione.

Questo comportamento si vede nella figura 9.8 dove \(d=0\) e \(c\ne0\). In questa figura, gli intervalli di previsione hanno quasi la stessa ampiezza per gli ultimi orizzonti di previsione, e le previsioni per gli ultimi punti sono vicine alla media dei dati.

Il valore di \(p\) è importante se i dati mostrano dei cicli. Per ottenere previsioni cicliche, è necessario avere \(p\ge2\), insieme ad alcune condizioni aggiuntive sui parametri. Per un modello AR(2), il comportamento ciclico si verifica se \(\phi_1^2+4\phi_2<0\) (come nel caso del modello delle esportazioni egiziane). In questo caso, il periodo medio dei cicli è16

\[ \frac{2\pi}{\text{arc cos}(-\phi_1(1-\phi_2)/(4\phi_2))}. \]

Grafici delle funzioni ACF e PACF

Di solito non è possibile dire, semplicemente osservando il grafico di una serie storica, quali valori di \(p\) e \(q\) siano appropriati per i dati. Tuttavia, a volte è possibile usare il grafico ACF, e il grafico PACF strettamente correlato, per determinare i valori appropriati di \(p\) e \(q\).

Ricordiamo che il grafico ACF, mostra le autocorrelazioni empiriche, che misurano la relazione tra \(y_t\) e \(y_{t-k}\) per diversi valori di \(k\). Se \(y_t\) e \(y_{t-1}\) sono correlati, allora anche \(y_{t-1}\) e \(y_{t-2}\) devono essere correlati. Tuttavia, \(y_t\) e \(y_{t-2}\) potrebbero essere correlati, semplicemente perché sono entrambi correlati a \(y_{t-1}\), anziché a causa di qualche nuova informazione contenuta in \(y_{t-2}\) che potrebbe essere usata per prevedere \(y_t\).

Per superare questo problema, possiamo usare le autocorrelazioni parziali. Queste misurano la relazione tra \(y_{t}\) e \(y_{t-k}\) dopo aver rimosso gli effetti dei ritardi \(1, 2, 3, \ldots, k - 1\). Quindi la prima autocorrelazione parziale è identica alla prima autocorrelazione, perché non c’è nulla da rimuovere. Ogni autocorrelazione parziale può essere stimata come l’ultimo coefficiente di un modello autoregressivo. In particolare, \(\alpha_k\), il \(k\)-esimo coefficiente di autocorrelazione parziale, è uguale alla stima di \(\phi_k\) in un modello AR(\(k\)). In pratica, esistono algoritmi più efficienti per calcolare \(\alpha_k\), che, tuttavia, danno gli stessi risultati.

Le figure 9.9 e 9.10 riportano i grafici delle autocorrelazioni empiriche globali (ACF) e parziali (PACF) per la serie delle esportazioni mostrata in figura 9.7. Le autocorrelazioni parziali hanno gli stessi valori critici, \(\pm 1.96/\sqrt{T}\), delle autocorrelazioni globali, e tali valori sono, usualmente, riportati nei grafici delle autocorrelazioni empiriche, come in figura 9.10.

global_economy %>%
  filter(Code == "EGY") %>%
  ACF(Exports) %>%
  autoplot()
ACF delle esportazioni egiziane.

Figura 9.9: ACF delle esportazioni egiziane.

global_economy %>%
  filter(Code == "EGY") %>%
  PACF(Exports) %>%
  autoplot()
PACF delle esportazioni egiziane.

Figura 9.10: PACF delle esportazioni egiziane.

Un modo conveniente per riprodurre allo stesso tempo il grafico della serie storica assieme alle sue funzioni di autocorrelazioni empiriche (ACF e PACF), è usare la funzione gg_tsdisplay() con plot_type = "partial".

Se i dati sono realizzazioni di un modello ARIMA(\(p\),\(d\),0) o ARIMA(0,\(d\),\(q\)), allora i grafici ACF e PACF possono essere molto utili nel determinare i valori di \(p\) or \(q\). Se, invece, \(p\) e \(q\) sono entrambi positivi, allora i grafici ACF e PACF non sono di molto aiuto nel determinare gli opportuni valori degli ordini \(p\) e \(q\).

La serie potrebbe essere modellata da un modello ARIMA(\(p\),\(d\),0), se i grafici delle funzioni ACF e PACF della serie, eventualmente differenziata, mostrano i seguenti comportamenti:

  • la funzione di autocorrelazione globale (ACF) decresce verso lo zero in modo esponenziale o sinusoidale;
  • nella PACF, l’autocorrelazione a ritardo \(p\) è significativa, ma tutte le altre autocorrelazioni oltre il ritardo \(p\) non lo sono.

Viceversa, la serie potrebbe essere modellata da un ARIMA(0,\(d\),\(q\)) se le funzioni di autocorrelazione empiriche ACF e PACF della serie, eventualmente differenziata, mostrano i seguenti comportamenti:

  • la funzione di autocorrelazione parziale (PACF) decresce verso lo zero in modo esponenziale o sinusoidale;
  • nella ACF, l’autocorrelazione a ritardo \(q\) è significativa, ma tutte le altre autocorrelazioni oltre il ritardo \(p\) non lo sono.

Dalla figura 9.9, possiamo osservare che la funzione ACF decresce verso lo zero seguento un comportamento sinusoidale, mentre dalla figura 9.10 vediamo che la PACF ha l’ultimo valore significativo al ritardo 4. Questo comportamento è coerente con un modello ARIMA(4,0,0).

fit2 <- global_economy %>%
  filter(Code == "EGY") %>%
  model(ARIMA(Exports ~ pdq(4,0,0)))
report(fit2)
#> Series: Exports 
#> Model: ARIMA(4,0,0) w/ mean 
#> 
#> Coefficients:
#>          ar1      ar2     ar3      ar4  constant
#>       0.9861  -0.1715  0.1807  -0.3283    6.6922
#> s.e.  0.1247   0.1865  0.1865   0.1273    0.3562
#> 
#> sigma^2 estimated as 7.885:  log likelihood=-140.5
#> AIC=293.1   AICc=294.7   BIC=305.4

Questo modello è solo leggermente peggiore del modello ARIMA(2,0,1) identificato dalla funzione ARIMA() (con un valore AICc di 294.70 confrontato a 294.29).

Possiamo anche specificare particolari valori di pdq() e quindi utilizzare la funzione ARIMA(). Per esempio, per trovare il miglior modello con \(p\in\{1,2,3\}\), \(q\in\{0,1,2\}\) e \(d=1\), si può usare ARIMA(y ~ pdq(p=1:3, d=1, q=0:2)).


  1. arc cos è la funzione inversa del coseno. Dovreste essere in grado di trovarla sulla vostra calcolatrice. Può essere chiamata acos o cos\(^{-1}\).↩︎