7.4 Modelli ARIMA: proprietà
In questa sezione discutiamo tre proprietà fondamentali dei modelli ARIMA, ottenendo condizioni sulla stazionarietà, una equazione ricorsiva per la funzione di autocovarianza (nel caso stazionario) e infine accennando al problema della stima dei parametri sulla base delle osservazioni, che include anche il problema della selezione del modello, ossia la scelta degli ordini \((p,d,q)\).
Consideriamo quindi un processo \((X_t)_{t \in \mathcal{T}}\) \(\operatorname{ARIMA}(p,d,q)\) di parametri \((\alpha_i)_{i=1}^p\) e \((\beta_j)_{j=1}^q\).
7.4.1 Stazionarietà
Il problema della stazionarietà è stato discusso nel caso \(\operatorname{ARIMA}(1,0,0)\), l’equazione lineare con smorzamento, in cui era stata ottenuta la condizione necessaria \(|\alpha_1| <1\) (e sufficiente, purché \(X_0\) fosse gaussiano di varianza opportuna). L’esempio della passeggiata aleatoria, pensato come \(\operatorname{ARIMA}(0,1,0)\) mostra che in tal caso la stazionarietà non vale.
Prima di presentare il risultato generale, osserviamo che i processi a media mobile, ossia \(\operatorname{ARIMA}(0,0,q)\) possono sempre essere stazionari (se si definiscono \(X_0\), \(X_1\), …, \(X_{q-1}\) opportunamente). Infatti l’equazione che li definisce, \[ X_t = q(L) W_t = W_t + \beta_1 W_{t-1} + \beta_2 W_{t-2} + \ldots + \beta_q W_{t-q}\] se estesa anche per \(t = 0\), \(1\), , \(q-1\) considerando il rumore bianco gaussiano definito anche per tempi negativi, è un caso particolare di convoluzione di un processo stazionario (il rumore bianco gaussiano \(W\)) con un filtro (dato dai coefficienti \(\beta_j\)), e quindi abbiamo già osservato che preserva la stazionarietà.
Per comprendere la stazionarietà nel caso generale, l’idea formale è di “risolvere” l’equazione del modello \[ p(L)(1-L)^d X_t = q(L) W_t,\] dividendo formalmente per \(p(L)(1-L)^d\). Si ottiene \[ X_ t =\frac{ q(L)}{ p(L)(1-L)^d} W_t,\] una scrittura che però non ha molto senso (sappiamo definire solo i polinomi nell’operatore ritardo \(L\), non certo le funzioni razionali). Tuttavia, sviluppando la funzione come serie di Taylor \[ \frac{ q(z)}{p(z)(1-z)^d} = \sum_{k=0}^\infty b_k z^k,\] possiamo almeno tentare di definire \(X\) nel seguente modo, \[ X_t = \sum_{k=0}^\infty b_k L^k W_t,\] avendo definito \(W_t\) anche per \(t\) negativi, che è una sorta di modello a media mobile con \(q=\infty\). La stazionarietà sarebbe allora un caso limite di quanto osservato prima, per \(q\) finito. Ovviamente tutto il problema sta nel mostrare che la serie effettivamente converge, fatto che dipende dalla crescita dei coefficienti \(b_k\) al tendere di \(k \to \infty\) e in ultima analisi agli zeri (complessi) del denominatore \(p(z)(1-z)^d\). Il risultato preciso che si può dimostrare è il seguente.
Teorema 7.1 Dati \((p,d,q)\) e coefficienti \((\alpha_i)_{i=1}^p\), \((\beta_j)_{j=1}^q\), posto \[ p(z) = 1- \sum_{i=1}^p \alpha_i z^i,\] allora esiste un modello \(\operatorname{ARIMA}(p,d,q)\) stazionario con tali coefficienti se \(d=0\) e tutte le radici complesse di \(p(z)\) hanno modulo \(|z|>1\), ossia \[ \text{se $z \in \mathbb{C}$ è tale che $p(z) = 0$, allora $|z|>1$.}\]
Verifichiamo che il teorema recupera la condizione trovata per l’equazione lineare con smorzamento. In tal caso vale \[ p(z) =1 -\alpha_1 z,\] la cui unica radice è \(z = 1/\alpha_1\). Essa ha modulo maggiore di \(1\) se e solo se \(|\alpha_1|<1\), che è appunto la condizione trovata.
7.4.2 Autocovarianza
Per costruzione i processi ARIMA hanno media nulla (nel caso fosse rilevante ammettere una media \(m\) non nulla basta modellizzare la differenza \(X_t - m\) come un ARIMA). L’equazione permette anche di ottenere una formula ricorsiva per la funzione di autocovarianza.
Vale infatti (supponiamo \(d=0\) per semplicità), per \(t\ge \max\cur{p,q}\), \[ \E{ X_s p(L) X_t} = \E{ X_s q(L) W_t} = \E{ X_s (W_t+\sum_{j=1}^q \beta_j W_{t-j})}\] Possiamo supporre che \(X_s\) sia indipendente dal rumore bianco \(W_r\), purché \(r>s\). In particolare, se \(t-q >s\), il membro a destra contiene solamente termini nulli, perché del tipo \[ \E{ X_s W_r}\] con \(r>s\). Ne segue che \[ 0 = \E{ X_s p(L) X_t } = \E{ X_s (X_t - \sum_{i=1}^p \alpha_i X_{t-i} )} = C(s,t)-\sum_{i=1}^p \alpha_i C(s,t-i).\] Riorganizzando i termini, troviamo che \[ C(s,t) = \sum_{i=1}^p \alpha_i C(s, t-i),\] purché \(t > s +q\). In particolare, \[ C(0,t) = \sum_{i=1}^p \alpha_i C(0, t-i), \quad \text{se $t>q$.}\] che è particolarmente utile nel caso in cui \(X\) sia stazionario.
Queste formule ricorsive sono dette equazioni di Yule-Walker e permettono di ricavare la funzione di autocovarianza per intervalli temporali (lag) grandi. In particolare, notiamo che nel caso \(p=1\), si riduce alla relazione già trovata \[ C(0,t) = \alpha_1 C(0,t-1)\] che iterando porta a \[ C(0,t) = \alpha_1^t C(0,0).\] Recuperiamo il fatto che la funzione diventa esponenzialmente piccola (nel caso stazionario \(|\alpha_1|<1\)) al crescere di \(t\). Questo fatto vale più in generale per processi \(\operatorname{ARIMA}\) stazionari. Un caso “limite” è quello dei processi a media mobile, ossia \(\operatorname{ARIMA}(0,0,q)\). In questo caso \(\alpha_i = 0\) e quindi \[ C(0,t) = 0\] è identicamente nulla se \(t>q\).
7.4.3 Stima dei parametri
A partire dall’osservazione di una serie storica \((x_t)_{t=0}^n\), come stimare i parametri di un processo ARIMA che la descrivono nel modo migliore? Abbiamo già osservato che la stima di massima verosimiglianza può fornire una risposta nel caso del rumore bianco gaussiano, della passeggiata aleatoria e dell’equazione lineare con smorzamento. In tutti e tre i casi il metodo si riduce alla minimizzazione dei residui quadratici (che appunto sono per ipotesi gaussiane indipendenti).
Si può quindi proporre lo stesso per un modello generale, dove tuttavia la nozione di residuo va chiarita, perché dall’equazione \[ p(L)(1-L)^d X_t = q(L)W_t\] è necessario ricavare il rumore bianco gaussiano \(W_t\), scrivendo \[ W_t = - \sum_{j=1}^q \beta_j W_{t-j} +p(L)(1-L)^dX_t\] Supponendo di osservare \(X_t = x_t\), questa equazione permette di definire ricorsivamente i residui \[ w_t = - \sum_{j=1}^q \beta_j w_{t-j} + p(L)(1-L)^d x_t,\] da cui infine la stima di massima verosimiglianza si ottiene minimizzando i residui quadratici (come funzione dei coefficienti \(\alpha = (\alpha_i)_{i=1}^p\) e \(\beta = (\beta_j)_{j=1}^q)\): \[ (\alpha_{\mle}, \beta_{\mle}) \in \operatorname{arg} \min_{\alpha, \beta} \sum_{t} w_t^2\] Per la risoluzione ci affidiamo a metodi numerici (in particolare se \(q\neq 0\)).
Remark. In R la funzione arima() oppure la funzione Arima() dalla libreria forecast permette di stimare i coefficienti di un modello arima (di ordine specificato) a partire dalle osservazioni. Avendo stimato i coefficienti la funzione forecast() permette anche di ottenere delle stime sui valori futuri come previsti dalle equazioni ricorsive del modello (con i coefficienti stimati) accompagnate da stime dell’incertezza (deviazione standard) dovute al termine di rumore bianco gaussiano.
La stima dei coefficienti non esaurisce tuttavia il problema, perché rimane da determinare l’ordine del modello, ossia la tripla di numeri \((p,d,q)\). È chiaro che, più grandi sono \(p\) e \(q\), più coefficienti avremo a disposizione e migliore sarà l’aderenza del modello ai dati osservati. Tuttavia, con una eccessiva aderenza si potrebbe incappare nel problema dell’overfit, e quindi non ottenere ad esempio previsioni ragionevoli. Per questo ci si serve di indicatori che tengano conto di tali fenomeni, come ad esempio gli indici AIC o BIC, utili per confrontare diversi modelli (sono da preferire i modelli con indici più piccolo). La libreria R forecast contiene il comando auto.arima() che restituisce automaticamente il miglior modello ARIMA a partire dai dati osservati e usando uno di questi criteri.
Esempio 7.3 Consideriamo il dataset precaricato in R AirPassengers che raccoglie una serie storica con i dati relativi al numero di passeggeri nelle linee aree internazionali dal 1949 al 1960. Usiamo i comandi descritti sopra sul modello.
plot(AirPassengers, lwd = 3, col = miei_colori[2])
Figura 7.3: Grafico della serie storica considerata
# Usiamo direttamente il comando
# auto.arima dalla libreria 'forecast'
library("forecast")
AP_auto <- auto.arima(AirPassengers)
# Con la funzione summary() possiamo
# vedere le informazioni principali
summary(AP_auto)## Series: AirPassengers
## ARIMA(2,1,1)(0,1,0)[12]
##
## Coefficients:
## ar1 ar2 ma1
## 0.5960 0.2143 -0.9819
## s.e. 0.0888 0.0880 0.0292
##
## sigma^2 estimated as 132.3: log likelihood=-504.92
## AIC=1017.85 AICc=1018.17 BIC=1029.35
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 1.342299 10.84619 7.86754 0.4206976 2.800458 0.245628 -0.001248475Vediamo che la funzione auto.arima() propone un modello ARIMA con stagionalità di periodo \(12\) (mesi) e ordine \((2,1,1)(0,1,0)\) (la seconda tripla si riferisce alla stagionalità). Precisamente, la funzione auto.arima() non determina automaticamente il periodo \(12\) e questo va indicato prima di applicarla ai dati, nel momento in cui si definisce un oggetto di tipo serie storica (in inglese time series) in R. Partendo da un vettore di dati osservati, il comando è ts(), che contiene l’opzione frequency (se non specificata è posta uguale ad \(1\)). In generale si può ricorrere alla funzione di autocorrelazione empirica (vedere la sezione successiva) o ad analisi spettrale per determinare eventuali stagionalità e il loro periodo. Un comando automatico è findfrequency().
# Con la funzione forecast() possiamo
# effettuare semplici previsioni per un
# numero di mesi futuri specificato
previsione <- forecast(AP_auto, 24)
# La funzione plot() si occupa di
# rappresentare sia i dati osservati
# che la previsione (e pure le bande di
# errore date dalle deviazioni standard
# stimate)
plot(previsione, xlab = "anno", ylab = "numero passeggeri",
col = miei_colori[2], lwd = 3)
Figura 7.4: dati osservati e previsione con un modello SARIMA
Come in molti altri problemi di stima, particolare attenzione va prestata alla gaussianità dei residui, pure forniti dalla funzione auto.arima().
par(mfrow = c(1, 2))
hist(AP_auto$residuals, col = miei_colori[1],
freq = FALSE, xlab = "residui", ylab = "frequenze",
main = "istogramma dei residui")
valori <- seq(min(AP_auto$residuals), max(AP_auto$residuals),
by = 0.1)
lines(valori, dnorm(valori, mean = mean(AP_auto$residuals),
sd = sd(AP_auto$residuals)), col = miei_colori[2],
lwd = 3)
qqnorm(AP_auto$residuals, col = miei_colori[1],
pch = 16)
qqline(AP_auto$residuals, col = miei_colori[2],
lwd = 3)
Figura 7.5: istogramma e QQ-plot dei residui