#  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   CAPÍTULO 6   %%%%%%%%%%%%%%%%%%%%%%%%%%%%% #
################################################################################
# Identificação de Modelos
################################################################################
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% # %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
               
# Identificando série simuladas

# Obs: O R calcula o IC baseado em uma série não correlacionada, por isso deve 
# ser interpretada com cuidado!!

# Sugiro usar a opção ci.type = "ma" para um processo em que as autocorrelações 
# são nulas para v > q (MA(q)).
ar.sim <- arima.sim(n=10, list(ar=c(0.3)))
par(mfrow=c(1,2))
acf.ar <- acf(ar.sim)
acf.ar <- acf(ar.sim, ci.type="ma")
pacf(ar.sim)
plot(acf.ar,ci.type="ma",ci=0.95)
plot(acf.ar,ci.type="ma",ci=0.90)
plot(acf.ar,ci.type="ma",ci=0.99)


ar.sim <- arima.sim(n=10000, list(ar=c(0.3)))
par(mfrow=c(1,2))
acf.ar <- acf(ar.sim, ci.type="ma")
pacf(ar.sim)

ar.sim <- arima.sim(n=10000, list(ar=c(-0.2, 0.5)))
acf.ar <- acf(ar.sim, ci.type="ma")
pacf(ar.sim)

ar.sim <- arima.sim(n=10000, list(ma=c(-0.2, 0.5)))
acf.ar <- acf(ar.sim, ci.type="ma")
pacf(ar.sim)

ar.sim <- arima.sim(n=10000, list(ma=c(-0.2)))
acf.ar <- acf(ar.sim, ci.type="ma")
pacf(ar.sim)

ar.sim <- arima.sim(n=10000, list(ar=c(0.7,0.1), ma=c(-0.2, 0.5)))
acf.ar <- acf(ar.sim, ci.type="ma")
pacf(ar.sim)

ar.sim <- arima.sim(n=10000, list(ar=c(0.99)))
acf.ar <- acf(ar.sim, ci.type="ma")
pacf(ar.sim)

ar.sim <- arima.sim(n=10000, list(ar=c(0.1, 0.3, -0.4)))
acf.ar <- acf(ar.sim, ci.type="ma")
pacf(ar.sim)

ar.sim <- arima.sim(n=10000, list(ar=c(-0.1, 0.3, -0.4, 0.2, 0.17)))
acf.ar <- acf(ar.sim, ci.type="ma")
pacf(ar.sim)

ar.sim <- arima.sim(n=10000, list(ma=c(0.1, -0.3, 0.4)))
acf.ar <- acf(ar.sim, ci.type="ma")
pacf(ar.sim)

ar.sim <- arima.sim(n=50000, list(ar=c(0.1, 0.3, -0.4), ma=c(0.1, -0.5, 0.3)))
acf.ar <- acf(ar.sim, ci.type="ma")
pacf(ar.sim)


################################################################################
# Exemplo 6.3 - série ICV
################################################################################
ICV <- read.csv2("C:\\dados\\ICV.csv", head=T, dec=".")
ICV
icv<-ICV[,2]
plot(icv, main="ICV de jan-70 a jun-79", lwd = "2", col = "51", xlab = "anos", 
ylim = c(10,850))
summary(icv)
length(icv)

# ------------------------------------------------------------------------------
# Verificação da Variabilidade para n pequeno(médias móveis)
# ------------------------------------------------------------------------------
n <- 114
j <- floor(n/13);j
abscissa <- NULL
ordenada <- NULL
for (i in 1: (n-j)){
abscissa[i] <- mean(icv[i:(i+j)])
ordenada[i] <- sd(icv[i:(i+j)])
}
plot(abscissa,ordenada, xlab='Média', ylab= 'Desvio-padrão', 
main='Verificação da Variabilidade')

