##------------------------------------------------------##
## INSERINDO OS DADOS ##
##------------------------------------------------------##
dados <- read.csv(file.choose(), sep = ";", dec = ",", header = T)
attach(dados)
##-------------------------------------------------------------##
## PLOTANDO A SÉRIE TEMPORAL ##
##-------------------------------------------------------------##
plot(ts(IPI, start = c(1985, 1), frequency = 12), type = "o",
pch = 21, bg = "yellow", las = 1, ylab = "Quantidade",
xlab = "Ano", main = "Série Completa")
abline(v = seq(1986,2000,1), lty = 2, col = "gray")
##----------------------------------------------------##
## PLOTANDO A SÉRIE SELECIONADA ##
##----------------------------------------------------##
serie <- IPI[1:127]
plot(ts(serie, start = c(1985, 1), frequency = 12), type = "o",
xlab = "Ano", ylab = "Quantidade", las = 1, pch = 21, bg = "yellow",
main = "Série Selecionada")
abline(v = seq(1986, 1995, 1), lty = 2, col = "gray")
##---------------------------------------------------------------------##
## MODELANDO A SÉRIE PELO MÉTODO DE HOLT-WINTERS ##
##---------------------------------------------------------------------##
## Equações de Recorrência ##
s = 12
a = 0.97
d = 0.2
c = 0.01
# Nível
Niv <- numeric(0)
Niv[s] = mean(serie[1:s])
# Tendência
Tend <- numeric(0)
Tend[s] = 0
# Sazonalidade
Saz <- numeric(0)
Saz[1:s] = serie[1:s]/Niv[s]
for (i in (s+1):length(serie)) {
Niv[i] <- a*(serie[i]-Saz[i-s])+(1-a)*(Niv[i-1]+Tend[i-1])
Tend[i] <- c*(Niv[i]-Niv[i-1])+(1-c)*Tend[i-1]
Saz[i] <- d*(serie[i]-Niv[i])+(1-d)*Saz[i-s]
}
##--------------------------------------------------------------------##
## PLOTANDO AS ESTIMATIVAS DA SÉRIE SELECIONADA ##
##--------------------------------------------------------------------##
plot(Niv, type = "l", las = 1, col = "blue",
xlab = "Nível", ylab = "Estimativa")
plot(Saz, type = "l", las = 1, col = "blue",
xlab = "Sazonalidade", ylab = "Estimativa")
plot(Tend, type = "l", las = 1, col = "blue",
xlab = "Tendência", ylab = "Estimativa")
##------------------------------------------------##
## OBTENDO O MODELO OBTIDO ##
##------------------------------------------------##
modelo <- numeric(0)
for(i in s:length(serie)){
modelo[i] = Niv[i] + Tend[i] + Saz[i]
}
##------------------------------------------------##
## PLOTANDO O MODELO OBTIDO ##
##------------------------------------------------##
plot(serie, type = "o", xlab = "Ano", ylab = "Quantidade",
las = 1, pch = 21, bg = "yellow",
main = "Série Selecionada", axes = F)
lines(modelo, col = "red")
axis(1, at = seq(1, length(serie), 12), labels = seq(1985, 1995))
axis(2, seq(65, 127, 5), las = 1)
box()
legend(1, 126.9, c("Série dos Dados", "Modelo"), lty = c(1,1),
col = c("black", "red"), pt.bg = c("yellow", NA),
pch = c(21,NA), bty = "n", y.intersp = 0.8)
##----------------------------------------------##
## PREVISÕES DO DO MODELO ##
##----------------------------------------------##
h = length(IPI)-length(serie)
Prev <- numeric(0)
Prev[length(serie)] <- modelo[length(serie)]
for (i in 1:h) {
Prev[length(serie) + i] <- ifelse(i <= s, Niv[length(serie)]+i*Tend[length(serie)]+Saz[length(serie)+i-s],
ifelse(i > s & i <= 2*s, Niv[length(serie)]+i*Tend[length(serie)]+Saz[length(serie)+i-2*s],
ifelse(i > 2*s & i <= 3*s, Niv[length(serie)]+i*Tend[length(serie)]+Saz[length(serie)+i-3*s],
ifelse(i > 3*s & i <= 4*s, Niv[length(serie)]+i*Tend[length(serie)]+Saz[length(serie)+i-4*s],
Niv[length(serie)]+i*Tend[length(serie)]+Saz[length(serie)+i-5*s]))))
}
##------------------------------------------------------------------##
## PLOTANDO AS PREVISÕES DA SÉRIE SELECIONADA ##
##------------------------------------------------------------------##
plot(IPI[(length(serie)+1):length(IPI)], type = "o", xlab = "Ano", ylab = "Quantidade",
las = 1, pch = 21, bg = "yellow", axes = F)
lines(Prev[(length(serie)+1):length(IPI)], type = "l", col = "red")
axis(1, at = seq(1,60,2), labels = c("Ago/95", "Out/95","Dez/95",
"Fev/96", "Abr/96", "Jun/96",
"Ago/96", "Out/96","Dez/96",
"Fev/97", "Abr/97", "Jun/97",
"Ago/97", "Out/97","Dez/97",
"Fev/98", "Abr/98", "Jun/98",
"Ago/98", "Out/98","Dez/98",
"Fev/99", "Abr/99", "Jun/99",
"Ago/99", "Out/99","Dez/99",
"Fev/00", "Abr/00", "Jun/00"), cex.axis = 0.5)
axis(2, seq(80,200,10), las = 1)
box()
legend(1,200, c("IPI Observado", "Previsão Direta"), lty = c(1,1), col = c("black", "red"),
pt.bg = c("yellow", NA), pch = c(21,NA), bty = "n", y.intersp = 0.8)
##------------------------------------------------------------------##
## PLOTANDO PREVISÃO DESATUALIZADA DO MODELO ##
##------------------------------------------------------------------##
plot(IPI, type = "o", pch = 21, bg = "yellow", las = 1,
ylab = "Quantidade", xlab = "Ano", axes = F)
lines(modelo, type = "l", col = "red")
lines(Prev, type = "l", col = "green4")
axis(1, at = seq(0,187, 12), labels = c("1985", "1986", "1987",
"1988", "1989", "1990",
"1991", "1992", "1993",
"1994", "1995", "1996",
"1997", "1998", "1999",
"2000"))
axis(2, seq(60, 200, 10), las = 1)
legend(1,200, c("IPI Observado", "Estimativa", "Previsão"),
col = c("black", "red", "green4"), pch = c(21, NA, NA),
pt.bg = c("yellow", NA, NA), bty = "n", y.intersp = 0.7,
lty = c(1,1,1))
box()
##----------------------------------------------------------------##
## PLOTANDO OS COEFICIENTES DO MODELO ##
##----------------------------------------------------------------##
Prev1 <- numeric(0)
Prev1[length(serie)] = modelo[length(serie)]
for (i in 1:h) {
Niv[length(serie)+i] <- a*(IPI[length(serie)+i]-Saz[length(serie)+i-s])+(1-a)*(Niv[length(serie)+i-1]+Tend[length(serie)+i-1])
Tend[length(serie)+i] <- c*(Niv[length(serie)+i]-Niv[length(serie)+i-1])+(1-c)*Tend[length(serie)+i-1]
Saz[length(serie)+i] <- d*(IPI[length(serie)+i]-Niv[length(serie)+i])+(1-d)*Saz[length(serie)+i-s]
Prev1[length(serie)+i] <- ifelse(i <= s, Niv[length(serie)+i-1]+Tend[length(serie)+i-1]+Saz[length(serie)+i-s],
ifelse(i > s & i <= 2*s, Niv[length(serie)+i-1]+Tend[length(serie)+i-1]+Saz[length(serie)+i-2*s],
ifelse(i > 2*s & i <= 3*s, Niv[length(serie)+i-1]+Tend[length(serie)+i-1]+Saz[length(serie)+i-3*s],
ifelse(i > 3*s & i <= 4*s, Niv[length(serie)+i-1]+Tend[length(serie)+i-1]+Saz[length(serie)+i-4*s],
Niv[length(serie)+i-1]+Tend[length(serie)+i-1]+Saz[length(serie)+i-5*s]))))
}
plot(Niv[(length(serie)+1):length(IPI)], type = "l", col = "blue", xlab = "Nível", ylab = "Estimativa")
plot(Saz[(length(serie)+1):length(IPI)], type = "l", col = "blue", xlab = "Sazonalidade", ylab = "Estimativa")
plot(Tend[(length(serie)+1):length(IPI)], type = "l", col = "blue", xlab = "Tendência", ylab = "Estimativa")
##---------------------------------------------------------------##
## PLOTANDO PREVISÃO ATUALIZADA DO MODELO ##
##---------------------------------------------------------------##
plot(IPI[(length(serie)+1):length(IPI)], type = "o", xlab = "Ano", ylab = "Quantidade",
las = 1, pch = 21, bg = "yellow", axes = F)
lines(Prev1[(length(serie)+1):length(IPI)], type = "l", col = "red")
axis(1, at = seq(1,60,2), labels = c("Ago/95", "Out/95","Dez/95",
"Fev/96", "Abr/96", "Jun/96",
"Ago/96", "Out/96","Dez/96",
"Fev/97", "Abr/97", "Jun/97",
"Ago/97", "Out/97","Dez/97",
"Fev/98", "Abr/98", "Jun/98",
"Ago/98", "Out/98","Dez/98",
"Fev/99", "Abr/99", "Jun/99",
"Ago/99", "Out/99","Dez/99",
"Fev/00", "Abr/00", "Jun/00"), cex.axis = 0.5)
axis(2, seq(80,150,10), las = 1)
box()
legend(1,150, c("IPI Observado", "Previsão Atualizada"), lty = c(1,1), col = c("black", "red"),
pt.bg = c("yellow", NA), pch = c(21,NA), bty = "n", y.intersp = 0.8)
##---------------------------------------------------------------##
## PLOTANDO PREVISÃO ATUALIZADA DO MODELO ##
##---------------------------------------------------------------##
plot(IPI, type = "o", pch = 21, bg = "yellow", las = 1,
ylab = "Quantidade", xlab = "Ano", axes = F)
lines(modelo, type = "l", col = "red")
lines(Prev1, type = "l", col = "green4")
axis(1, at = seq(0,187, 12), labels = c("1985", "1986", "1987",
"1988", "1989", "1990",
"1991", "1992", "1993",
"1994", "1995", "1996",
"1997", "1998", "1999",
"2000"))
axis(2, seq(60, 150, 10), las = 1)
legend(1,150, c("IPI Observado", "Estimativa", "Previsão"),
col = c("black", "red", "green4"), pch = c(21, NA, NA),
pt.bg = c("yellow", NA, NA), bty = "n", y.intersp = 0.7,
lty = c(1,1,1))
box()
##--------------------------------------------------------##
## MEDIDAS DE ACURÁCIA DO MODELO ##
##--------------------------------------------------------##
EQM <- sum((serie[s:length(serie)]-modelo[s:length(serie)])^2)/(length(serie)-s+1)
MAPE <- (sum(abs(serie[s:length(serie)]-modelo[s:length(serie)])/serie[s:length(serie)])*100)/(length(serie)-s+1)
MAD <- sum(abs(serie[s:length(serie)]-modelo[s:length(serie)]))/(length(serie)-s+1)
data.frame(EQM,MAPE,MAD)
## EQM MAPE MAD
## 1 0.1392587 0.2209183 0.2110252
##----------------------------------------------------------##
## MEDIDAS DE ACURÁCIA DE PREVISÁO ##
##----------------------------------------------------------##
EQM1 <- sum((IPI[(length(serie)+1):length(IPI)]-Prev[(length(serie)+1):length(IPI)])^2)/h
MAPE1 <- (sum(abs(IPI[(length(serie)+1):length(IPI)]-Prev[(length(serie)+1):length(IPI)])/IPI[(length(serie)+1):length(IPI)])*100)/h
MAD1 <- sum(abs(IPI[(length(serie)+1):length(IPI)]-Prev[(length(serie)+1):length(IPI)]))/h
EQM2 <- sum((IPI[(length(serie)+1):length(IPI)]-Prev1[(length(serie)+1):length(IPI)])^2)/h
MAPE2 <- (sum(abs(IPI[(length(serie)+1):length(IPI)]-Prev1[(length(serie)+1):length(IPI)])/IPI[(length(serie)+1):length(IPI)])*100)/h
MAD2 <- sum(abs(IPI[(length(serie)+1):length(IPI)]-Prev1[(length(serie)+1):length(IPI)]))/h
Previsão <- c("Desatualizado", "Atualizado")
EQM <- c(EQM1, EQM2)
MAPE <- c(MAPE1, MAPE2)
MAD <- c(MAD1, MAD2)
data.frame(Previsão, EQM, MAPE, MAD)
## Previsão EQM MAPE MAD
## 1 Desatualizado 286.9519 12.978538 14.204537
## 2 Atualizado 105.4109 7.436937 8.553676