Tópicos adicionais no domínio do tempo

Apresentamos aqui um material que pode ser considerado especial ou avançado no domínio do tempo. Dedicamos este espaço para estudar um dos tópicos no domínio do tempo mais úteis e interessantes, modelos de espaço de estados. Consequentemente, não abordamos agora modelos de espaço de estado ou tópicos relacionados, dos quais existem muitos. O conteúdo aqui é de tópicos independentes que podem ser lidos em qualquer ordem.

A maioria das seções depende de um conhecimento básico dos modelos ARMA, previsão e estimação. Algumas seções, por exemplo, a seção sobre modelos de memória longa, requerem algum conhecimento de análise espectral e tópicos relacionados. Além da memória longa, discutimos o teste de raiz unitária, modelos GARCH, modelos de limiar, funções de regressão com atraso ou de transferência e tópicos selecionados em modelos ARMAX multivariados.

1. Modelos ARMA de memória longa e diferenciação fracionária
2. Teste da raiz unitária
3. Modelos GARCH
4. Modelos de limiar
5. Regressão com retardo e modelagem da função de transferência
6. Modelo ARMAX: ARMA multivariado
7. Exercícios
8. Bibliografia

1. Modelos ARMA de memória longa e diferenciação fracionária

O processo ARMA(p,q) convencional é freqüentemente chamado de processo de memória curta porque os coeficientes na representação \[\begin{equation} X_t \, = \, \sum_{j=0}^\infty \psi_j W_{t-j}, \end{equation}\] obtidos resolvendo \[\begin{equation} \phi(z)\psi(z) \, = \, \theta(z), \end{equation}\] são dominados por decaimento exponencial.

Conforme apontado, este resultado implica que a função de autocorrelação (ACF) do processo de memória curta satisfaz \(\rho(h)\to 0\) exponencialmente rápido quando \(h\to \infty\). Quando o ACF amostral de uma série temporal decai lentamente então, diferenciar a série até parecer estacionária tem sido uma solução. Seguindo esse conselho com a série de variações glaciais, a primeira diferença dos logaritmos dos dados é representada como uma média móvel de primeira ordem.

Uma análise mais aprofundada dos resíduos leva ao ajuste de um modelo ARIMA(1,1,1), \[\begin{equation} \nabla X_t \, = \, \phi \nabla X_{t-1}+W_t+\theta W_{t-1}, \end{equation}\] onde entendemos \(X_t\) como a série transformada por logaritmo. Em particular, as estimativas dos parâmetros e os erros padrão, foram \(\widehat{\phi} = 0.23_{(0.05)}\), \(\widehat{\theta} = -0.89_{(0.03)}\) e \(\widehat{\sigma}_{_W}^2 = 0.23\).

O uso da primeira diferença \(\nabla X_t=(1-B)X_t\), no entanto, às vezes pode ser uma modificação muito severa no sentido de que o modelo não-estacionário pode representar uma superdiferença do processo original. Séries temporais de memória longa ou persistente foram consideradas em Hosking (1981) e Granger and Joyeux (1980) como compromissos intermediários entre os modelos do tipo ARMA de memória curta e os processos não estacionários totalmente integrados na classe Box–Jenkins.

A maneira mais fácil de gerar uma série de memória longa é pensar em usar o operador de diferença \((1-B)^d\) para valores fracionários de \(d\), por exemplo, \(|d|<\frac{1}{2}\), para que uma série básica de memória longa seja gerada como \[\begin{equation} (1-B)^d X_t \, = \, W_t, \end{equation}\] onde \(W_t\) ainda indica ruído branco com variância \(\sigma_{_W}^2\).

Definição 1:

O modelo de média móvel autoregressivo fraccionada de ordem (p,q) para uma série temporal estacionária \(X_t\) \((t\geq 1)\), com média \(\mu\) pode ser escrito, utilizando a notação de Box-Jenkins como \[\begin{equation} \tag{1} \phi(B)\Delta(B)(X_t-\mu) = \Theta(B)W_t, \end{equation}\] onde \[\begin{equation} \Phi(B)=1-\phi_1 B - \cdots -\phi_p B^p, \qquad \Theta(B)=1-\theta_1 B - \cdots - \theta_q B^q, \end{equation}\] e \[\begin{equation} \Delta(B)=(1-B)^d = \displaystyle \sum_{i=1}^\infty {d \choose i}(-1)^i B^i \end{equation}\] sendo \(W_t\) um ruído branco com variância \(\sigma_{_W}^2\).

A série fracionariamente diferenciada acima, para \(|d| <\frac{1}{2}\), é freqüentemente chamada de ruído fracionário, exceto quando \(d\) é zero. Agora, \(d\) se torna um parâmetro a ser estimado juntamente com \(\sigma_{_W}^2\).

Hosking (1984) mostrou que \(\gamma_\ell = \mbox{Cov}(X_t,X_{t-\ell}) = O(\ell^{2d-1})\) quando \(\ell\to\infty\). Assim, para \(0<d<\frac{1}{2}\), \(\sum_\ell \gamma_\ell\) diverge e pode dizer-se que a série temporal tem uma componente de memória longa.

Diferenciar o processo original, como na abordagem de Box-Jenkins, pode ser considerado simplesmente atribuindo um valor de \(d = 1\). Essa ideia foi estendida à classe de componentes fracionados integrados ARMA ou ARFIMA, em que \(-\frac{1}{2} \leq d \leq \frac{1}{2}\); quando \(d\) é negativo, o termo antipersistente é usado.

Processos de memória longa ocorrem na hidrologia, ver Hurst (1951) e McLeod and Hipel (1978) e em séries ambientais, como os vários dados que analisamos anteriormente, para citar alguns exemplos. Os dados de séries temporais de memória longa tendem a exibir autocorrelações amostrais que não são necessariamente grandes, como no caso de \(d = 1\), mas persistem por um longo tempo. A figura abaixo mostra o ACF amostral, para um retardo ou lag de 100, da série de variações transformada por logaritmo, que exibe o comportamento clássico de memória longa:

library(astsa)
par(mfrow = c(1,2), mar=c(3,3,3,1), mgp=c(1.6,.6,0), pch=19)
acf(log(varve), 100);grid()            # comparar este ACF com o de uma caminhada aleatória
acf(cumsum(rnorm(1000)), 100);grid()   # ACF de uma caminhada aleatória

Figura 1: Esquerda: gráfico da função de auocorrelação amostral do logaritmo da série de variações glaciais (varve). Direita: mostramos a função de auocorrelação amostral de uma caminhada aleatória, para comparar.

Para investigar suas propriedades, podemos usar a expansão binomial \((d>-1)\), e escrever \[\begin{equation} W_t \, = \, (1-B)^d X_t \, = \, \displaystyle \sum_{j=0}^\infty \pi_j B^j X_t \, = \, \sum_{j=0}^\infty \pi_j X_{t-j}, \end{equation}\] onde \[\begin{equation} \pi_j \, = \, \frac{\Gamma(j-d)}{\Gamma(j+1)\Gamma(-d)}, \end{equation}\] com \(\Gamma(x+1)=x\Gamma(x)\) sendo a função gama. Similarmente \((d < 1)\), podemos escrever \[\begin{equation} X_t \, = \, (1-B)^{-d} W_t \, = \, \displaystyle \sum_{j=0}^\infty \psi_j B^j W_t \, = \, \sum_{j=0}^\infty \psi_j W_{t-j}, \end{equation}\] onde \[\begin{equation} \psi_j \, = \, \frac{\Gamma(j+d)}{\Gamma(j+1)\Gamma(d)}\cdot \end{equation}\]

Quando \(|d| < 0.5\) os processos \(W_t\) e \(X_t\) anteriores são processos estacionários bem definidos (ver Brockwell and Davis 1991a). No caso de diferenciação fracionária, no entanto, os coeficientes satisfazem \[\begin{equation} \sum_j \pi_j^2 < \infty \qquad \mbox{e} \qquad \sum_j \psi_j^2 < \infty, \end{equation}\] em oposição à soma absoluta dos coeficientes nos processos ARMA.

Após algumas manipulações não triviais, pode-se mostrar que o ACF de \(X_t\) é \[\begin{equation} \rho(h) \, = \, \frac{\Gamma(h+d)\Gamma(1-d)}{\Gamma(h-d+1)\Gamma(d)} \approx h^{2d-1}, \end{equation}\] para \(h\) grande. A partir disso, vemos que para \(0 < d < \frac{1}{2}\) \[\begin{equation} \sum_{h=-\infty}^\infty |\rho(h)| \, = \, \infty, \end{equation}\] e, portanto, o termo é de memória longa.

Para examinar uma série como a das variações para um possível padrão de memória longa, pode ser conveniente procurar maneiras de estimar \(d\). Observando a expressão dos \(\pi_j\) acima, podemos derivar as recursões \[\begin{equation} \pi_{j+1}(d) \, = \, \frac{(j-d)\pi_j(d)}{j+1}, \end{equation}\] para \(j = 0,1,\cdots\), com \(\pi_0(d)=1\). Maximizar a verossimilhança conjunta dos erros sob normalidade, digamos, \(W_t(d)\), envolverá a minimização da soma dos erros ao quadrado \[\begin{equation} Q(d) \, = \, \displaystyle \sum_t W_t^2(d)\cdot \end{equation}\]

O método de Gauss-Newton usual, leva à expansão: \[\begin{equation} W_t(d) \, = \, W_t(d_0) + W_t'(d_0)(d-d_0), \end{equation}\] onde \[\begin{equation} W_t'(d_0) \, = \, \displaystyle \left.\dfrac{\partial W_t}{\partial d} \right|_{d=d_0} \end{equation}\] e \(d_0\) é uma estimativa inicial ou palpite no valor de \(d\). A configuração da regressão usual leva a \[\begin{equation} d \, = \, d_0 - \dfrac{\displaystyle \sum_t W_t'(d_0)W_t(d_0)}{\displaystyle \sum_t W_t'(d_0)^2}\cdot \end{equation}\]

As derivadas são calculadas recursivamente pela diferenciação de \(\pi_{j+1}(d)\) sucessivamente em relação a \(d\): \[\begin{equation} \pi_{j+1}'(d) \, = \, \dfrac{(j-d)\pi_j'(d)-\pi_j(d)}{j+1}, \end{equation}\] onde \(\pi_0'(d)=0\). Os erros são calculados a partir de uma aproximação a \(W_t=(1-B)^d X_t\), a saber \[\begin{equation} W_t(d) \, = \, \displaystyle\sum_{j=0}^t \pi_j(d)X_{j-t}\cdot \end{equation}\]

Aconselha-se omitir vários termos iniciais do cálculo e iniciar a somas \(\sum_t W_t'(d_0)W_t(d_0)\) e \(\sum_t W_t'(d_0)^2\) em algum valor razoavelmente grande de \(t\), tudo isso para ter uma aproximação razoável.

O modelo fracionário de ordem (p,q) pode ser aproximado pelo modelo \[\begin{equation} \tag{2} \phi(B)\Delta_r(B)(X_t-\mu) = \Theta(B)W_t, \end{equation}\] em que \[\begin{equation} \Delta_r(B)= \displaystyle \sum_{i=0}^r {d \choose i}(-1)^i B^i\cdot \end{equation}\]

Decorre do teorema de Kakeya-Enström (Henrici 1974) que \(\Delta_r (B)=0\) tem raízes fora do círculo unitário. Assim, (2) é estacionário para qualquer \(r\geq 0\). Além disso, para \(r\) suficientemente grande e \(d\) não maior do que \(\frac{1}{2}\) a diferença nos modelos (1) e (2) pode tornar-se desprezável.

Exemplo 1: Ajuste de memória longa da série de variações glaciais.

par(mfrow = c(1,1), mar=c(3,3,3,1), mgp=c(1.6,.6,0), pch=19)
plot(log(varve), type='l', ylab=expression(log(X[~t])), 
     xlab="Tempo", main = "Logaritmo da série de variações glaciais")
grid()

Consideramos analisar a série de variações glaciais discutida em vários exemplos. A figura acima mostra a série em escala logarítmica. Em exemplo anteriores foi observado que esta série poderia ser modelada como um processo ARIMA(1,1,1).

Ajustamos então o modelo diferenciado fracionariamente à série ajustada pela média, \(X_t - \overline{X}\). A aplicação do procedimento iterativo de Gauss-Newton, iniciando com \(d =0.1\) e omitindo as primeiras 30 iterações no cálculo, leva a um valor final de \(d =0.384\), o que implica o conjunto de coeficientes \(\pi_j (0.384)\), como dado na figura abaixo com \(\pi_0 (0.384) = 1\).

library(fracdiff)
lvarve = log(varve)-mean(log(varve))
varve.fd = fracdiff(lvarve, nar=0, nma=0, M=30)
varve.fd$d # = 0.3841688

## [1] 0.3841688

varve.fd$stderror.dpq # = 4.589514e-06 (resultado questionável!!)

## [1] 4.589514e-06

p = rep(1,31)
for (k in 1:30){ p[k+1] = (k-varve.fd$d)*p[k]/(k+1) }
par(mfrow = c(1,1), mar=c(3,3,3,1), mgp=c(1.6,.6,0), pch=19)
plot(1:30, p[-1], ylab=expression(pi(d)), xlab="Index", type="h");grid()

Figura 2: Coeficientes \(\pi_j (0.384)\), \(j = 1,2,\cdots,30\) na representação \(\pi_{j+1}(d)\).

R usa um procedimento de máxima verossimilhança truncada que foi discutido em Haslett and Raftery (1989), um procedimento um pouco mais elaborado do que simplesmente zerar os valores iniciais. O valor padrão de truncamento no R é \(M = 100\). No caso padrão, a estimativa é \(\widehat{d}=0.37\) com aproximadamente o mesmo erro padrão questionável.

res.fd = diffseries(log(varve), varve.fd$d) # frac diff resíduos
res.arima = resid(arima(log(varve), order=c(1,1,1))) # arima resíduos
par(mfrow = c(1,2), mar=c(3,3,1,1), mgp=c(1.6,.6,0), pch=19)
acf(res.arima, 30, xlim=c(1,30), ylim=c(-.1,.1), main=""); grid()
acf(res.fd, 30, xlim=c(1,30), ylim=c(-.1,.1), main="");grid()

Figura 3: À esquerda, função de autocorrelação amostral (ACF) dos resíduos do modelo ARIMA(1,1,1) ajustado à série de variações glaciais (varve), à direita mostramos o ACF dos resíduos do ajuste do modelo de memória longa, \[\begin{equation} (1-B)^d X_t = W_t, \end{equation}\] com \(d =0.384\).

library(forecast)
par(mfrow = c(1,1), mar=c(3,2,1,1), mgp=c(1.6,.6,0), pch=19)
ts.plot(lvarve, fitted(Arima(lvarve, order=c(1,1,1))), col = c(6,4), tyepe="l", xlab="Tempo");grid()
legend(0,2.0, legend = c("Observados","Estimados"), col = c(6,4), lty = 1)

A previsão de processos de memória longa é semelhante à previsão dos modelos ARIMA. Ou seja, escrevendo \(W_t=(1-B)^d X_t\) e \(\pi_{j+1}'(d) = \big((j-d)\pi_j'(d)-\pi_j(d)\big)/(j+1)\) podem ser usados para obter as previsões truncadas \[\begin{equation} \widetilde{X}^n_{n+m} \, = \, -\displaystyle \sum_{j=1}^n \pi_j(\widehat{d})\widetilde{X}^n_{n+m-j}, \end{equation}\] para \(m=1,2,\cdots\). Os limites de erro podem ser aproximados usando \[\begin{equation} P^n_{n+m} \, = \, \widehat{\sigma}^2_{_W} \Big( \displaystyle\sum_{j=0}^{m-1} \psi_j^2(\widehat{d})\Big) \end{equation}\] onde \[\begin{equation} \psi_j(\widehat{d}) \, = \, \dfrac{(j+\widehat{d})\psi_j(\widehat{d})}{j+1}, \end{equation}\] com \(\psi_0(\widehat{d})=1\).

Nenhum componente óbvio do tipo ARMA de memória curta pode ser visto na função de autocorrelação dos resíduos da série de variações glaciais fracionada e diferenciada mostrada na figura do Exemplo 1. É natural, no entanto, que existam casos em que componentes substanciais do tipo memória curta também estarão presentes nos dados que exibem memória longa.

Portanto, é natural definir o processo ARFIMA(p,d,q) geral, \(-0.5 < d < 0.5\) como \[\begin{equation} \phi(B)\nabla^d (X_t-\mu) \, = \, \theta(B)W_t\cdot \end{equation}\] Escrevendo o modelo na forma \[\begin{equation} \phi(B)\pi_d (B) (X_t-\mu) \, = \, \theta(B)W_t \end{equation}\] fica claro como calculamos os parámetros para o modelo mais geral.

Previsões para o processo ARFIMA(p,d,q) podem ser feitas notando que possamos igualar coeficientes em \[\begin{equation} \phi(z)\psi(z) \, = \, (1-z)^{-d}\theta(z) \end{equation}\] e \[\begin{equation} \theta(z)\pi(z) \, = \, (1-z)^d\phi(z) \end{equation}\] para obter as representações \[\begin{equation} X_t \, = \, \mu+\sum_{j=0}^\infty \phi_j W_{t-j} \qquad \mbox{e} \qquad W_t \, = \, \sum_{j=0}^\infty \pi_j (X_{t-j}-\mu)\cdot \end{equation}\] Podemos então proceder e obter \(\widetilde{X}^n_{n+m}\) e \(\widetilde{P}^n_{n+m}\).

Tratamentos abrangentes de modelos de séries temporais de memória longa são dados nos textos de Beran (1994), Palma (2007) e Robinson (1995) e deve-se notar que várias outras técnicas para estimar os parámetros, especialmente os parámetros de memória longa, pode ser desenvolvido no domínio da frequência.

Nesse caso, podemos pensar nas equações como geradas por uma série autoregressiva de ordem infinita com coeficientes \(\pi_j\) dados anteriormente. Usando a mesma abordagem de antes, obtemos \[\begin{array}{rcl} f_X(\omega) & = & \displaystyle \dfrac{\sigma^2_{_W}}{\displaystyle \left| \sum_{k=1}^\infty \pi_k e^{-2\pi \, i \, k\omega}\right|^2} \\ & = & \displaystyle \sigma_{_W}^2 \left|1-e^{-2\pi \, i \, \omega} \right|^{-2d} \, = \, \big( 4\mbox{sen}^2(\pi\omega)\big)^{-d}\sigma^2_{_W}, \end{array}\]

como representações equivalentes do espectro de um processo de memória longa. O longo espectro de memória se aproxima do infinito quando a frequência \(\omega\to 0\).

O principal motivo para definir a aproximação de Whittle à log-verossimilhança é propor seu uso para estimar o parámetro \(d\) no caso de memória longa como uma alternativa ao método do domínio do tempo mencionado anteriormente. A abordagem no domínio do tempo é étil devido à sua simplicidade e erros padrão facilmente calculáveis. Pode-se também usar uma abordagem de verossimilhança exata, desenvolvendo uma forma de inovação da verossimilhança, como em Brockwell and Davis (1991b).

Para a abordagem aproximada usando a verossimilhança de Whittle, consideramos usar a abordagem de Fox and Taqqu (1986) os quais mostraram que maximizar a log-verossimilhança de Whittle leva a um estimador consistente com a distribuição normal assintótica usual que seria obtida pelo tratamento da verossimilhança de Whittle como uma log-verossimilhança convencional, ver também Dahlhaus (1989), Robinson (1995) e Hurvich, Deo, and Brodsky (1998).

Infelizmente, as ordenadas do periodograma não são assintoticamente independentes (Hurvich and Beltrao 1993), embora uma quase verossimilhança na forma da aproximação Whittle funcione bem e tenha boas propriedades assintóticas.

Para ver como isso funcionaria no caso de memória puramente longa, escrevamos o espectro de memória longa como \[\begin{equation} f_X(\omega_k; d,\sigma^2_{_W}) \, = \, \sigma^2_{_W} g_k^{-d}, \end{equation}\] onde \[\begin{equation} g_k \, = \, 4\mbox{sen}^2(\pi\omega_k)\cdot \end{equation}\]

Então, diferenciando a log-verossimilhança, digamos \[\begin{equation} \ln \big(L(X;d,\sigma^2_{_W})\big) \, \approx \, -m\ln (\sigma^2_{_W})+d \sum_{k=1}^m \ln (g_k) -\dfrac{1}{\sigma^2_{_W}}\sum_{k=1}^m g^d_k I(\omega_k) \end{equation}\] em \(m = n/2-1\) frequências e resolvendo por \(\sigma^2_{_W}\) produz \[\begin{equation} \sigma^2_{_W}(d) \, = \, \dfrac{1}{m}\sum_{k=1}^m g^d_k I(\omega_k) \end{equation}\] como o estimador aproximado do máximo da verossimilhança para o parámetro de variância.

Para estimar \(d\), também podemos usar uma pesquisa em grade da log-verossimilhançaa perfilada \[\begin{equation} \ln \big( L(X;d)\big) \approx -m\ln\big(\sigma^2_{_W}(d)\big) +d\sum_{k=1}^m \ln(g_k)-m, \end{equation}\] no intervalo \((0;0.5)\), seguido de um procedimento de Newton–Raphson para convergência.

Exemplo 2: Espectros de memória longa para a série de variações glaciais.

No Exemplo 1, ajustamos um modelo de memória longa aos dados da variável de variações glaciais através de métodos no domínio do tempo. Ajustar o mesmo modelo usando métodos no domínio da frequência e a aproximação Whittle acima fornece \(\widehat{d}=0.380\), com um erro padrão estimado de \(0.028\).

O método anterior no domínio do tempo obtemos \(\widehat{d}=0.384\) com \(M = 30\) e \(\widehat{d}=0.370\) com \(M = 100\). Ambas as estimativas, obtidas através de métodos no domínio do tempo, apresentaram um erro padrão de cerca de \(4.6\times 10^{-6}\), o que parece implausível.

A estimativa da variância do erro é de \(\widehat{\sigma}_{_W}^2=0.2293\). No Exemplo 1, poderiamos ter usado var(res.fd) como uma estimativa; nesse caso, obtemos \(0.2298097\).

series = log(varve) # especificando a série a ser analisada
d0 = .1             # valor inicial de d
n.per = nextn(length(series))
m = (n.per)/2 - 1
per = Mod(fft(series-mean(series))[-1])^2 # removendo frequências zero
per = per/n.per                           # peridograma esccalado
g = 4*(sin(pi*((1:m)/n.per))^2)
# Função para calcular a log-verossimilhança
whit.like = function(d){
   g.d=g^d
   sig2 = (sum(g.d*per[1:m])/m)
   log.like = m*log(sig2) - d*sum(log(g)) + m
   return(log.like) }
# Estimação 
(est = optim(d0, whit.like, gr=NULL, method="L-BFGS-B", hessian=TRUE,
              lower=-.5, upper=.5, control=list(trace=1,REPORT=1)))

## iter    1 value -145.635158
## iter    2 value -150.428661
## iter    3 value -152.601414
## iter    4 value -152.942481
## iter    5 value -152.952753
## iter    6 value -152.952788
## iter    7 value -152.952788
## final  value -152.952788 
## converged

## $par
## [1] 0.3802707
## 
## $value
## [1] -152.9528
## 
## $counts
## function gradient 
##        8        8 
## 
## $convergence
## [1] 0
## 
## $message
## [1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"
## 
## $hessian
##          [,1]
## [1,] 1232.464

##-- Resultados: d.hat = .380, se(dhat) = .028, e sig2hat = .229 --##
cat("d.hat =", est$par, "se(dhat) = ",1/sqrt(est$hessian),"\n")

## d.hat = 0.3802707 se(dhat) =  0.02848478

g.dhat = g^est$par; sig2 = sum(g.dhat*per[1:m])/m
cat("sig2hat =",sig2,"\n")

## sig2hat = 0.2293286

sig2hat = 0.2293286 

Também pode-se considerar ajustar um modelo autoregressivo a esses dados. Seguindo essa abordagem, obteve-se um modelo auto-regressivo com \(p = 8\) e \[\begin{equation} \widehat{\phi}_{1:8}=\{0.34, 0.11, 0.04, 0.09, 0.08, 0.08, 0.02, 0.09\}, \end{equation}\] com \(\widehat{\sigma}_{_W}^2=0.23\) como a variância do erro.

Os dois log-espectros são mostrados na figura abaixo para \(\omega> 0\) e observamos que o espectro de memória longa acabará se tornando infinito, enquanto o espectro AR(8) á finito em \(\omega=0\).

u = spec.ar(log(varve), plot=FALSE) # produz AR(8)
g = 4*(sin(pi*((1:500)/2000))^2)
fhat = sig2*g^{-est$par} # estimador espectral de memória longa
ar.mle(log(varve)) # to get AR(8) estimates

## 
## Call:
## ar.mle(x = log(varve))
## 
## Coefficients:
##      1       2       3       4       5       6       7       8  
## 0.3408  0.1112  0.0398  0.0876  0.0775  0.0845  0.0226  0.0880  
## 
## Order selected 8  sigma^2 estimated as  0.2267

par(mfrow = c(1,1), mar=c(3,3,1,1), mgp=c(1.6,.6,0), pch=19)
plot(1:500/2000, log(fhat), type="l", ylab="log(spectrum)", xlab="frequency")
lines(u$freq[1:250], log(u$spec[1:250]), lty="dashed")
grid()

Figura 4: Estimadores espectrais de memória longa \(d =0.380\), linha sólida e AR(8), linha tracejada, para a série varve, série das variações glaciais paleoclimáticas.

Frequentemente, as séries temporais não têm uma memória puramente longa. Uma situação comum tem o componente de memória longa multiplicado por um componente de memória curta, levando a uma versão alternativa da forma \[\begin{equation} f_X(\omega; d,\theta) \, = \, g_k^{-d} f_0(\omega_k;\theta), \end{equation}\] onde \(f_0(\omega_k;\theta)\) pode ser o espectro de um processo autorregressivo de médias móveis com vetor de parâmetros \(\theta\)) ou pode não ser especificado. Se o espectro tiver uma forma paramétrica, a verossimilhança de Whittle poderá ser usada.