# ------------------------------------------------------------------------------
# Verificação da Variabilidade para n grande (intervalo de faixa)
# ------------------------------------------------------------------------------
n <- 114  # número de observações
g <- 8    # número de observações em cada grupo
n1 <- n-g
g1 <- g-1
j <- floor(n/g) # número de grupos
SEQ <- seq (1, n1, by=g)
abscissa1 <- NULL
ordenada1 <- NULL
for (i in SEQ){
abscissa1[i] <- mean(icv[i:(i+g1)])
ordenada1[i] <- sd(icv[i:(i+g1)])
}
plot(x=abscissa1,y=ordenada1, xlab='Média', ylab= 'Desvio-padrão', 
main='Verificação da Variabilidade')

# ------------------------------------------------------------------------------
# Transformando os dados
# ------------------------------------------------------------------------------
logicv <- log(icv)
par(mfrow=c(2,1))
ts.plot(icv, type = "l", main="série sem transformação")
ts.plot(logicv, type = "l", main="série com transformação")

# ------------------------------------------------------------------------------
# Identificando o modelo
# ------------------------------------------------------------------------------
par(mfrow=c(3,2))
acf(logicv, main = "log icv")
pacf(logicv, main = "log icv")

d1logicv <- diff(logicv)
acf(d1logicv, main = "1a diferença do log icv")
pacf(d1logicv, main = "1a diferença do log icv")

d2logicv <- diff(d1logicv)
acf(d2logicv, main = "2a diferença do log icv")
pacf(d2logicv, main = "2a diferença do log icv")

# ------------------------------------------------------------------------------
# Quantas diferenças?
# ------------------------------------------------------------------------------
# veja a FAC - 1o lag: ~ -0.5
# veja a FACP - mudança de sinais dos valores da FACP
# e quanto a variância?

# Variância da 1a diferença
vard1logicv <- var(d1logicv)
vard1logicv

# Variância da 2a diferença
vard2logicv <- var(d2logicv)
vard2logicv

# ------------------------------------------------------------------------------
# Teste de Raiz Unitária (Dickey-Fuller) 
# H0: há raiz unitária
# ------------------------------------------------------------------------------
require(urca)
?ur.df
logicv <- log(icv)
logicv1 <- ts(logicv)
# teste com intercepto, sem tendência linear e sem dummies sazonais
teste1 <- ur.df(logicv1, type=c("drift"), lags=0)
teste2 <- ur.df(d1logicv1, type=c("drift"), lags=0)
summary(teste1)
summary(teste2)

# H0 not rejected

# teste com intercepto, sem tendência linear e sem dummies sazonais
d1logicv <- diff(logicv)
d1logicv1 <- ts(d1logicv)
teste2 <- ADF.test (wts=d1logicv1,selectlags=list(Pmax=1),itsd=c(1,0,c(0)), regvar=0)
teste2
# H0 rejected at 1% significancy, thus d=1

# ------------------------------------------------------------------------------
# Testar se a média é estatisticamente diferente de zero (t-ratio)
# ------------------------------------------------------------------------------
var_d1logicv <- var(d1logicv)
r1 <- acf(d1logicv)
r1$acf[2]
varW_bar <- var_d1logicv*(1+r1$acf[2])/(113*(1-r1$acf[2]))
varW_bar
t_ratio <- mean(d1logicv)/(varW_bar)^0.5;t_ratio
# H0 é rejeitado

# ------------------------------------------------------------------------------
# UM MODELO PRELIMINAR É ARIMA(1,1,0) com uma constante
# (1 - B)(1 - phi_1*B)Yt = theta_0 + at, com Yt = log(icv)
# ------------------------------------------------------------------------------


