Trabalho de Séries Temporais
José Luiz de Carvalho Lisboa
25 de Junho
Descrição dos Dados
Para modelar o índice de custo de vida no Estado de São Paulo, utilizou-se o modelo de suavização exponencial de Holt. Onde os dados são referentes a 126 informações no período de Janeiro de 1970 a Junho de 1980. O modelo de suavização exponencial de Holt-Winters foi utilizado para modelar o índice de produto industrial, estes consistem em 187 observações período de Janeiro de 1985 a Julho de 2000.
##----------------------------------------------------##
## INSERINDO O BANCO DE DADOS ##
##----------------------------------------------------##
dados <- read.csv2(file.choose(), sep = ";", dec = ",", header = T)
attach(dados)
##----------------------------------------------------##
## PLOTANDO A SÉRIE DOS DADOS ##
##----------------------------------------------------##
plot(ts(ICV, start = c(1970, 1), frequency = 12), type = "l", col = "blue",
xlab = "Ano", ylab = "Quantidade", las = 1)
abline(v = seq(1971, 1980, 2), lty = 2)
##----------------------------------------------------##
## PLOTANDO A SÉRIE SELECIONADA ##
##----------------------------------------------------##
plot(ts(ICV[1:114], start = c(1970, 1), frequency = 12), type = "l", col = "blue", xlab = "Ano",
ylab = "Quantidade", las = 1)
abline(v = seq(1971, 1979, 2), lty = 2)
##-------------------------------------------------------------##
## MODELANDO A SÉRIE PELO MÉTODO DE HOLT ##
##-------------------------------------------------------------##
t = ICV[1:114]
a = 0.9
c = 0.3
Test <- t
Zest <- t
for (i in 3:114) {
Test[1] = NA
Zest[1] = NA
Test[2] = t[2]-t[1]
Zest[2] = t[2]
Zest[i] <- a*t[i]+(1-a)*(Zest[i-1]+Test[i-1])
Test[i] <- c*(Zest[i]-Zest[i-1])+(1-c)*Test[i-1]
}
##---------------------------------------------------##
## PLOTANDO A SÉRIE MODELADA ##
##---------------------------------------------------##
plot(ICV[1:114], type = "o", xlab = "Ano",
ylab = "Quantidade", las = 1, ylim = c(0.84, 779),
pch = 21, bg = "yellow", axes = F)
lines(Zest, col = "red", lwd = 2)
lines(Test, col = "green3", lwd = 2)
axis(1, at = seq(12,114,12), labels = c("1971", "1972", "1973",
"1974", "1975", "1976",
"1977", "1978", "1979"))
axis(2, seq(0,780,100), las = 1)
box()
legend(0,800, c("ICV Obsevado", "Estimativa do Nível", "Estimativa da Tendência"),
col = c("black", "red","green3"), x.intersp = 0.3, y.intersp = 0.7,
bty = "n", lty = c(1,1,1), pt.bg = c("yellow",NA, NA), pch = c(21,NA,NA))
##--------------------------------------------------------------##
## PREVISÃO DO MODELO COM BASE EM JUNHO ##
##--------------------------------------------------------------##
h = 1:12
Zprev <- Zest[114]+h*Test[114]
##------------------------------------------------------##
## PLOTANDO A PREVISÃO DA SÉRIE ##
##------------------------------------------------------##
plot(ICV[115:126], type = "o", xlab = "Ano",
ylab = "Quantidade", las = 1, pch = 21,
bg = "yellow", axes = F)
lines(Zprev, col = "red")
axis(1, at = seq(1,12), labels = c("Jul/1979", "Ago/1979", "Set/1979",
"Out/1979", "Nov/1979", "Dez/1979",
"Jan/1980", "Fev/1980", "Ma?/1980",
"Abr/1980", "Mai/1980", "Jun/1980"))
axis(2, seq(800,1375,50), las = 1)
box()
legend(1,1350, c("ICV Observado", "ICV Previsto"), col = c("black", "red"),
lty = c(1,1), pch = c(21, NA), pt.bg = c("yellow",NA),
bty = "n", x.intersp = 0.4, y.intersp = 0.8)
##---------------------------------------------------------##
## ESTIMANDO A PREVISÃO ATUALIZADA ##
##---------------------------------------------------------##
k <- ICV[115:126]
Azest <- k
Atest <- k
Aprev <- k
for (j in 2:12) {
Azest[1] = a*k[1]+(1-a)*(Zest[114]+Test[114])
Atest[1] = c*(Azest[1]-Zest[114])+(1-c)*Test[114]
Aprev[1] = Zest[114]+Test[114]
Azest[j] <- a*k[j]+(1-a)*(Azest[j-1]+Atest[j-1])
Atest[j] <- c*(Azest[j]-Azest[j-1])+(1-c)*Atest[j-1]
Aprev[j] <- Azest[j-1]+Atest[j-1]
}
##-----------------------------------------------------------------##
## PLOTANDO A PREVISÃO DA SÉRIE ATUALIZADA ##
##-----------------------------------------------------------------##
plot(ICV[115:126], type = "o", xlab = "Ano",
ylab = "Quantidade", las = 1, pch = 21,
bg = "yellow", axes = F)
lines(Aprev, col = "red")
axis(1, at = seq(1,12), labels = c("Jul/1979", "Ago/1979", "Set/1979",
"Out/1979", "Nov/1979", "Dez/1979",
"Jan/1980", "Fev/1980", "Mar/1980",
"Abr/1980", "Mai/1980", "Jun/1980"))
axis(2, seq(800, 1375, 50), las = 1)
box()
legend(1,1350, c("ICV Observado", "Previsão"), col = c("black", "red"),
lty = c(1,1), pt.bg = c("yellow", NA), pch = c(21,NA), bty = "n",
x.intersp = 0.4, y.intersp = 0.8)
##-----------------------------------------------------##
## PLOTANDO O MODELO DA SÉRIE ##
##-----------------------------------------------------##
Zmod <- Zest+Test
Zprev1 <- c(rep(NA,113), Zmod[114], Zprev)
plot(ICV[1:126], type = "o", xlab = "Ano",
ylab = "Quantidade", las = 1,
pch = 21, bg = "yellow", axes = F)
lines(Zmod, col = "red", lwd = 2)
lines(Zprev1, col = "green3", lwd = 2)
axis(1, at = seq(12,120,12), labels = c("1971", "1972", "1973",
"1974", "1975", "1976",
"1977", "1978", "1979",
"1980"))
axis(2, seq(0,1380,100), las = 1)
box()
legend(0,1300, c("ICV Obsevado", "Modelo Estimado", "Previsão Direta"),
col = c("black", "red","green3"), x.intersp = 0.3, y.intersp = 0.7,
bty = "n", lty = c(1,1,1), pt.bg = c("yellow",NA, NA), pch = c(21,NA,NA))
##----------------------------------------------------------------##
## PLOTANDO O MODELO DA SÉRIE ATUALIZADO ##
##----------------------------------------------------------------##
Zaprev <- c(rep(NA,113), Zmod[114], Aprev)
plot(ICV[1:126], type = "o", xlab = "Ano",
ylab = "Quantidade", las = 1,
pch = 21, bg = "yellow", axes = F)
lines(Zmod, col = "red", lwd = 2)
lines(Zaprev, col = "green3", lwd = 2)
axis(1, at = seq(12,120,12), labels = c("1971", "1972", "1973",
"1974", "1975", "1976",
"1977", "1978", "1979",
"1980"))
axis(2, seq(0,1380,100), las = 1)
box()
legend(0,1300, c("ICV Obsevado", "Modelo Estimado", "Previsão Atualizada"),
col = c("black", "red","green3"), x.intersp = 0.3, y.intersp = 0.7,
bty = "n", lty = c(1,1,1), pt.bg = c("yellow",NA, NA), pch = c(21,NA,NA))
##----------------------------------------------------##
## ERRO DE PREVISÃO DO MODELO ##
##----------------------------------------------------##
Erro <- ICV[115:126]-Zprev
Erroa <- ICV[115:126]-Aprev
par(mfrow = c(2,1))
plot(Erro, type = "o", xlab = "Ano",
ylab = "Erro Previsão Direta", las = 1, pch = 19,
col = "red", axes = F)
axis(1, at = seq(1,12), labels = c("Jul/1979", "Ago/1979", "Set/1979",
"Out/1979", "Nov/1979", "Dez/1979",
"Jan/1980", "Fev/1980", "Mar/1980",
"Abr/1980", "Mai/1980", "Jun/1980"))
axis(2, seq(-0.40,308,80), las = 1, cex.axis = 0.7)
box()
plot(Erroa, col = "green3", type = "o", las = 1,
axes = F, pch = 19, xlab = "Ano", ylab = "Erro Previsão Atualizado")
abline(h = 0, lty = 2, col = "red")
axis(1, at = seq(1,12), labels = c("Jul/1979", "Ago/1979", "Set/1979",
"Out/1979", "Nov/1979", "Dez/1979",
"Jan/1980", "Fev/1980", "Mar/1980",
"Abr/1980", "Mai/1980", "Jun/1980"))
axis(2, seq(-18,63,40), las = 1)
box()
##--------------------------------------------------------##
## MEDIDAS DE ACURÁCIA DO MODELO ##
##--------------------------------------------------------##
EQM <- sum((ICV[115:126]-Aprev)^2)/length(ICV[115:126])
MAPE <- (sum(abs(ICV[115:126]-Aprev)/ICV[115:126])*100)/length(ICV[115:126])
EMA <- sum(abs(ICV[115:126]-Aprev))/length(ICV[115:126])