In [1]:
options(OutDec = ",")
#==============================================================================
# Descricao
#==============================================================================
# Razao de crescimento trimestral do produto nacional bruto (PNB) real
# americano, sazonalmente ajustado, do segundo trimestre de 1947 ate'
# o primeiro trimestre de 1991.
#
#==============================================================================
# Carregando bibliotecas
#==============================================================================
silence <- suppressPackageStartupMessages # Para omitir mensagens de alertas
silence(library(tseries)) # Teste de Jarque-Bera: 'jarque.bera.test'
silence(library(forecast)) # Para previsao
silence(library(psych)) # Estatistica descritiva: 'describe'
In [2]:
#==============================================================================
# Lendo os dados do arquivo gnp.txt
#==============================================================================
gnp <- ts(scan("../dados/gnp.txt")*100,start=c(1947,2),frequency=4)
n <- length(gnp)
In [3]:
#==============================================================================
# Grafico da serie GNP (PNB)
#==============================================================================
par(mfrow=c(1,1),lwd=2.0,cex.lab=1.3,cex.axis=1.3,lab=c(10,10,5),
mar=c(5,5,2,2.5),cex.main=2.0,bty="n")
ts.plot(gnp,xlab="Trimestre",ylab="PNB",lwd=2,ylim=c(-3,5),main="",
xlim=c(1945,1995))
points(gnp,pch=15,cex=1.2)
In [4]:
#==============================================================================
# ACF
#==============================================================================
par(mfrow=c(1,1),lwd=2.0,cex.lab=1.5,cex.axis=1.5,lab=c(10,10,5),
mar=c(5,5,2,2.5),cex.main=2.0,bty="n")
acf(c(gnp),lwd=2,col="black",lag.max=16,xlab="Defasagem",main="",
ylim=c(-0.2,1),xlim=c(0,16))
In [5]:
#==============================================================================
# PACF
#==============================================================================
par(mfrow=c(1,1),lwd=2.0,cex.lab=1.5,cex.axis=1.5,lab=c(8,5,5),
mar=c(5,5,2,2.5),cex.main=2.0,bty="n")
pacf(c(gnp),lwd=2,col="black",lag.max=16,xlab="Defasagem",main="",
ylim=c(-0.2,1),xlim=c(0,16))
In [6]:
#==============================================================================
# Ajustando modelos AR(p)
#==============================================================================
# Por enquanto, utilizaremos a funcao "ar"
#==============================================================================
gnp.ar <- ar(gnp,order.max=15)
p <- gnp.ar$order
print(gnp.ar)
print("====================")
print(gnp.ar$aic) # Ver pagina 26 de aula_02.pdf
print("====================")
Call:
ar(x = gnp, order.max = 15)
Coefficients:
1 2 3
0,3463 0,1770 -0,1421
Order selected 3 sigma^2 estimated as 0,9676
[1] "===================="
0 1 2 3 4 5 6
27,5691310 2,6081086 1,5895550 0,0000000 0,2734771 2,2034466 4,0171066
7 8 9 10 11 12 13
5,9916210 5,8264833 7,5230025 7,8223499 9,5813222 7,3983158 8,9432291
14 15
10,9115648 12,8950445
[1] "===================="
In [7]:
#==============================================================================
# Calculando varios momentos amostrais dos residuos padronizados do
# modelo AR(3)
# Ver pagina 27 de aula_02.pdf; verificacao do modelo
#==============================================================================
options(digits=4)
a <- gnp.ar$resid[(p+1):n] # residuos
sigma2 <- gnp.ar$var # variancia do erro estimada
sigma <- sqrt(sigma2) # desvio padrao do erro estimado
resid <- a/sigma # residuos padronizados
gnp.describe <- describe(resid)
gnp.media <- gnp.describe[,3]
gnp.desvpad <- gnp.describe[,4]
gnp.mediana <- gnp.describe[,5]
gnp.assimetria <- gnp.describe[,11]
gnp.exc.curtose <- gnp.describe[,12]
STATS <- cbind(gnp.media,gnp.desvpad,gnp.mediana,gnp.assimetria,
gnp.exc.curtose)
print(STATS)
gnp.media gnp.desvpad gnp.mediana gnp.assimetria gnp.exc.curtose [1,] -0,003631 0,9971 -0,02014 0,04556 0,6847
In [8]:
# Se os erros forem normais, entao os residuos devem ter
# media 0 (zero), mediana 0 (zero), desvio padrao 1 (um),
# assimetria 0 (zero) e excesso de curtose 0 (zero).
In [9]:
#==============================================================================
# Verificando as hipoteses do modelo - ACF
#==============================================================================
par(mfrow=c(1,1),lwd=2.0,cex.lab=1.5,cex.axis=1.5,lab=c(10,10,5),
mar=c(5,5,2,2.5),cex.main=2.0,bty="n")
acf(resid,lwd=2,col="black",lag.max=16,xlab="Defasagem",main="",
xlim=c(0,16),ylim=c(-0.2,1))
In [10]:
#==============================================================================
# Teste de Ljung-Box
# Lembre-se da hipotese nula (H_0): independencia da serie
# Na verdade o teste e' para nao correlacao. Caso os dados sejam normais,
# correlacao nula implica independencia.
#==============================================================================
Box.test(resid,lag = floor(log(n)), type = "Ljung-Box",fitdf = p)
Box.test(resid,lag = 10, type = "Ljung-Box",fitdf = p)
Box-Ljung test data: resid X-squared = 2,1, df = 2, p-value = 0,3
Box-Ljung test data: resid X-squared = 7, df = 7, p-value = 0,4
In [11]:
# Conclusao: NAO rejeitamos a hipotese nula de correlacoes nulas da serie de
# residuos
In [12]:
#==============================================================================
# Teste de Jarque-Bera
# Lembre-se da hipotese nula (H_0): normalidade da serie
#==============================================================================
jarque.bera.test(resid)
Jarque Bera Test data: resid X-squared = 3,9, df = 2, p-value = 0,1
In [13]:
# Conclusao: NAO rejeitamos a hipotese nula de normalidade da serie
# residual
In [14]:
#==============================================================================
# Teste de significancia dos parametros do modelo AR(3)
#==============================================================================
# qnorm(1-0.1/2)
IC1 <- cbind(gnp.ar$ar - 1.64*sqrt(diag(gnp.ar$asy.var.coef)),
gnp.ar$ar + 1.64*sqrt(diag(gnp.ar$asy.var.coef))) # IC de aprox. 90%
colnames(IC1) <- c("5%","95%")
rownames(IC1) <- c("phi[1]","phi[2]","phi[3]")
print(IC1)
print("====================")
# qnorm(1-0.05/2)
IC2 <- cbind(gnp.ar$ar - 1.96*sqrt(diag(gnp.ar$asy.var.coef)),
gnp.ar$ar + 1.96*sqrt(diag(gnp.ar$asy.var.coef))) # IC de aprox. 95%
colnames(IC2) <- c("2,5%","97,5%")
rownames(IC2) <- c("phi[1]","phi[2]","phi[3]")
print(IC2)
print("====================")
5% 95%
phi[1] 0,22247 0,47003
phi[2] 0,04771 0,30622
phi[3] -0,26587 -0,01831
[1] "===================="
2,5% 97,5%
phi[1] 0,19832 0,494186
phi[2] 0,02249 0,331441
phi[3] -0,29002 0,005846
[1] "===================="
In [15]:
#==============================================================================
# Previsao
# Ver paginas 28 a 32 de aula_02.pdf
#==============================================================================
objprev <- forecast(gnp.ar,h=12)
gnp.prev <- ts(c(objprev[4]$mean),start=c(1991,2),frequency=4)
gnp.ipli <- ts(c(objprev[5]$lower[,2]),start=c(1991,2),frequency=4)
gnp.ipls <- ts(c(objprev[6]$upper[,2]),start=c(1991,2),frequency=4)
objprev
Point Forecast Lo 80 Hi 80 Lo 95 Hi 95 1991 Q2 0,1325 -1,1281 1,393 -1,795 2,060 1991 Q3 0,4666 -0,8674 1,801 -1,574 2,507 1991 Q4 0,7565 -0,6290 2,142 -1,362 2,875 1992 Q1 0,8048 -0,5810 2,191 -1,315 2,924 1992 Q2 0,8253 -0,5606 2,211 -1,294 2,945 1992 Q3 0,7998 -0,5868 2,186 -1,321 2,920 1992 Q4 0,7877 -0,5989 2,174 -1,333 2,908 1993 Q1 0,7761 -0,6106 2,163 -1,345 2,897 1993 Q2 0,7736 -0,6132 2,160 -1,347 2,894 1993 Q3 0,7724 -0,6144 2,159 -1,348 2,893 1993 Q4 0,7731 -0,6136 2,160 -1,348 2,894 1994 Q1 0,7735 -0,6132 2,160 -1,347 2,894
In [16]:
#==============================================================================
# Grafico da serie temporal com previsao
#==============================================================================
par(mfrow=c(1,1),lwd=2.0,cex.lab=1.3,cex.axis=1.3,lab=c(10,10,5),
mar=c(5,5,2,2.5),xpd=T,cex.main=2.0,bty="n")
ts.plot(gnp,xlab="Trimestre",ylab="PNB",xlim=c(1945,1995),ylim=c(-3,5),main="")
points(gnp,pch=15)
points(gnp.prev,pch=15,col="red")
points(gnp.ipli,pch=15,col="blue")
points(gnp.ipls,pch=15,col="blue")
lines(gnp.ipli,pch=15,col="blue")
lines(gnp.ipls,pch=15,col="blue")
