12.8 Covariâncias e correlações
12.8.1 Definições
Stationary time series are unchanging in respect to their general structure. (Wold 1938, 1)
Estacionariedade estrita: ocorre quando uma série temporal é invariante a deslocamentos (shifts), i.e., para qualquer subconjunto de tamanho \(n\) e qualquer inteiro \(\tau\), \((Y_{t_1},Y_{t_2},\ldots,Y_{t_n})\) tem a mesma distribuição conjunta de \((Y_{t_1+\tau},Y_{t_2+\tau},\ldots,Y_{t_n+\tau})\) (Peng 2020).
Estacionariedade fraca ou de segunda ordem: ocorre se a média for constante e a covariância entre quaisquer dois valores depender apenas da diferença de tempo entre esses dois valores (e não do valor de \(t\) em si), i.e., \(E(Y_t)=\mu\) e \(Cov(Y_t,Y_{t+\tau})=\gamma_{\tau}\) (Peng 2020).
Ruído branco: “The simplest possible example of a stationary random sequence is white noise. This consists of a sequence of mutually independent random variables, each with mean zero and finite variance \(\sigma^2\). Its autocovariance function is \(\gamma(t,s)=\sigma^2\)” if \(t=s\), 0 if \(t \ne s\). (Diggle 1990, 13–14)
Ruído branco gaussiano: É um ruído branco com distribuição normal.
Exemplo 12.9 A função stats::arima.sim() permite simular séries de ruído branco. Note que não se rejeita a normalidade de 94.3% das séries simuladas com \(\alpha=0.05\), uma vez que rand.gen = rnorm. Veja ?stats::arima.sim() para mais detalhes.
M <- 1000
wn <- vector(length = M)
for(i in 1:M){
set.seed(i); sim <- arima.sim(model = list(order = c(0, 0, 0)), n = 200)
wn[i] <- shapiro.test(sim)$p.value
}
sum(wn > 0.05)/M## [1] 0.943
12.8.2 Lags
plot_lag <- function(x, max.lag = 12){
par(mfrow=c(3,4))
n <- length(x)
for(i in 1:max.lag){
plot(x[1:(n-i)],x[(i+1):n], main = paste0('k = ',i))
}
par(mfrow=c(1,1))
}
12.8.3 Autocovariância
Considerando a abordagem de (Diggle 1990, 34–39), a (função de) autocovariância amostral de uma variável aleatória estacionária \(Y_{t}\) para \(k\) inteiro é \[\begin{equation} \hat{\gamma}(k) = cov(Y_{t},Y_{t-k}) = \frac{1}{n} \sum_{i=k+1}^{n} (y_{i}-\bar{y}) (y_{i-k}-\bar{y}) \tag{12.7} \end{equation}\]
Conforme (Venables and Ripley 2002, 390), o divisor \(n\) é utilizado mesmo que haja \(n-|t|\) termos, de forma a garantir que a sequência \(\hat{\gamma}(k)\) é a sequência de covariâncias de uma série temporal estacionária de segunda ordem.
12.8.4 Autocorrelação
Os lags que extrapolam os limites de
acf()indicam a ordem \(q\) do \(ARIMA(p,d,q)\).
Sendo \(\gamma(0)\) a variância de cada \(Y_{t}\), a função de autocorrelação amostral é \[\begin{equation} \hat{\rho}(k) = cor(Y_{t},Y_{t-k}) = \frac{\hat{\gamma}(k)}{\hat{\gamma}(0)} \tag{12.8} \end{equation}\]
Propriedades
- \(\rho(k)=\rho(-k)\)
- \(-1 \le \rho(k) \le +1\)
- Se \(Y_{t}\) e \(Y_{t-k}\) são independentes, então \(\rho(k)=0\)
A função stats::acf() calcula (e por padrão plota) estimativas da função de autocorrelação e autocovariância amostrais. Apresenta por padrão \(\lceil 10\log_{10}n \rceil\) lag.max para séries univariadas, onde \(n\) é o número de observações.
# simulando ruído branco gausssiano
set.seed(42); z <- rnorm(100)
# autocovariância
acf(z, type = 'covariance', main = 'stats::acf(type = "covariance")')
Os valores podem ser obtidos a partir as Eq. (12.7) e (12.8).
# autocovariância
acov <- function(x,k){
n <- length(x)
m <- mean(x)
g <- 0
for(i in (k+1):n){
g <- g + (1/n)*(x[i]-m)*(x[i-k]-m)
}
return(g)
}
round(sapply(0:16, acov, x = lh), 5)## [1] 0.29792 0.17146 0.05417 -0.04313 -0.05208 -0.04458 -0.00625 -0.00604 -0.00125 -0.04042 -0.04583 -0.02896 0.01458
## [14] 0.03562 0.02583 0.03542 0.04500##
## Autocovariances of series 'lh', by lag
##
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0.29792 0.17146 0.05417 -0.04312 -0.05208 -0.04458 -0.00625 -0.00604 -0.00125 -0.04042 -0.04583 -0.02896 0.01458 0.03562
## 14 15 16
## 0.02583 0.03542 0.04500# autocorelação
acor <- function(x,k){
r <- acov(x,k)/acov(x,0)
return(r)
}
round(sapply(0:16, acor, x = lh), 3)## [1] 1.000 0.576 0.182 -0.145 -0.175 -0.150 -0.021 -0.020 -0.004 -0.136 -0.154 -0.097 0.049 0.120 0.087 0.119 0.151##
## Autocorrelations of series 'lh', by lag
##
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
## 1.000 0.576 0.182 -0.145 -0.175 -0.150 -0.021 -0.020 -0.004 -0.136 -0.154 -0.097 0.049 0.120 0.087 0.119 0.151A função stats::acf() também permite tratar os NA. na.pass retorna o objeto inalterado, veja exemplos aqui.
# autocovariância
acf(forecast::gold, type = 'covariance',
na.action = na.pass,
main = 'stats::acf(type = "covariance")')
12.8.5 Autocorrelação parcial
Os lags que extrapolam os limites de
pacf()indicam a ordem \(p\) do \(ARIMA(p,d,q)\).
Seguindo a abordagem de (Enders 2014, 64–65), a autocorrelação parcial entre \(Y_{t}\) e \(Y_{t-k}\) elimina os efeitos dos valores intermediários \(Y_{t-1}\) e \(Y_{t-(k-1)}\). Sua estimativa pode ser obtida por \[\begin{equation} \hat{\phi}(1,1) = \hat{\rho}(1) \tag{12.9} \end{equation}\] \[\begin{equation} \hat{\phi}(k,k) = \frac{\hat{\rho}(k) - \sum_{i=1}^{k-1} \hat{\phi}(k-1,i) \; \hat{\rho}(k-i)}{1 - \sum_{i=1}^{k-1} \hat{\phi}(k-1,i) \; \hat{\rho}(i)}, \;\;\;\; k = 2,3,\ldots \tag{12.10} \end{equation}\] \[\begin{equation} \hat{\phi}(s,j) = \hat{\phi}(s-1,j) - \hat{\phi}(s,s) \hat{\phi}(s-1,s-j) \tag{12.11} \end{equation}\] onde \(j = 1,2,\ldots,s-1\). Note que \(\hat{\phi}(1,1) = corr(Y_{t},Y_{t-1})\) e \(\hat{\phi}(k,k) = corr(Y_{t}-\hat{Y}_{t},Y_{t-k}-\hat{Y}_{t-k})\), onde \(\hat{Y}_{t}\) e \(\hat{Y}_{t-k}\) são combinações lineares de \(\{Y_{t-1},\ldots,Y_{t-(k-1)}\}\) que minimizam o EQM de \(Y_{t}\) e \(Y_{t-k}\), respectivamente. Para mais detalhes veja (Durbin 1960).
As funções stats::pacf() e stats::acf(type = 'partial') calculam autocorrelações parciais.
# simulando ruído branco gausssiano
set.seed(42); z <- rnorm(100)
pacf(z, na.action = na.pass, main = 'stats::pacf()')
Exercício 12.16 Considere os conceitos de autocovariância, autocorrelação e autocorrelação parcial.
- Mostre que \(\gamma(0)\) é a variância de \(Y_{t}\).
- Verifique que \(\hat{\phi}(2,2) = \frac{\hat{\rho}(2) - \hat{\rho}(1)^2}{1 - \hat{\rho}(1)^2}\).
- Crie uma função em Python para calcular a autocovariância. Compare com
stats::acf()estatsmodels.tsa.stattools.acovf(). - Crie uma função em Python para calcular a autocorrelação. Compare com
stats::acf()estatsmodels.tsa.stattools.acf(). - Crie uma função em Python para calcular a autocorrelação. Compare com
stats::pacf()estatsmodels.tsa.stattools.pacf().
12.8.6 Covariância e correlação cruzadas
De acordo com (Diggle 1990, 202–3), a função de covariância cruzada de um processo aleatório estacionário bivariado \(\{(X_t,Y_t)\}\) é o conjunto de números \(\gamma_{xy}(k)\) definido, para cada inteiro \(k\), por \[\begin{equation} \gamma_{xy}(k) = cov(X_{t},Y_{t-k}) \tag{12.12} \end{equation}\]
Note que \(\gamma_{xx}(k)\) e \(\gamma_{yy}(k)\) são as funções de autocovariância de \(\{X_t\}\) e \(\{Y_t\}\) respectivamente, conforme Eq. (12.7). Da mesma forma, \[\gamma_{xy}(-k) = cov(X_t,Y_{t+k}) \\ \hspace{34pt} = cov(Y_{t+k},X_t) \\ \hspace{74pt} = cov(Y_t,X_{t-k}) = \gamma_{yx}(k)\]
A função de correlação cruzada de \(\{(X_t,Y_t)\}\) é \[\begin{equation} \rho_{xy}(k) = \frac{\gamma_{xy}(k)}{\sqrt{\gamma_{xx}(0) \gamma_{yy}(0)}} \tag{12.13} \end{equation}\]
O estimador da função de correlação cruzada é o correlograma cruzado tal que \[\begin{equation} \hat{\rho}_{xy}(k) = \frac{\hat{\gamma}_{xy}(k)}{\sqrt{\hat{\gamma}_{xx}(0) \hat{\gamma}_{yy}(0)}} \tag{12.14} \end{equation}\] onde \[\hat{\gamma}_{xx}(0) = \frac{1}{n} \sum_{t=1}^{n} (x_t-\bar{x})^2, \;\;\;\;\; \hat{\gamma}_{yy}(0) = \frac{1}{n} \sum_{t=1}^{n} (y_t-\bar{y})^2,\] \[\bar{x} = \frac{1}{n} \sum_{t=1}^{n} x_t, \;\;\;\;\; \bar{y} = \frac{1}{n} \sum_{t=1}^{n} y_t,\] \[\hat{\gamma}_{xy}(k) = \left\{ \begin{array}{l l} n^{-1} \sum_{t=k+1}^{n} (x_{t}-\bar{x}) (x_{t-k}-\bar{y}), & \quad k \ge 0\\ n^{-1} \sum_{t=1}^{n+k} (x_{t}-\bar{x}) (x_{t-k}-\bar{y}), & \quad k < 0\\ \end{array} \right. \]
A função stats::ccf calcula as estimativas da correlação cruzada e da covariância cruzada de duas séries univariadas.
Exemplo 12.10 Considere as séries temporais que das mortes mensais por bronquite, enfisema e asma no Reino Unido entre 1974–1979, ambos os sexos (ldeaths), homens (mdeaths) e mulheres (fdeaths), apresentada por (Diggle 1990, 238).
# covariância cruzada
stats::ccf(mdeaths, fdeaths, type = 'covariance',
ylab = 'cross-covariance', main = 'stats::ccf(type = "covariance")')
# correlação cruzada
stats::ccf(mdeaths, fdeaths,
ylab = 'cross-correlation', main = 'stats::ccf(type = "correlation")')
12.8.7 Periodograma
The periodogram is a summary description based on a representation of an observed time series as a superposition of sinusoidal waves of various frequencies. (Diggle 1990, 47)
A função stats::spectrum() estima a densidade espectral de uma série temporal. É uma função wrapper (empacotada) que chama os métodos spec.pgram() (periodograma via FFT com suavizador opcional) e spec.ar() (periodograma via modelo AR). Segundo a documentação, autores como (Thomson 1990) alertam fortemente que os espectros de modelos AR podem ser enganosos.
# janela gráfica 2x2
par(mfrow = c(2,2))
# Ruído branco gaussiano
stats::spectrum(z)
stats::spectrum(z, spans = 3)
stats::spectrum(z, spans = c(3,3))
stats::spectrum(z, spans = c(3,5))
# Luteinizing Hormone in Blood Samples
stats::spectrum(lh)
stats::spectrum(lh, spans = 3)
stats::spectrum(lh, spans = c(3,3))
stats::spectrum(lh, spans = c(3,5))
# Monthly corticosteroid drug subsidy in Australia from 1991 to 2008
stats::spectrum(h02)
stats::spectrum(h02, spans = 3)
stats::spectrum(h02, spans = c(3,3))
stats::spectrum(h02, spans = c(3,5))
# bivariate example - plots marginal spectra: now plot coherency and phase
mfdeaths <- ts.union(mdeaths, fdeaths)
mfdeaths.spc <- spec.pgram(mfdeaths, spans = c(3,3))
plot(mfdeaths.spc, plot.type = 'coherency')
plot(mfdeaths.spc, plot.type = 'phase')
Exercício 12.17 Considere o periodograma.
- Leia a lição 6.1 The Periodogram do curso STAT 510 - Applied Time Series Analysis da Penn State.
- Leia a discussão What is the confidence interval calculated in a spectral density periodogram in R?.
O pacote beyondWhittle (Meier et al. 2024) fornece implementações de procedimentos paramétricos, não paramétricos e semiparamétricos bayesianos para séries temporais univariadas e multivariadas. O pacote é baseado nos métodos apresentados em (Kirch et al. 2019), (Meier 2018) e (Tang et al. 2023).
# specify grid-points at which the tv-PSD is evaluated
res_time <- seq(0, 1, by = 0.005)
freq <- pi * seq(0, 1, by = 0.01)
# Ruído branco gaussiano
result <- beyondWhittle::gibbs_bdp_dw(z, m = 3, Ntotal = 10000, burnin = 2000,
res_time = res_time, freq = freq)
par(mfrow = c(2,2))
plot(result)
##
## Bayes factor 1.398142
## upper bound 27.280819# Luteinizing Hormone in Blood Samples
result <- beyondWhittle::gibbs_bdp_dw(lh, m = 3, Ntotal = 10000, burnin = 2000,
res_time = res_time, freq = freq)
par(mfrow = c(2,2))
plot(result)
##
## Bayes factor 6.206386
## upper bound 27.280819# Monthly corticosteroid drug subsidy in Australia from 1991 to 2008
result <- beyondWhittle::gibbs_bdp_dw(h02, m = 3, Ntotal = 10000, burnin = 2000,
res_time = res_time, freq = freq)
par(mfrow = c(2,2))
plot(result)
##
## Bayes factor 0.1023031
## upper bound 27.2808192