################################################################################
# Exemplo 6.4 - série Umidade
################################################################################
UMID <- read.csv2("C:\\dados\\atmosfera.csv", head=T, dec=".")
UMID
umid <- UMID [1:358,3];umid
temp <- UMID [,2]
summary(umid)
ts.plot(umid, lwd = "2", col = "1", xlab = "dias")
length(umid)
# ------------------------------------------------------------------------------
# Verificação da Variabilidade para n grande (intervalo de faixa)
# ------------------------------------------------------------------------------
n <- 358  # total number of observations
g <- 7    # number of observations in each group
g1 <- g-1
j <- floor(n/g) # number of groups
SEQ <- seq (1, n, by=g)
abscissa1 <- NULL
ordenada1 <- NULL
for (i in SEQ){
abscissa1[i] <- mean(umid[i:(i+g1)])
ordenada1[i] <- sd(umid[i:(i+g1)])
}
abscissa1
ordenada1

abscissa11 <- abscissa1[1:358]
ordenada11 <- ordenada1[1:358]

plot(x=abscissa11,y=ordenada11, xlab='Média', ylab= 'Desvio-padrão',
main='Verificação da Variabilidade')


# funcao BoxCox.ar e BoxCox
help.search("boxcox")
require(forecast)
require(TSA)
?BoxCox
?BoxCox.ar
umid1 <- as.ts(umid)
umid1
plot(umid1)

BoxCox.ar(umid1)
par(mfrow=c(2,1))
plot(BoxCox(umid1,2))
plot(umid1)

# ------------------------------------------------------------------------------
# Aparentemente sem necessidade de transformar os dados
# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
# Identificando o modelo
# ------------------------------------------------------------------------------
par(mfrow=c(3,2))
acf(umid)
pacf(umid)

d1umid <- diff(umid)
acf(d1umid, main = "1a diferença")
pacf(d1umid, main = "1a diferença")

d2umid <- diff(d1umid)
acf(d2umid, main = "2a diferença")
pacf(d2umid, main = "2a diferença")

# ------------------------------------------------------------------------------
# Quantas diferenças?
# ------------------------------------------------------------------------------
# veja a FAC - 1o lag: ~ -0.5
# veja a FACP - mudança de sinais dos valores da FACP

# ------------------------------------------------------------------------------
# Teste de Raiz Unitária (Dickey-Fuller) 
# H0: há raiz unitária
# ------------------------------------------------------------------------------
require(uroot)
ADF.test
umid1 <- ts(umid)
# teste com intercepto, sem tendência linear e sem dummies sazonais
teste <- ADF.test (wts=umid1,selectlags=list(Pmax=1),itsd=c(0,0,c(0)), regvar=0)
teste
# H0 rejected

# ------------------------------------------------------------------------------
# Testar se a média é estatisticamente diferente de zero (t-ratio)
# H0: E(W) = 0.
# ------------------------------------------------------------------------------
var_umid <- var(umid)
r1 <- acf(umid)
r1$acf[2]
varW_bar <- var_umid*(1+r1$acf[2])/(357*(1-r1$acf[2]))
varW_bar
t_ratio <- mean(umid)/(varW_bar)^0.5;t_ratio
# H0 é rejeitado

# ------------------------------------------------------------------------------
# UM MODELO PRELIMINAR É ARIMA(1,0,0) com uma constante
# (1 - phi_1*B)Zt = theta_0 + at
# ------------------------------------------------------------------------------


################################################################################
# Exemplo 6.5 - série CO
################################################################################
CO <- read.csv2("C:\\dados\\poluicao.csv", head=T, dec=".")
CO
names(CO)
co <- CO [1:358,3];co
summary(co)
ts.plot(co, lwd = "2", col = "1", xlab = "dias")
length(co)

# ------------------------------------------------------------------------------
# Verificação da Variabilidade para n grande (intervalo de faixa)
# ------------------------------------------------------------------------------
n <- 358  # total number of observations
g <- 7    # number of observations in each group
g1 <- g-1
j <- floor(n/g) # number of groups
SEQ <- seq (1, n, by=g)
abscissa1 <- NULL
ordenada1 <- NULL
for (i in SEQ){
abscissa1[i] <- mean(co[i:(i+g1)])
ordenada1[i] <- sd(co[i:(i+g1)])
}
abscissa1
ordenada1

