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'

#==============================================================================
# Lendo os dados do arquivo gnp.txt
#==============================================================================
gnp <- ts(scan("../dados/gnp.txt")*100,start=c(1947,2),frequency=4)
n   <- length(gnp)

#==============================================================================
# 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)

#==============================================================================
# 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))

#==============================================================================
# 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))

#==============================================================================
# Ajustando varios modelos AR(p), MA(q) e ARMA(p,q)
#==============================================================================
ajuste0  <- arima(gnp,c(0,0,0)) # Ruido branco

ajuste1  <- arima(gnp,c(1,0,0)) # AR(1)
ajuste2  <- arima(gnp,c(2,0,0)) # AR(2)
ajuste3  <- arima(gnp,c(3,0,0)) # AR(3)
ajuste4  <- arima(gnp,c(4,0,0)) # AR(4)
ajuste5  <- arima(gnp,c(5,0,0)) # AR(5)

ajuste6  <- arima(gnp,c(0,0,1)) # MA(1)
ajuste7  <- arima(gnp,c(0,0,2)) # MA(2)
ajuste8  <- arima(gnp,c(0,0,3)) # MA(3)
ajuste9  <- arima(gnp,c(0,0,4)) # MA(4)
ajuste10 <- arima(gnp,c(0,0,5)) # MA(5)

ajuste11 <- arima(gnp,c(1,0,1)) # ARMA(1,1)
ajuste12 <- arima(gnp,c(1,0,2)) # ARMA(1,2)
ajuste13 <- arima(gnp,c(2,0,1)) # ARMA(2,1)
ajuste14 <- arima(gnp,c(2,0,2)) # ARMA(2,2)
ajuste15 <- arima(gnp,c(1,0,3)) # ARMA(1,3)
ajuste16 <- arima(gnp,c(3,0,1)) # ARMA(3,1)
ajuste17 <- arima(gnp,c(3,0,2)) # ARMA(3,2)
ajuste18 <- arima(gnp,c(2,0,3)) # ARMA(2,3)
ajuste19 <- arima(gnp,c(3,0,3)) # ARMA(3,3)

#==============================================================================
# Comparando os modelos via AIC
#==============================================================================
options(digits=10)
aic <- as.matrix(c(
ajuste1$aic,
ajuste2$aic,
ajuste3$aic,
ajuste4$aic,
ajuste5$aic,
ajuste6$aic,
ajuste7$aic,
ajuste8$aic,
ajuste9$aic,
ajuste10$aic,
ajuste11$aic,
ajuste12$aic,
ajuste13$aic,
ajuste14$aic,
ajuste15$aic,
ajuste16$aic,
ajuste17$aic,
ajuste18$aic,
ajuste19$aic),19,1)
print(aic)

             [,1]
 [1,] 502,0772288
 [2,] 500,9401382
 [3,] 499,3350539
 [4,] 499,6572820
 [5,] 501,5917081
 [6,] 510,1888111
 [7,] 498,7314890
 [8,] 498,5249296
 [9,] 500,5151940
[10,] 499,2796709
[11,] 502,4087924
[12,] 499,2244183
[13,] 501,4580904
[14,] 498,0275216
[15,] 500,5218788
[16,] 499,8966825
[17,] 498,3545760
[18,] 498,8064330
[19,] 500,2774912

#==============================================================================
# Resultado da comparacao via AIC
# O ajuste14, ARMA(2,2), teve melhor ajuste segundo o AIC
#==============================================================================
ajuste14

Call:
arima(x = gnp, order = c(2, 0, 2))

Coefficients:
            ar1         ar2         ma1        ma2  intercept
      0,6086893  -0,4539728  -0,2985213  0,5992405  0,7700619
s.e.  0,1622047   0,1683798   0,1378623  0,1711154  0,1113143

sigma^2 estimated as 0,9242543:  log likelihood = -243,01,  aic = 498,03

#==============================================================================
# Teste de significancia dos parametros do modelo AR(2,2)
#==============================================================================
# qnorm(1-0.05/2)
IC <- cbind(ajuste14$coef - 1.96*sqrt(diag(ajuste14$var.coef)),
      ajuste14$coef + 1.96*sqrt(diag(ajuste14$var.coef))) # IC de aprox. 95%
colnames(IC) <- c("2,5%","97,5%")
rownames(IC) <- c("AR1","AR2","MA1","MA2","CST")
print(IC)
print("====================")

             2,5%          97,5%
AR1  0,2907681064  0,92661046261
AR2 -0,7839971390 -0,12394847304
MA1 -0,5687313562 -0,02831121789
MA2  0,2638542913  0,93462676591
CST  0,5518859785  0,98823790495
[1] "===================="

#==============================================================================
# Calculando varios momentos amostrais dos residuos padronizados do modelo 
# ARMA(2,2)
#==============================================================================
options(digits=4)
a               <- ajuste14$res    # residuos
sigma2          <- ajuste14$sigma2 # 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,] -7,857e-05       1,003     -0,0336        0,02503          0,5297

# 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).

#==============================================================================
# Analisando as hipoteses do modelo
#==============================================================================
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,1,2.5),cex.main=2.0,bty="n")
acf(c(resid),xlab="Defasagem",main="",lag.max=20,ylim=c(-0.2,1))

#==============================================================================
# Analisando as hipoteses do modelo
#==============================================================================
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,1,2.5),cex.main=2.0,bty="n")
pacf(c(resid),xlab="Defasagem",main="",lag.max=20,xlim=c(0,20))

#==============================================================================
# Analisando as hipoteses do modelo
#==============================================================================
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")
hist(resid,xlab="Resíduos",main="",ylab="Frequência",col="lightblue")

#==============================================================================
# Analisando as hipoteses do modelo
#==============================================================================
par(mfrow=c(1,1),lwd=2.0,cex.lab=1.5,cex.axis=1.5,lab=c(8,10,5),
     mar=c(5,5,2,2.5),cex.main=2.0,bty="n")
boxplot(resid,ylab="Resíduos",col="lightblue",ylim=c(-3.5,3.5),main="")

#==============================================================================
# Analisando as hipoteses do modelo
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")
qqnorm(resid,main="",xlab="Percentis teóricos",ylab="Percentis amostrais",
       pch=15,col="red",xlim=c(-4,4),ylim=c(-4,4))
qqline(resid,lwd=2)

#==============================================================================
# Teste de Ljung-Box
# Lembre-se da hipotese nula (H_0): independencia da serie
#==============================================================================
pq <- ajuste14$arma[1] + ajuste14$arma[2] # ARMA(2,2), pq <- 2+2
Box.test(resid,lag = floor(log(n)), type = "Ljung-Box",fitdf = pq)
print("====================")
Box.test(resid,lag = 10, type = "Ljung-Box",fitdf = pq) 
print("====================")

	Box-Ljung test

data:  resid
X-squared = 0,98, df = 1, p-value = 0,3

[1] "===================="

	Box-Ljung test

data:  resid
X-squared = 5, df = 6, p-value = 0,5

[1] "===================="

# Conclusao: NAO rejeitamos a hipotese nula de correlacoes nulas
# da serie residual

#==============================================================================
# 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 = 2,4, df = 2, p-value = 0,3

# Conclusao: NAO rejeitamos a hipotese nula de normalidade da serie de 
# residuos

#==============================================================================
# Previsao com o modelo ARMA(2,2)
#==============================================================================
objprev  <- forecast(ajuste14,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)

#==============================================================================
# Grafico da serie temporal com previsao via modelo ARMA(2,2)
#==============================================================================
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")