7.5 Stima della funzione di autocovarianza

In questa sezione affrontiamo il problema generale di stimare la funzione di media e di autocovarianza di un processo \(X\) a partire dall’osservazione dei valori per \(n\) tempi consecutivi \((X_t)_{t=0}^{n-1} = (x_t)_{t=0}^{n-1}\), anche detta nel linguaggio statistico una serie storica.

Si può pensare a questo problema come ad una generalizzazione del problema di stimare valor medio e varianza di una famiglia di variabili aleatorie indipendenti, tutte con la stessa legge. In particolare, abbiamo affrontato il caso gaussiano nella Sezione 5.4 e ottenuto come stime di massima verosimiglianza la media e la covarianza campionarie.

In questo caso introduciamo, invece dell’ipotesi di indipendenza, la stazionarietà del processo, ossia la matrice di covarianza delle variabili \((X_t)_{t =0}^{n-1}\) è costante sulle diagonali (oltre ad essere simmetrica): \[ (C(s,t))_{s,t=0}^{n-1} = (C(|t-s|))_{s,t=0}^{n-1}\] e supponiamo pure che il processo \(X\) sia gaussiano e centrato (ossia la funzione di media sia nota e costantemente nulla): in questo modo possiamo scrivere esplicitamente la funzione di verosimiglianza (supponendo che la matrice di covarianza sia invertibile) \[\begin{split} L( C ; x) & = p((X_t)_{t=0}^{n-1} = (x_t)_{t=0}^{n-1} | C) \\ & \propto \exp\bra{ - \frac 1 2 x^T C^{-1} x} \frac{1}{\sqrt{ \det C }},\end{split}\] dove abbiamo posto, per alleggerire la notazione, \(x = (x_t)_{t=0}^{n-1}\).

Possiamo determinare la stima di massima verosimiglianza per \(C\) con i soliti passaggi, ossia passando al logaritmo e cambiando di segno: si tratta di minimizzare la funzione \[ C \mapsto x^T C^{-1} x + \log \det C.\] Tuttavia è comunque difficile calcolare esplicitamente \(C_{\mle}\) (ma si può ricorrere a metodi numerici).

Per proseguire analiticamente e ottenere delle espressioni elementari conviene introdurre una ulteriore ipotesi matematica nella struttura della matrice di covarianza: non solo supponiamo che sia costante sulle diagonali, ma anche che sia circolante, ossia che valga l’identità, per ogni \(k =1,2,\ldots, n-1\), \[ C(k) = C(n-k).\] Questa ipotesi è giustificabile solo per semplificare i calcoli, non vi è una ragione particolare per ritenere che la funzione di autocovarianza di un processo stazionario la soddisfi, eccetto al più nel caso in cui il processo sia periodico di periodo \(n\), ossia valga \(X_{t+n} = X_t\) per ogni \(t\). Ma ricordiamo che \(n\) è solamente il numero di osservazioni, e di solito se vi è una periodicità, anche approssimata (si parla in tal caso di una stagionalità) essa è di periodo molto minore di \(n\). In ogni caso, una volta trovata \(C_{\mle}\) con questa ipotesi, possiamo proporre una modifica per il caso generale. Un’altra possibilità sarebbe di ragionare nel limite \(n \to \infty\), ma quento introdurrebbe ulteriori problemi tecnici.

Il vantaggio di supporre che la matrice \(C\) sia circolante è che, passando alla trasformata di Fourier a tempi finiti, essa diventa diagonale. Precisamente, usando la notazione della Sezione B.1, introduciamo la matrice \(F \in \mathbb{C}^{n\times n}\), \[ F_{\xi t} = e^{-2 \pi i \xi t/n},\] per \(\xi = 0,1, \ldots, (n-1)\), in modo che la trasformata di Fourier a tempi finiti di \((x_t)_{t =0}^{n-1}\) sia \[\hat x (\xi ) = \sum_{t=0}^{n-1} x_t e^{- 2 \pi i \xi t/n} = \sum_{t=0}^{n-1} F_{\xi t} x_t,\] ossia \(\hat x = F x\). Ricordando che \(x\) è l’osservazione del processo \(X\), la matrice di covarianza del vettore aleatorio \(FX\) si trasforma come al solito (formula per le trasformazioni affini) \[ \Sigma_{ FX} = F \Sigma_X \bar{F}^T = FC \bar{F}^T,\] dove l’unico accorgimento è che, essendo la matrice \(F\) complessa, la formula va modificata introducendo il coniugato del trasposto \(\bar{F}^T\) (invece del semplice trasposto).

