################################################################################
# Livro "ANÁLISE DE SÉRIES TEMPORAIS" - Morettin e Toloi (2004)
################################################################################

                              --- & ---

#  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   CAPÍTULO 1   %%%%%%%%%%%%%%%%%%%%%%%%%%%%% #
################################################################################
# INTRODUÇÃO
################################################################################
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% # %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

# PACOTES
require(lmtest)
require(tseries)
require(forecast)

################################################################################
# Explorando a Série ICV
################################################################################
dados5 <- read.table('D-BANESPA.txt');dados5
ICV <- read.table('D-BANESPA.txt', head=F)
dim(icv)
icv <- ICV[,1]
summary(icv)
ts.plot(icv)

icv<-ICV[,1]
plot(icv)
plot(icv, type="l")
plot(icv, type="l", lwd=2)
plot(icv, type="l", lwd=1.5, col="blue")
plot(icv, type="l", lwd=2, col="blue", xlab="Anos")
plot(icv, type="l", lwd=2, col="blue", xlab="Anos", ylab="ICV", ylim = c(10,800))

################################################################################
# Testing Normality
################################################################################
dev.off()
x11()
hist(icv, prob=T, col="red")
lines(density(icv), lwd=2)

################################################################################
# Fitting a regression model
################################################################################
t<-seq(1970,1979,length=length(icv))
t
t2<-t^2
reg<-lm(icv~t+t2)
reg
names(reg)
reg$fit
anova(reg)
summary(reg)
plot(icv, type="l")
lines(reg$fit,col=2,lwd=2)
resicv<-resid(reg)
plot(resicv)

################################################################################
# Checking the variability
################################################################################
icv1 <- (icv[1:8])
icv1
m1 <- mean(icv1)
d1 <- var(icv1)^.5
w1 <- (max(icv1)-min(icv1))

icv2 <- (icv[9:16])
icv2
m2 <- mean(icv2)
d2 <- var(icv2)^.5
w2 <- (max(icv2)-min(icv2))

icv3 <- (icv[17:24])
icv3
m3 <- mean(icv3)
d3 <- var(icv3)^.5
w3 <- (max(icv3)-min(icv3))

icv4 <- (icv[25:32])
icv4
m4 <- mean(icv4)
d4 <- var(icv4)^.5
w4 <- (max(icv4)-min(icv4))

icv5 <- (icv[33:40])
icv5
m5 <- mean(icv5)
d5 <- var(icv5)^.5
w5 <- (max(icv5)-min(icv5))

icv6 <- (icv[41:48])
icv6
m6 <- mean(icv6)
d6 <- var(icv6)^.5
w6 <- (max(icv6)-min(icv6))

icv7 <- (icv[49:56])
icv7
m7 <- mean(icv7)
d7 <- var(icv7)^.5
w7 <- (max(icv7)-min(icv7))

icv8 <- (icv[57:64])
icv8
m8 <- mean(icv8)
d8 <- var(icv8)^.5
w8 <- (max(icv8)-min(icv8))

icv9 <- (icv[65:72])
icv9
m9 <- mean(icv9)
d9 <- var(icv9)^.5
w9 <- (max(icv9)-min(icv9))

icv91 <- (icv[73:80])
icv91
m91 <- mean(icv91)
d91 <- var(icv91)^.5
w91 <- (max(icv91)-min(icv91))

icv92 <- (icv[81:88])
icv92
m92 <- mean(icv92)
d92 <- var(icv92)^.5
w92 <- (max(icv92)-min(icv92))

icv93 <- (icv[89:96])
icv93
m93 <- mean(icv93)
d93 <- var(icv93)^.5
w93 <- (max(icv93)-min(icv93))

icv94 <- (icv[97:104])
icv94
m94 <- mean(icv94)
d94 <- var(icv94)^.5
w94 <- (max(icv94)-min(icv94))

icv95 <- (icv[105:114])
icv95
m95 <- mean(icv95)
d95 <- var(icv95)^.5
w95 <- (max(icv95)-min(icv95))