Entretanto, existe uma quantidade substancial de literatura semiparamétrica que desenvolve os estimadores quando o espectro subjacente \(f_0(\omega;\theta)\) é desconhecido. Uma classe de estimadores semi-paramétricos gaussianos simplesmente usa a mesma verossimilhança de Whittle avaliada em uma sub-banda de baixas frequências, digamos \(m'=\sqrt{n}\). Existe alguma latitude na seleção de uma banda relativamente livre de interferências de baixa frequência devido ao componente de memória curta.

Se o espectro for altamente parametrizado, pode-se estimar usando a log-verossimilhança de Whittle e estimar em conjunto os parámetros \(d\) e \(\theta\) usando o método de Newton-Raphson. Se estivermos interessados em estimares não paramétricos, usando o estimador espectral suavizado convencional para o periodograma ajustado para o componente de memória longa, digamos que \(g_k^d I(\omega_k)\) pode ser uma abordagem possível.

Geweke and Porter-Hudak (1983) desenvolveram um método aproximado para estimar \(d\) com base em um modelo de regressão. Observe que podemos escrever uma equação simples para o logaritmo do espectro como \[\begin{equation} \ln \big(f_X(\omega_k;d) \big) \, = \, \ln \big(f_0(\omega_k;\theta) \big)-d \ln \Big( 4\mbox{sen}^2(\pi\omega_k)\Big), \end{equation}\] com as frequências \(\omega_k = k/n\) restritas a uma faixa \(k = 1,2,\cdots,m\) perto da frequência zero com \(m =\sqrt{n}\) como o valor recomendado.

A relação acima sugere o uso de um modelo de regressão linear simples da forma, \[\begin{equation} \ln \big( I(\omega_k)\big) \, = \, \beta_0 -d \ln\big(4\mbox{sen}^2(\pi\omega_k) \big)+\epsilon_k, \end{equation}\] para o periodograma para estimar os parámetros \(\sigma^2_{_W}\) e \(d\). Nesse caso, executa-se mínimos quadrados usando \(\ln \big( I(\omega_k)\big)\) como variável dependente e \(\ln\big(4\mbox{sen}^2(\pi\omega_k) \big)\) como variável independente para \(k = 1,2,\cdots,m\). A estimativa da inclinação resultante é então usada como uma estimativa de \(-d\). Para uma boa discussão de vários métodos alternativos para selecionar \(m\), (veja Hurvich and Deo 1999).

O pacote R fracdiff também fornece esse método através do comando fdGPH(); consulte o arquivo de ajuda para obter mais informações.

fdGPH(log(varve), bandw=.9) # m = n^bandw dhat = 0.383 se(dhat) = 0.041

## $d
## [1] 0.3833677
## 
## $sd.as
## [1] 0.04085055
## 
## $sd.reg
## [1] 0.04163713

2. Teste da raiz unitária

Conforme discutido na seção anterior, o uso da primeira diferença \(\nabla X_t =(1-B)X_t\) pode ser uma modificação muito severa, no sentido de que o modelo não-estacionário pode representar uma superdiferença do processo original. Por exemplo, considere um processo causal de AR(1), assumimos ao longo desta seção que o ruído é gaussiano, \[\begin{equation} X_t \, = \, \phi X_{t-1}+W_t\cdot \end{equation}\]

A aplicação de \((1-B)\) em ambos os lados mostra que a diferença, \(\nabla X_t = \phi\nabla X_{t-1}+\nabla W_t\) ou \[\begin{equation} Y_t \, = \, \phi Y_{t-1}+W_t-W_{t-1}, \end{equation}\] onde \(Y_t=\nabla X_t\), introduz problemas estranhos de correlação e invertibilidade. Significa que, enquanto \(X_t\) é um processo causal AR(1), trabalhar com o processo diferenciado \(Y_t\) será problemático porque é um processo ARMA(1,1) não invertível.

O teste da raiz unitária fornece uma maneira de testar se \(X_t = \phi X_{t-1}+W_t\) é uma caminhada aleatória, o caso nulo, em oposição a um processo causal, a alternativa. Ou seja, fornece um procedimento para testar \[\begin{equation} H_0 \, : \, \phi\, = \, 1 \qquad \mbox{contra} \qquad H_1 \, : \, |\phi| < 1\cdot \end{equation}\]

Uma estatística de teste óbvia seria considerar \(\widehat{\phi}-1\), adequadamente normalizado, na esperança de desenvolver uma estatística de teste assintoticamente normal, sendo \(\widehat{\phi}\) um dos estimadores ideais de \(\phi\). Infelizmente, a teoria não funciona no caso nulo porque o processo não é estacionário. Além disso, a estimativa próxima ao limite da estacionariedade produz distribuições amostrais altamente distorcidas e isso é uma boa indicação de que o problema será atípico.

Para examinar o comportamento de \(\widehat{\phi}-1\) sob a hipótese nula de que \(\phi= 1\) ou mais precisamente, que o modelo é um passeio aleatório, \(\displaystyle X_t = \sum_{j=1}^t W_j\) ou \(X_t = X_{t-1}+W_t\) com \(X_0 = 0\), considere o estimador de mínimos quadrados de \(\phi\). Observando que \(\mu_X = 0\), o estimador de mínimos quadrados pode ser escrito como \[\begin{equation} \widehat{\phi} \, = \, \frac{\displaystyle \sum_{t=1}^n X_t X_{t-1}}{\displaystyle \sum_{t=1}^n X_{t-1}^2} \, = \, 1+ \frac{\displaystyle \frac{1}{n}\sum_{t=1}^n W_t X_{t-1}}{\displaystyle \frac{1}{n}\sum_{t=1}^n X_{t-1}^2}, \end{equation}\] onde escrevemos \(X_t=X_{t-1}+W_t\) no numerador: lembre-se de que \(X_0 = 0\) e na configuração de mínimos quadrados, estamos regredindo \(X_t\) em \(X_{t-1}\) para \(t = 1,,\cdots,n\). Portanto, sob \(H_0\) verdadeira, temos que \[\begin{equation} \widehat{\phi} -1 \, = \, \frac{\displaystyle \frac{1}{n\sigma_{_W}^2}\sum_{t=1}^n W_t X_{t-1}}{\displaystyle \frac{1}{n\sigma_{_W}^2}\sum_{t=1}^n X_{t-1}^2}\cdot \end{equation}\]

Considere o numerador da expressão acima de \(\widehat{\phi}\). Observe primeiro que elevando ao quadrado os dos dois lados da \(X_t=X_{t-1}+W_t\), obtemos \(X_t^2=X_{t-1}^2+2X_{t-1}W_t+W_t^2\), de modo que \[\begin{equation} X_{t-1}W_t \, = \, \frac{1}{2}\big( X_t^2-X_{t-1}^2-W_t^2\big), \end{equation}\] e somando \[\begin{equation} \frac{1}{n\sigma_{_W}^2}\sum_{t=1}^n X_{t-1}W_t \, = \, \displaystyle \frac{1}{2}\left( \frac{X_n^2}{n\sigma_{_W}^2}-\frac{\sum_{t=1}^n W_t^2}{n\sigma_{_W}^2}\right)\cdot \end{equation}\] Porque \(\displaystyle X_n=\sum_{t=1}^n W_t\), temos \(X_n\sim N(0,n\sigma_{_W}^2)\), de maneira que \[\begin{equation} \frac{1}{n\sigma_{_W}^2} X_n \sim \chi^2_1, \end{equation}\] tem distribuição qui-quadrado com um grau de liberdade. Além disso, porque \(W_t\) é ruído branco gaussiano \[\begin{equation} \frac{1}{n}\sum_{t=1}^n W_t^2 \overset{P}{\longrightarrow} \sigma_{_W}^2 \qquad \mbox{ou} \qquad \frac{1}{n\sigma_{_W}^2}\sum_{t=1}^n W_t^2 \overset{P}{\longrightarrow} 1\cdot \end{equation}\] Consequentemente, quando \(n\to\infty\) \[\begin{equation} \frac{1}{n\sigma_{_W}^2}\sum_{t=1}^n X_{t-1} W_t \overset{D}{\longrightarrow} \displaystyle \frac{1}{2}\big( \chi^2_1-1\big)\cdot \end{equation}\]

A seguir, focaremos no denominador da expressão acima de \(\widehat{\phi}\). Primeiro, apresentamos o movimento browniano padrão.

Definição 2. Movimento Browniano.

Um processo estocástico de tempo contínuo \(\{W(t) \, : \, t\geq 0\}\) é chamado de movimento Browniano padrão se satisfizer as seguintes condições:

1. \(W(0)=0\);

2. \(\{ W(t_2)-W(t_1), \, W(t_3)-W(t_2), \, \cdots, \, W(t_n)-W(t_{n-1})\}\) forem independentes para qualquer coleção de pontos, \(0\leq t_1 < t_2 < \cdots < t_n\) e inteiros \(n<2\);

3. \(W(t+\Delta t)-W(t)\sim N(0,\Delta t)\) para \(\Delta t>0\).

Além de (i) - (iii) acima, supõe-se que quase todos os caminhos amostrais de \(W(t)\) sejam contínuos em \(t\). O resultado para o denominador usa o Teorema do Limite Central Funcional, que pode ser encontrado em Billingsley (1999). Em particular, se \(\xi_1,\cdots,\xi_n\) for uma sequência de variáveis aleatórias independentes idênticamente distribuídas com média 0 e variância 1; então, para \(0\leq t\leq 1\), o processo contínuo no tempo \[\begin{equation} S_n(t) \, = \, \frac{1}{\sqrt{n}}\sum_{j=1}^{[[nt]]} \xi_j \overset{D}{\longrightarrow} W(t), \end{equation}\] quando \(n\to\infty\), onde \([[\cdot]]\) é o maior inteiro e \(W(t)\) é um movimento Browniano padrão em \([0,1]\).

A intuição aqui é para \(k=[[nt]]\) e fixado \(t\), o Torema do Limite Central implica que \[\begin{equation} \sqrt{t}\frac{1}{\sqrt{k}}\sum_{j=1}^k \xi_j \sim N(0,t), \qquad \mbox{quando} \qquad n\to\infty\cdot \end{equation}\]

Observe que o sob a hipótese nula \(X_s=W_1+\cdots+W_s\sim N(0,s\sigma_{_W}^2)\) e, pelo resultado acima, temos que \[\begin{equation} \frac{X_s}{\sigma_{_W} \sqrt{n}} \overset{D}{\longrightarrow} W(s)\cdot \end{equation}\] Deste fato, podemos mostrar que, quando \(n\to\infty\), \[\begin{equation} \frac{1}{n}\sum_{t=1}^n \left( \frac{X_{t-1}}{\sigma_{_W} \sqrt{n}}\right)^2 \, \overset{D}{\longrightarrow} \, \displaystyle \int_0^1 W^2(t)\mbox{d}t\cdot \end{equation}\]

Obtemos então que, quando \(n\to\infty\), \[\begin{equation} n(\widehat{\phi} -1) \, = \, \frac{\displaystyle \frac{1}{n\sigma_{_W}^2}\sum_{t=1}^n W_t X_{t-1}}{\displaystyle \frac{1}{n^2\sigma_{_W}^2}\sum_{t=1}^n X_{t-1}^2} \, \overset{D}{\longrightarrow} \, \frac{\displaystyle \frac{1}{2}\big( \chi_1^2 -1\big)}{\displaystyle \int_0^1 W^2(t)\mbox{d}t}\cdot \end{equation}\]

A estatística \(n(\widehat{\phi}-1)\) é conhecida como raiz unitária ou estatística Dickey-Fuller (DF) (Dickey and Fuller 1979), embora a estatística real do teste DF seja normalizada de forma um pouco diferente. Derivações relacionadas foram discutidas em Rao (n.d.) e em Evans and Savin (1981). Como a distribuição da estatística de teste não tem uma forma fechada, os quantis da distribuição devem ser calculados por aproximação numérica ou por simulação. O pacote R tseries fornecem este teste junto com testes mais gerais que mencionamos brevemente.

Para um modelo mais geral, notamos que o teste DF foi estabelecido observando que se \(X_t=\phi X_{t-1}+W_t\), então \(\nabla X_t=(\phi-1)X_{t-1}+W_t=\gamma X_{t-1}+W_t\) e poderiamos testar \(H_):\gamma=0\) regredindo \(\nabla X_t\) em \(X_{t-1}\). Formando assim uma estatística de Wald e derivar sua distribuição limite, a derivação anterior baseada no movimento browniano é devida a Phillips (1987).

O teste foi estendido para acomodar modelos AR(p), \(X_t=\displaystyle \sum_{j=1}^p \phi_j X_{t-j}+W_t\), como segue. Subtraindo \(X_{t-1}\) de ambos os lados, se obtém \[\begin{equation} \nabla X_t \, = \, \gamma X_{t-1} +\sum_{j=1}^{p-1} \psi_j \nabla X_{t-j} +W_t, \end{equation}\] onde \(\displaystyle \gamma = \sum_{j=1}^p \phi_j-1\) e \(\displaystyle \psi_j=-\sum_{i=j}^p \phi_i\), pata \(j=2,\cdots,p\).

Para uma verificação rápida da expressão acima, quando \(p = 2\), observe que \[\begin{equation} X_t \, = \, (\phi_1+\phi_2)X_{t-1} - \phi_2(X_{t-1}-X_{t-2})+W_t, \end{equation}\] subtraindo \(X_{t-1}\) de ambos os lados. Para testar a hipótese de que o processo tem uma raiz unitária em 1, ou seja, o AR polinomial \(\phi(z)=0\) quando \(z = 1\), podemos testar \(H_0: \gamma = 0\) estimando \(\gamma\) na regressão de \(\nabla X_t\) em \(X_{t-1},\nabla X_{t-1},\cdots,\nabla X_{t-p+1}\) e formando um teste de Wald baseado em \[\begin{equation} t_{_\gamma} \, = \, \frac{\widehat{\gamma}}{\sqrt{\mbox{Var}(\widehat{\gamma})}}\cdot \end{equation}\]

Este teste leva ao chamado teste aumentado de Dickey-Fuller (ADF). Enquanto os cálculos para obter a distribuição nula assintótica mudam, as ideias básicas e o maquinário permanecem os mesmos que no caso simples. A escolha de \(p\) é crucial e discutiremos algumas sugestões no exemplo. Para modelos ARMA(p,q), o teste ADF pode ser usado assumindo que \(p\) é grande o suficiente para capturar a estrutura de correlação essencial. Outra alternativa é o teste de Phillips-Perron (PP), que difere dos testes do ADF principalmente no modo como eles lidam com a correlação serial e a heteroscedasticidade nos erros.

Pode-se estender o modelo para incluir uma tendência constante ou mesmo não-estocástica. Por exemplo, considere o modelo \[\begin{equation} X_t \, = \, \beta_0+\beta_1 t+\phi X_{t-1}+W_t\cdot \end{equation}\] Se assumirmos \(\beta_1 = 0\), então, sob a hipótese nula, \(\phi= 1\), o processo é um passeio aleatório com tendência constante \(\beta_0\). Sob a hipótese alternativa, o processo é um AR(1) causal, com média \(\mu_X = \beta_0(1-\phi)\). Se não podemos assumir \(\beta_1 = 0\), então o interesse aqui é testar a hipótese nula que \((\beta_1,\phi)=(0,1)\), simultaneamente, versus a alternativa que \(\beta_1\neq 0\) e \(|\phi|< 1\). Neste caso, a hipótese nula seria que o processo é um passeio aleatório com tendência, versus a hipótese alternativa de que o processo é de tendência estacionária.

Exemplo 3: Teste da raíz unitária na série de variações glaciais.

Neste exemplo, usamos o pacote R tseries para testar a hipótese nula de que o logaritmo da série de variações glaciais tem uma raiz unitária, versus a hipótese alternativa de que o processo é estacionário.

Testamos a hipótese nula usando os testes DF, ADF e PP disponíveis. Observe que, em cada caso, a equação de regressão geral incorpora uma tendência constante e linear. No teste ADF, o número padrão de componentes AR incluídos no modelo, digamos \(k\), é \([[(n-1)^\frac{1}{3}]]\), que corresponde ao limite superior sugerido na taxa na qual o número de lags, \(k\), deve ser feito para crescer com o tamanho da amostra para a configuração geral ARMA(p,q). Para o teste PP, o valor padrão de \(k\) é \([[0.04n^\frac{1}{4}]]\).

library(tseries)
adf.test(log(varve), k=0) # DF test

## 
##  Augmented Dickey-Fuller Test
## 
## data:  log(varve)
## Dickey-Fuller = -12.857, Lag order = 0, p-value = 0.01
## alternative hypothesis: stationary

adf.test(log(varve)) # ADF test

## 
##  Augmented Dickey-Fuller Test
## 
## data:  log(varve)
## Dickey-Fuller = -3.5166, Lag order = 8, p-value = 0.04071
## alternative hypothesis: stationary

pp.test(log(varve)) # PP test

## 
##  Phillips-Perron Unit Root Test
## 
## data:  log(varve)
## Dickey-Fuller Z(alpha) = -304.54, Truncation lag parameter = 6, p-value
## = 0.01
## alternative hypothesis: stationary

Cada teste acima rejeita a hipótese nula de que a série de variações glaciais possui uma raiz unitária. A conclusão desses testes suporta a conclusão da seção anterior de que o logaritmo desta série é de memória longa ao invês de integrada.

3. Modelos GARCH

Vários problemas, como a precificação de opções em finanças, motivaram o estudo da volatilidade ou variabilidade de uma série temporal. Os modelos ARMA foram usados para modelar a média condicional de um processo quando a variância condicional era constante. Usando um AR(1) como exemplo, assumimos \[\begin{equation} \mbox{E}(X_t \, | \, X_{t-1},X_{t-2},\cdots) \, = \, \phi X_{t-1} \qquad \mbox{Var}(X_t \, | \, X_{t-1},X_{t-2},\cdots) \, = \, \sigma_{_W}^2\cdot \end{equation}\]

Em muitos problemas, no entanto, a suposição de uma variância condicional constante será violada. Modelos como o auto-regressivo condicionalmente heterocedástico ou o modelo ARCH, introduzidos pela primeira vez por Engle (1982), foram desenvolvidos para modelar mudanças na volatilidade. Esses modelos foram posteriormente estendidos aos modelos ARCH ou GARCH generalizados por Bollerslev (1986).

Nesses problemas, estamos preocupados em modelar a taxa de retorno ou crescimento de uma série. Por exemplo, se \(X_t\) for o valor de um ativo no momento \(t\), então o retorno ou ganho relativo, \(r_t\), do ativo no momento \(t\) será \[\begin{equation} r_t \, = \, \frac{X_t-X_{t-1}}{X_{t-1}}\cdot \end{equation}\] Esta definição implica que \(X_t = (1-r_t)X_{t-1}\). Assim, se o retorno representar uma pequena porcentagem, em magnitude, de mudança percentual então \[\begin{equation} \nabla \log(X_t)\approx r_t\cdot \end{equation}\]