abscissa11 <- abscissa1[1:358]
ordenada11 <- ordenada1[1:358]

plot(x=abscissa11,y=ordenada11, xlab='Média', ylab= 'Desvio-padrão',
main='Verificação da Variabilidade')

# ------------------------------------------------------------------------------
# Transformação dos dados
# ------------------------------------------------------------------------------
logco <- log(co)
ts.plot(logco)

# ------------------------------------------------------------------------------
# Identificando o modelo
# ------------------------------------------------------------------------------
par(mfrow=c(3,2))
acf(logco)
pacf(logco)

d1logco <- diff(logco)
acf(d1logco, main = "1a diferença")
pacf(d1logco, main = "1a diferença")

d2logco <- diff(d1logco)
acf(d2logco, main = "2a diferença")
pacf(d2logco, main = "2a diferença")

# ------------------------------------------------------------------------------
# Quantas diferenças?
# ------------------------------------------------------------------------------

# ------------------------------------------------------------------------------
# UM MODELO PRELIMINAR É ARIMA(3,0,0) com uma constante, com Yt = log(icv)
# ------------------------------------------------------------------------------



#  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   CAPÍTULO 7   %%%%%%%%%%%%%%%%%%%%%%%%%%%%% #
################################################################################
# ESTIMAÇÃO DE MODELOS ARIMA(p,d,q)
################################################################################
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% # %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%                                                                                


################################################################################
# Exemplo 7.1 EMV-Procedimento Condicional
################################################################################
Zt <- c(150,147,143,148,153,149,155,162,170,172)
Zt
# Ajuste de um modelo IMA(1,1)
ima <- arima(Zt,order=c(0,1,1));ima
# Obs: o mínimo de S*(theta) ocorre ao redor de theta = 0.4.
# Variância do parâmetro
var_theta <- (1-0.4^2)/9
ima$var.coef


################################################################################
# Exemplo 7.3 Estimativa dos parâmetros das séries simuladas
################################################################################

# ------------------------------------------------------------------------------
# Processo AR(1)
# ------------------------------------------------------------------------------
ar.sim <- arima.sim(n=10000, list(ar=c(0.8)))
ar_fit <- arima(ar.sim, order=c(1,0,0), include.mean=F);ar_fit
#significância do parâmetro - altamente significativo
0.7936/0.0061

# ------------------------------------------------------------------------------
# Processo MA(1)
# ------------------------------------------------------------------------------
ma.sim <- arima.sim(n=10000, list(ma=c(-0.8)))
ma_fit <- arima(ma.sim, order=c(0,0,1), include.mean=F);ma_fit
#significância do parâmetro - altamente significativo
-0.7919/0.0059

# ------------------------------------------------------------------------------
# Processo ARMA(1,1)
# ------------------------------------------------------------------------------
arma.sim <- arima.sim(n=10000, list(ar=c(0.8), ma=c(-0.3)))
arma_fit <- arima(arma.sim, order=c(1,0,1), include.mean=F);arma_fit
#significância dos parâmetros - altamente significativos
0.8043/0.009
-0.2954/0.0146


################################################################################
# Exemplo 7.4 Retornando ao exemplo 6.3 - série ICV
################################################################################

# ------------------------------------------------------------------------------
# UM MODELO PRELIMINAR É ARIMA(1,1,0) com uma constante
# (1 - B)(1 - phi_1*B)Yt = theta_0 + at, com Yt = log(icv) 
# ------------------------------------------------------------------------------

# arima_icv <- arima(logicv, order=c(1,1,0)); arima_icv
# A função "arima" não permite incluir uma constante quando d=1! Optar pela 
# opção "xreg"!!