Questo cambio di coordinate dalla base dei “tempi” a quella delle “frequenze” ha l’effetto di diagonalizzare la matrice delle covarianze. Infatti, scrivendo \[ \hat{C}(\xi) = \sum_{k=0}^{n-1} e^{-2 \pi i \xi k/n} C(k),\] troviamo che \[ \begin{split} (\Sigma_{FX})_{\xi \ell} &= ( F C \bar{F}^T)_{\xi \ell} \\ & =\sum_{s,t=0}^{n-1} F_{\xi s } C(s,t) \bar{F}_{t \ell}\\ & =\sum_{s,t=0}^{n-1} e^{-2\pi i \xi s/n} C(s-t)e^{2\pi i \ell t/n}\\ & = \sum_{t=0}^{n-1} e^{2\pi i \ell t/n} e^{-2\pi i \xi t/ n} \sum_{s=0}^{n-1} e^{-2\pi i \xi (s-t)/n} C(s-t)\\ & = \hat{C}(\xi) \sum_{t=0}^{n-1} e^{-2\pi i (\xi-\ell)t/n} \\ & = \hat{C}(\xi) n \delta_0(\xi-\ell) \end{split} \] dove ricordiamo la notazione \(\delta_0(x)\) per la funzione che vale \(1\) se \(x=0\), e \(0\) altrimenti. Abbiamo usato l’ipotesi che \(C\) sia circolante per dedurre che, per ogni \(t\), vale \[ \sum_{s=0}^{n-1} e^{-2\pi i \xi (s-t)/n} C(s-t) = \sum_{k=0}^{n-1} e^{-2\pi i \xi k/n} C(k) = \hat{C}(\xi).\] Notiamo tra l’altro che i numeri \(n \hat{C}(\xi)\) sono quindi (multipli degli) autovalori della matrice di covarianza \(C\) e quindi sono tutti positivi (e non nulli avendo supposto che \(C\) sia invertibile). In queste nuove coordinate, le componenti del vettore \(FX\) sono non correlate e quindi, essendo gaussiane, indipendenti. La verosimiglianza assume quindi una espressione molto più trattabile: \[ L( C; \hat x) = p( FX = \hat x | C ) \propto \exp\bra{ - \frac 1 {2n} \sum_{\xi=0}^{n-1} \frac{| \hat x (\xi)|^2 }{\hat C(\xi)} } \frac{1}{\sqrt{\prod_{\xi=0}^{n-1} \hat C(\xi)}} \] Di conseguenza, passando ai logaritmi e moltiplicando tutto per \(-2\), la stima di massima verosimiglianza si ottiene minimizzando la funzione \[ C \mapsto \sum_{\xi=1}^{n-1} \sqa{ \frac 1 n\frac{ | \hat x(\xi)|^2}{\hat{C}(\xi)} + \log \hat{C}(\xi)}\] A questo punto si può trattare formalmente le \(\hat{C}(\xi)\) come i parametri da stimare e ottenere le stime di massima verosimiglianza \[ \hat{C}_{\mle} (\xi) = \frac{ \abs{ \hat x(\xi)}^2}{n} = \frac{1} n \abs{ \sum_{t=0}^{n-1} x_t e^{-2\pi i t \xi /n}}^2.\] Invertendo la trasformata di Fourier, possiamo infine ottenere le stime di massima verosimiglianza cercate \[\begin{split} C_{\mle} (k) & = \frac 1 n \sum_{\xi=0}^{n-1} \hat{C}_{\mle}(\xi) e^{2 \pi i k \xi/n} \\ & =\frac 1 n \sum_{\xi = 0}^{n-1} \frac 1 n \abs{ \sum_{t=0}^{n-1} x_t e^{-2\pi i t \xi/n}}^2 e^{2 \pi i k \xi/n}\\ & = \frac 1 {n^2} \sum_{\xi=0}^{n-1} \sum_{s,t=0}^{n-1} x_t x_s e^{-2\pi i t \xi/n}e^{2\pi i s \xi/n} e^{2 \pi i k \xi/n}\\ & = \frac 1 {n^2} \sum_{s,t=0}^{n-1} x_t x_s \sum_{\xi=0}^{n-1} e^{-2 \pi i (t-s-k)\xi/n}\end{split}\] Discutiamo l’ultima epressione che abbiamo trovato: sicuramente, se \(t=s+k\), allora il termine \[ \sum_{\xi=0}^{n-1} e^{-2 \pi i (t-s-k)\xi/n} = \sum_{\xi=0}^{n-1} 1 = n,\] tuttavia questo non è l’unico caso, perché potrebbe anche accadere che \(t = s+k-n\) e allora ugualmente si avrebbe \[ \sum_{\xi=0}^{n-1} e^{-2 \pi i (t-s-k)\xi/n} = \sum_{\xi=0}^{n-1} e^{2 \pi i \xi } = \sum_{\xi=0}^{n-1} 1 = n.\] In tutti gli altri casi possibili, si trova invece \[ \sum_{\xi=0}^{n-1} e^{-2 \pi i (t-s-k)\xi/n} = 0,\] e quindi concludiamo che \[ C_{\mle} (k) = \frac 1 n \sum_{s=0}^{n-1-k}x_s x_{s+k} + \frac 1 n \sum_{s=n-k+1}^{n-1}x_s x_{s+k-n}.\]