Qualquer valor, \(\nabla \log(X_t)\) ou \((X_t-X_{t-1})/X_{t-1}\), será chamado de retorno e indicado por \(r_t\). Lembre-se que, se \(r_t = (X_t-X_{t-1}/X_{t-1}\) for uma pequena porcentagem, então \(\log(1 + r_t)\approx r_t\). É mais fácil programar \(\nabla \log(X_t)\), portanto, isso geralmente será usado em vez de calcular \(r_t\) diretamente. Embora seja um nome impróprio, \(\nabla \log(X_t)\) é frequentemente chamado de log-retorno; mas os retornos não estão sendo registrados. Uma alternativa ao modelo GARCH é o modelo de volatilidade estocástica; porque eles são modelos de espaço de estados.

Normalmente, para séries financeiras, o retorno \(r_t\), não tem uma variância condicional constante e períodos altamente voláteis tendem a ser agrupados. Em outras palavras, há uma forte dependência de explosões súbitas de variabilidade em um retorno sobre o passado da própria série. Por exemplo, a figura abaixo mostra os retornos diários do Dow Jones Industrial Average (DJIA) de 20 de abril de 2006 a 20 de abril de 2016. Nesse caso, como é típico, o retorno \(r_t\) é razoavelmente estável, exceto para rajadas de curto prazo de alta volatilidade.

par(mar=c(4,4,4,1))
plot(nyse, xlab="Tempo", ylab="Retorno da Bolsa de Valores de Nova Iorque (NYSE)",
     main="Retornos de mercado ponderados pelo valor diário\n 
     entre 2 de fevereiro de 1984 e 31 de dezembro de 1991")
grid()

O modelo ARCH mais simples o ARCH(1), modela o retorno como \[\begin{equation} r_t \, = \, \sigma_t \epsilon_t \qquad \mbox{e} \qquad \sigma_t^2 \, = \, \alpha_0 +\alpha_1 r_{t-1}^2, \end{equation}\] onde \(\epsilon_t\) é o ruído branco gaussiano padrão, \(\epsilon_t \sim N(0,1)\) independentes. A suposição de normalidade pode ser relaxada. Vamos discutir isso mais tarde. Assim como nos modelos ARMA, devemos impor algumas restrições aos parámetros do modelo para obter propriedades desejáveis. Uma restrição óbvia é que \(\alpha_0,\alpha_1\geq 0\), porque \(\sigma_t^2\) é a variância.

Como veremos, os modelos ARCH(1) retornam como um processo de ruído branco com variância condicional não constante e essa variância condicional depende do retorno anterior. Primeiro, observe que a distribuição condicional de \(r_t\) dado \(r_{t-1}\) é gaussiana: \[\begin{equation} r_t \,| \, r_{t-1} \, \sim \, N\big( 0,\alpha_0+\alpha_1 r_{t-1}^2\big)\cdot \end{equation}\]

Além disso, é possível escrever o modelo ARCH(1) como um modelo AR(1) não-Gaussiano no quadrado dos retornos \(r_t^2\). Para isso, reescrevemos \[\begin{equation} r_t^2 \, = \, \sigma_t^2 \epsilon_t^2 \qquad \mbox{e} \qquad \alpha_0+\alpha_1 r_{t-1}^2 \, = \, \sigma_t^2, \end{equation}\] e subtraindo as duas equações obtemos \[\begin{equation} r_t^2 \, - \, (\alpha_0+\alpha_1 r_{t-1}^2) \, = \, \sigma_t^2 \epsilon_t^2 -\sigma_t^2\cdot \end{equation}\] Agora, escreva esta equação como \[\begin{equation} r_t^2 \, = \, \alpha_0+\alpha_1 r_{t-1}^2+\nu_t, \end{equation}\] onde \(\nu_t=\sigma_t^2(\epsilon_t^2-1)\). Devido a \(\epsilon_t^2\) ser o quadrado de uma variável aleatória \(N(0,1)\), \(\epsilon_t^2-1\) é uma variável aleatória \(\chi_1^2\) deslocada, para ter média zero.

Para explorar as propriedades do mosdelo ARCH, definimos \(\mathcal{R}_s = \{r_s, r_{s-1},\cdots \}\). Então, usando que \(r_t | r_{t-1} \sim N\big( 0,\alpha_0+\alpha_1 r_{t-1}^2\big)\), vemos imediatamente que \(r_t\) tem média zero: \[\begin{equation} \mbox{E}(r_t) \, = \, \mbox{E}\big(\mbox{E}(r_t \, | \, \mathcal{R}_{t-1})\big) \, = \, \mbox{E}\big(\mbox{E}(r_t \, | \, r_{t-1})\big) \, = \, 0\cdot \end{equation}\] Porque \(\mbox{E}\big(\mbox{E}(r_t|\mathcal{R}_{t-1})\big) = 0\), o processo \(r_t\) é dito ser uma diferença martingal.

Por ser \(r_t\) uma diferença martingal, também é uma sequência não correlacionada. Por exemplo, com \(h> 0\), \[\begin{array}{rcl} \mbox{Cov}(r_{t+h}, r_t) & = & \mbox{E}(r_t \, r_{t+h}) \, = \, \mbox{E}\big(\mbox{E}(r_t \, r_{t+h}|\mathcal{R}_{t+h-1})\big) \\ & = & \mbox{E}\big(r_t \, \mbox{E}(r_{t+h} \, | \, \mathcal{R}_{t+h-1})\big) \, = \, 0\cdot \end{array}\]

A última linha acima segue porque \(r_t\) pertence ao conjunto de informações \(R_{t+h-1}\) para \(h> 0\) e \[\begin{equation} \mbox{E}(r_{t+h} | \mathcal{R}_{t+h-1})=0\cdot \end{equation}\]

Um argumento similar estabelecerá o fato de que o processo de erro \(\nu_t\) é também uma diferença martingal e, consequentemente, uma sequência não correlacionada. Se a variância de \(\nu_t\) for finita e constante em relação ao tempo e \(0\leq \alpha_1 < 1\), então com base no Teorema 1, \(r_t^2 = \alpha_0+\alpha_1 r_{t-1}^2+\nu_t\) especifica um processo causal \(AR(1)\) para \(r^2_t\). Portanto, \(\mbox{E}(r^2_t)\) e \(\mbox{Var}(r^2_t)\) devem ser constantes em relação ao tempo \(t\). Isso implica que \[\begin{equation} \mbox{E}(r_t^2) \, = \, \mbox{Var}(r_t) \, = \, \frac{\alpha_0}{1-\alpha_1} \end{equation}\] e, depois de algumas manipulações, \[\begin{equation} \mbox{E}(r_t^4) \, = \, \frac{3\alpha_0^2}{(1-\alpha_1)^2}\frac{1-\alpha_1^2}{1-3\alpha_1^2}, \end{equation}\] desde que \(3\alpha_1^2 < 1\). Observe que \[\begin{equation} \mbox{Var}(r_t^2) \, = \, \mbox{E}(r_t^4) \, - \, \big(\mbox{E}(r_t^2)\big)^2, \end{equation}\] que existe apenas se \(0 < \alpha_1 < 1/\sqrt{3}\approx 0.58\). Além disso, estes resultados implicam que a curtose, \(\kappa\) de \(r_t\) é \[\begin{equation} \kappa \, = \, \frac{\mbox{E}(r_t^4)}{\big(\mbox{E}(r_t^2)\big)^2} \, = \, 3\frac{1-\alpha_1^2}{1-3\alpha_1^2}, \end{equation}\] a qual nunca é menor que 3, a curtose da distribuição normal.

Assim, a distribuição marginal dos retornos \(r_t\), é leptocúrtica ou tem “caudas grossas”. Resumindo, se \(0\leq \alpha_1 < 1\), o processo \(r_t\) em si é ruído branco e sua distribuição não condicional é distribuída simetricamente em torno de zero; essa distribuição é leptocúrtica. Se, além disso, \(3\alpha_1^2 < 1\), o quadrado do processo, \(r^2_t\), segue um modelo causal AR(1) com ACF dado por \(\rho_{_{r_t^2}}(h)=\alpha_1^h\geq 0\), para todo \(h>0\). Se \(3\alpha_1^2 \geq 1\), mas \(\alpha_1 < 1\), pode-se mostrar que \(r^2_t\) é estritamente estacionário com variância infinita (ver Douc, Moulines, and Stoffer 2014).

A estimação dos parámetros \(\alpha_0\) e \(\alpha_1\) do modelo ARCH(1) é tipicamente realizada por máxima verossimilhança condicional. A verossimilhança condicional dos dados \(r_2,\cdots,r_n\) dado \(r_1\), é dada por \[\begin{equation} L(\alpha_0,\alpha_1|r_1) \, = \, \prod_{t=2}^n f_{\alpha_0,\alpha_1}(r_t|r_{t-1}), \end{equation}\] onde a densidade \(f_{\alpha_0,\alpha_1}(r_t|r_{t-1})\) é a densidade normal \(N\big( 0,\alpha_0+\alpha_1 r_{t-1}^2\big)\). Portanto, a função critério a ser minimizada, \(\ell(\alpha_0,\alpha_1)=-\ln\big(L(\alpha_0,\alpha_1|r_1)\big)\) é dada por \[\begin{equation} \ell(\alpha_0,\alpha_1) \, = \, \displaystyle \frac{1}{2}\sum_{t=2}^n \ln \big( \alpha_0+\alpha_1 r_{t-1}^2\big) \, + \, \frac{1}{2}\sum_{t=2}^n \left( \frac{r_t^2}{\alpha_0 +\alpha_1 r_{t-1}^2}\right)\cdot \end{equation}\]

A estimação é realizada por métodos numéricos, nesse caso, expressões analíticas para o vetor gradiente, \(\ell^{(1)}(\alpha_0,\alpha_1)\) e a matriz Hessiana, \(\ell^{(2)}(\alpha_0,\alpha_1)\), podem ser obtidas por cálculos diretos. Por exemplo, o vetor gradiente \(2\times 1\), \(\ell^{(1)}(\alpha_0,\alpha_1)\), é dado por \[\begin{equation} \begin{pmatrix} \displaystyle \partial \ell\big/ \partial\alpha_0 \\ \displaystyle \partial \ell\big/ \partial\alpha_1 \end{pmatrix} \, = \, \sum_{t=2}^n \begin{pmatrix} 1 \\ r_{t-1}^2 \end{pmatrix} \times \frac{\displaystyle \alpha_0+\alpha_1 r_{t-1}^2-r_t^2} {\displaystyle \big( \alpha_0+\alpha_1 r_{t-1}^2\big)^2}\cdot \end{equation}\]

O cálculo da matriz Hessiana é deixado como um exercício. A verossimilhança do modelo ARCH tende a ser plana, a menos que \(n\) seja muito grande. Uma discussão sobre esse problema pode ser encontrada em Shephard (1996).

Também é possível combinar uma regressão ou um modelo ARMA para a média com um modelo ARCH para os erros. Por exemplo, uma regressão com o modelo de erros ARCH(1) teria as observações \(X_t\) como função linear dos \(p\) regressores, \(z_p = (z_{t1},\cdots,z_{tp})^\top\) e o ruído \(Y_t\) ARCH(1), digamos, \[\begin{equation} X_t \, = \, \beta^\top z_t + Y_t, \end{equation}\] onde \(Y_t\) satisfaz \(r_t = \sigma_t\epsilon_t\) e \(\sigma_t^2 = \alpha_0 +\alpha_1 r_{t-1}^2\) mas, neste caso, não é observado. Da mesma forma, por exemplo, um modelo AR(1) para os dados \(X_t\) exibindo erros ARCH(1) seria \[\begin{equation} X_t \, = \, \phi_0 +\phi_1 X_{t-1}+Y_t\cdot \end{equation}\] Esses tipos de modelos foram explorados por Weiss (1984).

Exemplo 4: Análise do PIB dos EUA.

No Exemplo 39, na Seção Modelos ARIMA, ajustamos um modelo MA(2) e um modelo AR(1) à série do PIB dos EUA e concluímos que os resíduos de ambos os ajustes pareciam se comportar como um processo de ruído branco. No Exemplo 43, concluímos que o AR(1) é provavelmente o melhor modelo nesse caso.

Foi sugerido que a série do PIB dos EUA tem erros ARCH e, neste exemplo, investigaremos essa reivindicação. Se o termo de ruído do PIB for ARCH, os quadrados dos resíduos do ajuste devem se comportar como um processo AR(1) não gaussiano.

library(astsa)
u = sarima(diff(log(gnp)), 1, 0, 0)

## initial  value -4.589567 
## iter   2 value -4.654150
## iter   3 value -4.654150
## iter   4 value -4.654151
## iter   4 value -4.654151
## iter   4 value -4.654151
## final  value -4.654151 
## converged
## initial  value -4.655919 
## iter   2 value -4.655921
## iter   3 value -4.655922
## iter   4 value -4.655922
## iter   5 value -4.655922
## iter   5 value -4.655922
## iter   5 value -4.655922
## final  value -4.655922 
## converged
## <><><><><><><><><><><><><><>
##  
## Coefficients: 
##       Estimate     SE t.value p.value
## ar1     0.3467 0.0627  5.5255       0
## xmean   0.0083 0.0010  8.5398       0
## 
## sigma^2 estimated as 9.029569e-05 on 220 degrees of freedom 
##  
## AIC = -6.44694  AICc = -6.446693  BIC = -6.400958 
## 

A figura abaixo à direita mostra o ACF e o PACF dos resíduos ao quadrado, parece que pode haver alguma dependência, embora pequena, nos resíduos.

acf2(resid(u$fit)^2, 20)

##      [,1] [,2] [,3] [,4]  [,5] [,6]  [,7] [,8] [,9] [,10] [,11] [,12] [,13]
## ACF  0.12 0.13 0.03 0.13  0.01 0.05 -0.03 0.06 0.08 -0.08  0.09  0.10  0.01
## PACF 0.12 0.12 0.00 0.12 -0.02 0.02 -0.04 0.05 0.08 -0.12  0.11  0.09 -0.05
##      [,14] [,15] [,16] [,17] [,18] [,19] [,20]
## ACF   0.04  0.15 -0.02  0.04 -0.05  0.01  0.05
## PACF  0.05  0.13 -0.09  0.01 -0.05  0.00  0.04

Figura 5: Funções de autocorrelação e autocorrelação parcial amostrais dos quadrados dos resíduos do ajuste AR(1) no PNB dos EUA.

Usamos o pacote R do fGarch para ajustar um modelo AR(1)-ARCH(1) aos retornos do PIB dos EUA com os seguintes resultados. Notamos que o modelo GARCH(1,0) especifica um ARCH(1) no código abaixo.

library(fGarch)
summary(garchFit(~arma(1,0)+garch(1,0), diff(log(gnp))))

## 
## Series Initialization:
##  ARMA Model:                arma
##  Formula Mean:              ~ arma(1, 0)
##  GARCH Model:               garch
##  Formula Variance:          ~ garch(1, 0)
##  ARMA Order:                1 0
##  Max ARMA Order:            1
##  GARCH Order:               1 0
##  Max GARCH Order:           1
##  Maximum Order:             1
##  Conditional Dist:          norm
##  h.start:                   2
##  llh.start:                 1
##  Length of Series:          222
##  Recursion Init:            mci
##  Series Scale:              0.01015924
## 
## Parameter Initialization:
##  Initial Parameters:          $params
##  Limits of Transformations:   $U, $V
##  Which Parameters are Fixed?  $includes
##  Parameter Matrix:
##                      U          V    params includes
##     mu     -8.20681904   8.206819 0.8205354     TRUE
##     ar1    -0.99999999   1.000000 0.3466459     TRUE
##     omega   0.00000100 100.000000 0.1000000     TRUE
##     alpha1  0.00000001   1.000000 0.1000000     TRUE
##     gamma1 -0.99999999   1.000000 0.1000000    FALSE
##     delta   0.00000000   2.000000 2.0000000    FALSE
##     skew    0.10000000  10.000000 1.0000000    FALSE
##     shape   1.00000000  10.000000 4.0000000    FALSE
##  Index List of Parameters to be Optimized:
##     mu    ar1  omega alpha1 
##      1      2      3      4 
##  Persistence:                  0.1 
## 
## 
## --- START OF TRACE ---
## Selected Algorithm: nlminb 
## 
## R coded nlminb Solver: 
## 
##   0:     682.89527: 0.820535 0.346646 0.100000 0.100000
##   1:     308.43148: 0.763492 0.258112  1.06104 0.352453
##   2:     306.07332: 0.681276 0.195897  1.04763 0.304072
##   3:     301.00807: 0.561958 0.448458 0.825277 0.0402737
##   4:     298.88361: 0.383716 0.465477 0.632947 0.385969
##   5:     296.74288: 0.504144 0.389445 0.683635 0.247795
##   6:     296.67703: 0.497724 0.366843 0.688130 0.229496
##   7:     296.60039: 0.500011 0.385702 0.703145 0.211105
##   8:     296.59692: 0.515645 0.374174 0.690079 0.194961
##   9:     296.56381: 0.513570 0.367018 0.702272 0.200013
##  10:     296.55723: 0.523440 0.363126 0.708406 0.194151
##  11:     296.55632: 0.522578 0.364913 0.710104 0.194839
##  12:     296.55598: 0.520871 0.364956 0.710924 0.193212
##  13:     296.55568: 0.519486 0.366571 0.710213 0.194510
##  14:     296.55568: 0.519509 0.366597 0.710266 0.194512
##  15:     296.55568: 0.519511 0.366586 0.710290 0.194452
##  16:     296.55568: 0.519505 0.366562 0.710299 0.194464
##  17:     296.55568: 0.519526 0.366560 0.710295 0.194472
##  18:     296.55568: 0.519522 0.366563 0.710295 0.194471
## 
## Final Estimate of the Negative LLH:
##  LLH:  -722.2849    norm LLH:  -3.253536 
##           mu          ar1        omega       alpha1 
## 0.0052779470 0.3665625611 0.0000733096 0.1944713364 
## 
## R-optimhess Difference Approximated Hessian Matrix:
##                  mu          ar1         omega        alpha1
## mu     -2749495.406 -24170.12492  4.546826e+06   -1586.69157
## ar1      -24170.125   -390.26682  1.253879e+04      -6.73379
## omega   4546825.892  12538.78540 -1.590043e+10 -706934.13922
## alpha1    -1586.692     -6.73379 -7.069341e+05    -142.53953
## attr(,"time")
## Time difference of 0.04847479 secs
## 
## --- END OF TRACE ---
## 
## 
## Time to Estimate Parameters:
##  Time difference of 0.2405105 secs
## 
## Title:
##  GARCH Modelling 
## 
## Call:
##  garchFit(formula = ~arma(1, 0) + garch(1, 0), data = diff(log(gnp))) 
## 
## Mean and Variance Equation:
##  data ~ arma(1, 0) + garch(1, 0)
## <environment: 0x610c13e12e70>
##  [data = diff(log(gnp))]
## 
## Conditional Distribution:
##  norm 
## 
## Coefficient(s):
##         mu         ar1       omega      alpha1  
## 0.00527795  0.36656256  0.00007331  0.19447134  
## 
## Std. Errors:
##  based on Hessian 
## 
## Error Analysis:
##         Estimate  Std. Error  t value Pr(>|t|)    
## mu     5.278e-03   8.996e-04    5.867 4.44e-09 ***
## ar1    3.666e-01   7.514e-02    4.878 1.07e-06 ***
## omega  7.331e-05   9.011e-06    8.135 4.44e-16 ***
## alpha1 1.945e-01   9.554e-02    2.035   0.0418 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log Likelihood:
##  722.2849    normalized:  3.253536 
## 
## Description:
##  Fri Oct  4 17:05:28 2024 by user:  
## 
## 
## Standardised Residuals Tests:
##                                  Statistic     p-Value
##  Jarque-Bera Test   R    Chi^2   9.1180361 0.010472337
##  Shapiro-Wilk Test  R    W       0.9842406 0.014336496
##  Ljung-Box Test     R    Q(10)   9.8743260 0.451587524
##  Ljung-Box Test     R    Q(15)  17.5585456 0.286584403
##  Ljung-Box Test     R    Q(20)  23.4136291 0.268943681
##  Ljung-Box Test     R^2  Q(10)  19.2821016 0.036822454
##  Ljung-Box Test     R^2  Q(15)  33.2364836 0.004352735
##  Ljung-Box Test     R^2  Q(20)  37.7425920 0.009518988
##  LM Arch Test       R    TR^2   25.4162474 0.012969006
## 
## Information Criterion Statistics:
##       AIC       BIC       SIC      HQIC 
## -6.471035 -6.409726 -6.471669 -6.446282

Observe que os \(p\)-valores no parágrafo das estimativas são bilaterais, portanto devem ser reduzidos pela metade ao considerar os parâmetros do modelo ARCH. Neste exemplo, obtemos \(\widehat{\phi}_0 =0.005\), chamado \(\mu\) na saída e \(\widehat{\phi}_1 =0.367\), chamado \(ar1\) para as estimativas dos parámetros AR(1).

As estimativas dos parámetros ARCH(1) são \(\alpha_0 = 0\), chamado ômega e para a constante \(\alpha_1 =0.194\), que é significativo com um \(p\)-valor de cerca de \(0.0418/2=0.02\). Existem vários testes realizados nos resíduos \([R]\) ou nos resíduos ao quadrado \([R^2]\).

Por exemplo, a estatística Jarque – Bera testa os resíduos do ajuste para normalidade com base na assimetria e a curtose observadas e parece que os resíduos têm alguma assimetria e curtose não normais. A estatística Shapiro-Wilk testa os resíduos do ajuste para normalidade com base nas estatísticas de ordem empíricas. Os outros testes, principalmente baseados na estatística \(Q\), são usados nos resíduos e seus quadrados.

O modelo ARCH(1) pode ser estendido ao modelo geral ARCH(p) da seguinte maneira óbvia, retemos \(r_t=\sigma_t \epsilon_t\) mas estendemos \[\begin{equation} \sigma_t^2 \, = \, \alpha_0 +\alpha_1 r_{t-1}^2+\cdots+\alpha_p r_{t-p}^2\cdot \end{equation}\] A estimação para o modelo ARCH(p) também segue de maneira similar a partir da discussão da estimação para o modelo do ARCH(1). Ou seja, a verossimilhança condicional dos dados \(r_{p+1},\cdots,r_n\) dados \(r_{1},\cdots,r_p\), é dada por \[\begin{equation} L(\alpha \, | \, r_1,\cdots,r_n) \, = \, \displaystyle \prod_{t=p+1}^n f_\alpha(r_t \, | \, r_{t-1},\cdots,r_{t-m}), \end{equation}\] onde \(\alpha=(\alpha_0,\alpha_1,\cdots,\alpha_p)\) e sob a suposição de normalidade, as densidades condicionais \(f_\alpha(\cdot|\cdot)\) são, para \(t>p\), dadas por \[\begin{equation} r_t \, | \, r_{t-1},\cdots,r_{t-p} \, \sim \, N\big( 0,\alpha_0 +\alpha_1 r_{t-1}^2+\cdots+\alpha_p r_{t-p}^2\big)\cdot \end{equation}\]

Outra extensão do ARCH é o modelo ARCH generalizado ou GARCH, desenvolvido por Bollerslev (1986). Por exemplo, um modelo GARCH(1,1) mantém \(r_t=\sigma_t\epsilon_t\) mas se estende da seguinte maneira: \[\begin{equation} \sigma_t^2 \, = \, \alpha_0+\alpha_1 r_{t-1}^2+\beta_1 \sigma_{t-1}^2\cdot \end{equation}\] Sob a condição de \(\alpha_1 + \beta_1 < 1\), usando manipulações semelhantes às anteriores, o modelo GARCH(1,1) admite um modelo ARMA(1,1) não gaussiano para o processo quadrado \[\begin{equation} r_t^2 \, = \, \alpha_0 +(\alpha_0+\beta_1)r_{t-1}^2 +\nu -\beta_1 v_{t-1}, \end{equation}\] onde \(v_t=\sigma_t^2(\epsilon_t^2-1)\). A representação acima segue escrevendo \[\begin{array}{rcl} r_t^2-\sigma_t^2 & = & \sigma_t^2 (\epsilon_t^2-1) \\ \beta_1 (r_{t-1}^2-\sigma_{t-1}^2) & = & \beta_1\sigma_{t-1}^2(\epsilon_{t-1}^2-1), \end{array}\]

subtraindo a segunda equação da primeira e usando o fato de que, \(\sigma_t^2-\beta_1 \sigma_{t-1}^2=\alpha_0+\alpha_1 r_{t-1}^2\), no lado esquerdo do resultado. O modelo GARCH(p,q) retém \(\sigma_t^2 = \alpha_0+\alpha_1 r_{t-1}^2+\beta_1 \sigma_{t-1}^2\) e se estende até \[\begin{equation} \sigma_t^2 \, = \, \alpha_0+\sum_{j=1}^p \alpha_j r_{t-j}^2 +\sum_{j=1}^q \beta_j \sigma_{t=j}^2\cdot \end{equation}\]

A estimação por máxima verossimilhançaa condicional dos parâmetros do modelo GARCH(m,r) é semelhante ao caso ARCH(m), em que a verossimilhança condicional é produto de densidade \(N(0,\sigma_t^2)\) com \(\sigma^2_t\) dadas acima e onde o condicionamento está nas primeiras \(\max(m,r)\) observações, com \(\sigma^2_1 = \cdots=\sigma^2_r = 0\). Uma vez obtidas as estimativas dos parámetros, o modelo pode ser usado para obter previsões um pouco à frente da volatilidade, digamos \(\sigma^2_{t+1}\), dadas por \[\begin{equation} \widehat{\sigma}^2_{t+1} \, = \, \widehat{\alpha}_0+\sum_{j=1}^p \widehat{\alpha}_j r^2_{t+1-j}+\sum_{j=1}^q \widehat{\beta}_j\widehat{\sigma}^2_{t+1-j}\cdot \end{equation}\] Exploramos esses conceitos no exemplo a seguir.

Exemplo 5: Análise GARCH dos retornos diários.

Como mencionado anteriormente, os retornos diários do Dow Jones Industrial Average (DJIA) exibem recursos clássicos do modelo GARCH. Além disso, observa-se a existência de alguma autocorrelação de baixo nível e, para incluir esse comportamento, usamos o pacote R fGarch para ajustar um modelo AR(1)-GARCH(1,1) à série:

library(xts)
djiar = diff(log(djia$Close))[-1]
acf2(djiar)    # exibe alguma autocorrelação

##      [,1]  [,2] [,3]  [,4]  [,5] [,6]  [,7] [,8]  [,9] [,10] [,11] [,12] [,13]
## ACF  -0.1 -0.06 0.05 -0.02 -0.06 0.01 -0.02 0.02 -0.01  0.04 -0.01  0.04  0.01
## PACF -0.1 -0.07 0.04 -0.02 -0.06 0.00 -0.02 0.03 -0.01  0.04 -0.01  0.04  0.02
##      [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25]
## ACF  -0.04 -0.06  0.06  0.00 -0.07  0.02  0.05 -0.06  0.04  0.01 -0.01     0
## PACF -0.04 -0.06  0.04  0.02 -0.06  0.00  0.04 -0.04  0.03  0.00  0.00     0
##      [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35] [,36] [,37]
## ACF      0  0.03 -0.03  0.01  0.02 -0.02  0.01  0.01 -0.09  0.03  0.03 -0.02
## PACF     0  0.04 -0.03  0.00  0.02  0.00  0.00  0.01 -0.07  0.01  0.02  0.01
##      [,38] [,39] [,40] [,41] [,42] [,43] [,44] [,45] [,46] [,47] [,48] [,49]
## ACF      0  0.03  0.02  0.01  0.00 -0.06     0     0 -0.01  0.02  0.00 -0.04
## PACF     0  0.01  0.04  0.02  0.01 -0.06     0     0 -0.01  0.01 -0.01 -0.05
##      [,50] [,51] [,52] [,53] [,54] [,55] [,56] [,57] [,58] [,59] [,60] [,61]
## ACF  -0.04 -0.02  0.03  0.00  0.00  0.01 -0.03  0.02  0.02 -0.07  0.02  0.02
## PACF -0.04 -0.03  0.01  0.01  0.02  0.01 -0.03  0.01  0.02 -0.05  0.01  0.02

acf2(djiar^2)  # exala autocorrelação

##      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]
## ACF   0.2 0.41 0.19 0.31 0.34 0.31 0.32 0.22 0.32  0.24  0.43  0.28  0.25  0.13
## PACF  0.2 0.39 0.08 0.15 0.25 0.13 0.11 0.01 0.11  0.05  0.23  0.08 -0.07 -0.16
##      [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25] [,26]
## ACF   0.22  0.26  0.27  0.27  0.17  0.23  0.25  0.19  0.29  0.15  0.18  0.15
## PACF -0.03  0.04  0.03  0.05 -0.01  0.00  0.09 -0.09  0.07 -0.01  0.00  0.00
##      [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35] [,36] [,37] [,38]
## ACF   0.28  0.22  0.22  0.14  0.16  0.20  0.14  0.25  0.10  0.17  0.12  0.16
## PACF  0.12  0.02 -0.03 -0.06  0.00  0.01 -0.06  0.06 -0.03 -0.05  0.02 -0.06
##      [,39] [,40] [,41] [,42] [,43] [,44] [,45] [,46] [,47] [,48] [,49] [,50]
## ACF   0.15  0.10  0.08  0.08  0.11  0.13  0.13  0.08  0.09  0.11  0.07  0.08
## PACF -0.06 -0.06 -0.03  0.01  0.00  0.03  0.02 -0.01  0.00  0.06 -0.02 -0.05
##      [,51] [,52] [,53] [,54] [,55] [,56] [,57] [,58] [,59] [,60] [,61]
## ACF   0.07  0.08  0.06  0.09  0.06  0.12  0.09  0.06  0.08  0.04  0.08
## PACF  0.03  0.04  0.02 -0.02 -0.09  0.08  0.06  0.00  0.01 -0.02  0.00

Figura 6:Funções de autocorrelação (ACF) e autocorrelação parcial (PACF) da série \(\nabla \log(r_t)\), sendo \(r_t\) os retornos diários do valor de fechamento da Bolsa de Valores de Nova York (DJIA) a esquerda e as mesmas funções de \(\nabla \log(r_t^2)\).

A seguir o ajuste do modelo AR(1)-GARCH(1,1):

library(fGarch)
summary(djia.g <- garchFit(~arma(1,0)+garch(1,1), data=djiar, cond.dist='std'))

## 
## Series Initialization:
##  ARMA Model:                arma
##  Formula Mean:              ~ arma(1, 0)
##  GARCH Model:               garch
##  Formula Variance:          ~ garch(1, 1)
##  ARMA Order:                1 0
##  Max ARMA Order:            1
##  GARCH Order:               1 1
##  Max GARCH Order:           1
##  Maximum Order:             1
##  Conditional Dist:          std
##  h.start:                   2
##  llh.start:                 1
##  Length of Series:          2517
##  Recursion Init:            mci
##  Series Scale:              0.01210097
## 
## Parameter Initialization:
##  Initial Parameters:          $params
##  Limits of Transformations:   $U, $V
##  Which Parameters are Fixed?  $includes
##  Parameter Matrix:
##                      U           V      params includes
##     mu     -0.15336279   0.1533628  0.01533395     TRUE
##     ar1    -0.99999999   1.0000000 -0.10129752     TRUE
##     omega   0.00000100 100.0000000  0.10000000     TRUE
##     alpha1  0.00000001   1.0000000  0.10000000     TRUE
##     gamma1 -0.99999999   1.0000000  0.10000000    FALSE
##     beta1   0.00000001   1.0000000  0.80000000     TRUE
##     delta   0.00000000   2.0000000  2.00000000    FALSE
##     skew    0.10000000  10.0000000  1.00000000    FALSE
##     shape   1.00000000  10.0000000  4.00000000     TRUE
##  Index List of Parameters to be Optimized:
##     mu    ar1  omega alpha1  beta1  shape 
##      1      2      3      4      6      9 
##  Persistence:                  0.9 
## 
## 
## --- START OF TRACE ---
## Selected Algorithm: nlminb 
## 
## R coded nlminb Solver: 
## 
##   0:     2966.5649: 0.0153339 -0.101298 0.100000 0.100000 0.800000  4.00000
##   1:     2944.7775: 0.0153351 -0.0991718 0.0803853 0.105688 0.794226  3.99986
##   2:     2910.9452: 0.0153425 -0.0864392 0.0198842 0.162519 0.809277  4.00053
##   3:     2891.9991: 0.0153434 -0.0855699 0.0350005 0.168390 0.817706  4.00077
##   4:     2882.6364: 0.0153798 -0.0489500 0.0211392 0.164590 0.840540  4.00286
##   5:     2881.9308: 0.0153823 -0.0484866 0.0204254 0.166791 0.845389  4.00314
##   6:     2881.5676: 0.0153882 -0.0491693 0.0164091 0.164254 0.847788  4.00365
##   7:     2880.9567: 0.0153959 -0.0493841 0.0174461 0.161534 0.852254  4.00434
##   8:     2880.6300: 0.0154048 -0.0492056 0.0154785 0.158065 0.855753  4.00516
##   9:     2880.3573: 0.0154165 -0.0487285 0.0159462 0.154690 0.859699  4.00626
##  10:     2880.1251: 0.0154338 -0.0478634 0.0146375 0.150753 0.862375  4.00793
##  11:     2879.9538: 0.0154636 -0.0476574 0.0149098 0.147973 0.865376  4.01086
##  12:     2879.8208: 0.0155046 -0.0483136 0.0138769 0.146462 0.866629  4.01495
##  13:     2879.7187: 0.0155485 -0.0477907 0.0141952 0.145783 0.867557  4.01934
##  14:     2879.4147: 0.0157804 -0.0423446 0.0141397 0.141896 0.868276  4.04239
##  15:     2878.3899: 0.0174942 -0.0569599 0.00849084 0.140326 0.874697  4.21043
##  16:     2878.2901: 0.0193411 -0.0502775 0.0142182 0.160120 0.862580  4.37136
##  17:     2876.3443: 0.0203758 -0.0460842 0.0117796 0.167860 0.852936  4.44569
##  18:     2874.7360: 0.0206345 -0.0556032 0.0129810 0.138773 0.869446  4.42715
##  19:     2874.6798: 0.0206349 -0.0554790 0.0122662 0.138551 0.869226  4.42717
##  20:     2874.6640: 0.0206433 -0.0553401 0.0124782 0.138563 0.869396  4.42765
##  21:     2874.6497: 0.0206511 -0.0550613 0.0122159 0.138402 0.869458  4.42809
##  22:     2874.6350: 0.0206596 -0.0549252 0.0124189 0.138406 0.869616  4.42858
##  23:     2874.6215: 0.0206675 -0.0546669 0.0121782 0.138244 0.869660  4.42903
##  24:     2874.6075: 0.0206760 -0.0545335 0.0123731 0.138241 0.869807  4.42952
##  25:     2874.5945: 0.0206841 -0.0542918 0.0121481 0.138082 0.869840  4.42998
##  26:     2874.5812: 0.0206927 -0.0541616 0.0123363 0.138075 0.869979  4.43047
##  27:     2874.5688: 0.0207009 -0.0539341 0.0121235 0.137918 0.870004  4.43094
##  28:     2874.5560: 0.0207095 -0.0538073 0.0123055 0.137908 0.870137  4.43143
##  29:     2874.5440: 0.0207178 -0.0535925 0.0121029 0.137755 0.870156  4.43190
##  30:     2874.5317: 0.0207265 -0.0534690 0.0122790 0.137743 0.870284  4.43239
##  31:     2874.5201: 0.0207349 -0.0532657 0.0120850 0.137594 0.870299  4.43287
##  32:     2874.5082: 0.0207436 -0.0531457 0.0122556 0.137580 0.870422  4.43336
##  33:     2874.4969: 0.0207520 -0.0529528 0.0120691 0.137435 0.870434  4.43384
##  34:     2874.4854: 0.0207608 -0.0528362 0.0122343 0.137420 0.870552  4.43433
##  35:     2874.4745: 0.0207693 -0.0526528 0.0120547 0.137279 0.870563  4.43482
##  36:     2874.4633: 0.0207780 -0.0525396 0.0122148 0.137264 0.870677  4.43531
##  37:     2874.4526: 0.0207866 -0.0523652 0.0120413 0.137127 0.870686  4.43580
##  38:     2874.4418: 0.0207954 -0.0522552 0.0121965 0.137111 0.870797  4.43629
##  39:     2874.4314: 0.0208040 -0.0520890 0.0120287 0.136979 0.870805  4.43678
##  40:     2874.4209: 0.0208129 -0.0519822 0.0121792 0.136962 0.870913  4.43728
##  41:     2874.4108: 0.0208215 -0.0518237 0.0120168 0.136833 0.870919  4.43776
##  42:     2874.4005: 0.0208304 -0.0517201 0.0121627 0.136816 0.871024  4.43826
##  43:     2874.3906: 0.0208391 -0.0515687 0.0120052 0.136692 0.871029  4.43875
##  44:     2874.3806: 0.0208480 -0.0514682 0.0121469 0.136674 0.871132  4.43925
##  45:     2874.3709: 0.0208568 -0.0513235 0.0119941 0.136554 0.871136  4.43974
##  46:     2874.3611: 0.0208657 -0.0512261 0.0121316 0.136536 0.871236  4.44024
##  47:     2874.3516: 0.0208745 -0.0510876 0.0119833 0.136419 0.871240  4.44073
##  48:     2874.3421: 0.0208835 -0.0509931 0.0121169 0.136401 0.871337  4.44123
##  49:     2874.3328: 0.0208923 -0.0508605 0.0119728 0.136287 0.871340  4.44172
##  50:     2874.3234: 0.0209013 -0.0507689 0.0121027 0.136269 0.871435  4.44222
##  51:     2874.3143: 0.0209101 -0.0506418 0.0119626 0.136159 0.871438  4.44271
##  52:     2874.3051: 0.0209191 -0.0505530 0.0120889 0.136141 0.871530  4.44321
##  53:     2874.2962: 0.0209280 -0.0504312 0.0119525 0.136034 0.871532  4.44370
##  54:     2874.2872: 0.0209370 -0.0503451 0.0120754 0.136016 0.871622  4.44420
##  55:     2874.2784: 0.0209460 -0.0502282 0.0119427 0.135911 0.871624  4.44470
##  56:     2874.2695: 0.0209550 -0.0501448 0.0120624 0.135894 0.871712  4.44520
##  57:     2874.2609: 0.0209640 -0.0500326 0.0119332 0.135792 0.871714  4.44569
##  58:     2874.2522: 0.0209730 -0.0499518 0.0120497 0.135775 0.871799  4.44619
##  59:     2874.2437: 0.0209820 -0.0498440 0.0119238 0.135676 0.871801  4.44669
##  60:     2874.2351: 0.0209911 -0.0497658 0.0120373 0.135658 0.871884  4.44718
##  61:     2867.5226: 0.0359780 -0.0470911 0.0137990 0.144361 0.852357  5.25825
##  62:     2865.5718: 0.0404097 -0.0578345 0.0109641 0.131220 0.867066  5.40064
##  63:     2864.1663: 0.0502410 -0.0565564 0.0123810 0.135121 0.866818  5.32198
##  64:     2861.9246: 0.0702717 -0.0514217 0.0120283 0.121547 0.870490  5.60050
##  65:     2861.7513: 0.0725452 -0.0527447 0.0111041 0.127470 0.868333  5.61948
##  66:     2861.6463: 0.0716706 -0.0559045 0.0109502 0.125504 0.869885  5.75772
##  67:     2861.6053: 0.0710293 -0.0558188 0.0109748 0.124511 0.870269  5.90135
##  68:     2861.6001: 0.0709467 -0.0554294 0.0109826 0.124373 0.870101  5.97019
##  69:     2861.6000: 0.0709462 -0.0553237 0.0110033 0.124471 0.869990  5.97824
##  70:     2861.6000: 0.0709470 -0.0553083 0.0109935 0.124438 0.870019  5.97922
##  71:     2861.6000: 0.0709479 -0.0553160 0.0109957 0.124445 0.870013  5.97877
##  72:     2861.6000: 0.0709477 -0.0553152 0.0109957 0.124445 0.870013  5.97877
## 
## Final Estimate of the Negative LLH:
##  LLH:  -8249.619    norm LLH:  -3.27756 
##            mu           ar1         omega        alpha1         beta1 
##  8.585366e-04 -5.531522e-02  1.610141e-06  1.244448e-01  8.700129e-01 
##         shape 
##  5.978768e+00 
## 
## R-optimhess Difference Approximated Hessian Matrix:
##                   mu           ar1         omega        alpha1         beta1
## mu     -4.791650e+07 -4.661860e+04 -1.205637e+09 -3.466876e+04 -9.149130e+04
## ar1    -4.661860e+04 -2.490882e+03 -1.234076e+06 -8.202990e+00 -9.744044e+01
## omega  -1.205637e+09 -1.234076e+06 -1.703966e+13 -5.515893e+08 -8.452241e+08
## alpha1 -3.466876e+04 -8.202990e+00 -5.515893e+08 -3.611230e+04 -4.469030e+04
## beta1  -9.149130e+04 -9.744044e+01 -8.452241e+08 -4.469030e+04 -6.270080e+04
## shape  -9.970368e+02  8.557384e-02 -3.050371e+06 -1.832760e+02 -2.340769e+02
##                shape
## mu     -9.970368e+02
## ar1     8.557384e-02
## omega  -3.050371e+06
## alpha1 -1.832760e+02
## beta1  -2.340769e+02
## shape  -2.547432e+00
## attr(,"time")
## Time difference of 0.6388414 secs
## 
## --- END OF TRACE ---
## 
## 
## Time to Estimate Parameters:
##  Time difference of 3.733712 secs
## 
## Title:
##  GARCH Modelling 
## 
## Call:
##  garchFit(formula = ~arma(1, 0) + garch(1, 1), data = djiar, cond.dist = "std") 
## 
## Mean and Variance Equation:
##  data ~ arma(1, 0) + garch(1, 1)
## <environment: 0x610c15efdb50>
##  [data = djiar]
## 
## Conditional Distribution:
##  std 
## 
## Coefficient(s):
##          mu          ar1        omega       alpha1        beta1        shape  
##  8.5854e-04  -5.5315e-02   1.6101e-06   1.2444e-01   8.7001e-01   5.9788e+00  
## 
## Std. Errors:
##  based on Hessian 
## 
## Error Analysis:
##          Estimate  Std. Error  t value Pr(>|t|)    
## mu      8.585e-04   1.470e-04    5.842 5.16e-09 ***
## ar1    -5.532e-02   2.023e-02   -2.735 0.006238 ** 
## omega   1.610e-06   4.459e-07    3.611 0.000305 ***
## alpha1  1.244e-01   1.660e-02    7.496 6.55e-14 ***
## beta1   8.700e-01   1.526e-02   57.022  < 2e-16 ***
## shape   5.979e+00   7.917e-01    7.551 4.31e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log Likelihood:
##  8249.619    normalized:  3.27756 
## 
## Description:
##  Fri Oct  4 17:05:34 2024 by user:  
## 
## 
## Standardised Residuals Tests:
##                                   Statistic    p-Value
##  Jarque-Bera Test   R    Chi^2  310.0086428 0.00000000
##  Shapiro-Wilk Test  R    W        0.9820293 0.00000000
##  Ljung-Box Test     R    Q(10)   16.8224625 0.07838591
##  Ljung-Box Test     R    Q(15)   26.4481300 0.03356807
##  Ljung-Box Test     R    Q(20)   28.7109940 0.09360789
##  Ljung-Box Test     R^2  Q(10)   15.3676068 0.11922325
##  Ljung-Box Test     R^2  Q(15)   19.1365028 0.20761476
##  Ljung-Box Test     R^2  Q(20)   22.9288301 0.29230264
##  LM Arch Test       R    TR^2    15.0397629 0.23926913
## 
## Information Criterion Statistics:
##       AIC       BIC       SIC      HQIC 
## -6.550353 -6.536453 -6.550364 -6.545309

O comando plot permite aceder a diversos gráficos descritivos, escolhendo which = 1 obtemos a série original. Uma outra possibilidade seria especificar quais destes gráficos queremos visualizar, para isso devemos escrever:

par(mfrow=c(2,2))
plot(djia.g, which = 1) # para ver as opções de gráficos de séries temporais
plot(djia.g, which = c(7,9,13))

Para explorar as previsões da volatilidade obtidas pelo modelo GARCH, calculamos e mostramos graficamente parte dos dados em torno das crises financeiras de 2008, juntamente com as previsões um passo à frente da volatilidade correspondente, \(\widehat{\sigma}_t\) como uma linha sólida na figura abaixo.

par(mfrow=c(1,2))
djiar1 = djiar[time(djiar)>= "2007-11-20" & time(djiar) <= "2009-11-02"]
plot(djiar1, col="lightblue")
lines(xts(djia.g@sigma.t[match(djiar["2007-11-20"][[1]],djiar):
         match(djiar["2009-11-02"][[1]],djiar)], time(djiar1)), lwd=2, col="red")
grid()

Figura 7: Previsões GARCH um passo à frente da volatilidade \(\widehat{\sigma}_t\) do DJIA, sobreposta em parte dos dados, incluindo a crise financeira de 2008.

Outro modelo que mencionamos brevemente é o modelo ARCH de potência assimétrica. O modelo mantêm \(r_t=\sigma_t \epsilon_t\), mas a variância condicional é modelada como \[\begin{equation} \sigma_t^\delta \, = \, \alpha_0 + \sum_{j=1}^p \alpha_j\big(|r_{t-j}|-\gamma_j r_{t-j} \big)^\delta + \sum_{j=1}^q \beta_j \sigma_{t-j}^\delta\cdot \end{equation}\]

Observe que o modelo é GARCH quando \(\delta= 2\) e \(\gamma_j = 0\), para \(j\in\{1,\cdots,p\}\). Os parámetros \(\gamma_j\), tais que \(|\gamma_j|\leq 1\) são os parámetros de alavancagem, que são uma medida de assimetria e \(\delta>0\) é o parámetro para o termo de potência. Um valor positivo ou negativo de \(\gamma_j\) significa que choques positivos ou negativos passados têm um impacto mais profundo na volatilidade condicional atual do que choques negativos ou positivos passados. Esse modelo combina a flexibilidade de um expoente variável com o coeficiente de assimetria para levar em consideração o efeito de alavancagem. Além disso, para garantir que \(\sigma_t> 0\), assumimos que \(\alpha_0> 0\), \(\alpha_j\geq 0\) com pelo menos um \(\alpha_j> 0\) e \(\beta_j\geq 0\). Continuamos a análise dos retornos do DJIA no exemplo a seguir.

Exemplo 6: Análise APARCH dos retornos do DJIA.

O pacote R fGarch foi usado para ajustar um modelo AR-APARCH aos retornos do DJIA discutidos no Exemplo 5. Como no exemplo anterior, incluímos um AR(1) no modelo para contabilizar a média condicional.

Nesse caso, podemos pensar no modelo como \(r_t = \mu_t + Y_t\), onde \(\mu_t\) é um AR(1) e \(Y_t\) é um ruído APARCH com variância condicional modelada como acima com erros \(t-Student\). A volatilidade prevista é, obviamente, diferente dos valores mostrados na figura do Exemplo 5, mas parece semelhante quando representada graficamente.

summary(djia.ap <- garchFit(~arma(1,0)+aparch(1,1), data=djiar, cond.dist='std'))

## 
## Series Initialization:
##  ARMA Model:                arma
##  Formula Mean:              ~ arma(1, 0)
##  GARCH Model:               aparch
##  Formula Variance:          ~ aparch(1, 1)
##  ARMA Order:                1 0
##  Max ARMA Order:            1
##  GARCH Order:               1 1
##  Max GARCH Order:           1
##  Maximum Order:             1
##  Conditional Dist:          std
##  h.start:                   2
##  llh.start:                 1
##  Length of Series:          2517
##  Recursion Init:            mci
##  Series Scale:              0.01210097
## 
## Parameter Initialization:
##  Initial Parameters:          $params
##  Limits of Transformations:   $U, $V
##  Which Parameters are Fixed?  $includes
##  Parameter Matrix:
##                      U           V      params includes
##     mu     -0.15336279   0.1533628  0.01533395     TRUE
##     ar1    -0.99999999   1.0000000 -0.10129752     TRUE
##     omega   0.00000100 100.0000000  0.10000000     TRUE
##     alpha1  0.00000001   1.0000000  0.10000000     TRUE
##     gamma1 -0.99999999   1.0000000  0.10000000     TRUE
##     beta1   0.00000001   1.0000000  0.80000000     TRUE
##     delta   0.00000000   2.0000000  2.00000000     TRUE
##     skew    0.10000000  10.0000000  1.00000000    FALSE
##     shape   1.00000000  10.0000000  4.00000000     TRUE
##  Index List of Parameters to be Optimized:
##     mu    ar1  omega alpha1 gamma1  beta1  delta  shape 
##      1      2      3      4      5      6      7      9 
##  Persistence:                  0.901 
## 
## 
## --- START OF TRACE ---
## Selected Algorithm: nlminb 
## 
## R coded nlminb Solver: 
## 
##   0:     2956.2216: 0.0153339 -0.101298 0.100000 0.100000 0.100000 0.800000  2.00000  4.00000
##   1:     2933.7103: 0.0153349 -0.0994951 0.0816703 0.105086 0.101441 0.794372  1.99989  3.99987
##   2:     2915.7667: 0.0153364 -0.0969033 0.0659020 0.116672 0.103734 0.795785  2.00000  3.99997
##   3:     2868.4756: 0.0153445 -0.0848310 0.0240069 0.174769 0.115713 0.826262  2.00000  4.00127
##   4:     2866.0082: 0.0153484 -0.0806167 0.0223779 0.179100 0.121372 0.835425  2.00000  4.00179
##   5:     2865.7462: 0.0153525 -0.0768065 0.0131299 0.176106 0.127204 0.838913  1.99992  4.00228
##   6:     2863.1438: 0.0153538 -0.0756661 0.0180386 0.176968 0.129009 0.841881  2.00000  4.00245
##   7:     2856.2559: 0.0153874 -0.0476933 0.00966761 0.170084 0.175631 0.860615  1.99995  4.00607
##   8:     2853.4276: 0.0154247 -0.0296623 0.0252149 0.165491 0.224425 0.839248  1.99938  4.00935
##   9:     2853.1938: 0.0154256 -0.0300866 0.0153713 0.163947 0.225567 0.836621  1.99923  4.00940
##  10:     2850.4660: 0.0154260 -0.0302179 0.0197141 0.165283 0.226119 0.839057  1.99934  4.00949
##  11:     2849.7425: 0.0154283 -0.0310641 0.0167839 0.166273 0.228959 0.841965  1.99939  4.00980
##  12:     2849.5692: 0.0154315 -0.0322119 0.0190914 0.167177 0.232819 0.844052  1.99942  4.01021
##  13:     2848.6409: 0.0154353 -0.0333690 0.0169749 0.166916 0.237263 0.843022  1.99933  4.01062
##  14:     2848.3152: 0.0154386 -0.0342576 0.0174318 0.167026 0.241137 0.846288  1.99937  4.01105
##  15:     2847.6373: 0.0154426 -0.0351956 0.0157169 0.166737 0.245766 0.845141  1.99926  4.01149
##  16:     2839.1179: 0.0155788 -0.0588461 0.0224608 0.179509 0.397542 0.822415  1.99659  4.02659
##  17:     2838.8052: 0.0155793 -0.0585367 0.0207575 0.178821 0.397822 0.822264  1.99652  4.02663
##  18:     2838.1726: 0.0155842 -0.0554703 0.0209085 0.175805 0.400777 0.827695  1.99629  4.02724
##  19:     2837.7162: 0.0155992 -0.0507970 0.0191626 0.176094 0.405520 0.825669  1.99483  4.02886
##  20:     2837.4475: 0.0156202 -0.0464274 0.0206223 0.176945 0.410239 0.825177  1.99250  4.03120
##  21:     2837.1785: 0.0156567 -0.0447657 0.0194430 0.176216 0.412875 0.824921  1.98763  4.03523
##  22:     2836.9040: 0.0156936 -0.0461077 0.0202590 0.174871 0.412771 0.826868  1.98246  4.03931
##  23:     2836.6516: 0.0157301 -0.0471688 0.0190315 0.172895 0.412978 0.828105  1.97731  4.04342
##  24:     2836.3797: 0.0157682 -0.0466784 0.0200044 0.172404 0.414501 0.829093  1.97209  4.04780
##  25:     2836.1207: 0.0158048 -0.0455259 0.0190642 0.171430 0.416707 0.829267  1.96710  4.05211
##  26:     2835.8676: 0.0158426 -0.0452688 0.0199533 0.170875 0.418023 0.830310  1.96188  4.05658
##  27:     2835.6267: 0.0158798 -0.0455789 0.0189661 0.169464 0.418888 0.831065  1.95666  4.06102
##  28:     2835.3886: 0.0159172 -0.0456461 0.0197800 0.168803 0.419888 0.832210  1.95143  4.06555
##  29:     2835.1543: 0.0159542 -0.0451712 0.0189081 0.167816 0.421426 0.832533  1.94629  4.07010
##  30:     2834.9241: 0.0159914 -0.0447683 0.0197123 0.167375 0.422845 0.833416  1.94116  4.07471
##  31:     2834.7006: 0.0160283 -0.0446969 0.0188494 0.166312 0.424006 0.833876  1.93599  4.07932
##  32:     2834.4814: 0.0160655 -0.0447235 0.0195903 0.165737 0.425022 0.834911  1.93082  4.08398
##  33:     2834.2655: 0.0161023 -0.0445152 0.0187773 0.164800 0.426283 0.835268  1.92569  4.08866
##  34:     2834.0526: 0.0161393 -0.0442157 0.0195129 0.164393 0.427582 0.836108  1.92059  4.09338
##  35:     2833.8441: 0.0161762 -0.0440328 0.0187322 0.163510 0.428810 0.836443  1.91548  4.09810
##  36:     2833.6395: 0.0162132 -0.0439777 0.0194243 0.163042 0.429877 0.837354  1.91037  4.10286
##  37:     2833.4380: 0.0162501 -0.0438482 0.0186721 0.162195 0.431041 0.837682  1.90527  4.10763
##  38:     2833.2392: 0.0162871 -0.0436345 0.0193522 0.161825 0.432241 0.838481  1.90021  4.11243
##  39:     2833.0436: 0.0163240 -0.0434436 0.0186346 0.161057 0.433449 0.838752  1.89514  4.11724
##  40:     2832.8513: 0.0163612 -0.0433460 0.0192843 0.160671 0.434536 0.839569  1.89008  4.12207
##  41:     2832.6617: 0.0163983 -0.0432372 0.0185881 0.159919 0.435658 0.839852  1.88502  4.12692
##  42:     2832.4746: 0.0164355 -0.0430720 0.0192227 0.159591 0.436797 0.840603  1.87999  4.13178
##  43:     2822.2407: 0.0240943 -0.0195785 0.0280234 0.146120 0.600308 0.860703 0.840968  5.13460
##  44:     2819.9890: 0.0240944 -0.0196461 0.0297368 0.146988 0.600356 0.862147 0.841026  5.13461
##  45:     2819.7948: 0.0240980 -0.0205470 0.0281911 0.146928 0.600991 0.863159 0.841962  5.13508
##  46:     2819.5859: 0.0241091 -0.0209043 0.0285820 0.146900 0.600826 0.863833 0.843477  5.13656
##  47:     2819.4916: 0.0241103 -0.0211540 0.0281429 0.146655 0.600817 0.863952 0.845800  5.13679
##  48:     2819.3963: 0.0241121 -0.0213389 0.0282730 0.146627 0.600722 0.864336 0.848138  5.13709
##  49:     2819.3040: 0.0241142 -0.0216029 0.0278552 0.146362 0.600700 0.864453 0.850446  5.13742
##  50:     2819.2151: 0.0241151 -0.0217604 0.0279896 0.146344 0.600616 0.864809 0.852807  5.13759
##  51:     2819.1308: 0.0241156 -0.0219735 0.0276007 0.146108 0.600609 0.864889 0.855155  5.13768
##  52:     2819.0525: 0.0241140 -0.0220674 0.0277400 0.146128 0.600558 0.865197 0.857525  5.13751
##  53:     2818.9797: 0.0241115 -0.0221969 0.0273818 0.145946 0.600585 0.865223 0.859866  5.13720
##  54:     2818.9068: 0.0241092 -0.0222721 0.0275194 0.145971 0.600548 0.865509 0.862227  5.13692
##  55:     2818.8339: 0.0241077 -0.0224172 0.0271744 0.145770 0.600561 0.865539 0.864583  5.13672
##  56:     2818.7607: 0.0241066 -0.0225111 0.0273113 0.145773 0.600513 0.865827 0.866957  5.13655
##  57:     2818.6876: 0.0241062 -0.0226714 0.0269775 0.145555 0.600514 0.865862 0.869321  5.13646
##  58:     2818.6145: 0.0241061 -0.0227809 0.0271137 0.145538 0.600456 0.866153 0.871699  5.13639
##  59:     2818.5416: 0.0241067 -0.0229523 0.0267898 0.145306 0.600447 0.866191 0.874064  5.13639
##  60:     2818.4692: 0.0241076 -0.0230722 0.0269257 0.145276 0.600384 0.866481 0.876442  5.13640
##  61:     2818.4041: 0.0241063 -0.0231831 0.0266160 0.145084 0.600402 0.866484 0.878798  5.13612
##  62:     2818.3364: 0.0241058 -0.0232726 0.0267534 0.145071 0.600356 0.866746 0.881177  5.13596
##  63:     2818.2573: 0.0241116 -0.0235511 0.0264330 0.144755 0.600302 0.866843 0.883403  5.13659
##  64:     2818.1796: 0.0241171 -0.0237603 0.0265664 0.144647 0.600196 0.867175 0.885669  5.13715
##  65:     2818.1059: 0.0241218 -0.0239882 0.0262598 0.144364 0.600156 0.867238 0.887969  5.13756
##  66:     2818.0350: 0.0241255 -0.0241360 0.0263963 0.144298 0.600079 0.867528 0.890320  5.13781
##  67:     2817.9652: 0.0241298 -0.0243298 0.0261004 0.144038 0.600053 0.867570 0.892653  5.13807
##  68:     2817.8963: 0.0241342 -0.0244711 0.0262362 0.143973 0.599983 0.867852 0.895001  5.13831
##  69:     2817.8274: 0.0241398 -0.0246664 0.0259464 0.143709 0.599956 0.867894 0.897314  5.13864
##  70:     2817.7592: 0.0241456 -0.0248119 0.0260812 0.143636 0.599889 0.868176 0.899643  5.13893
##  71:     2817.6912: 0.0241526 -0.0250069 0.0257969 0.143369 0.599865 0.868218 0.901932  5.13931
##  72:     2817.6243: 0.0241597 -0.0251510 0.0259308 0.143293 0.599806 0.868496 0.904241  5.13964
##  73:     2817.5574: 0.0241681 -0.0253415 0.0256518 0.143026 0.599789 0.868536 0.906505  5.14005
##  74:     2817.4916: 0.0241767 -0.0254823 0.0257846 0.142948 0.599740 0.868811 0.908791  5.14040
##  75:     2817.4257: 0.0241867 -0.0256684 0.0255103 0.142683 0.599733 0.868850 0.911022  5.14084
##  76:     2817.3609: 0.0241969 -0.0258055 0.0256420 0.142603 0.599699 0.869122 0.913277  5.14122
##  77:     2817.2959: 0.0242087 -0.0259862 0.0253719 0.142339 0.599705 0.869159 0.915465  5.14169
##  78:     2817.2321: 0.0242206 -0.0261181 0.0255022 0.142260 0.599689 0.869427 0.917681  5.14208
##  79:     2817.1680: 0.0242344 -0.0262923 0.0252360 0.141999 0.599713 0.869463 0.919816  5.14257
##  80:     2817.1051: 0.0242482 -0.0264181 0.0253649 0.141921 0.599718 0.869728 0.921984  5.14298
##  81:     2817.0420: 0.0242641 -0.0265851 0.0251023 0.141663 0.599764 0.869762 0.924055  5.14349
##  82:     2816.9799: 0.0242799 -0.0267039 0.0252296 0.141587 0.599794 0.870024 0.926163  5.14390
##  83:     2816.9163: 0.0243022 -0.0269415 0.0249643 0.141289 0.599857 0.870104 0.927736  5.14488
##  84:     2816.8692: 0.0243071 -0.0268958 0.0250993 0.141350 0.599953 0.870265 0.930024  5.14426
##  85:     2816.8305: 0.0243063 -0.0268219 0.0248626 0.141301 0.600094 0.870153 0.932199  5.14328
##  86:     2816.7711: 0.0243245 -0.0269168 0.0249870 0.141234 0.600176 0.870398 0.934250  5.14357
##  87:     2816.7144: 0.0243514 -0.0272256 0.0247187 0.140902 0.600283 0.870526 0.934891  5.14498
##  88:     2816.6679: 0.0243583 -0.0271709 0.0248546 0.140956 0.600399 0.870674 0.937155  5.14435
##  89:     2816.6202: 0.0243672 -0.0271585 0.0246204 0.140837 0.600550 0.870604 0.939405  5.14382
##  90:     2816.5605: 0.0243913 -0.0272849 0.0247434 0.140738 0.600663 0.870870 0.941100  5.14439
##  91:     2816.5046: 0.0244199 -0.0275626 0.0244856 0.140427 0.600820 0.870987 0.941792  5.14563
##  92:     2816.4498: 0.0244391 -0.0275902 0.0246160 0.140398 0.600959 0.871188 0.943828  5.14563
##  93:     2816.4016: 0.0244519 -0.0275781 0.0243891 0.140275 0.601140 0.871122 0.946012  5.14517
##  94:     2816.3445: 0.0244777 -0.0276762 0.0245110 0.140198 0.601301 0.871367 0.947627  5.14562
##  95:     2816.2924: 0.0245078 -0.0279438 0.0242608 0.139905 0.601498 0.871479 0.948081  5.14684
##  96:     2816.2375: 0.0245327 -0.0279990 0.0243856 0.139858 0.601678 0.871694 0.949812  5.14705
##  97:     2816.1902: 0.0245484 -0.0279792 0.0241669 0.139744 0.601890 0.871623 0.951911  5.14660
##  98:     2816.1360: 0.0245757 -0.0280508 0.0242884 0.139689 0.602087 0.871846 0.953460  5.14692
##  99:     2816.0877: 0.0246066 -0.0283119 0.0240468 0.139415 0.602310 0.871951 0.953700  5.14810
## 100:     2816.0355: 0.0246322 -0.0283406 0.0241717 0.139389 0.602524 0.872143 0.955401  5.14818
## 101:     2815.9920: 0.0246464 -0.0282868 0.0239647 0.139305 0.602754 0.872048 0.957498  5.14752
## 102:     2815.9397: 0.0246757 -0.0283625 0.0240840 0.139250 0.602982 0.872267 0.958866  5.14788
## 103:     2815.8957: 0.0247063 -0.0286352 0.0238502 0.138983 0.603218 0.872375 0.958800  5.14912
## 104:     2815.8452: 0.0247331 -0.0286561 0.0239745 0.138963 0.603458 0.872556 0.960419  5.14916
## 105:     2815.8047: 0.0247459 -0.0285791 0.0237787 0.138900 0.603697 0.872442 0.962505  5.14837
## 106:     2813.4253: 0.0276797 -0.0247133 0.0241936 0.130265 0.616571 0.873922  1.09058  5.17431
## 107:     2812.5865: 0.0304454 -0.0420449 0.0218254 0.114614 0.626011 0.887083  1.03010  5.30445
## 108:     2810.5922: 0.0338870 -0.0386060 0.0201970 0.115783 0.658178 0.891626  1.03486  5.35422
## 109:     2808.9409: 0.0373444 -0.0312035 0.0153824 0.119147 0.692089 0.890359  1.05908  5.38998
## 110:     2806.7660: 0.0404647 -0.0385326 0.0162136 0.120838 0.702312 0.889608  1.08720  5.49759
## 111:     2805.0914: 0.0436369 -0.0524132 0.0169129 0.113035 0.720936 0.892846  1.06748  5.59831
## 112:     2803.4267: 0.0544190 -0.0413269 0.0141070 0.0950836 0.921342 0.910146 0.960211  5.56332
## 113:     2799.8675: 0.0559346 -0.0549879 0.0164249 0.106965 0.979412 0.894811  1.06810  5.74935
## 114:     2799.2320: 0.0522843 -0.0519325 0.0200920 0.103246  1.00000 0.895146  1.06053  5.77393
## 115:     2798.1942: 0.0487045 -0.0512814 0.0216800 0.104878  1.00000 0.891831  1.02628  5.83936
## 116:     2797.4598: 0.0439573 -0.0559564 0.0188027 0.0982294  1.00000 0.896196  1.08059  6.42234
## 117:     2796.9661: 0.0437284 -0.0488382 0.0203989 0.0984266  1.00000 0.892755  1.09220  6.83454
## 118:     2796.8946: 0.0447421 -0.0484723 0.0201060 0.0973110  1.00000 0.893595  1.09086  7.06828
## 119:     2796.8587: 0.0434068 -0.0485694 0.0201879 0.0976983  1.00000 0.893534  1.08956  7.27075
## 120:     2796.8497: 0.0431628 -0.0482850 0.0201765 0.0974632  1.00000 0.894567  1.07653  7.31167
## 121:     2796.8474: 0.0433511 -0.0481651 0.0202061 0.0979972  1.00000 0.894505  1.07120  7.30015
## 122:     2796.8472: 0.0432635 -0.0482294 0.0202626 0.0980864  1.00000 0.894460  1.07011  7.28789
## 123:     2796.8472: 0.0432502 -0.0481709 0.0202595 0.0980905  1.00000 0.894456  1.07024  7.28554
## 124:     2796.8472: 0.0432527 -0.0481837 0.0202595 0.0980896  1.00000 0.894456  1.07023  7.28577
## 
## Final Estimate of the Negative LLH:
##  LLH:  -8311.583    norm LLH:  -3.302178 
##            mu           ar1         omega        alpha1        gamma1 
##  0.0005233999 -0.0481836532  0.0001798036  0.0980895866  0.9999999900 
##         beta1         delta         shape 
##  0.8944562433  1.0702336395  7.2857681269 
## 
## R-optimhess Difference Approximated Hessian Matrix:
##                   mu           ar1         omega        alpha1        gamma1
## mu     -9.257238e+07 -1.304435e+05 -6.612961e+08 -2.435863e+06 -1.277034e+05
## ar1    -1.304435e+05 -2.960676e+03 -7.707977e+05 -1.492688e+03 -7.809030e+01
## omega  -6.612961e+08 -7.707977e+05 -9.932961e+09 -3.284934e+07 -1.721812e+06
## alpha1 -2.435863e+06 -1.492688e+03 -3.284934e+07 -1.662165e+05 -8.715904e+03
## gamma1 -1.277034e+05 -7.809030e+01 -1.721812e+06 -8.715904e+03 -9.621191e+03
## beta1  -3.642694e+06 -3.791562e+03 -5.191557e+07 -2.185954e+05 -1.146445e+04
## delta  -2.481011e+05 -2.643141e+02 -3.568303e+06 -1.497204e+04 -7.830997e+02
## shape  -4.448795e+03 -2.850967e+00 -4.321021e+04 -2.197835e+02 -1.153887e+01
##                beta1         delta         shape
## mu     -3.642694e+06 -2.481011e+05  -4448.794711
## ar1    -3.791562e+03 -2.643141e+02     -2.850967
## omega  -5.191557e+07 -3.568303e+06 -43210.212959
## alpha1 -2.185954e+05 -1.497204e+04   -219.783487
## gamma1 -1.146445e+04 -7.830997e+02    -11.538867
## beta1  -3.201300e+05 -2.153576e+04   -285.455733
## delta  -2.153576e+04 -1.520674e+03    -19.344024
## shape  -2.854557e+02 -1.934402e+01     -1.118535
## attr(,"time")
## Time difference of 1.192186 secs
## 
## --- END OF TRACE ---
## 
## 
## Time to Estimate Parameters:
##  Time difference of 7.492474 secs
## 
## Title:
##  GARCH Modelling 
## 
## Call:
##  garchFit(formula = ~arma(1, 0) + aparch(1, 1), data = djiar, 
##     cond.dist = "std") 
## 
## Mean and Variance Equation:
##  data ~ arma(1, 0) + aparch(1, 1)
## <environment: 0x610c0d763da8>
##  [data = djiar]
## 
## Conditional Distribution:
##  std 
## 
## Coefficient(s):
##         mu         ar1       omega      alpha1      gamma1       beta1  
##  0.0005234  -0.0481837   0.0001798   0.0980896   1.0000000   0.8944562  
##      delta       shape  
##  1.0702336   7.2857681  
## 
## Std. Errors:
##  based on Hessian 
## 
## Error Analysis:
##          Estimate  Std. Error  t value Pr(>|t|)    
## mu      5.234e-04   1.525e-04    3.432 0.000598 ***
## ar1    -4.818e-02   1.934e-02   -2.491 0.012728 *  
## omega   1.798e-04   3.443e-05    5.222 1.77e-07 ***
## alpha1  9.809e-02   1.030e-02    9.525  < 2e-16 ***
## gamma1  1.000e+00   1.045e-02   95.729  < 2e-16 ***
## beta1   8.945e-01   1.049e-02   85.279  < 2e-16 ***
## delta   1.070e+00   1.350e-01    7.928 2.22e-15 ***
## shape   7.286e+00   1.123e+00    6.489 8.61e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log Likelihood:
##  8311.583    normalized:  3.302178 
## 
## Description:
##  Fri Oct  4 17:05:47 2024 by user:  
## 
## 
## Standardised Residuals Tests:
##                                  Statistic     p-Value
##  Jarque-Bera Test   R    Chi^2  245.155735 0.000000000
##  Shapiro-Wilk Test  R    W        0.983058 0.000000000
##  Ljung-Box Test     R    Q(10)   15.595834 0.111801446
##  Ljung-Box Test     R    Q(15)   26.450945 0.033541649
##  Ljung-Box Test     R    Q(20)   30.170733 0.067133999
##  Ljung-Box Test     R^2  Q(10)   19.176877 0.038072747
##  Ljung-Box Test     R^2  Q(15)   30.466634 0.010345747
##  Ljung-Box Test     R^2  Q(20)   35.364488 0.018246766
##  LM Arch Test       R    TR^2    29.577347 0.003231656
## 
## Information Criterion Statistics:
##       AIC       BIC       SIC      HQIC 
## -6.598000 -6.579468 -6.598020 -6.591274

plot(djiar1, col="lightblue")

lines(xts(djia.ap@sigma.t[match(djiar["2007-11-20"][[1]],djiar):
    match(djiar["2009-11-02"][[1]],djiar)], time(djiar1)), lwd=2, col="red")

Na maioria das aplicações, a distribuição do ruído \(\epsilon_t\), raramente é normal. O pacote R fGarch permite que diversas distribuições sejam ajustadas aos dados; consulte o arquivo de ajuda para obter informações. Algumas desvantagens do GARCH e modelos relacionados são as seguintes.

1. O modelo GARCH assume que retornos positivos e negativos têm o mesmo efeito porque a volatilidade depende de retornos quadrados; os modelos assimétricos ajudam a aliviar esse problema.

2. Esses modelos geralmente são restritivos devido às restrições rígidas nos parámetros do modelo, por exemplo, para um ARCH(1), \(0\leq \alpha_1^2 < \frac{1}{3}\).

3. A verossimilhança é plana, a menos que \(n\) seja muito grande.

4. Os modelos tendem a predizer a volatilidade porque respondem lentamente a grandes retornos isolados.

Várias extensõs ao modelo original foram propostas para superar algumas das deficiências que acabamos de mencionar. Por exemplo, discutimos o fato de que o pacote de funções fGarch permite dinâmica de retorno assimétrica. No caso de persistência na volatilidade, o modelo GARCH integrado ou IGARCH pode ser usado. Lembre-se, mostramos que o modelo GARCH(1,1) pode ser escrito como \[\begin{equation} r_t^2 \, = \, \alpha_0 +(\alpha_1+\beta_1)r_{t-1}^2+\nu_t-\beta_1 \nu_{t-1}, \end{equation}\] e \(r_t^2\) é estacionário se \(\alpha_1+\beta_1 < 1\). No modelo IGARCH fazemos \(\alpha_1+\beta_1=1\), nesse caso, o modelo IGARCH(1,1) é \[\begin{equation} r_t \, = \, \sigma_t\epsilon_t \qquad \mbox{e} \qquad \sigma_t^2 \, = \, \alpha_0+(1-\beta_1)r_{t-1}^2+\beta_1\sigma_{t-1}^2\cdot \end{equation}\]

Existem muitas extensões diferentes no modelo ARCH básico que foram desenvolvidas para lidar com as várias situações observadas na prática. Os leitores interessados podem encontrar as discussões gerais em Engle, Nelson, and Bollerslev (1994) e Shephard (1996) valem a pena ler as referidas referências. Além disso, Gouriéroux (1997) faz uma apresentação detalhada do ARCH e modelos relacionados com aplicações financeiras e contêm uma extensa bibliografia. Dois excelentes textos sobre análise financeira de séries temporais são Chan (2002) e Tsay (2002).

Por fim, discutimos brevemente modelos de volatilidade estocástica. O componente de volatilidade, \(\sigma^2_t\), no GARCH e os modelos relacionados são condicionalmente não estocásticos. Por exemplo, no modelo ARCH(1), sempre que o retorno anterior for avaliado em, digamos \(c\), isto é, \(r_{t-1}=c\), deve ser o caso \(\sigma^2_t =\alpha_0 + \alpha_1 c^2\). Essa suposição parece um pouco irrealista. O modelo de volatilidade estocástica adiciona um componente estocástico à volatilidade da seguinte maneira. No modelo GARCH, um retorno, digamos, \(r_t\), é \[\begin{equation} r_t=\sigma_t \epsilon_t, \qquad \mbox{o qual implica que} \qquad \log(r_t^2) \, = \, \log(\sigma^2_t)+\log(\epsilon_t^2)\cdot \end{equation}\] Assim, as observações \(\log(r^2_t)\) são geradas por dois componentes, a volatilidade não observada, \(\log(\sigma^2_t)\) e o ruído não observado, \(\log(\epsilon_t^2)\).

Enquanto, por exemplo, o modelo GARCH(1,1) modela a volatilidade sem erro, \(\sigma^2_{t + 1} = \alpha_0+ \alpha_1 r_t^2 + \beta_1\sigma^2_t\), o modelo básico de volatilidade estocástica assume que a variável latente registrada é um processo autoregressivo, \[\begin{equation} \log(\sigma^2_{t+1}) \, = \, \phi_0+\phi_1 \log(\sigma^2_t)+W_t, \end{equation}\] em que \(W_t\sim N(0,\sigma_{_W}^2)\) são independentes.

A introdução do termo de ruído torna o processo de volatilidade latente estocástico. As duas expressões acima juntas compõem o modelo de volatilidade estocástica. Dadas \(n\) observações, o objetivo é estimar os parâmetros \(\phi_0\), \(\phi_1\) e \(\sigma^2_{_W}\) e prever a volatilidade futura.

4. Modelos de limiar

Discutimos anteriormente o fato de que, para uma série temporal estacionária, a melhor predição linear para a frente no tempo é igual à melhor predição linear para trás no tempo. Este resultado decorreu do fato de que a matriz de variâncias e covariâncias de \(X_{1:n} = \{X_1,X_2,\cdots,X_n\}\), digamos, \(\Gamma = \{\gamma(i-j)\}_{ij,=1}^n\), é a mesma do que a matriz de variâncias e covariâncias de \(X_{n:1} = \{X_n,X_{n-1},X_{n-2},\cdots,X_1\}\). Além disso, se o processo for gaussiano, as distribuições de \(X_{1:n}\) e \(X_{n:1}\) são idênticas. Nesse caso, um gráfico de tempo de \(X_{1:n}\), ou seja, um gráfico dos dados traçados para frente no tempo, deve ser semelhante a um gráfico de tempo de \(X_{n:1}\), ou seja, um gráfico dos dados traçados para trás no tempo.

Existem, no entanto, muitas séries que não se enquadram nesta categoria. Por exemplo, a figura abaixo mostra o gráfico das mortes mensais por pneumonia e influenza para cada 10.000 habitantes nos EUA por 11 anos, de 1968 a 1978.

library(astsa)
par(mar=c(4,4,1,1))
plot(flu, type="b", xlab="Anos", pch=19, cex=0.7)
grid()

Normalmente, o número de mortes tende a aumentar ↑ mais rápido do que diminuir ↘, especialmente durante epidemias. Assim, se os dados fossem mostrados para trás no tempo, essa série tenderia a aumentar mais lentamente do que a diminuir.

Além disso, se as mortes mensais por pneumonia e influenza seguirem um processo linear gaussiano, não esperariamos ver tantos surtos de mudanças positivas e negativas que ocorrem periodicamente nesta série. Também, embora o número de mortes seja normalmente maior durante os meses de inverno, os dados não são perfeitamente sazonais, ou seja, embora o pico da série muitas vezes ocorra em janeiro, nos demais anos o pico ocorre em fevereiro ou março, portanto, os modelos ARMA sazonais não captariam esse comportamento.

Existem muitas abordagens para modelar séries não lineares que podem ser usadas (ver Priestley 1988); aqui, nos concentramos na classe de modelos de limiar TARMA, apresentada em Tong (1983), Tong (1990). A ideia básica desses modelos é ajustar modelos ARMA lineares locais e seu apelo seria que podemos usar a intuição de ajustar modelos ARMA lineares globais.

Por exemplo, um modelo de limiar de \(k\)-regimes SETARMA tem a forma \[\begin{equation*} X_t \, = \, \left\{ \begin{array}{rcl} \displaystyle \phi_0^{(1)}+\sum_{i=1}^{p_1} \phi_i^{(1)}X_{t-i}+W_t^{(1)}+\sum_{j=1}^{q_1} \theta_j^{(1)} W_{t-j}^{(1)}, & \mbox{caso} & X_{t-d}\leq r_1, \\ \displaystyle \phi_0^{(2)}+\sum_{i=1}^{p_2} \phi_i^{(2)}X_{t-i}+W_t^{(2)}+\sum_{j=1}^{q_2} \theta_j^{(2)} W_{t-j}^{(2)}, & \mbox{caso} & r_1< X_{t-d}\leq r_2, \\ \vdots & & \vdots \\ \displaystyle \phi_0^{(k)}+\sum_{i=1}^{p_k} \phi_i^{(k)}X_{t-i}+W_t^{(k)}+\sum_{j=1}^{q_k} \theta_j^{(k)} W_{t-j}^{(k)}, & \mbox{caso} & r_{k-1}< X_{t-d}, \\ \end{array} \right. \end{equation*}\] onde \(W_t^{(j)}\sim N(0,\sigma_{_j}^2)\) independentes indeticamente distribuídos, para \(j = 1,\cdots,k\), o inteiro positivo \(d\) é um atraso ou delay específico e \(-\infty < r_1 < \cdots < r_{k-1} < \infty\) consitui uma partição dos números reais.

Esses modelos permitem mudanças nos coeficientes ARMA ao longo do tempo e essas mudanças são determinadas pela comparação de valores anteriores, retrocedidos por um intervalo de tempo igual a \(d\), a valores de limite fixos. Cada modelo ARMA diferente é conhecido como um regime. Na definição acima, os valores \((p_j, q_j)\) da ordem dos modelos ARMA podem diferir em cada regime, embora em muitas aplicações, eles sejam iguais.

Estacionaridade e a invertibilidade são preocupações óbvias ao ajustar modelos de série temporal. Para os modelos de séries temporais de limiar, como modelos TAR, TMA e TARMA, no entanto, as condições de estacionaridade e invertibilidade na literatura são menos conhecidas em geral e frequentemente são modelos restritos de ordem um.

O modelo pode ser generalizado para incluir a possibilidade de que os regimes dependam de uma coleção de valores passados do processo ou de que os regimes dependam de uma variável exôgena, como casos presa-predador. Por exemplo, o Lince Canadense (Canadian lynx) foi exaustivamente estudado, consulte o conjunto de dados R lynx e esta série é normalmente usada para demonstrar o ajuste de modelos de limiar.

A presa do Lince varia de pequenos roedores a veados, com a Lebre sendo sua presa predominantemente favorita. Na verdade, em certas áreas, o Lince está tão intimamente ligado à Lebre que sua população aumenta e diminui conforme muda a população da Lebre, embora outras fontes de alimento possam ser abundantes. Nesse caso, parece razoável substituir \(X_{t-d}\) por \(Y_{t-d}\), onde \(Y_t\) é o tamanho da população de Lebre.

A popularidade dos modelos TAR se deve ao fato de serem relativamente simples de especificar, estimar e interpretar, em comparação com muitos outros modelos de série temporais não lineares. Além disso, apesar de sua aparente simplicidade, a classe de modelos TAR pode reproduzir muitos fenômenos não lineares. No exemplo a seguir, usamos esse métodos para ajustar um modelo de limiar para a série mensal de mortes por pneumonia e influenza mencionadas anteriormente.

Exemplo 7: Modelagem de limiar da série de influenza

library(astsa)
# Mostrando os dados com as iniciais do mês como pontos
par(mar=c(4,4,2,1))
plot(flu, type="c", xlab="Anos", cex=0.7)
Months = c("J","F","M","A","M","J","J","A","S","O","N","D")
points(flu, pch=Months, cex=.8, font=2)
grid()

Como discutido anteriormente, o exame da figura acima nos leva a acreditar que a série temporal de mortes mensais por pneumonia e influenza, digamos \(flu_t\), não é linear. Também é evidente da figura que há uma ligeira tendência negativa nos dados.

Descobrimos que a maneira mais conveniente de ajustar um modelo de limiar a esses dados, enquanto removemos a tendência, é trabalhar com as primeiras diferenças. Os dados diferenciados, digamos \[\begin{equation*} X_t \, = \, flu_t \, - \, flu_{t-1} \end{equation*}\] são exibidos na figura abaixo como pontos (+) representando as observações.

# Iniciar análise
dflu = diff(flu)
par(mar=c(4,4,2,1))
lag1.plot(dflu, corr=FALSE, pch=19, cex=0.7) # scatterplot com ajuste lowess
abline(v=.5, col=4, lty=6)

Figura 8: Gráfico de dispersão de \(dflu_t = flu_t - flu_{t-1}\) versus \(dflu_{t-1}\) com um ajuste lowess sobreposto (linha). Uma linha tracejada vertical indica \(dflu_{t-1} = 0.05\).

A não linearidade dos dados é mais pronunciada no gráfico das primeiras diferenças \(X_t\). Claramente, \(X_t\) sobe lentamente por alguns meses e, então, em algum momento no inverno, tem a possibilidade de pular para um grande número quando \(X_t\) exceder cerca de \(0.05\). Se o processo der um grande salto, ocorrerá uma diminuição significativa subsequente em \(X_t\).

Outro gráfico revelador é o gráfico de defasagem de \(X_t\) versus \(X_{t-1}\) mostrado na figura acima, este sugere a possibilidade de dois regimes lineares baseados se \(X_{t-1}\) excede ou não 0.05.

Como uma análise inicial, ajustamos o seguinte modelo de limiar \[\begin{array}{rcl} X_t & = & \displaystyle \alpha^{(1)} \, + \, \sum_{j=1}^p \phi_j^{(1)} X_{t-j} \, + \, W_t^{(1)}, \qquad X_{t-1}< 0.05; \\ X_t & = & \displaystyle \alpha^{(2)} \, + \, \sum_{j=1}^p \phi_j^{(2)} X_{t-j} \, + \, W_t^{(2)}, \qquad X_{t-1}\geq 0.05, \end{array}\]

com \(p = 6\), assumindo que seria maior do que o necessário. O modelo acima é fácil de ajustar usando duas execuções de regressão linear, uma quando \(X_t <0.05\) e a outra quando \(X_t\geq 0.05\). Os detalhes são fornecidos no código R a continuação.

Usando o modelo final, as previsões com um mês de antecedência podem ser feitas e são mostradas na figura abaixo como uma linha. O modelo se sai extremamente bem em prever uma epidemia de gripe; o pico em 1976, no entanto, foi perdido por este modelo. Quando ajustamos um modelo com um limiar menor de: 04, as epidemias de gripe foram um tanto subestimadas, mas a epidemia de gripe no oitavo ano foi prevista um mês antes. Escolhemos o modelo com um limite de 0.05 porque o diagnóstico residual não mostrou nenhum desvio óbvio da suposição do modelo (exceto para um outlier em 1976); o modelo com um limite de 0.04 ainda tinha alguma correlação restante nos resíduos e havia mais de um outlier. Finalmente, a previsão para além de um mês de antecedência para este modelo é complicada, mas existem algumas técnicas aproximadas (ver Tong 1983).

Os seguintes comandos podem ser usados para realizar esta análise no R.

thrsh = 0.05 # limiar
dflu1 = as.numeric(dflu)
Z = ts.intersect(dflu, dplyr::lag(dflu1,1), dplyr::lag(dflu1,2), 
                 dplyr::lag(dflu1,3), dplyr::lag(dflu1,4))[-c(1:4),]
ind1 = ifelse(Z[,2] < thrsh, 1, NA) # indicador < thrsh
ind2 = ifelse(Z[,2] < thrsh, NA, 1) # indicador >= thrsh
X1 = Z[,1]*ind1
X2 = Z[,1]*ind2
summary(fit1 <- lm(X1~ Z[,2:5]) ) # case 1

## 
## Call:
## lm(formula = X1 ~ Z[, 2:5])
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.13312 -0.02049  0.00218  0.01667  0.26666 
## 
## Coefficients:
##                               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                   0.004471   0.004894   0.914 0.363032    
## Z[, 2:5]dplyr::lag(dflu1, 1)  0.506650   0.078319   6.469  3.2e-09 ***
## Z[, 2:5]dplyr::lag(dflu1, 2) -0.200086   0.056573  -3.537 0.000604 ***
## Z[, 2:5]dplyr::lag(dflu1, 3)  0.121047   0.054463   2.223 0.028389 *  
## Z[, 2:5]dplyr::lag(dflu1, 4) -0.110938   0.045979  -2.413 0.017564 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.04578 on 105 degrees of freedom
##   (17 observations deleted due to missingness)
## Multiple R-squared:  0.3763, Adjusted R-squared:  0.3526 
## F-statistic: 15.84 on 4 and 105 DF,  p-value: 3.568e-10

summary(fit2 <- lm(X2~ Z[,2:5]) ) # case 2

## 
## Call:
## lm(formula = X2 ~ Z[, 2:5])
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.089975 -0.036825 -0.006328  0.040765  0.129509 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                   0.40794    0.04675   8.726 1.53e-06 ***
## Z[, 2:5]dplyr::lag(dflu1, 1) -0.74833    0.16644  -4.496 0.000732 ***
## Z[, 2:5]dplyr::lag(dflu1, 2) -1.03231    0.21137  -4.884 0.000376 ***
## Z[, 2:5]dplyr::lag(dflu1, 3) -2.04504    1.05000  -1.948 0.075235 .  
## Z[, 2:5]dplyr::lag(dflu1, 4) -6.71178    1.24538  -5.389 0.000163 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.0721 on 12 degrees of freedom
##   (110 observations deleted due to missingness)
## Multiple R-squared:  0.9207, Adjusted R-squared:  0.8943 
## F-statistic: 34.85 on 4 and 12 DF,  p-value: 1.618e-06

D = cbind(rep(1, nrow(Z)), Z[,2:5]) # matriz de planejamento
p1 = D %*% coef(fit1) # obtendo previsões
p2 = D %*% coef(fit2)
prd = ifelse(Z[,2] < thrsh, p1, p2)
par(mar=c(4,4,2,1))
plot(dflu, ylim=c(-.5,.5), type='p', pch=3)
lines(prd)
prde1 = sqrt(sum(resid(fit1)^2)/df.residual(fit1) )
prde2 = sqrt(sum(resid(fit2)^2)/df.residual(fit2) )
prde = ifelse(Z[,2] < thrsh, prde1, prde2)
tx = time(dflu)[-(1:4)]
xx = c(tx, rev(tx))
yy = c(prd-2*prde, rev(prd+2*prde))
polygon(xx, yy, border=8, col=gray(.6, alpha=.25) )
abline(h=.05, col=4, lty=6)
grid()

Figura 9: Primeiras diferenças do número de mortes mensais por pneumonia e influenza nos EUA (+); previsões com um mês de antecedência (linha contínua) com \(\pm 2\) os limites de erro de previsão. A linha horizontal é o limiar.

Finalmente, notamos que existe um pacote R denominado tsDyn que pode ser usado para ajustar esses modelos; presumimos que o dflu já existe.

library(tsDyn) # carga do pacote
# vignette("tsDyn") # para detalhes do pacote
(u = setar(dflu, m=4, thDelay=0, th=.05)) # ajustando o modelo e visualizando resultados

## 
## Non linear autoregressive model
## 
## SETAR model ( 2 regimes)
## Coefficients:
## Low regime:
##      const.L       phiL.1       phiL.2       phiL.3       phiL.4 
##  0.004471044  0.506649694 -0.200086031  0.121047354 -0.110938271 
## 
## High regime:
##    const.H     phiH.1     phiH.2     phiH.3     phiH.4 
##  0.4079353 -0.7483325 -1.0323129 -2.0450407 -6.7117769 
## 
## Threshold:
## -Variable: Z(t) = + (1) X(t)+ (0)X(t-1)+ (0)X(t-2)+ (0)X(t-3)
## -Value: 0.05 (fixed)
## Proportion of points in low regime: 86.61%    High regime: 13.39%

(u = setar(dflu, m=4, thDelay=0)) # ajustando o modelo limiar (=.036)

## 
## Non linear autoregressive model
## 
## SETAR model ( 2 regimes)
## Coefficients:
## Low regime:
##       const.L        phiL.1        phiL.2        phiL.3        phiL.4 
##  0.0006269563  0.4608089284 -0.2243720404  0.1100931813 -0.1307031988 
## 
## High regime:
##    const.H     phiH.1     phiH.2     phiH.3     phiH.4 
##  0.2035231 -0.4071318 -1.4686776  0.3768388 -0.8298225 
## 
## Threshold:
## -Variable: Z(t) = + (1) X(t)+ (0)X(t-1)+ (0)X(t-2)+ (0)X(t-3)
## -Value: 0.03646
## Proportion of points in low regime: 84.25%    High regime: 15.75%

BIC(u); AIC(u) # se você quiser experimentar outros modelos; m = 3 funciona bem também

## [1] -678.5372

## [1] -710.1644

plot(u) # gráficos

grid()

O limite encontrado aqui é 0.036, que inclui algumas observações a mais do que usando 0.04, mas sofre das mesmas desvantagens observadas anteriormente.

5. Regressão com retardo e modelagem da função de transferência

Consideramos anteriormente a regressão defasada em uma abordagem no domínio da frequência com base na coerência. Por exemplo, considere as séries de SOI e Recrutamento; estas séries são exibidas na figura a continuação

library(astsa)
par(mar=c(2,2,2,1), mfrow=c(1,2)) 
soi = astsa::soi 
plot(soi, ylab="", xlab="", main="Índice de Oscilação do Sul")
grid()
rec = astsa::rec 
plot(rec, ylab="", xlab="", main="Recrutamento")
grid()

Nesse exemplo, o interesse estava em prever a série de Recrutamento como a saída, digamos \(Y_t\), a partir do SOI como entrada, digamos \(X_t\). Consideramos o modelo de regressão defasado como \[\begin{equation*} Y_t \, = \, \sum_{j=0}^\infty \alpha_j X_{t-j} \, + \, \eta_t \, = \, \alpha(B)X_t \, + \, \eta_t, \end{equation*}\] onde \(\sum_j |\alpha_j|<\infty\). Assumimos que o processo de entrada \(X_t\) e o processo de ruído \(\eta_t\) acima são estacionários e mutuamente independentes. Os coeficientes \(\alpha_0,\alpha_1,\cdots\) descrevem os pesos atribuídos aos valores anteriores de \(X_t\) usados na previsão de \(Y_t\) e usamos a notação \[\begin{equation*} \alpha(B) \, = \, \sum_{j=0}^\infty \alpha_j B^j\cdot \end{equation*}\]

Na formulação de Box and Jenkins (1970), atribuímos modelos ARIMA(p,d,q) e ARIMA(\(p_\eta,d_\eta,q_\eta\), às séries \(X_t\) e \(\eta_t\), respectivamente. Também assumimos que o ruído \(\eta_t\), era ruído branco. Os componentes em notação de retrocesso, para o caso simple ARMA(p,q) modelando da entrada e no ruído, teria a representação \[\begin{equation*} \phi(B)X_t \, = \, \theta(B)W_t \qquad \mbox{e} \qquad \phi_\eta(B)\eta_t \, = \, \theta_\eta(B)Z_t, \end{equation*}\] onde \(W_t\) e \(Z_t\) são processos de ruído branco independentes com variâncias \(\sigma_{_W}^2\) e \(\sigma_{_Z}^2\), respectivamente.

Box and Jenkins (1970) propuseram que padrões sistemáticos frequentemente observados nos coeficientes \(\alpha_j\); para \(j = 1,2,\cdots\), muitas vezes pode ser expressos como uma proporção de polinômios envolvendo um pequeno número de coeficientes, juntamente com um atraso especificado \(d\), então \[\begin{equation*} \alpha(B) \, = \, \dfrac{\delta(B)B^d}{\omega(B)}, \end{equation*}\] onde \[\begin{equation*} \omega(B) \, = \, 1-\omega_1 B -\omega_2 B^2 - \cdots - \omega_r B^r \end{equation*}\] e \[\begin{equation*} \delta(B) \, = \, \delta_0+\delta_1 B + \cdots + \delta_s B^s \end{equation*}\] são os operadores indicados; nesta seção, achamos conveniente representar o operador inverso, digamos \(\omega(B)^{-1}\), como \(1/\omega(B)\).

Determinar um modelo parcimonioso envolvendo uma forma simples para \(\alpha(B)\) e estimar todos os parámetros no modelo acima são as principais tarefas na metodologia da função de transferência. Devido ao grande número de parámetros, é necessário desenvolver uma metodologia sequencial.

Suponha que nos concentremos primeiro em encontrar o modelo ARIMA para a entrada \(X_t\) e aplicando esse operador a ambos os lados na primeira equação nesta seção, obtendo o novo modelo \[\begin{equation*} \widetilde{Y}_t \, = \, \dfrac{\phi(B)}{\theta(B)}Y_t \, = \, \alpha(B)\dfrac{\phi(B)}{\theta(B)}X_t \, + \, \dfrac{\phi(B)}{\theta(B)}\eta_t \, = \, \alpha(B)W_t \, \widetilde{\eta}_t, \end{equation*}\] onde \(W_t\) e o ruído transformado \(\widetilde{\eta}_t\) são independentes.

A série \(W_t\) é uma versão pré-branqueada da série de entrada e sua correlação cruzada com a série de saída transformada \(\widetilde{Y}_t\) será \[\begin{equation*} \gamma_{_{\widetilde{Y},W}}(h) \, = \, \mbox{E}(\widetilde{Y}_{t+h}W_t) \, = \, \mbox{E}\left( \sum_{j=0}^\infty \alpha_j W_{t+h-j}W_t\right) \, = \, \sigma_{_W}^2 \alpha_h, \end{equation*}\] porque a função de autocovariância do ruído branco será zero, exceto quando \(j = h\) acima. Portanto, ao calcular a correlação cruzada entre a série de entrada pré-branqueada e a série de saída transformada, deve-se obter uma estimativa grosseira do comportamento de \(\alpha(B)\).

Exemplo 8: Relacionando o SOI com a série de Recrutamento

Damos um exemplo simples do procedimento sugerido para o SOI e a série de Recrutamento. A Figura 10 abaixo mostra as funções de autocorrelação (ACF) e autocorrelação parcial (PACF) amostrais do SOI destendido e fica claro, a partir do PACF, que uma série autoregressiva com \(p = 1\) fará um trabalho razoável.

library(astsa)
soi.d = resid(lm(soi~time(soi), na.action=NULL)) # SOI destendida (sem tendência)
par(mar=c(4,4,2,1), pch=19)
acf2(soi.d)

##      [,1] [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8] [,9] [,10] [,11] [,12] [,13]
## ACF  0.59 0.35  0.18  0.01 -0.15 -0.23 -0.22 -0.14 0.01  0.19  0.33  0.38  0.28
## PACF 0.59 0.00 -0.03 -0.12 -0.16 -0.08  0.01  0.07 0.15  0.18  0.16  0.06 -0.11
##      [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25]
## ACF   0.07 -0.10 -0.22 -0.35 -0.43 -0.38 -0.24 -0.08  0.11  0.28  0.32  0.22
## PACF -0.25 -0.15 -0.06 -0.06 -0.07  0.01  0.02  0.01  0.06  0.10  0.04 -0.05
##      [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35] [,36] [,37]
## ACF   0.06 -0.07 -0.20 -0.32 -0.42 -0.36 -0.20 -0.05  0.15  0.32  0.38  0.29
## PACF -0.11 -0.04 -0.04 -0.05 -0.11 -0.01  0.04 -0.03  0.07  0.10  0.08 -0.01
##      [,38] [,39] [,40] [,41] [,42] [,43] [,44] [,45] [,46] [,47] [,48]
## ACF   0.15 -0.02 -0.16 -0.27 -0.30 -0.28 -0.16  0.04  0.20  0.37  0.40
## PACF -0.05 -0.10 -0.04  0.00  0.07 -0.01 -0.02  0.02 -0.02  0.09  0.01

Figura 10: Funções ACF e PACF da série SOI destendida.

Ajustando a série obtemos \(\widehat{\phi} = 0.588\) com \(\widehat{\sigma}_{_W}^2 = 0.092\), aplicamos o operador \(1-0.588 B\) a ambos, \(X_t\) e \(Y_t\) e calculamos a função de correlação cruzada.

fit = arima(soi.d, order=c(1,0,0))
fit

## 
## Call:
## arima(x = soi.d, order = c(1, 0, 0))
## 
## Coefficients:
##          ar1  intercept
##       0.5875     0.0008
## s.e.  0.0379     0.0344
## 
## sigma^2 estimated as 0.09181:  log likelihood = -102.1,  aic = 210.19

ar1 = as.numeric(coef(fit)[1]) # = 0.5875
soi.pw = resid(fit)
rec.fil = stats::filter(rec, filter=c(1, -ar1), sides=1)
par(mar=c(4,4,3,1), pch=19)
ccf(soi.pw, rec.fil, ylab="CCF", na.action=na.omit, panel.first=grid())
grid()

Figura 11: Função de correlação cruzada (CCF) amostral da série SOI pré-branqueada, destendida e da série de Recrutamento transformada de forma semelhante; atrasos negativos indicam que o SOI lidera o Recrutamento.

Observando a mudança aparente de \(d = 5\) meses e a diminuição daí em diante, parece plausível levantar a hipótese de um modelo da forma \[\begin{equation*} \alpha(B) \, = \, \delta_0 B^5 (1+\omega_1 B+\omega_2 B^2+ \cdots ) \, = \, \dfrac{\delta_0 B^5}{1-\omega_1 B} \end{equation*}\] para a função de transferência. Nesse caso, esperariamos que \(\omega_1\) fosse negativo.

No código acima soi.pw representa a série SOI destendida pré-branqueada e rec.fil representa a série de Recrutamento filtrada.

Em alguns casos, podemos postular a forma dos componentes separados \(\delta(B)\) e \(\omega(B)\), então podemos escrever a equação \[\begin{equation*} Y_t \, = \, \dfrac{\delta(B)B^d}{\omega(B)} X_t \, + \, \eta_t \end{equation*}\] como \[\begin{equation*} \omega(B)Y_t \, = \, \delta(B)B^d X_t \, + \, \omega(B)\eta_t, \end{equation*}\] ou em forma de modelo de regressão \[\begin{equation*} Y_t \, = \, \sum_{k=1}^r \omega_k Y_{t-k} \, + \, \sum_{k=0}^s \delta_k X_{t-d+k} \, + \, U_t, \end{equation*}\] onde \(U_t = \omega(B)\eta_t\).

Assim que tivermos \(Y_t\), podemos ajustar o modelo se esquecermos \(\eta_t\) e permitirmos que \(U_t\) tenha qualquer comportamento ARMA. Ilustramos essa técnica no exemplo a seguir.

Exemplo 9: Modelo de função de transferência para SOI e Recrutamento

Ilustramos o procedimento para ajustar um modelo de regressão defasado da forma sugerida no Exemplo 8 para a série SOI destendida \(X_t\) e a série de Recrutamento \(Y_t\).

Os resultados relatados aqui são praticamente os mesmos que os resultados obtidos a partir da abordagem no domínio da frequência. Com base no Exemplo 8, determinamos que \[\begin{equation*} Y_t \, = \, \alpha + \omega_1 Y_{t-1} +\delta_0 X_{t-5} + U_t \end{equation*}\] é um modelo razoável.

soi.d = resid(lm(soi~time(soi), na.action=NULL))
fish = ts.intersect(rec, RL1=dplyr::lag(as.numeric(rec),1), SL5=dplyr::lag(as.numeric(soi.d),5))
fish = fish[-c(1:5),]
(u = lm(fish[,1]~fish[,2:3], na.action=NULL))

## 
## Call:
## lm(formula = fish[, 1] ~ fish[, 2:3], na.action = NULL)
## 
## Coefficients:
##    (Intercept)  fish[, 2:3]RL1  fish[, 2:3]SL5  
##         8.8971          0.8556        -20.3771

acf2(resid(u)) # sugere AR1

##      [,1]  [,2]  [,3]  [,4]  [,5] [,6]  [,7] [,8]  [,9] [,10] [,11] [,12] [,13]
## ACF   0.4  0.09 -0.04 -0.15 -0.02 0.07 -0.01 0.02 -0.06 -0.10 -0.09 -0.09 -0.09
## PACF  0.4 -0.08 -0.06 -0.12  0.11 0.05 -0.09 0.04 -0.07 -0.04 -0.06 -0.03 -0.06
##      [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25]
## ACF  -0.06 -0.07 -0.07  0.01  0.01 -0.02 -0.02  0.01  0.03  0.04  0.02  0.00
## PACF -0.03 -0.04 -0.04  0.05 -0.03 -0.04 -0.03  0.05 -0.01 -0.01 -0.01 -0.02
##      [,26] [,27] [,28] [,29] [,30] [,31] [,32]
## ACF   0.03  0.08  0.06  0.08  0.05 -0.01 -0.12
## PACF  0.03  0.06  0.00  0.05  0.01 -0.01 -0.14

Figura 12: CCF amostral da série SOI pré-branqueada, destendida e da série de Recrutamento transformada de forma semelhante; atrasos negativos indicam que o SOI lidera o Recrutamento.

(arx = sarima(fish[,1], 1, 0, 0, xreg=fish[,2:3])) # final model

## initial  value 2.050589 
## iter   2 value 1.963560
## iter   3 value 1.962035
## iter   4 value 1.956727
## iter   5 value 1.956486
## iter   6 value 1.956230
## iter   7 value 1.956056
## iter   8 value 1.956027
## iter   9 value 1.956024
## iter  10 value 1.956024
## iter  10 value 1.956024
## final  value 1.956024 
## converged
## initial  value 1.955587 
## iter   2 value 1.955586
## iter   3 value 1.955585
## iter   4 value 1.955584
## iter   5 value 1.955584
## iter   6 value 1.955584
## iter   7 value 1.955584
## iter   8 value 1.955584
## iter   8 value 1.955584
## iter   8 value 1.955584
## final  value 1.955584 
## converged
## <><><><><><><><><><><><><><>
##  
## Coefficients: 
##           Estimate     SE  t.value p.value
## ar1         0.4487 0.0503   8.9183       0
## intercept  12.3323 1.5746   7.8321       0
## RL1         0.8005 0.0234  34.2778       0
## SL5       -21.0307 1.0915 -19.2674       0
## 
## sigma^2 estimated as 49.93217 on 444 degrees of freedom 
##  
## AIC = 6.771366  AICc = 6.771567  BIC = 6.817178 
## 

## $fit
## 
## Call:
## arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D, Q), period = S), 
##     xreg = xreg, transform.pars = trans, fixed = fixed, optim.control = list(trace = trc, 
##         REPORT = 1, reltol = tol))
## 
## Coefficients:
##          ar1  intercept     RL1       SL5
##       0.4487    12.3323  0.8005  -21.0307
## s.e.  0.0503     1.5746  0.0234    1.0915
## 
## sigma^2 estimated as 49.93:  log likelihood = -1511.79,  aic = 3033.57
## 
## $degrees_of_freedom
## [1] 444
## 
## $ttable
##           Estimate     SE  t.value p.value
## ar1         0.4487 0.0503   8.9183       0
## intercept  12.3323 1.5746   7.8321       0
## RL1         0.8005 0.0234  34.2778       0
## SL5       -21.0307 1.0915 -19.2674       0
## 
## $ICs
##      AIC     AICc      BIC 
## 6.771366 6.771567 6.817178

Com base nessas técnicas, o modelo ajustado é o mesmo obtido no dom[inio da frequência, a saber, \[\begin{equation*} Y_t \, = \, 12.3323 + 0.8005 Y_{t-1}-21.0307 X_{t-5} +U_t \qquad \mbox{e} \qquad U_t \, = \, 0.4487 U_{t-1}+W_t, \end{equation*}\] onde \(W_t\) é um ruído branco com \(\widehat{\sigma}_{_W}^2=49.93\).

A Figura 12 mostra as funções ACF e PACF amostrais do ruído \(U_t\) estimado, mostrando que um AR(1) é apropriado. Além disso, a Figura 13 exibe a série de Recrutamento e as previsõ de um passo à frente com base no modelo final.

pred = fish[,1] + resid(arx$fit) # previsões um-passo-à-frente
ts.plot(pred, fish[,1], col=c('gray90',1), lwd=c(7,1))
grid()

Figura 13: A série de recrutamento (linha) e as previsões de um passo à frente (amostra cinza) com base no modelo de função de transferência final.

Para completar, concluímos a discussão do método mais complicado de Box-Jenkins para ajustar modelos de função de transferência. Observamos, entretanto, que o método não tem otimização geral reconhecível e geralmente não é melhor ou pior do que o método discutido anteriormente.

A forma de \(Y_t\) sugere fazer uma regressão nas versões defasadas de ambas as séries de entrada e saída para obter \(\widehat{\beta}\), a estimativa do vetor de regressão \((r + s + 1)\times 1\) \[\begin{equation*} \beta \, = \, (\omega_1,\cdots,\omega_r,\delta_0,\delta_1,\cdots,\delta_s)^\top \cdot \end{equation*}\] Os resíduos da regressão, digamos, \[\begin{equation*} \widehat{U}_t \, = \, Y_t - \widehat{\beta}^\top Z_t, \end{equation*}\] onde \[\begin{equation*} Z_t \, = \, (Y_{t-1},\cdots,Y_{t-r},X_{t-d},\cdots,X_{t-d-s})^\top \end{equation*}\] denota o vetor usual de variáveis independentes, poderia ser usado para aproximar o melhor modelo ARMA para o processo de ruído \(\eta_t\), porque podemos calcular um estimador para esse processo usando \(\widehat{U}_t\) e \(\widehat{\omega}(B)\) e aplicando o operador de média móvel para obter \(\widehat{\eta}_t\).

Ajustando um modelo ARMA(\(p_\eta,q_\eta\)) para este ruído estimado então completa a especificação. O anterior sugere o seguinte procedimento sequencial para ajustar o modelo da função de transferência aos dados.

1. Ajuste um modelo ARMA à série de entrada \(X_t\) para estimar os parámetros \(\phi_1,\cdots,\phi_p,\theta_1\cdots,\theta_q,\sigma_{_W}^2\) na utilizando a expressão \[\begin{equation*} \phi(B)X_t \, = \, \theta(B)W_t\cdot \end{equation*}\] Retenha os coeficientes ARMA para uso na etapa (ii) e os resíduos ajustados \(\widehat{W}_t\) para uso na etapa (iii).

2. Aplicar o operador determinado na etapa (i), ou seja, \[\begin{equation*} \widehat{\phi}(B)Y_t \, = \, \widehat{\theta}(B)\widetilde{Y}_t, \end{equation*}\] para determinar a série de saída transformada \(\widetilde{Y}_t\).

3. Use a função de correlação cruzada entre \(\widetilde{Y}_t\) e \(\widehat{W}_t\) nas etapas (i) e (ii) para sugerir uma forma para os componentes do polinômio \[\begin{equation*} \alpha(B) \, = \, \dfrac{\delta(B)B^d}{\omega(B)} \end{equation*}\] e estimar o atraso de tempo \(d\).

4. Obter \(\widehat{\beta}=(\widehat{\omega}_1,\cdots,\widehat{\omega}_r,\widehat{\delta}_1,\cdots,\widehat{\delta}_s)\) ajustando uma regressão linear da forma \[\begin{equation*} Y_t \, = \, \sum_{k=1}^r \omega_k Y_{t-k} \, + \, \sum_{k=0}^s \delta_k X_{t-d+k} \, + \, U_t\cdot \end{equation*}\] Retenha os resíduos \(\widehat{U}_t\) para uso na etapa (v).

5. Aplique a transformação de médias móveis \(U_t = \omega(B)\eta_t\) aos resíduos \(\widehat{U}_t\) para encontrar a série de ruído \(\widehat{\eta}_t\), e ajuste um modelo ARMA ao ruído, obtendo os coeficientes estimados em \(\widehat{\phi}_\eta\) e \(\widehat{\theta}_\eta\).

O procedimento acima é bastante razoável mas, como mencionado anteriormente, não é ideal em nenhum sentido. A estimação simultânea de mínimos quadrados, com base em \(X_t\) e \(Y_t\) observados, pode ser realizada observando que o modelo de função de transferência pode ser escrito como \[\begin{equation*} Y_t \, = \, \dfrac{\delta(B)B^d}{\omega(B)} X_t \, + \, \dfrac{\theta_\eta(B)}{\phi_\eta(B)}Z_t, \end{equation*}\] que pode ser colocado na forma \[\begin{equation*} \omega(B)\phi_\eta(B)Y_t \, = \, \phi_\eta(B)\delta(B)B^d X_t \, + \, \omega(B)\theta_\eta(B)Z_t \end{equation*}\] e é claro, podemos usar mínimos quadrados para minimizar \(\sum_t Z_t^2\), como nas seções anteriores.

No Exemplo 9, simplesmente permitimos que \[\begin{equation*} U_t\dfrac{\theta_\eta(B)}{\phi_\eta(B)}Z_t \end{equation*}\] acima tivesse qualquer estrutura ARMA. Finalmente, mencionamos que também podemos expressar a função de transferência na forma de espaço de estados como um modelo ARMAX; consulte a Seção 6 a seguir.

6. Modelo ARMAX: modelo ARMA multivariado

Para compreender os modelos de séries temporais multivariadas e suas capacidades, apresentamos primeiro uma introdução às técnicas de regressão de séries temporais multivariadas. Como todos os processos são vetoriais, suspendemos o uso de negrito para vetores. Uma extensão útil do modelo básico de regressão univariada seria o caso em que temos mais de uma série de resultados, isto é, análise de regressão multivariada.

Suponha, em vez de uma única variável resposta \(Y_t\), existe uma coleção de \(k\) variáveis resposta \(Y_{t1},Y_{t2},\cdots,Y_{tk}\) as quais estão relacionados às entradas como \[\begin{equation*} Y_{ti} \, = \, \beta_{i1}Z_{t1} + \beta_{i2} Z_{t2}+\cdots + \beta_{ir}Z_{tr}+W_{ti} \end{equation*}\] para cada uma das \(i = 1,2,\cdots,k\) variáveis de saída.

Assumimos que as variáveis \(W_{ti}\) são correlacionadas sobre o identificador de variável \(i\), mas ainda são independentes ao longo do tempo. Formalmente, assumimos \(\mbox{Cov}(W_{si},W_{tj}) = \sigma_{ij}\) para \(s = t\) e é zero caso contrário. Então, escrevendo a expressão acima em notação matricial, com \(Y_t = (Y_{t1},Y_{t2},\cdots,Y_{tk})^\top\) sendo o vetor de respostas ou saídas e \(B = \{\beta_{ij}\}\), \(i=1,\cdots,k\) e \(j=1,\cdots,r\) sendo uma matriz \(k\times r\) contendo os coeficientes de regressão, leva à forma de aparência simples \[\begin{equation*} Y_t \, = \, BZ_t \, + \, W_t\cdot \end{equation*}\]

Aqui, o processo \(k\times 1\) vetorial \(W_t\) é assumido como uma coleção de vetores independentes com matriz de covariâncias comum \(\mbox{E}(W_t W_t^\top)=\Sigma_{_W}\), a matriz \(k\times k\) contendo as covariâncias \(\sigma_{ij}\). Sob a suposição de normalidade, o estimador de máxima verossimilhança para a matriz de regressão é \[\begin{equation*} \widehat{B} \, = \, \left(\sum_{t=1}^n y_tz_t^\top \right)\left(\sum_{t=1}^n z_tz_t^\top \right)^{-1}\cdot \end{equation*}\]

A matriz de covariâncias do erro \(\Sigma_{_W}\) é estimada por \[\begin{equation*} \widehat{\Sigma}_{_W} \, = \, \dfrac{1}{n-r}\sum_{t=1}^n (y_t-\widehat{B}z_t)(y_t-\widehat{B}z_t)^\top\cdot \end{equation*}\]

A incerteza nos estimadores pode ser avaliada a partir de \[\begin{equation*} se\big( \widehat{\beta}_{ij}\big) \, = \, \sqrt{c_{ii}\widehat{\sigma}_{jj}}, \end{equation*}\] para \(i=1,\cdots,r\) e \(j=1,\cdots,k\), onde se denota o erro padrão estimado \(\widehat{\sigma}_{jj}\), o \(j\)-ésimo elemento da diagonal de \(\widehat{\Sigma}_{_W}\) e \(c_{ii}\) é o \(i\)-ésimo elemento da diagonal de \((\sum_{t=1}^n z_tz_t^\top)^{-1}\).

Além disso, o Critéio de Informação e Akaike muda para \[\begin{equation*} AIC \, = \, \ln\Big(|\widehat{\Sigma}_{_W}|\Big) \, + \, \dfrac{2}{n}\left(kr+\dfrac{k(k+1)}{2}\right) \end{equation*}\] e o \(BIC\) substitui o segundo termo da expressão acima por \(K\ln(n)/n\) onde \(K = kr + k(k + 1)/2\). Bedrick and Tsai (1994) forneceram uma expressão para o AIC corrigido no caso multivariado como \[\begin{equation*} AICc \, = \, \ln\Big(|\widehat{\Sigma}_{_W}|\Big) \, + \, \dfrac{k(r+n)}{n-k-r-1}\cdot \end{equation*}\]

Muitos conjuntos de dados envolvem mais de uma série temporal e muitas vezes estamos interessados nas possíveis dinâmicas relacionadas a todas as séries. Nesta situação, estamos interessados em modelar e prever \(k\times 1\) séries temporais com valor vetorial \(X_t=(X_{t1},\cdots,X_{tk})^\top\), \(t=0,\pm 1, \pm 2,\cdots\).

Infelizmente, estender modelos ARMA univariados para o caso multivariado não é tão simples. O modelo multivariado autorregressivo, entretanto, constitui uma extensão direta do modelo AR univariado.

Para o modelo autoregressivo vetorial de primeira ordem VAR(1), temos \[\begin{equation*} X_t \, = \, \alpha +\Phi X_{t-1}+W_t, \end{equation*}\] onde \(\Phi\) é uma matriz \(k\times k\) de transição que expressa a dependência de \(X_t\) em \(X_{t-1}\). O processo de ruído branco vetorial \(W_t\) é considerado normal multivariado com média zero e matriz de covariâncias \[\begin{equation*} \mbox{E}\big(W_t W_t^\top \big) \, = \, \Sigma_{_W}\cdot \end{equation*}\]

O vetor \(\alpha=(\alpha_1,\alpha_2,\cdots,\alpha_k)^\top\) aparece como a constante na configuração da regressão. Se \(\mbox{E}(X_t)=\mu\), então \(\alpha=(\mbox{I}-\Phi)\mu\).

Observe a semelhança entre o modelo VAR e o modelo de regressão linear multivariado anteriormente apresentado, nesta mesma seção. As fórmulas de regressão continuam e podemos, ao observar \(x_1,\cdots,x_n\), configurar o modelo VAR(1) com \(y_t = x_t\), \(B = (\alpha,\Phi)\) e \(z_t = (1,x_{t-1}^\top)^\top\). Então, a estimativa de máxima verossimilhança condicional para a matriz de covariâncias é dada por \[\begin{equation*} \widehat{\Sigma}_{_W} \, = \, \dfrac{1}{n-1}\sum_{t=2}^n (x_t-\widehat{\alpha}-\widehat{\Phi}x_{t-1})(x_t-\widehat{\alpha}-\widehat{\Phi}x_{t-1})^\top\cdot \end{equation*}\]

A forma especial assumida para a componente constante \(\alpha\), do modelo vetorial AR pode ser generalizada para incluir um vetor \(r\times 1\) fixo de entradas \(U_t\). Ou seja, poderíamos ter proposto o modelo vetorial ARX, \[\begin{equation*} X_t \, = \, \Gamma U_t \, + \, \sum_{j=1}^P \Phi_j X_{t-j} \, + \, W_t, \end{equation*}\] onde \(\Gamma\) seria a matriz \(p\times r\) de parâmetros.

O X em ARX refere-se ao processo vetorial exôgeno que denotamos aqui por \(U_t\). A introdução de variáveis exôgenas através de substituir \(\alpha\) por \(\Gamma U_t\) não apresenta nenhum problema especial em fazer inferências e frequentemente deixaremos o X de lado por ser supérfluo.

Exemplo 10. Poluição, clima e mortalidade.

Por exemplo, para a série tridimensional composta de mortalidade cardiovascular \(X_{t1}\), temperatura \(X_{t2}\) e níveis de partículas \(X_{t3}\), definamos \(X_t = (X_{t1}, X_{t2}, X_{t3})\) como um vetor de dimensão \(k = 3\). Podemos imaginar relações dinâmicas entre as três séries definidas como a relação de primeira ordem, \[\begin{equation*} X_{t1} \, = \, \alpha_1 +\beta_1 t +\phi_{11} X_{t-1,1} + \phi_{12}X_{t-1,2} + \phi_{13} X_{t-1,3} +W_{t1}, \end{equation*}\] que expressa o valor atual da mortalidade como uma combinação linear de tendência e seu valor passado imediato e os valores anteriores de temperatura e níveis de partículas.

Similarmente, \[\begin{equation*} X_{t2} \, = \, \alpha_2 +\beta_2 t +\phi_{21} X_{t-1,1} + \phi_{22}X_{t-1,2} + \phi_{23} X_{t-1,3} +W_{t2} \end{equation*}\] e \[\begin{equation*} X_{t3} \, = \, \alpha_3 +\beta_3 t +\phi_{31} X_{t-1,1} + \phi_{32}X_{t-1,2} + \phi_{33} X_{t-1,3} +W_{t3}, \end{equation*}\] expressam a dependência da temperatura e dos níveis de partículas nas outras séries. Claro, existem métodos para a identificação preliminar desses modelos e discutiremos esses métodos em breve.

O modelo ARX acima é \[\begin{equation*} X_t \, = \, \Gamma U_t + \Phi X_{t-1} + W_t, \end{equation*}\] onde \(\Gamma = (\alpha | \beta)\) é \(3\times 2\) e \(U_t = (1,t)^\top\) é \(2\times 1\).

Ao longo de grande parte desta seção, usaremos o pacote R vars para ajustar modelos AR vetoriais por meio de mínimos quadrados.

library(astsa)
library(vars)
x = cbind(cmort, tempr, part)
ajuste = VAR(x, p=1, type='both') # 'both' ajusta constante + tendência
summary(ajuste)

## 
## VAR Estimation Results:
## ========================= 
## Endogenous variables: cmort, tempr, part 
## Deterministic variables: both 
## Sample size: 507 
## Log Likelihood: -5116.02 
## Roots of the characteristic polynomial:
## 0.8931 0.4953 0.1444
## Call:
## VAR(y = x, p = 1, type = "both")
## 
## 
## Estimation results for equation cmort: 
## ====================================== 
## cmort = cmort.l1 + tempr.l1 + part.l1 + const + trend 
## 
##           Estimate Std. Error t value Pr(>|t|)    
## cmort.l1  0.464824   0.036729  12.656  < 2e-16 ***
## tempr.l1 -0.360888   0.032188 -11.212  < 2e-16 ***
## part.l1   0.099415   0.019178   5.184 3.16e-07 ***
## const    73.227292   4.834004  15.148  < 2e-16 ***
## trend    -0.014459   0.001978  -7.308 1.07e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 5.583 on 502 degrees of freedom
## Multiple R-Squared: 0.6908,  Adjusted R-squared: 0.6883 
## F-statistic: 280.3 on 4 and 502 DF,  p-value: < 2.2e-16 
## 
## 
## Estimation results for equation tempr: 
## ====================================== 
## tempr = cmort.l1 + tempr.l1 + part.l1 + const + trend 
## 
##           Estimate Std. Error t value Pr(>|t|)    
## cmort.l1 -0.244046   0.042105  -5.796 1.20e-08 ***
## tempr.l1  0.486596   0.036899  13.187  < 2e-16 ***
## part.l1  -0.127661   0.021985  -5.807 1.13e-08 ***
## const    67.585598   5.541550  12.196  < 2e-16 ***
## trend    -0.006912   0.002268  -3.048  0.00243 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 6.4 on 502 degrees of freedom
## Multiple R-Squared: 0.5007,  Adjusted R-squared: 0.4967 
## F-statistic: 125.9 on 4 and 502 DF,  p-value: < 2.2e-16 
## 
## 
## Estimation results for equation part: 
## ===================================== 
## part = cmort.l1 + tempr.l1 + part.l1 + const + trend 
## 
##           Estimate Std. Error t value Pr(>|t|)    
## cmort.l1 -0.124775   0.079013  -1.579    0.115    
## tempr.l1 -0.476526   0.069245  -6.882 1.77e-11 ***
## part.l1   0.581308   0.041257  14.090  < 2e-16 ***
## const    67.463501  10.399163   6.487 2.10e-10 ***
## trend    -0.004650   0.004256  -1.093    0.275    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 12.01 on 502 degrees of freedom
## Multiple R-Squared: 0.3732,  Adjusted R-squared: 0.3683 
## F-statistic: 74.74 on 4 and 502 DF,  p-value: < 2.2e-16 
## 
## 
## 
## Covariance matrix of residuals:
##        cmort  tempr   part
## cmort 31.172  5.975  16.65
## tempr  5.975 40.965  42.32
## part  16.654 42.323 144.26
## 
## Correlation matrix of residuals:
##        cmort  tempr   part
## cmort 1.0000 0.1672 0.2484
## tempr 0.1672 1.0000 0.5506
## part  0.2484 0.5506 1.0000

Para este caso particular, obtemos \[\begin{equation} \widehat{\alpha} \, = \, (73.23, 67.59, 67.46)^\top, \qquad \widehat{\beta} \, = \, (-0.014, -0.007, -0.005)^\top, \end{equation}\] \[\begin{equation*} \widehat{\Phi} \, = \, \begin{pmatrix} 0.46_{(0.04)} & -0.36_{(0.03)} & 0.10_{(0.02)} \\ -0.24_{(0.04)} & 0.49_{(0.04)} & -0.13_{(0.02)} \\ -0.12_{(0.08)} & -0.48_{(0.07)} & 0.58_{(0.04)} \end{pmatrix} \qquad \mbox{e} \qquad \widehat{\Sigma} \, = \, \begin{pmatrix} 31.17 & 5.98 & 16.65 \\ 5.98 & 40.96 & 42.32 \\ 16.65 & 52.32 & 144.26 \end{pmatrix}, \end{equation*}\] onde os erro padrão são dados entre parênteses.

plot(ajuste)

Para o vetor \(X_{t1},X_{t2},X_{t3}) = (M_t,T_t,P_t)\), com \(M_t\), \(T_t\) e \(P_t\) denotando mortalidade, temperatura e nível de particulado, respectivamente, obtemos a equação de predição para mortalidade, \[\begin{equation*} \widehat{M}_t \, = \, 73.23 -0.014 t + 0.46 M_{t-1} -0.36 T_{t-1}+0.10 P_{t-1}\cdot \end{equation*}\] Comparando a mortalidade observada e prevista com este modelo leva a um \(R^2\) de cerca de 0.69.

Pode-se estender o processo VAR(1) para ordens superiores, VAR(p). Para fazer isso, escrevemos o vetor de regressores como \[\begin{equation*} Z_t \, = \, (1, X_{t-1}^\top, X_{t-2}^\top,,\cdots,X_{t-p}^\top)^\top \end{equation*}\] e a matriz de regressão como \(B = (\alpha,\Phi_1,\Phi_2,\cdots,\Phi_p)\). Então, este modelo de regressão pode ser escrito como \[\begin{equation*} X_t \, = \, \alpha + \sum_{j=1}^p \Phi_j X_{t-j}+W_t, \end{equation*}\] para \(t=p+1,\cdots,n\).

A matriz \(k\times k\) de somas de produtos do erro torna-se \[\begin{equation*} SSE \, = \, \sum_{t=p+1}^n (x_t-Bz_t)(x_t-Bz_t)^\top, \end{equation*}\] de modo que o estimativa de máxima verossimilhança condicional para a matriz de covariâncias de erro \(\Sigma_{_W}\) é \[\begin{equation*} \widehat{\Sigma}_{_W} \, = \, \dfrac{SSE}{n-p}, \end{equation*}\] como no caso da regressão multivariada, exceto que agora apenas existem \(n-p\) resíduos.

Para o caso multivariado, descobrimos que o Critério de Schwarz assume a forma \[\begin{equation*} BIC \, = \, \ln\Big(|\widehat{\Sigma}_{_W}|\Big) \, + \, \dfrac{k^2p}{n}\ln(n), \end{equation*}\] fornece classificações mais razoáveis do que o AIC ou AICc, a versão corrigida. O resultado é consistente com os relatados nas simulações de Lütkepohl (1993).

Exemplo 11: Poluição, clima e mortalidade (continuação)

Usamos o pacote R vars primeiro para selecionar um modelo VAR(p) e depois ajustar o modelo. Os critérios de seleção usados no pacote são AIC, Hannan-Quinn (HQ; Hannan and Quinn (1979), BIC (SC) e Erro de Previsão Final (FPE).

O procedimento de Hannan-Quinn é semelhante ao BIC, mas com \(\ln(n)\) substituído por \(2\ln\big(\ln(n)\big)\) no termo de penalização. O FPE encontra o modelo que minimiza o erro de predição um passo à frente do quadrado médio aproximado (ver Akaike 1969 para detalhes); raramente é usado.

VARselect(x, lag.max=10, type="both")

## $selection
## AIC(n)  HQ(n)  SC(n) FPE(n) 
##      9      5      2      9 
## 
## $criteria
##                   1           2           3           4           5           6
## AIC(n)     11.73780    11.30185    11.26788    11.23030    11.17634    11.15266
## HQ(n)      11.78758    11.38149    11.37738    11.36967    11.34557    11.35176
## SC(n)      11.86463    11.50477    11.54689    11.58541    11.60755    11.65996
## FPE(n) 125216.91717 80972.28678 78268.19568 75383.73647 71426.10041 69758.25113
##                  7           8           9          10
## AIC(n)    11.15247    11.12878    11.11915    11.12019
## HQ(n)     11.38144    11.38760    11.40784    11.43874
## SC(n)     11.73587    11.78827    11.85473    11.93187
## FPE(n) 69749.89175 68122.40518 67476.96374 67556.45243

Observe que o BIC escolhe o modelo de ordem \(p = 2\), enquanto os critérios AIC e FPE selecionam um modelo de ordem \(p = 9\) e o critério Hannan-Quinn seleciona um modelo de ordem \(p = 5\).

Ajustando o modelo selecionado pelo BIC, obtemos \[\begin{equation} \widehat{\alpha}=(56.10,49.88,59.58), \qquad \widehat{\beta}=(-0.011,-0.005,-0.008), \end{equation}\] \[\begin{array}{c} \widehat{\Phi}_1 \, = \, \begin{pmatrix} 0.30_{(0.04)} & -0.20_{(0.04)} & 0.04_{(0.02)} \\ -0.11_{(0.05)} & 0.26_{(0.05)} & -0.05_{(0.03)} \\ 0.08_{(0.09)} & -0.39_{(0.09)} & 0.39_{(0.05)} \end{pmatrix}, \\ \widehat{\Phi}_2 \, = \, \begin{pmatrix} 0.28_{(0.04)} & -0.08_{(0.03)} & 0.07_{(0.03)} \\ -0.04_{(0.05)} & 0.36_{(0.05)} & -0.10_{(0.03)} \\ -0.33_{(0.09)} & 0.05_{(0.09)} & 0.38_{(0.05)} \end{pmatrix}, \end{array}\]

sendo que os erros padrão são fornecidos entre parênteses. A estimative de \(\Sigma_{_W}\) é \[\begin{equation*} \widehat{\Sigma}_{_W} \, = \, \begin{pmatrix} 28.03 & 7.08 & 16.33 \\ 7.08 & 37.63 & 40.88 \\ 16.33 & 40.88 & 123.45 \end{pmatrix}\cdot \end{equation*}\]

summary(ajuste.VAR <- VAR(x, p=2, type="both")) 

## 
## VAR Estimation Results:
## ========================= 
## Endogenous variables: cmort, tempr, part 
## Deterministic variables: both 
## Sample size: 506 
## Log Likelihood: -4987.186 
## Roots of the characteristic polynomial:
## 0.8807 0.8807 0.5466 0.4746 0.4746 0.4498
## Call:
## VAR(y = x, p = 2, type = "both")
## 
## 
## Estimation results for equation cmort: 
## ====================================== 
## cmort = cmort.l1 + tempr.l1 + part.l1 + cmort.l2 + tempr.l2 + part.l2 + const + trend 
## 
##           Estimate Std. Error t value Pr(>|t|)    
## cmort.l1  0.297059   0.043734   6.792 3.15e-11 ***
## tempr.l1 -0.199510   0.044274  -4.506 8.23e-06 ***
## part.l1   0.042523   0.024034   1.769  0.07745 .  
## cmort.l2  0.276194   0.041938   6.586 1.15e-10 ***
## tempr.l2 -0.079337   0.044679  -1.776  0.07639 .  
## part.l2   0.068082   0.025286   2.692  0.00733 ** 
## const    56.098652   5.916618   9.482  < 2e-16 ***
## trend    -0.011042   0.001992  -5.543 4.84e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 5.295 on 498 degrees of freedom
## Multiple R-Squared: 0.7227,  Adjusted R-squared: 0.7188 
## F-statistic: 185.4 on 7 and 498 DF,  p-value: < 2.2e-16 
## 
## 
## Estimation results for equation tempr: 
## ====================================== 
## tempr = cmort.l1 + tempr.l1 + part.l1 + cmort.l2 + tempr.l2 + part.l2 + const + trend 
## 
##           Estimate Std. Error t value Pr(>|t|)    
## cmort.l1 -0.108889   0.050667  -2.149  0.03211 *  
## tempr.l1  0.260963   0.051292   5.088 5.14e-07 ***
## part.l1  -0.050542   0.027844  -1.815  0.07010 .  
## cmort.l2 -0.040870   0.048587  -0.841  0.40065    
## tempr.l2  0.355592   0.051762   6.870 1.93e-11 ***
## part.l2  -0.095114   0.029295  -3.247  0.00125 ** 
## const    49.880485   6.854540   7.277 1.34e-12 ***
## trend    -0.004754   0.002308  -2.060  0.03993 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 6.134 on 498 degrees of freedom
## Multiple R-Squared: 0.5445,  Adjusted R-squared: 0.5381 
## F-statistic: 85.04 on 7 and 498 DF,  p-value: < 2.2e-16 
## 
## 
## Estimation results for equation part: 
## ===================================== 
## part = cmort.l1 + tempr.l1 + part.l1 + cmort.l2 + tempr.l2 + part.l2 + const + trend 
## 
##           Estimate Std. Error t value Pr(>|t|)    
## cmort.l1  0.078934   0.091773   0.860 0.390153    
## tempr.l1 -0.388808   0.092906  -4.185 3.37e-05 ***
## part.l1   0.388814   0.050433   7.709 6.92e-14 ***
## cmort.l2 -0.325112   0.088005  -3.694 0.000245 ***
## tempr.l2  0.052780   0.093756   0.563 0.573724    
## part.l2   0.382193   0.053062   7.203 2.19e-12 ***
## const    59.586169  12.415669   4.799 2.11e-06 ***
## trend    -0.007582   0.004180  -1.814 0.070328 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 11.11 on 498 degrees of freedom
## Multiple R-Squared: 0.4679,  Adjusted R-squared: 0.4604 
## F-statistic: 62.57 on 7 and 498 DF,  p-value: < 2.2e-16 
## 
## 
## 
## Covariance matrix of residuals:
##        cmort  tempr   part
## cmort 28.034  7.076  16.33
## tempr  7.076 37.627  40.88
## part  16.325 40.880 123.45
## 
## Correlation matrix of residuals:
##        cmort  tempr   part
## cmort 1.0000 0.2179 0.2775
## tempr 0.2179 1.0000 0.5998
## part  0.2775 0.5998 1.0000

Usando a notação do exemplo anterior, o modelo de predição para a mortalidade cardiovascular é estimado como \[\begin{equation*} \widehat{M}_t \, = \, 56-0.11t +0.3 M_{t-1}-0.2 T_{t-1} +0.04 P_{t-1} +0.28 M_{t-2} -0.08 T_{t-2} +0.07 P_{t-2}\cdot \end{equation*}\]

Para examinar os resíduos, podemos representar graficamente as correlações cruzadas dos resíduos e examinar a versão multivariada do teste Q da seguinte forma:

acf(resid(ajuste.VAR), 52)

serial.test(ajuste.VAR, lags.pt=12, type="PT.adjusted")

## 
##  Portmanteau Test (adjusted)
## 
## data:  Residuals of VAR object ajuste.VAR
## Chi-squared = 162.35, df = 90, p-value = 4.602e-06

Figura 14: Na diagonal, funções de autocorrelação (ACFs) e fora das diagonais, funções de correlação cruzada para os resíduos do VAR(2) tridimensional, ajustado ao conjunto de dados de poluição, clima e mortalidade em LA. Fora da diagonal, a série com o segundo nome é aquela que lidera.

A matriz de correlação cruzada é mostrada numa saída acima. A Figura 14 mostra os ACFs das séries de resíduos individuais ao longo da diagonal. Por exemplo, o primeiro gráfico diagonal é o ACF de \(M_t -\widehat{M}_t\) e assim por diante. As diagonais externas exibem os CCFs entre pares de séries de resíduos.

Se o título do gráfico fora da diagonal for \(x \, \& \, y\), então \(y\) lidera o gráfico; isto é, na diagonal superior, o gráfico mostra \(\mbox{Corr}\big(x(t+Lag), y(t)\big)\), enquanto na diagonal inferior, se o título for \(x \, \& \, y\), você obtêm um gráfico de \(\mbox{Corr}\big(x( t + Lag), y(t)\big)\) sim! é a mesma coisa, mas as defasagens são negativas na diagonal inferior. O gráfico está rotulado de forma estranha, basta lembrar que a segunda série com nome é a que lidera.

Na Figura 14, notamos que a maioria das correlações nas séries de resíduos são desprezíveis, no entanto, a correlação de ordem zero da mortalidade com os resíduos de temperatura é de cerca de 0.22 e a mortalidade com resíduos do material particulado é de cerca de 0.28.

Utilize o comando a seguir para visualizar os valores exatos.

head(acf(resid(ajuste.VAR), 52, plot = FALSE)$acf)

## , , 1
## 
##             [,1]        [,2]        [,3]
## [1,] 1.000000000 0.217879438 0.277509279
## [2,] 0.019047975 0.005389606 0.005820489
## [3,] 0.019666570 0.029266506 0.057572004
## [4,] 0.001875098 0.113509001 0.072168281
## [5,] 0.037066723 0.089215447 0.089481560
## [6,] 0.023449770 0.006534372 0.041550540
## 
## , , 2
## 
##              [,1]         [,2]         [,3]
## [1,]  0.217879438  1.000000000  0.599821673
## [2,] -0.006913463 -0.027315346 -0.020280582
## [3,] -0.041739638 -0.073264190 -0.057460803
## [4,] -0.036525178 -0.019103187 -0.005234533
## [5,]  0.122311098  0.099738803 -0.037019996
## [6,] -0.042626752  0.001307769  0.072799848
## 
## , , 3
## 
##              [,1]         [,2]        [,3]
## [1,]  0.277509279  0.599821673  1.00000000
## [2,] -0.005606983 -0.006272803 -0.07188693
## [3,] -0.050471197 -0.029277742 -0.10440933
## [4,] -0.051476471  0.024355781  0.08185922
## [5,]  0.108826391  0.115183087  0.10143172
## [6,] -0.022393526 -0.062856200  0.03997762

Isso significa que o modelo AR não está capturando o efeito simultâneo da temperatura e da poluição sobre a mortalidade, lembre-se de que os dados evoluem ao longo de uma semana. É possível ajustar modelos simultâneos; veja Reinsel (1997) para mais detalhes. Assim, não inesperadamente, o teste Q rejeita a hipótese nula de que o ruído seja branco.

A estatística do teste Q é dada por \[\begin{equation*} Q \, = \, n^2\sum_{h=1}^H \dfrac{1}{n-h} \mbox{tr}\left( \widehat{\Gamma}_{_W}(h){\widehat{\Gamma}_{_W}(0)}^{-1}\widehat{\Gamma}_{_W}(h){\widehat{\Gamma}_{_W}(0)}^{-1} \right), \end{equation*}\] onde \[\begin{equation*} \widehat{\Gamma}_{_W}(h) \, = \, \dfrac{1}{n}\sum_{t=1}^{n-h} \widehat{W}_{t+h} \widehat{W}_t^\top, \end{equation*}\] e \(\widehat{W}_{t}\) é o processo residual. Sob a hipóteses nula de que \(W_t\) seja ruído branco, a estatístia Q tem distribuição assintótica \(\chi^2\) com \(k^2(H-p)\) graus de liberdade.

Finalmente, a previsão segue de uma maneira direta do caso univariado. Usando o pacote R vars, use o comando predict e o comando fanchart, que produz um bom gráfico:

(fit.pr = predict(ajuste.VAR, n.ahead = 24, ci = 0.95)) # 24 semanas à frente

## $cmort
##           fcst    lower     upper       CI
##  [1,] 87.26921 76.89173  97.64668 10.37748
##  [2,] 87.02842 76.09665  97.96020 10.93178
##  [3,] 87.29573 75.40528  99.18618 11.89045
##  [4,] 87.32528 74.90938  99.74119 12.41591
##  [5,] 87.28237 74.33199 100.23275 12.95038
##  [6,] 87.17621 73.74964 100.60277 13.42657
##  [7,] 87.00116 73.11597 100.88634 13.88519
##  [8,] 86.78430 72.46755 101.10104 14.31674
##  [9,] 86.53258 71.81185 101.25330 14.72073
## [10,] 86.25716 71.16321 101.35111 15.09395
## [11,] 85.96792 70.53387 101.40197 15.43405
## [12,] 85.67068 69.93014 101.41121 15.74053
## [13,] 85.37231 69.35882 101.38580 16.01349
## [14,] 85.07662 68.82233 101.33091 16.25429
## [15,] 84.78756 68.32281 101.25230 16.46474
## [16,] 84.50751 67.86030 101.15472 16.64721
## [17,] 84.23850 67.43429 101.04272 16.80422
## [18,] 83.98179 67.04340 100.92019 16.93840
## [19,] 83.73823 66.68589 100.79057 17.05234
## [20,] 83.50824 66.35971 100.65678 17.14853
## [21,] 83.29196 66.06268 100.52125 17.22929
## [22,] 83.08928 65.79254 100.38601 17.29673
## [23,] 82.89988 65.54709 100.25266 17.35279
## [24,] 82.72333 65.32418 100.12248 17.39915
## 
## $tempr
##           fcst    lower    upper       CI
##  [1,] 70.33437 58.31182 82.35692 12.02255
##  [2,] 69.09470 56.77089 81.41851 12.32381
##  [3,] 69.22457 56.09460 82.35454 13.12997
##  [4,] 68.83591 55.32996 82.34186 13.50595
##  [5,] 69.01016 55.04606 82.97425 13.96409
##  [6,] 69.02784 54.66842 83.38726 14.35942
##  [7,] 69.24058 54.48088 84.00028 14.75970
##  [8,] 69.42757 54.29875 84.55638 15.12882
##  [9,] 69.67548 54.19959 85.15137 15.47589
## [10,] 69.92480 54.13266 85.71694 15.79214
## [11,] 70.18945 54.11108 86.26782 16.07837
## [12,] 70.45350 54.11980 86.78720 16.33370
## [13,] 70.71606 54.15682 87.27531 16.55925
## [14,] 70.97206 54.21547 87.72865 16.75659
## [15,] 71.21930 54.29151 88.14710 16.92780
## [16,] 71.45577 54.38057 88.53098 17.07521
## [17,] 71.68018 54.47895 88.88141 17.20123
## [18,] 71.89173 54.58344 89.20001 17.30829
## [19,] 72.08997 54.69128 89.48865 17.39868
## [20,] 72.27479 54.80021 89.74938 17.47459
## [21,] 72.44632 54.90833 89.98430 17.53798
## [22,] 72.60484 55.01416 90.19551 17.59068
## [23,] 72.75079 55.11652 90.38505 17.63426
## [24,] 72.88470 55.21454 90.55486 17.67016
## 
## $part
##           fcst    lower    upper       CI
##  [1,] 56.20206 34.42555 77.97858 21.77651
##  [2,] 56.97066 34.09865 79.84267 22.87201
##  [3,] 54.68760 29.56147 79.81373 25.12613
##  [4,] 54.06954 27.96958 80.16950 26.09996
##  [5,] 53.02247 26.03038 80.01457 26.99209
##  [6,] 52.27031 24.64432 79.89629 27.62598
##  [7,] 51.57798 23.45803 79.69794 28.11995
##  [8,] 50.95265 22.44679 79.45851 28.50586
##  [9,] 50.41566 21.61341 79.21791 28.80225
## [10,] 49.92440 20.88946 78.95934 29.03494
## [11,] 49.49682 20.28223 78.71141 29.21459
## [12,] 49.11220 19.75756 78.46685 29.35465
## [13,] 48.77353 19.31057 78.23649 29.46296
## [14,] 48.47221 18.92535 78.01907 29.54686
## [15,] 48.20602 18.59443 77.81761 29.61159
## [16,] 47.97048 18.30903 77.63193 29.66145
## [17,] 47.76256 18.06281 77.46231 29.69975
## [18,] 47.57915 17.85006 77.30824 29.72909
## [19,] 47.41758 17.66608 77.16908 29.75150
## [20,] 47.27540 17.50683 77.04398 29.76857
## [21,] 47.15042 17.36888 76.93196 29.78154
## [22,] 47.04067 17.24932 76.83202 29.79135
## [23,] 46.94438 17.14563 76.74314 29.79875
## [24,] 46.85998 17.05566 76.66431 29.80432

par(mar=c(3,3,1,1))
fanchart(fit.pr, colors = "red") # gráfico da previsão + erro

Figura 15: As previsões, em vermelho, de um modelo VAR(2) ajustado aos dados de mortalidade - poluição em LA.

Obviamente, a estimação via Yule-Walker, mínimos quadrados incondicionais e máxima verossimilhança (MLE) seguem diretamente das contrapartes univariadas. Para modelos VAR(p) puros, a estrutura de autocovariância leva à versão multivariada das equações de Yule-Walker: \[\begin{array}{c} \Gamma(h) \, = \, \displaystyle \sum_{j=1}^P \Phi_j \Gamma(h-j), \qquad h=1,2,\cdots, \\ \Gamma(0) \, = \, \displaystyle \sum_{j=1}^P \Phi_j \Gamma(-j)+\Sigma_{_W}, \end{array}\]

onde \(\Gamma(h)=\mbox{Cov}(X_{t+h},X_t)\) é uma matriz \(k\times k\) e \(\Gamma(-h)=\Gamma(h)^\top\).

A estimativa da matriz de autocovariâncias é semelhante ao caso univariado, ou seja, com \(\overline{x}=\frac{1}{n}\sum_{t=1}^n x_t\) como estimativa de \(\mu=\mbox{E}(X_t)\), \[\begin{equation*} \widehat{\Gamma}(h) \, = \, \dfrac{1}{n}\sum_{t=1}^{n-h} (x_{t+h}-\overline{x})(x_{t+h}-\overline{x})^\top, \qquad h=0,1,2,\cdots,n-1, \end{equation*}\] e \(\widehat{\Gamma}(-h)=\widehat{\Gamma}(h)^\top\).

Se \(\widehat{\gamma}_{i,j}(h)\) denota o elemento na \(i\)-ésima linha e \(j\)-ésima coluna de \(\widehat{\Gamma}(h)\), as funções de correlação cruzada (CCF), são estimadas por \[\begin{equation*} \widehat{\rho}_{i,j}(h) \, = \, \dfrac{\widehat{\gamma}_{i,j}(h)}{\sqrt{\widehat{\gamma}_{i,i}(0)\widehat{\gamma}_{j,j}(0)}}, \qquad h=0,1,2,\cdots,n-1\cdot \end{equation*}\] Quando \(i = j\) na expressão acima, obtemos a função de autocorrelação estimada (ACF) das séries individuais.

Embora a estimação por mínimos quadrados tenha sido usada em exemplos anteriores, poderíamos também ter usado a estimação de Yule-Walker, estimação por máxima verossimilhança condicional ou incondicional. Como no caso univariado, os estimadores de Yule-Walker, os estimadores de máxima verossimilhança e os estimadores de mínimos quadrados são assintoticamente equivalentes.

Para exibir a distribuição assintótica dos estimadores dos parámetros de autorregressão, escrevemos \[\begin{equation*} \phi \, = \, \mbox{vec}(\Phi_1,\cdots,\Phi_P), \end{equation*}\] onde o operador \(\mbox{vec}\) empilha as colunas de uma matriz em um vetor. Por exemplo, para um modelo bivariado AR(2), \[\begin{equation*} \phi \, = \, \mbox{vec}(\Phi_1,\phi_2) \, = \, (\Phi_{1_{11}},\Phi_{1_{21}},\Phi_{1_{12}},\Phi_{1_{22}},\Phi_{2_{11}},\Phi_{2_{21}},\Phi_{2_{12}},\Phi_{2_{22}})^\top, \end{equation*}\] onde \(\Phi_{l_{ij}}\) é o \(ij\)-ésimo elemento de \(\Phi_l\), \(l=1,2,\cdots\).

---

Teorema 1. Distribuição em grandes amostras dos estimadores VAR

Denotemos por \(\widehat{\phi}\) o vetor de estimadores dos parámetros obtidos via Yule – Walker, mínimos quadrados ou máxima verossimilhança para um modelo AR(p) \(k\)-dimensional. Então, \[\begin{equation*} \sqrt{n}\big(\widehat{\phi}-\phi\big) \, \sim \, \displaystyle AN\big(0,\Sigma_{_W} \otimes \Gamma_{PP}^{-1}\big), \end{equation*}\] onde \(\Gamma_{PP}=\{\Gamma(i-j)\}_{i,j=1}^P\) é uma matriz \(kp\times kp\) e \(\Sigma_{_W} \otimes \Gamma_{PP}^{-1}=\{\sigma_{ij}\Gamma_{PP}^{-1}\}_{i,l=1}^k\) é uma matriz \(k^2p\times k^2p\) com \(\sigma_{i,j}\) denotando o \(ij\)-ésimo elemento de \(\Sigma_{_W}\).

Demonstração Ver Lütkepohl (1993).

---

A matriz de variâncias e covariâncias do estimador \(\widehat{\phi}\) é aproximada substituindo \(\Sigma_{_W}\) por \(\widehat{\Sigma}_{_W}\) e substituindo \(\Gamma(h)\) por \(\widehat{\Gamma}(h)\) em \(\Gamma_{PP}\). A raiz quadrada dos elementos diagonais de \(\Sigma_{_W} \otimes \Gamma_{PP}^{-1}\) dividida por \(\sqrt{n}\) fornece os erros padrão individuais.

Para o exemplo de dados de mortalidade, os erros padrão estimados para o ajuste VAR(2) estão listados no Exemplo 11; embora esses erros padrão tenham sido obtidos de uma execução de regressão, eles também poderiam ter sido calculados usando o Teorema 1.

A série temporal vetorial \(X_t\) de dimensão \(k\times 1\), para \(t = 0,\pm 1,\pm 2,\cdots\) é dita ser VARMA(p,q) se \(X_t\) for estacionária e \[\begin{equation*} X_t \, = \, \alpha+\Phi_1X_{t-1}+\cdots+\Phi_p X_{t-p}+W_t+\Theta_1 W_{t-1}+\cdots+\Theta_q W_{t-q}, \end{equation*}\] com \(\Phi_p\neq 0\), \(\Theta_q\neq 0\) e \(\Sigma_{_W}> 0\), ou seja, \(\Sigma_{_W}\) é definida positivo.

As matrizes de coeficientes \(\Phi_j\), \(j = 1,\cdots,p\) e \(\Theta_j\), \(j = 1,\cdots,q\) são matrizes \(k\times k\). Se \(X_t\) tiver média \(\mu\), então \(\alpha=(\mbox{I}-\Phi_1-\cdots-\Phi_p)\mu\). Como no caso univariado, teremos que colocar uma série de condições no modelo ARMA multivariado para garantir que o modelo seja único e tenha propriedades desejáveis, como causalidade. Essas condições serão discutidas em breve.

Como no modelo VAR, a forma especial assumida para o componente constante pode ser generalizada para incluir um vetor de entradas \(r\times 1\) fixo, \(U_t\). Ou seja, poderiamos ter proposto o modelo vetorial ARMAX, \[\begin{equation*} X_t \, = \, \Gamma U_t \, + \, \sum_{j=1}^p \Phi_j X_{t-j} \, + \ \sum_{j=1}^q \Theta_j W_{t-j} \, + \, W_t, \end{equation*}\] onde \(\Gamma\) é uma matriz \(p\times r\) de parámetros.

Embora estender modelos AR univariados ou MA puros para o caso vetorial seja bastante fácil, estender modelos ARMA univariados para o caso multivariado não é uma questão simples. Nossa discussão será breve, mas os leitores interessados podem obter mais detalhes em Lütkepohl (1993), Reinsel (1997) e Tiao and Tsay (1989).

No caso multivariado, o operador autoregressivo assume a forma \[\begin{equation*} \Phi(B) \, = \, \mbox{I}-\Phi_1 B-\cdots-\Phi_p B^p, \end{equation*}\] e o operador de médias móveis é \[\begin{equation*} \Theta(B) \, = \, \mbox{I}+\Theta_1 B + \cdots +\Theta_q B^q\cdot \end{equation*}\]

O modelo VARMA(p,q) de média zero é então escrito de forma concisa como \[\begin{equation*} \Phi(B)X_t \, = \, \Theta(B)W_t\cdot \end{equation*}\]

O modelo é considerado causal se as raízes de \(|\Phi(z)|\), onde \(|\cdot|\) denota o determinante, estão fora do círculo unitário \(|z|> 1\); isto é, \(|\Phi(z)|\neq 0\) para qualquer valor \(z\) tal que \(|z|\leq 1\). Neste caso, podemos escrever \[\begin{equation*} X_t \, = \, \Psi(B)W_t, \end{equation*}\] onde \(\displaystyle \Psi=\sum_{j=0}^\infty \Psi_j B^j\), \(\Psi_0=\mbox{I}\) e \(\displaystyle \sum_{j=0}^\infty ||\Psi_j ||<\infty\).

O modelo será considerado invertível se as raízes de \(|\Theta(z)|\) estiverem fora do círculo unitário. Então, podemos escrever \[\begin{equation*} W_t \, = \, \Pi (B)X_t, \end{equation*}\] onde \(\displaystyle \Pi(B)=\sum_{j=0}^\infty \Pi_j B^j\), \(\Pi_0=\mbox{I}\) e \(\displaystyle \sum_{j=0}^\infty ||\Pi_j||<\infty\).

Analogamente ao caso univariado, podemos determinar as matrizes \(\Psi_j\) resolvendo \(\Psi(z) = \Phi(z)^{-1}\Theta(z)\), \(|z|\leq 1\) e as matrizes \(\Pi_j\) resolvendo \(\Pi(z) = \Theta(z)^{-1}\Phi(z)\), \(|z|\leq 1\).

Para um modelo causal, podemos escrever \(X_t = \Psi(B)W_t\) de forma que a estrutura geral de autocovariância de um modelo ARMA(p,q) seja, para \(h\geq 0\), \[\begin{equation*} \Gamma(h) \, = \, \mbox{Cov}(X_{t+h},X_t) \, = \, \sum_{j=0}^\infty \Psi_{j+h}\Sigma_{_W}\Psi_j^\top, \end{equation*}\] e \(\Gamma(-h)=\Gamma(h)^\top\). Para processos MA(q) puros, a expressão acima torna-se \[\begin{equation*} \Gamma(h) \, = \, \sum_{j=0}^{q-h} \Theta_{j+h}\Sigma_{_W}\Theta_j^\top, \end{equation*}\] onde \(\Theta_0=\mbox{I}\). Claro que, \(\Gamma(h)=0\) para \(h>q\).

Como no caso univariado, precisaremos de condições para a unicidade do modelo. Essas condições são semelhantes às condições no caso univariado em que os polinômios autoregressivos e de médias móveis não têm fatores comuns. Para explorar os problemas de unicidade que encontramos com modelos ARMA multivariados, considere um processo AR(1) bivariado, \(X_t = (X_{t,1},X_{t,2})^\top\), dado por \[\begin{array}{rcl} X_{t,1} & = & \phi X_{t-1,2}+W_{t,1},\\ X_{t,2} & = & W_{t,2}, \end{array}\]

onde \(W_{t,1}\) e \(W_{t,2}\) são processos de ruído branco independentes e \(|\phi|<1\). Ambos os processos, \(X_{t,1}\) e \(X_{t,2}\) são causais e invertíveis. Além disso, os processos são estacionários conjuntamente porque \[\begin{equation*} \mbox{Cov}(X_{t+h,1},X_{t,2}) \, = \, \phi\mbox{Cov}(X_{t+h-1,2},X_{t,2}) \, = \, \phi \gamma_{2,2}(h-1) \, = \, \phi\sigma_{_{W_2}}^2\delta_1^h, \end{equation*}\] que não depende de \(t\). Observe, \(\delta_1^h = 1\) quando \(h = 1\), caso contrário, \(\delta_1^h=0\).

Em notação matricial, podemos escrever este modelo como \[\begin{equation*} X_t \, = \, \Phi X_{t-1}+W_t, \qquad \mbox{onde} \qquad \Phi=\begin{pmatrix} 0 & \phi \\ 0 & 0 \end{pmatrix}\cdot \end{equation*}\] Podemos escrever a expressão anterior em notação de operador como \[\begin{equation*} \Phi(B)X_t \, = \, W_t, \qquad \mbox{onde} \qquad \Phi(z)=\begin{pmatrix} 1 & -\phi z \\ 0 & 1 \end{pmatrix}\cdot \end{equation*}\] Além disso, este modelo pode ser escrito como um modelo ARMA(1,1) bivariado \[\begin{equation*} X_t \, = \, \Phi X_{t-1}+\Theta W_{t-1}+W_t, \end{equation*}\] onde \[\begin{equation*} \Phi_1=\begin{pmatrix} 0 & \phi+\theta \\ 0 & 0 \end{pmatrix} \qquad \mbox{e} \qquad \Theta_1=\begin{pmatrix} 0 & -\theta \\ 0 & 0 \end{pmatrix}, \end{equation*}\] e \(\theta\) arbitrário.

Para verificar isto, podemos escrever o modelo ARMA(1,1) como \[\begin{equation*} \Phi_1(B)X_t \, = \, \Theta_1(B)W_t \end{equation*}\] ou \[\begin{equation*} \Theta_1(B)^{-1}\Phi_1(B)X_t \, = \, W_t, \end{equation*}\] onde \[\begin{equation*} \Phi_1(z)=\begin{pmatrix} 1 & -(\phi+\theta)z \\ 0 & 1 \end{pmatrix} \qquad \mbox{e} \qquad \Theta_1(z)=\begin{pmatrix} 1 & -\theta z \\ 0 & 1 \end{pmatrix}\cdot \end{equation*}\] Então \[\begin{equation*} \Theta_1(z)^{-1}\Phi_1(z) \, = \, \begin{pmatrix} 1 & \theta z \\ 0 & 1 \end{pmatrix}\begin{pmatrix} 1 & -(\phi+\theta)z \\ 0 & 1 \end{pmatrix} \, = \, \begin{pmatrix} 1 & -\phi z \\ 0 & 1 \end{pmatrix} \, = \, \Phi(z), \end{equation*}\] onde \(\Phi(z)\) é o polinômio associado ao modelo bivariado AR(1). Por ser \(\theta\) arbitrário, os parámetros do modelo ARMA(1,1) não são identificáveis. Não existe nenhum problema, entretanto, no ajuste do modelo AR(1) dado anteriormente.

O problema na discussão anterior foi causado pelo fato de que tanto \(\Theta(B)\) quanto \(\Theta(B)^{-1}\) são finitos; esse operador de matriz é chamado de unimodular. Se \(U(B)\) for unimodular, \(|U(z)|\) será constante. Também é possível que dois modelos ARMA(p,q) multivariados aparentemente diferentes, digamos \[\begin{equation*} \Phi (B)X_t \, = \, \Theta (B)W_t \qquad \mbox{e} \qquad \Phi_*(B)X_t \, = \, \Theta_*(B)W_t, \end{equation*}\] sejam relacionados por meio de um operador unimodular, \(U(B)\) como \[\begin{equation*} \Phi_*(B) \, = \, U(B)\Phi(B) \qquad \mbox{e} \qquad \Theta_*(B) \, = \, U(B)\Theta(B), \end{equation*}\] de forma que as ordens de \(\Phi(B)\) e \(\Theta(B)\) sejam iguais à ordens de \(\Phi_*(B)\) e \(\Theta_*(B)\), respectivamente. Por exemplo, considere os modelos bivariados ARMA(1,1) dados por \[\begin{equation*} \Phi X_t = \begin{pmatrix} 1 & -\phi B \\ 0 & 1 \end{pmatrix} X_t \, = \, \begin{pmatrix} 1 & \theta B \\ 0 & 1 \end{pmatrix}W_t \, = \, \Theta W_t \end{equation*}\] e \[\begin{equation*} \Phi_*(B) X_t = \begin{pmatrix} 1 & (\alpha-\phi) B \\ 0 & 1 \end{pmatrix} X_t \, = \, \begin{pmatrix} 1 & (\alpha+\theta) B \\ 0 & 1 \end{pmatrix}W_t \, = \, \Theta_*(B) W_t, \end{equation*}\] onde \(\alpha\), \(\phi\) e \(\theta\) são constantes arbitrárias.

Este resultado implica que os dois modelos possuem a mesma função de autocovariância \(\Gamma(h)\). Dois desses modelos ARMA(p,q) são considerados equivalentes observacionalmente.

Como mencionado anteriormente, além de exigir causalidade e invertibilidade, precisaremos de algumas suposições adicionais no caso multivariado para garantir que o modelo seja único. Para garantir a identificabilidade dos parámetros do modelo multivariado ARMA(p,q), precisamos das seguintes duas condições adicionais:

1. os operadores matriciais \(\Phi(B)\) e \(\Theta(B)\) não têm fatores à esquerda comuns além dos unimodulares, isto é, se \(\Phi(B) = U(B)\Phi_*(B)\) e \(\Theta(B) = U(B)\Theta_*(B)\), o fator comum deve ser unimodular e

2. com \(q\) tão pequeno quanto possível e \(p\) tão pequeno quanto possível para aquele \(q\), a matriz \((\Phi_p,\Theta_q)\) deve ser de posto completo \(k\).

Uma sugestão para evitar a maioria dos problemas mencionados é ajustar apenas modelos vetoriais AR(p) em situações multivariadas. Embora essa sugestão possa ser razoável para muitas situações, essa filosofia não está de acordo com a lei da parcimônia porque podemos ter que ajustar um grande número de parámetros para descrever a dinâmica de um processo.

As inferências assintóticas para o caso geral de modelos ARMA vetoriais são mais complicadas do que os modelos AR puros; detalhes podem ser encontrados em Reinsel (1997) ou Lütkepohl (1993), por exemplo. Também observamos que a estimação para modelos VARMA podem ser reformuladas no problema de estimação para modelos de espaço de estados.

Exemplo 12: O algoritmo Spliid para ajustar ARMA vetorial

Vale a pena mencionar um algoritmo simples para ajustar modelos vetoriais ARMA devido a Spliid (1983); porque usa repetidamente as equações de regressão multivariadas.

Considere um modelo ARMA(p,q) geral para uma série temporal com uma média diferente de zero \[\begin{equation*} X_t \, = \, \alpha+\Phi_1 X_{t-1}+\cdots+\Phi_p X_{t-p}+W_t+\Theta_1 X_{t-1}+\cdots+\Theta_q W_{t-q}\cdot \end{equation*}\] Se \(\mu=\mbox{E}(X_t)\), então \(\alpha=(\mbox{I}-\Phi_1-\cdots-\Phi_p)\mu\).

Se \(W_{t-1},\cdots,W_{t-q}\) foram observados, poderiamos reorganizar a equação acima como um modelo de regressão multivariado \[\begin{equation*} X_t \, = \, \mathcal{B}Z_t + W_t, \end{equation*}\] com \[\begin{equation*} Z_t \, = \, \big(1,X_{t-1}^\top,\cdots,X_{t-p}^\top,W_{t-1}^\top,\cdots,W_{t-q}^\top \big)^\top \end{equation*}\] e \[\begin{equation*} \mathcal{B} \, = \, \big( \alpha,\Phi_1,\cdots,\Phi_p,\Theta_1,\cdots,\Theta_q\big), \end{equation*}\] para \(t=p+1,\cdots,n\).

Dado um estimador inicial \(\mathcal{B}_0\) de \(\mathcal{B}\), podemos reconstruir \(\{W_{t-1},\cdots,W_{t-q}\}\) fazendo \[\begin{equation*} W_{t-j} \, = \, X_{t-j}-\mathcal{B}_0Z_{t-j}, \qquad t=p+1,\cdots,n, \quad j=1,\cdots,q, \end{equation*}\] onde, se \(q>p\), colocamos \(W_{t-j}=0\) para \(t-j\leq 0\).

Os novos valores de \(\{W_{t-1},\cdots,W_{t-q}\}\) são então colocados nos regressores \(Z_t\) e uma nova estimativa, digamos \(\mathcal{B}_1\), é obtida. O valor inicial, \(\mathcal{B}_0\), pode ser calculado ajustando uma autorregressão pura de ordem \(p\) ou superior e tomando \(\Theta_1 = \cdots = \Theta_q = 0\). O procedimento é então iterado até que as estimativas dos parámetros se estabilizam.

O algoritmo frequentemente converge, mas não para os estimadores de máxima verossimilhança. A experiência sugere que os estimadores podem estar razoavelmente próximos dos estimadores de máxima verossimilhança. O algoritmo pode ser considerado uma maneira rápida e fácil de ajustar um modelo VARMA inicial como ponto de partida para usar a estimação de máxima verossimilhança, que é melhor feita por meio de modelos de espaço de estado.

Usamos o pacote R marima para ajustar o modelo ARMA(2;1) vetorial ao conjunto de dados de mortalidade - poluição. Notamos que a mortalidade é destendida antes da análise. As previsões de mortalidade um passo à frente são mostradas na figura a continuação.

library(astsa)
library(marima)
model = define.model(kvar=3, ar=c(1,2), ma=c(1))
arp = model$ar.pattern; map = model$ma.pattern
cmort.d = resid(detr <- lm(cmort~ time(cmort), na.action=NULL))
xdata = matrix(cbind(cmort.d, tempr, part), ncol=3) # removendo atributo ts
fit = marima(xdata, ar.pattern=arp, ma.pattern=map, means=c(0,1,1), penalty=1)

## All cases in data,  1  to  508  accepted for completeness.
## 508 3  = MARIMA - dimension of data

innov = t(resid(fit)); plot.ts(innov); acf(innov, na.action=na.pass)

# valores ajustados para cmort
pred = ts(t(fitted(fit))[,1], start=start(cmort), freq=frequency(cmort)) + 
  detr$coef[1] + detr$coef[2]*time(cmort)
plot(pred, ylab="Cardiovascular Mortality", lwd=2, col=4, cex=0.7); points(cmort)
# imprimindo as estimativas e correspondentes estatísticas t^2
short.form(fit$ar.estimates, leading=FALSE)

## , , Lag=1
## 
##        x1=y1     x2=y2     x3=y3
## y1 -0.311106  0.000000 -0.113831
## y2  0.000000 -0.656331  0.048208
## y3 -0.108892  0.000000 -0.860911
## 
## , , Lag=2
## 
##        x1=y1     x2=y2     x3=y3
## y1 -0.333368  0.132523 -0.047171
## y2  0.000000 -0.200135  0.054616
## y3  0.179117 -0.102144 -0.151250

short.form(fit$ar.fvalues, leading=FALSE)

## , , Lag=1
## 
##        x1=y1    x2=y2      x3=y3
## y1 51.211590  0.00000   7.909997
## y2  0.000000 41.73927   3.141522
## y3  1.571085  0.00000 113.299823
## 
## , , Lag=2
## 
##        x1=y1     x2=y2    x3=y3
## y1 67.241600 11.890684 2.515931
## y2  0.000000  8.101381 2.903037
## y3  4.860958  1.768893 6.477311

short.form(fit$ma.estimates, leading=FALSE)

## , , Lag=1
## 
##        x1=y1     x2=y2     x3=y3
## y1  0.000000 -0.186710 -0.106470
## y2 -0.114446 -0.446431  0.000000
## y3  0.000000 -0.278378 -0.672962

short.form(fit$ma.fvalues, leading=FALSE)

## , , Lag=1
## 
##       x1=y1     x2=y2    x3=y3
## y1 0.000000 14.514812  4.75441
## y2 4.683436 16.375957  0.00000
## y3 0.000000  8.079861 47.56415

grid()

fit$resid.cov # estimativa da matriz cov de ruído

##           u1        u2        u3
## u1 27.346796  6.503133  13.79950
## u2  6.503133 36.202990  38.12202
## u3 13.799495 38.122018 109.20702

Figura 16: As previsões de um VAR(2) ajustado aos dados de mortalidade - poluição em LA.

7. Exercícios

1- O conjunto de dados arf, no pacote astsa contêm 1000 observações simuladas de um modelo ARFIMA(1,1,0) com \(\phi=0.75\) e \(d=0.4\).

1. Mostre um gráfico dos dados e comentar.

2. Faça os gráficos ACF e PACF dos dados e comente.

3. Estime os parámetros e teste a significância dos estimadore \(\widehat{\phi}\) e \(\widehat{d}\).

4. Explique por que, usando os resultados das partes (a) e (b), parece razoável diferenciar os dados antes da análise. Ou seja, se \(X_t\) representa os dados, explique por que podemos optar por ajustar um modelo ARMA ao \(\nabla X_t\).

5. Mostre os gráficos do ACF e o PACF do \(\nabla X_t\) e comente.

6. Ajuste um modelo ARMA para \(\nabla X_t\) e comentar.

2- Calcule a função de autocorrelação (ACF) amostral dos valores absolutos dos retornos da NYSE, ou seja, dos retornos ou mudança percentual da Bolsa de Valores de Nova York (NYSE) entre 2 de Fevereiro de 1984 a 31 de Dezembro de 1991 exibidos na figura abaixo, até o atraso 200 e comente se o ACF indica memória longa. Ajuste um modelo ARFIMA aos valores absolutos e comente.

par(mar=c(4,4,4,2))
plot(nyse, xlab="Tempo", ylab="Retorno da Bolsa de Valores de Nova Iorque (NYSE)", 
     main="Retornos de mercado ponderados pelo valor diário\n 
     entre 2 de fevereiro de 1984 e 31 de dezembro de 1991")
grid()

3- Mostra graficamene a série de temperaturas globais globtemp (pacote astsa) e, em seguida, teste se existe uma raiz unitária versus a alternativa de que o processo é estacionário usando os três testes, DF, ADF e PP, discutidos. Comente.

4- Mostre graficamente a série do Produto Interno Bruto (PIB), de sigla GNP en inglês (gnp, pacote astsa) e teste a existência de raiz unitária contra a alternativa de que o processo é explosivo. Comente.

5- Verifique que \[\begin{equation*} \nabla X_t = \gamma X_{t-1} +\sum_{j=1}^{p-1} \psi_j \nabla X_{t-j} +W_t, \end{equation*}\] onde \(\gamma = \sum_{j=1}^p \phi_j-1\) e \(\psi_j=-\sum_{i=j}^p \phi_i\), para \(j=2,\cdots,p\).

6- Os preços semanais do petróleo bruto em dólares por barril estão no arquivo oil (pacote astsa). Investigue se a taxa de crescimento do preço semanal do petróleo apresenta comportamento GARCH. Nesse caso, ajuste um modelo apropriado à taxa de crescimento.

7- O pacote stats contêm os preços de fechamento diário de quatro dos principais índices de ações europeus; digite help (EuStockMarkets) para obter detalhes. Ajuste um modelo GARCH aos retornos de uma dessas séries e discuta suas descobertas.
Observação: o conjunto de dados contêm valores reais e não retornos. Portanto, os dados devem ser transformados antes do ajuste do modelo.

8- O vetor gradiente \(\ell^{(1)}(\alpha_0,\alpha_1)\), para um modelo ARCH(1) foi mostrado ser \[\begin{equation} \ell^{(1)}(\alpha_0,\alpha_1) = \begin{pmatrix} \displaystyle \partial \ell \big/ \partial\alpha_0 \\ \displaystyle \partial \ell \big/ \partial\alpha_1 \end{pmatrix} \, = \, \sum_{t=2}^n \begin{pmatrix} 1 \\ r_{t-1}^2 \end{pmatrix} \times \frac{\displaystyle \alpha_0+\alpha_1 r_{t-1}^2-r_t^2}{\displaystyle \big( \alpha_0+\alpha_1 r_{t-1}^2\big)^2}\cdot \end{equation}\]
Prove este resultado e depois use-o para calcular a matriz Hessiana \[\begin{equation} \ell^{(2)}(\alpha_0,\alpha_1) \, = \, \begin{pmatrix} \displaystyle \partial^2 \ell \big/ \partial\alpha_0^2 & \displaystyle \partial^2 \ell \big/ \partial\alpha_0\partial\alpha_1 \\ \displaystyle \partial^2 \ell \big/ \partial\alpha_0\partial\alpha_1 & \displaystyle \partial^2 \ell \big/ \partial\alpha_1^2 \end{pmatrix}\cdot \end{equation}\]

9- Os dados de manchas solares sunspotz (pacote astsa) foram mostrados anteriormente e são apresentados na figura abaixo. A partir do gráfico discuta por que é razoável ajustar um modelo de limiar aos dados e, em seguida, ajuste um modelo de limiar.

library(astsa)
par(mar=c(3,3,1,1), mgp=c(1.6,.6,0))
plot(soi) # gráfico dos dados
grid()

10- Os dados no arquivo climhyd contém 454 meses de valores medidos para as variáveis climáticas: temperatura do ar (Temp), ponto de orvalho (DewPt), cobertura de nuvens (CldCvr), velocidade do vento (WndSpd), precipitação (Precip) e influxo (Inflow), no Lago Shasta; os dados são exibidos na Figura VII.3. Gostaríamos de examinar as possíveis relações entre os fatores climáticos e o influxo para o Lago Shasta.

1. Ajuste um modelo \(ARIMA(0,0,0)\times (0,1,1)_{12}\) para (i) a precipitação transformada \(P_t = \sqrt{\mbox{Precip}}\) e (ii) influxo transformado \(I_t = \ln(\mbox{Inflow})\).

2. Ajuste um modelo ARIMA na parte (a) para a precipitação transformada para a série de influxo para gerar os resíduos de influxo pré-branqueados assumindo o modelo para a precipitação. Calcule a correlação cruzada entre os resíduos de influxo usando o modelo ARIMA de precipitação e os resíduos de precipitação usando o modelo de precipitação e interprete. Use os coeficientes do modelo ARIMA para construir os resíduos de influxo transformados.

11- Para o conjunto de dados climhyd (pacote astsa), considere prever os influxos transformados \(I_t = \ln(\mbox{Inflow})\) a partir dos valores da precipitação transformados \(P_t = \sqrt{\mbox{Precip}}\) usando um modelo de função de transferência da forma \[\begin{equation*} (1-B^{12})I_t \, = \, \alpha(B)(1-B^{12})P_t+\eta_t, \end{equation*}\] onde assumimos que a diferenciação sazonal é uma coisa razoável a se fazer. Você pode pensar nisso como adequado ajustando \[\begin{equation*} Y_t \, = \, \alpha(B)X_t + eta_t, \end{equation*}\] onde \(Y_t\) e \(X_t\) são os influxos e precipitações transformados sazonalmente diferenciados.

1. Argumente que \(X_t\) pode ser ajustado por médias móveis sazonais de primeira ordem e use a transformação obtida para pré-branquear a série \(X_t\).

2. Aplique a transformação aplicada em (a) à série \(Y_t\) e calcule a função de correlação cruzada relacionando a série pré-branqueada à série transformada. Defenda uma função de transferência da forma \[\begin{equation*} \alpha(B) \, = \,\dfrac{\delta_0}{1-\omega_1 B}\cdot \end{equation*}\]

3. Escreva o modelo geral obtido na forma de regressão para estimar \(\delta_0\) e \(\omega_1\). Observe que você estará minimizando as somas ao quadrado dos resíduos para a série de ruído transformado \((1-\widehat{\omega}_1 B)\eta_t\). Retenha os resíduos para modelagem adicional envolvendo o ruído \(\eta_t\). O resíduo observado é \(U_t=(1-\widehat{\omega}_1 B)\eta_t\).

4. Ajuste os resíduos de ruído obtidos em (c) com um modelo ARMA e forneça a forma final sugerida por sua análise nas partes anteriores.

5. Discuta o problema de prever \(Y_{t + m}\) usando o passado infinito de \(Y_t\) e o presente e o passado infinito de \(X_t\). Determine o valor previsto e a variação da previsão.

12- Considere o conjunto de dados econ5 (pacote astsa) contendo o desemprego trimestral nos EUA (unemp), PIB (gnp), consumo (consum) e investimento privado (prinv) e governamental (govinv) de 1948-III a 1988-II. O componente sazonal foi removido dos dados. Concentrando-se no desemprego (\(U_t\)), no PIB (\(G_t\)) e no consumo (\(C_t\)), ajuste um modelo ARMA vetorial aos dados depois de logaritmizar cada série e, em seguida, remover a tendência linear. Ou seja, ajustar um modelo ARMA vetorial para \(X_t = (X_{1t}, X_{2t}, X_{3t})^\top\) onde, por exemplo, \[\begin{equation} X_{1t} = \log(U_t)-\widehat{\beta}_0-\widehat{\beta}_1t, \end{equation}\] sendo \(\widehat{\beta}_0\) e \(\widehat{\beta}_1\) as estimativas de mínimos quadrados para a regressão de \(\log(U_t)\) no tempo \(t\). Execute um conjunto completo de diagnósticos nos resíduos.

8. Bibliografia

Akaike, H. 1969. “Fitting Autoregressive Models for Prediction.” Annals of the Institute of Statistical Mathematical, no. 21: 243–47.

Bedrick, E. J., and C.-L. Tsai. 1994. “Model Selection for Multivariate Regression in Small Samples.” Biometrics, no. 50: 226–31.

Beran, J. 1994. Statistics for Long Memory Processes. New York: Chapman; Hall.

Billingsley, P. 1999. Convergence of Probability Measures. 2nd edition. New York: Wiley.

Bollerslev, T. 1986. “Generalized Autoregressive Conditional Heteroscedasticity.” Journal of Econometrics, no. 31: 307–27.

Box, G. E. P., and G. M. Jenkins. 1970. Time Series Analysis, Forecasting, and Control. Oakland, CA: Holden-Day.

Brockwell, P. J., and R. A. Davis. 1991a. Time Series: Theory and Methods. 2nd ed. New York: Springer-Verlag.

———. 1991b. Time Series: Theory and Methods. 2nd ed. New York: Springer-Verlag.

Chan, N. H. 2002. Time Series: Applications to Finance. New York: Wiley.

Dahlhaus, R. 1989. “Efficient Parameter Estimation for Self-Similar Processes.” Annals of Statistics, no. 17: 1749–66.

Dickey, D., and W. Fuller. 1979. “Distribution of the Estimators for Autoregressive Time Series with a Unit Root.” Journal of the American Statistical Association, no. 74: 427–31.

Douc, R., E. Moulines, and D. S. Stoffer. 2014. Nonlinear Time Series: Theory, Methods, and Applications with r Examples. Boca Raton: CRC Press.

Engle, R. F. 1982. “Autoregressive Conditional Heteroscedasticity with Estimates of the Variance of United Kingdom Inflation.” Econometrica, no. 50: 987–1007.

Engle, R. F., D. Nelson, and T. Bollerslev. 1994. ARCHModels. In Handbook of Econometrics. Edited by R. Engle and D. McFadden. Vol. IV. Amsterdam: North Holland.

Evans, G. B. A., and N. E. Savin. 1981. “The Calculation of the Limiting Distribution of the Least Squares Estimator of the Parameter in a Random Walk Model.” Annals of Statistics, no. 9: 1114–18.

Fox, R., and M. S. Taqqu. 1986. “Large Sample Properties of Parameter Estimates for Strongly Dependent Stationary Gaussian Time Series.” Annals of Statistics, no. 14: 517–32.

Geweke, J. F., and S. Porter-Hudak. 1983. “The Estimation and Application of Long-Memory Time Series Models.” Journal of the Time Series Analysis, no. 4: 221–38.

Gouriéroux, C. 1997. ARCH Models and Financial Applications. New York: Springer-Verlag.

Granger, C. W. J., and R. Joyeux. 1980. “An Introduction to Long-Memory Time Series Models and Fractional Differencing.” Jounal of the Time Series Analysis, no. 1: 15–29.

Hannan, E. J., and B. G. Quinn. 1979. “The Determination of the Order of an Autoregression.” Journal of the Royal Statistical Society, B, no. 41: 190–95.

Haslett, J., and A. E. Raftery. 1989. “Space-Time Modelling with Long-Memory Dependence: Assessing Ireland’s Wind Power Resource.” Applied Statistics, no. 38: 1–21.

Henrici, P. 1974. Applied and Computational Complex Analysis. 1st ed. New York: Wiley.

Hosking, J. R. M. 1981. “Fractional Differencing.” Biometrika, no. 68: 165–76.

———. 1984. “Modelling Persistence in Hydrological Time Series Using Fractional Differencing.” Water Resources Research, no. 20: 1898–908.

Hurst, H. 1951. “Long Term Storage Capacity of Reservoirs.” Transactions of the American Society of Civil Engineers, no. 116: 778–808.

Hurvich, C. M., and K. I. Beltrao. 1993. “Asymptotics for the Low-Requency Oridnates of the Periodogram for a Long-Memory Time Series.” Journal of the Time Series Analysis, no. 14: 455–72.

Hurvich, C. M., and R. S. Deo. 1999. “Plug-in Selection of the Number of Frequencies in Regression Estimates of the Memory Parameter of a Long-Memory Time Series.” Journal of the Time Series Analysis, no. 20: 331–41.

Hurvich, C. M., R. S. Deo, and J. Brodsky. 1998. “The Mean Squared Error of Geweke and Porter-Hudak’s Estimator of the Memory Parameter of a Long-Memory Time Series.” Journal of the Time Series Analysis, no. 19: 19–46.

Lütkepohl, H. 1993. Introduction to Multiple Time Series Analysis. 2nd edition. Berlin: Springer-Verlag.

McLeod, A. I., and K. W. Hipel. 1978. “Preservation of the Rescales Adjusted Range.” Water Resources Research, no. 14: 491–518.

Palma, W. 2007. Long-Memory Time Series: Theory and Methods. New York: Wiley.

Phillips, P. C. B. 1987. “Time Series Regression with a Unit Root.” Econometrica, no. 55: 227–301.

Priestley, M. B. 1988. Nonlinear and Nonstationary Time Series Analysis. London: Academic Press.

Rao, M. M. n.d. “Asymptotic Distribution of an Estimator of the Boundary Parameter of an Unstable Process.” Annals of Statistics, no. 6: 185–90.

Reinsel, G. C. 1997. Elements of Multivariate Time Series Analysis. 2nd edition. New York: Springer-Verlag.

Robinson, P. M. 1995. “Gaussian Semiparametric Estimation of Long Range Dependence.” Annals of Statistics, no. 23: 1630–61.

Shephard, N. 1996. Statistical Aspects of ARCH and Stochastic Volatility. In Time Series Models in Econometrics, Finance and Other Fields. Edited by D. R. Cox D. V. Hinkley and O. E. Barndorff-Nielson. London: Chapman; Hall.

Spliid, H. 1983. “A Fast Estimation Method for the Vector Autoregressive Moving Average Model with Exogenous Variables.” Journal of the American Statistical Associations, no. 78: 843–49.

Tiao, G. C., and R. S. Tsay. 1989. “Model Specification in Multivariate Time Series.” Journal of the Royal Statistics Society, B, no. 51: 157–213.

Tong, H. 1983. Threshold Models in Nonlinear Time Series Analysis. Springer Lecture Notes in Statistics 21. New York: Springer-Verlag.

———. 1990. Nonlinear Time Series: A Dynamical System Approach. Oxford: Oxford University Press.

Tsay, R. S. 2002. Analysis of Financial Time Series. New York: Wiley.

Weiss, A. A. 1984. “ARMA Models with ARCH Errors.” Journal of the Time Series Analysis, no. 5: 129–43.