Laboratório 3 - Séries Temporais 1

Neste Laboratório iremos reproduzir o script 6, que foi visto em sala de aula, e descrever e interpretar os resultados obtidos.

Simulate AR(1)

Nesta primeira atividade, iremos simular um processo auto regressivo de ordem p = 1, com 30 observações e com paramêtro igual a -0.8264. Também estaremos descartando as 10 primeiras observações desta série e utilizando apenas as 20 restantes. Ao calcular a acf e a pacf, vemos que de fato, o processo se apresenta como um AR(1), pois temos uma queda no gráfico da ACF e no PACF vemos que para lags > p= 1, a função assume valores iguais a 0 (dentro do intervalo). O código para executar tais procedimentos, assim como o output, podem ser vistos a seguir.

ar.sim<-arima.sim( n=30, list(ar=c(-.8264)))
plot(seq(11,30),ar.sim[11:30])
lines(seq(11,30),ar.sim[11:30])

acf(ar.sim)

pacf(ar.sim)

Forecast 1-step

Nesta etapa iremos fazer previsões de 1 passo a frente do modelo simulado AR(1). A previsão feita a um passo da série temporal acima pode ser feita no R através do código abaixo:

ar.forecast<-vector()
ar.forecast.se<-vector()
for(n in 1:21){
AR_fit <-arima(ar.sim[1:(9+n)], order=c(1,0,0))
#print(AR_fit)
ar.forecast[n]<-predict(AR_fit, n.ahead = 1)$pred
ar.forecast.se[n]<-predict(AR_fit, n.ahead = 1)$se
}

plot(seq(10,30),ar.forecast)
lines(seq(10,30),ar.forecast)

True values vs forecasts / Confidence interval

Agora iremos comparar os valores obtidos através do forecast 1-step com os valores originais do processo AR(1) simulado. Ou seja, estamos comparando os valores originais da série com os valores previstos através da previsão feita para um passo a frente. Como podemos ver abaixo, a linha em vermelho representa as previsão para n.ahead=1 e as linhas em preto são os valores observados do processo AR(1) que simulamos anteriormente.
Também podemos construir intervalos de confiança para as previsões que realizamos aqui, ou seja, fazer previsões intervalar.

plot(seq(11,30),ar.sim[11:30])
lines(seq(11,30),ar.sim[11:30])
#plot(seq(10,30),ar.forecast,col = 2, lty = 2)
lines(seq(10,30),ar.forecast,col = 2, lty = 1)

# confidence interval:

library(graphics)
plot(seq(10,30),ar.sim[10:30])
lines(seq(10,30),ar.sim[10:30])
points(seq(10,30),ar.forecast,col = 2, lty = 1)
lines(seq(10,30),ar.forecast,col = 2, lty = 1)
points(seq(10,30),ar.forecast - 2*ar.forecast.se, type = "l", col = 2, lty = 2)
points(seq(10,30),ar.forecast + 2*ar.forecast.se, type = "l", col = 2, lty = 2)

AR(2)

A seguir podemos ver a ACF e a PACF para o desfasamento 24 de um modelo auto regressivo de ordem p=2. Os valores dos parâmetros deste modelo AR(2) são 1.5 e -0.75. Iremos discutir o que esses gráficos nos informam a respeito do modelo da série temporal. Pelos gráficos, conseguimos identificar que o processo se trata de uma série que pode ser bem modelada pelo AR(2), uma vez que na ACF vemos uma queda rápida e na PACF vemos que para os lags > p = 2, a função de autocorrelação parcial assume valor igual a 0. Os resultados aqui discutidos facilitam na identificação do modelo e de sua ordem adequado para modelar uma série temporal.

# O ACF e o PACF, para o desfasamento 24, de um modelo AR(2): x_t= phi_1 X_{t-1}+phi_2 X_{t-2}+w_t,
# com phi_1 = 1,5 e phi_2 = -.75.
acf = ARMAacf(ar=c(1.5,-.75), ma=0, 24)
pacf = ARMAacf(ar=c(1.5,-.75), ma=0, 24, pacf=T)
 par(mfrow=c(1,2))
 plot(acf, type="h", xlab="lag")
 abline(h=0)
 plot(pacf, type="h", xlab="lag")
 abline(h=0)

Exemplo 3.16 do livro S&S - Dados recruit.dat

Aqui iremos usar o que já foi discutido neste laboratório para um conjunto de dados do livro do Shumway e Stoffer. O conjunto de dados em questão consiste em 453 observações (1 para cada mês) de recrutamento ao longo dos anos. A princípio, iremos plotar o gráfico da série temporal para observar o seu comportamento.
Pelo gráfico abaixo, vemos que a série aparentemente é estacionária e apresenta algum componente sazonal.
Pelo teste ADF, verificamos a estacionaridade, uma vez que rejeitamos a hipótese nula de que a série é não estacionária, pois o p-valor é muito baixo e menor que 0.05. Calculando a ACF e a PACF, conseguimos confirmar a estacionaridade por meio da queda de valores na ACF que é rápida. Pela ACF também vemos indicios da sazonalidade, através da correlação existente para lags múltiplos e distantes, o que da ideia da existência de um ciclo.
Agora iremos identificar um modelo para esta série temporal. Observando a PACF, vemos que para lags > 2, esta função assume valor igual a 0, ou seja, até o lag p=2, temos valores significativos para a PACF. Logo, através do comportamento descrito da ACF e PACF, o modelo adequado é o AR(2). Por fim, ajustamos este modelo aos dados para obter as estimativas de seus 2 parâmetros.

library(tseries)

rec = scan("rec.dat")
rec = ts(rec,frequency=12,start=1950)

plot(rec)

adf.test(rec)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  rec
## Dickey-Fuller = -5.4298, Lag order = 7, p-value = 0.01
## alternative hypothesis: stationary

#par(mfrow=c(2,1))
 acf(rec, 48)

 pacf(rec, 48)

 # O ACF e o PACF  dado nas Figuras acima s?o consistentes com o comportamento de um AR (2)
 # o  ACF tem ciclos correspondendo aproximadamente a um per?odo de 12 meses, e o
 # PACF tem grandes valores para h = 1, 2 e ent?o ? essencialmente zero para maior ordem de defasagens.
 ##  Assim: x_t=phi_0 + phi_1 X_{t-1}+phi_2 X_{t-2}+w_t
 fit=ar.ols(rec,aic=F,order.max=2,demean=F,intercept=T)
 fit # estimates

## 
## Call:
## ar.ols(x = rec, aic = F, order.max = 2, demean = F, intercept = T)
## 
## Coefficients:
##       1        2  
##  1.3541  -0.4632  
## 
## Intercept: 6.737 (1.111) 
## 
## Order selected 2  sigma^2 estimated as  89.72

 fit$asy.se # standard errors

## $x.mean
## [1] 1.110599
## 
## $ar
## [1] 0.04178901 0.04187942