La formula trovata è la somma due contributi, i quali ricordano rispettivamente il primo la covarianza campionaria tra \(X_s\) e il processo “traslato” avanti nel tempo di \(k\) istanti, \(X_{s+k}\) e il secondo la covarianza campionaria tra \(X_s\) e il traslato indietro di \(n-k\) istanti, \(X_{s+k-n}\). Questo riflette la condizione ulteriore che abbiamo imposto nella funzione di autocovarianza, ossia che \(C(k)= C(n-k)\). Osserviamo che la prima somma consiste di \(n-k\) termini, mentre la seconda di \(k\) termini perciò per \(k\) molto più piccolo di \(n\), possiamo supporre che la prima dia un contributo più rilevante.

Tornando al caso generale, possiamo proporre come stima per \(C\) semplicemente il primo termine, ossia la covarianza campionaria tra \(X_s\) e il traslato \(X_{s+k}\). Ovviamente questo pone dei problemi, perché le osservazioni disponibili sono solo fino al tempo \(n-1\) e quindi dovremo sommare solo \(n-k\) termini. In generale definiamo allora la funzione di autocovarianza campionaria (o empirica) come \[ c(k) = \frac{1}{n-k}\sum_{s=0}^{n-1-k} (x_s - \bar{x}_0)(x_{s+k} - \bar{x}_{k}),\] dove le medie campionarie \(\bar{x}_0\), \(\bar{x}_k\) sono rispettivamente sui primi \(n-k\) e sugli ultimi \(n-k\) valori. Equivalentemente \(c(k)\) può essere pensato come la covarianza tra due variabili aleatorie definite nel seguente modo: si sceglie \(S \in \cur{0,1, \ldots, n-1-k}\) casuale uniforme e si considerano i valori \(x_{S}\) (prima variabile) e \(x_{S+k}\) (seconda variabile). Basandoci su questa idea possiamo allora definire anche le varianze campionarie \[ \sigma^2_0 = \frac{1}{n-k}\sum_{s=0}^{n-1-k} (x_s - \bar{x}_0)^2\] e \[ \sigma^2_k = \frac{1}{n-k}\sum_{s=0}^{n-1-k} (x_{s+k} - \bar{x}_k)^2\] e quindi la funzione di autocorrelazione campionaria, data da \[ \operatorname{acf}(k) = \frac{c(k)}{\sigma_0 \sigma_k},\] che assume sempre valori tra \([-1,1]\) (è il coefficiente di correlazione tra le due variabili \(x_{S}\) e \(x_{S+k}\) definite sopra). Questa funzione è la più utilizzata in pratica per stimare la funzione di autocorrelazione di un processo a partire dalle osservazioni. In R è disponibile tramite il comando acf().