theta0 <- seq(1,length(logicv),1)
arima_icv <- arima(logicv,order=c(1,1,0), xreg = theta0);arima_icv
# O que o R chama de "intercept" nos resultados é a média do processo e não a 
# constante (intercepto).

arima_icv1 <- arima(diff(logicv), order=c(1,0,0), include.mean=T); arima_icv1
# Podemos optar por x = diff(logicv), mas a previsão resultará em valores
# diferençados!!

# ------------------------------------------------------------------------------
# Cálculo da constante
# ------------------------------------------------------------------------------
theta0 <- 0.0212*(1-0.5073); theta0

# ------------------------------------------------------------------------------
# Cálculo do erro padrão
# ------------------------------------------------------------------------------
ep <- (arima_icv1$var.coef)^0.5;ep

# ------------------------------------------------------------------------------
# Significância dos parâmetros
# ------------------------------------------------------------------------------
# phi1 
0.5073/0.0822
# intercept
0.0212/0.0018

# ------------------------------------------------------------------------------
# Outra forma de estimar modelos ARIMA(p,d,q)
# ------------------------------------------------------------------------------
require(tseries)
arma_icv <- arma(d1logicv,order=c(1,0), include.intercept=T); arma_icv
# nessa função o resultado "intercept" é de fato a constante!


################################################################################
# Exemplo 7.5 Retornando ao exemplo 6.4 - série UMIDADE
################################################################################

# ------------------------------------------------------------------------------
# UM MODELO PRELIMINAR É ARIMA(1,0,0) com uma constante
# (1 - phi_1*B)Zt = theta_0 + at
# ------------------------------------------------------------------------------

arima_umid <- arima(umid, order=c(1,0,0), include.mean=T); arima_umid

# ------------------------------------------------------------------------------
# Cálculo da constante
# ------------------------------------------------------------------------------
theta0 <- arima_umid$coef[2]*(1-arima_umid$coef[1]); theta0

# ------------------------------------------------------------------------------
# Cálculo do erro padrão
# ------------------------------------------------------------------------------
ep <- (arima_umid$var.coef)^0.5;ep

# ------------------------------------------------------------------------------
# Significância dos parâmetros
# ------------------------------------------------------------------------------
# phi1 
arima_umid$coef[1]/ep[1,1]
# intercept
arima_umid$coef[2]/ep[2,2]

# ------------------------------------------------------------------------------
# Outra forma de estimar modelos ARIMA(p,d,q)
# ------------------------------------------------------------------------------
arma_umid <- arma(umid,order=c(1,0), include.intercept=T); arma_umid


################################################################################
# Exemplo 7.6 Retornando ao exemplo 6.5 - série CO
################################################################################

# ------------------------------------------------------------------------------
# UM MODELO PRELIMINAR É ARIMA(3,0,0) com uma constante
# ------------------------------------------------------------------------------
arima_co <- arima(logco, order=c(3,0,0), include.mean=T); arima_co

# ------------------------------------------------------------------------------
# Cálculo da constante
# ------------------------------------------------------------------------------
theta0 <- arima_co$coef[4]*(1-arima_co$coef[1]-arima_co$coef[2]-arima_co$coef[3])
theta0

# ------------------------------------------------------------------------------
# Cálculo do erro padrão
# ------------------------------------------------------------------------------
ep <- (arima_co$var.coef)^0.5;ep

# ------------------------------------------------------------------------------
# Significância dos parâmetros
# ------------------------------------------------------------------------------
# phi1 
arima_co$coef[1]/ep[1,1]
# phi2 
arima_co$coef[2]/ep[2,2]
# phi3 
arima_co$coef[3]/ep[3,3]
# intercept
arima_co$coef[4]/ep[4,4]

# ------------------------------------------------------------------------------
# Outra forma de estimar modelos ARIMA(p,d,q)
# ------------------------------------------------------------------------------
arma_co <- arma(logco,order=c(3,0), include.intercept=T); arma_co