ZbarDP <- cbind(rbind(m1, m2, m3, m4, m5, m6, m7, m8, m9, m91, m92, m93, m94, m95),rbind(d1, d2, d3, d4, d5, d6, d7, d8, d9, d91, d92, d93, d94, d95))
ZbarDP

ZbarW <- cbind(rbind(m1, m2, m3, m4, m5, m6, m7, m8, m9, m91, m92, m93, m94, m95),rbind(w1, w2, w3, w4, w5, w6, w7, w8, w9, w91, w92, w93, w94, w95))
ZbarW

par(mfrow = c(1,2))
plot(ZbarDP)
plot(ZbarW)

################################################################################
# Checking the variability - moving average
################################################################################
dev.off()
x11()
j <- floor(114/12);j
n <- 114
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(x=abscissa,y=ordenada, xlab='Média', ylab= 'Desvio-padrão', 
main='Verificação da Variabilidade')

################################################################################
# Checking the variability - grouping (for large samples)
################################################################################
n <- 114  # total number of observations
g <- 12   # number of observations in each group
g1 <- g-1
j <- floor(n/g);j # number of groups
SEQ <- seq (1, n, by=g)
abscissa1 <- NULL
ordenada1 <- NULL
for (i in SEQ){
abscissa1[i] <- mean(icv[i:(i+g1)])
ordenada1[i] <- sd(icv[i:(i+g1)])
}
abscissa1
ordenada1

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

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


#  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   CAPÍTULO 5   %%%%%%%%%%%%%%%%%%%%%%%%%%%%% #
################################################################################
# MODELOS ARIMA (p,d,q)
################################################################################
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% # %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


# ------------------------------------------------------------------------------
# Modelos Auto-regressivos - AR(p)
# ------------------------------------------------------------------------------

# AR(1) com phi = 0.6
ar.sim <- arima.sim(n=10000, list(ar=c(0.6)))
ar.sim
plot(ar.sim, main="processo AR(1) simulado com phi=0.6", type="l")
par(mfrow=c(1,2))
acf.ar <- acf(ar.sim)
pacf(ar.sim)
acf.ar$acf[2]

# AR(1) com phi = -0.6
simulacao <- arima.sim(n=10000, list(ar=c(-0.6)))
acf(simulacao)
pacf(simulacao)
acf.ar$acf[2]

# AR(1) com phi = 0.1
simulacao <- arima.sim(n=10000, list(ar=c(0.1)))
acf(simulacao)
pacf(simulacao)
acf.ar$acf[2]

# AR(1) com phi = 0.01
simulacao <- arima.sim(n=10000, list(ar=c(0.01)))
ar.sim
acf(simulacao)
pacf(simulacao)
acf.ar$acf[2]

# AR(1) com phi = 0.9
simulacao <- arima.sim(n=10000, list(ar=c(0.9)))
acf(simulacao)
pacf(simulacao)
acf.ar$acf[2]

# AR(1) com phi = 0.99
simulacao <- arima.sim(n=10000, list(ar=c(0.99)))
acf(simulacao)
pacf(simulacao)
acf.ar$acf[2]

# AR(2) com phi_1 = 0.8 e phi_2 = -0.5
ar.sim <- arima.sim(n=10000, list(ar=c(0.3,-0.5)))
acf(simulacao)
pacf(simulacao)

# AR(2) com phi_1 = -0.8 e phi_2 = 0.5
simulacao <- arima.sim(n=10000, list(ar=c(-0.3,0.5)))
acf(simulacao)
pacf(simulacao)

# AR(2) com phi_1 = -0.8 e phi_2 = 0.5
# verificar condicoes de estacionariedade
ar.sim <- arima.sim(n=10000, list(ar=c(-0.8,0.5)))
acf(simulacao)
pacf(simulacao)

dev.off()
x11()