Esempio 7.4 Consideriamo la serie dei residui ottenuti dalla stima di un modello ARIMA con stagionalità sulla serie AirPassengers. Ci aspettiamo che sia rappresentabile come un rumore bianco gaussiano, di cui la funzione di autocovarianza è molto semplice (è \(\delta_0\)).

ACF_residui <- acf(AP_auto$residuals)

Figura 7.6: autocorrelazione empirica dei residui

# Osserviamo che nelle ascisse
# l'intervallo temporale (lag) è
# espresso in multipli del periodo di
# 12 mesi, per via della stagionalità
# della serie di partenza.

La funzione di autocorrelazione campionaria è uno strumento utile per determinare eventuali stagionalità e il loro periodo \(k\), che si ottiene in corrispondenza di “picchi” della funzione (più precisamente, massimi locali). Bisogna tuttavia osservare che in presenza di una componente lineare (detto anche trend) della serie storica, la funzone di autocorrelazione campionaria tende ad essere uniformemente a valori grandi (e quindi maschera le stagionalità).

acf(AirPassengers)

Figura 7.7: la funzione di autocorrelazione è uniformemente grande per via del trend

Un modo efficace per rimuovere questo effetto è passare ad una derivata discreta della serie osservata, tramite la funzione diff().

acf(diff(AirPassengers))

Figura 7.8: la derivata discreta rimuove il trend e la funzione di autocorrelazione evidenzia la stagionalità a 12 mesi

Oltre alla funzione di autocorrelazione campionaria, ci possiamo chiedere cosa accada dei calcoli svolti passando alla trasformata di Fourier, nel caso in cui l’ipotesi semplificativa \(C(k) = C(n-k)\) non sia valida. Anche in questo caso, i passaggi sono approssimativamente validi purché \(n\) diventi molto grande. In tal caso la stima che abbiamo trovato \[ \hat{C}(\xi) \approx \frac{ \abs{ \hat x(\xi)}^2}{n} = \frac{1} n \abs{ \sum_{t=0}^{n-1} x_t e^{-2\pi i t \xi /n}}^2\] diventa esatta nel limite \(n \to \infty\) e passando al valor medio. Precisamente, vale il seguente teorema.

Teorema 7.2 (di Wiener-Khinchin) Sia \((X_t)_{t =0}^\infty\) un processo a valori reali, stazionario in senso lato, con media nulla \(\E{X_t} = 0\), e tale che \[ \sum_{k=0}^\infty |C(k)| < \infty.\] Per ogni \(\xi \in [0,1]\), si ponga \[ \hat{C}(\xi) = \sum_{k\in \mathbb{Z}} C(|k|)e^{-2\pi i k \xi}.\] Allora vale il limite \[ \lim_{ n \to \infty} \frac{1} n \E{ \abs{ \sum_{t=0}^{n-1} X_t e^{-2\pi i t \xi }}^2} = \hat{C}(\xi).\]

In virtù della formula sopra, la trasformata di Fourier \(\hat{C}\) della funzione di autocovarianza \(C\) è detta anche densità spettrale di potenza (in inglese power spectral density), perché rappresenta il valore medio (sia nel tempo che nel senso della probabilità) dell’energia del processo associata alla frequenza \(\xi\).

Esempio 7.5 Lo spettrogramma di una serie storica (ossia il modulo della trasformata di Fourier, eventualmente in scala logaritmica) è uno strumento utile per determinare eventuali periodicità. In R si può utilizzare direttamente il comando spectrum(), che più precisamente fonisce una stima della densità spettrale di potenza.

spectrum(AirPassengers)

Figura 7.9: Stima della densità spettrale di potenza della serie AirPassengers

Lo spettro indica chiaramente la stagionalità (ricordiamo che è già indicato che un periodo corrisponde a \(12\) mesi). Vediamo un esempio diverso nel caso dei residui (che ricordiamo sono modellizzabili come un rumore bianco gaussiano).

spectrum(AP_auto$residuals)

Figura 7.10: Stima della densità spettrale di potenza dei residui della serie AirPassengers dopo un fit con un modello ARIMA

In questo caso lo spettro non ha picchi particolari, segno in particolare dell’assenza di stagionalità. Abbiamo anche calcolato che nel caso di rumore bianco la densità spettrale di potenza (teorica) è costante.