# ------------------------------------------------------------------------------
# Modelos de Médias Móveis - MA(q)
# ------------------------------------------------------------------------------
?arima.sim
# MA(1) com theta = 0.8
# Obs: o R considera o modelo com -theta
ma.sim <- arima.sim(n=10000, list(ma=c(-0.8)))
ma.sim
plot(ma.sim, main="processo MA(1) simulado com theta=-0.8", type="l")
par(mfrow=c(1,2))
acf.ma <- acf(ma.sim)
pacf(ma.sim)
acf.ma$acf[2]

# lembrando que ro_j = -theta/(1+theta^2)
ro_j <- -0.8/(1+0.8^2)
ro_j

# MA(1) com theta = 0.1
ma.sim <- arima.sim(n=10000, list(ma=c(-0.1)))
ma.sim
acf.ma <- acf(ma.sim)
pacf(ma.sim)

# MA(1) com theta = 0.9
ma.sim <- arima.sim(n=10000, list(ma=c(-0.9)))
ma.sim
acf.ma <- acf(ma.sim)
pacf(ma.sim)

# MA(1) com theta = 0.99
ma.sim <- arima.sim(n=10000, list(ma=c(-0.99)))
ma.sim
acf.ma <- acf(ma.sim)
pacf(ma.sim)

# MA(2) com theta_1 = 0.5 e theta_2 = -0.3
ma.sim <- arima.sim(n=10000, list(ma=c(-0.5,0.3)))
ma.sim
acf.ma <- acf(ma.sim)
pacf(ma.sim)

# MA(2) com theta_1 = -0.5 e theta_2 = 0.3
ma.sim <- arima.sim(n=10000, list(ma=c(0.5,-0.3)))
ma.sim
acf.ma <- acf(ma.sim)
pacf(ma.sim)

# MA(2) com theta_1 = -0.5 e theta_2 = 0.8
# verificar condicoes de invertibilidade
# Cuidado!! O R não acusa processos MA(q) NÃO invertíveis.
ma.sim <- arima.sim(n=10000, list(ma=c(0.5,-0.8)))
ma.sim
acf.ma <- acf(ma.sim)
pacf(ma.sim)

# ------------------------------------------------------------------------------
# Modelos Auto-regressivos e de Médias Móveis - ARMA(p,q)
# ------------------------------------------------------------------------------
# ARMA(1,1) com phi = 0.8 e theta = 0.3
arma.sim <- arima.sim(n=10000, list(ar=c(0.8), ma=c(-0.3)))
arma.sim
ts.plot(arma.sim)
acf.arma <- acf(arma.sim)
pacf(arma.sim)
acf.arma

# ARMA(2,1) com phi_1 = 0.7, phi_2 = 0.2 e theta = 0.3
arma.sim <- arima.sim(n=10000, list(ar=c(0.7,0.2), ma=c(0.3)))
arma.sim
ts.plot(arma.sim)
acf.arma <- acf(arma.sim)
pacf(arma.sim)

# ------------------------------------------------------------------------------
# Modelos Auto-regressivos Integrados e de Médias Móveis - ARIMA(p,d,q)
# ------------------------------------------------------------------------------
# ARIMA(1,1,0) com phi = 0.8 
ts.sim <- arima.sim(list(order = c(1,1,0), ar = 0.4), n = 1000)
ts.plot(ts.sim)
acf.arima <- acf(ts.sim)
pacf(ts.sim)
acf(diff(ts.sim))
pacf(diff(ts.sim))

# ARIMA(0,1,1) com theta = 0.3 
ts.sim <- arima.sim(list(order = c(0,1,1), ma = -0.3), n = 1000)
ts.plot(ts.sim)
acf.arima <- acf(ts.sim)
pacf(diff(ts.sim))
acf(diff(ts.sim))
pacf(diff(ts.sim))

# ARIMA(1,1,1) com phi = 0.8 e theta = 0.3
ts.sim <- arima.sim(list(order = c(1,1,1), ar = 0.8, ma = -0.3), n = 1000)
ts.plot(ts.sim)
acf.arima <- acf(ts.sim)
pacf(diff(ts.sim))
acf(diff(ts.sim))
pacf(diff(ts.sim))