Time Series: Polynomial Regression Model

Para a aplicação do método de regressão polinomial, considerou-se 187 amostras referentes ao Índice de Produto Industrial do Brasil (IPI), no período de Janeiro de 1985 a Julho de 2000. Os dados foram obtidos no repositório de Pedro A. Morettin, no link https://www.ime.usp.br/~pam/IPI.XLS.

################################################################
##                                                            ##
##                MÉTODO DE DECOMPOSIÇÃO                      ##
##          Modelo de Regressão - Dados Não Sazonais          ##
##                                                            ##
################################################################

st <- read.csv2(file.choose())  ; attach(st) #INSERINDO OS DADOS
                                             # rm(st) apaga o vetor st      

################################################################
##                    PLOTANDO A SÉRIE                        ##
################################################################

plot(ts(IPI, start = c(1985,1), frequency = 12), xlab = "Ano", 
     ylab = "Índice de Produto Industrial", las = 1, type = 'o', 
     pch = 21, bg = "yellow", cex.lab = 1.1) ; abline(v = 
                            c(1986:2000), lty = 2, col = "gray")

################################################################
##          SELECIONANDO A SÉRIE PARA MODELAGEM               ##
################################################################

dt <- IPI[1:175]                 # Série para modelagem

plot(ts(dt, start = c(1985,1), frequency = 12), xlab = "Ano", 
     ylab = "Índice de Produto Industrial", las = 1, type = 'o', 
     pch = 21, bg = "yellow", cex.lab = 1.1) ; abline(v = 
                            c(1986:1999), lty = 2, col = "gray")

################################################################
##                      MATRIZ DO MODELO                      ##
################################################################

P = 4
X = matrix(1:length(dt), length(dt), P+1)   # Matriz do modelo
for (i in 1:(P+1)) {
  X[, i] = X[, i]^(i-1)
}

################################################################
##              OBTENDO COEFICIENTES DO MODELO                ##
################################################################
B = solve(t(X)%*%X, tol = 1.23565e-225)%*%t(X)%*%dt
Z = X%*%B

plot(dt, xlab = "Period", ylab = "Value", las = 1, type = 'o', 
     pch = 21, bg = "yellow", cex.lab = 1.1)
lines(Z, col = "red", pch = 18, type = "o") ; abline(v = seq(13,
                                  length(dt),12), lty = 2, col = 
                                    "gray")
legend(x = length(dt)/2, y = min(dt)-21,
       xpd = TRUE,
       xjust = .5,
       yjust = 1,
       horiz = TRUE,
       legend = c("Samples",
                  "Model"),
       lty = c(1, 1),
       col = c("black", "red"),
       pch = c(21, 18), pt.bg = c("yellow", NA),
       bty = "n")

################################################################
##                OBTENDO PREVISÕES DO MODELO                 ##
################################################################
##------------------------------------------------------------##
##          MODELO DE PRIMEIRA ORDEM (MODELO LINEAR)          ##
##------------------------------------------------------------##

t = 12                                     # Número de Previsões
Xp = matrix(1:(length(dt)+t), length(dt)+t, P+1)
for (i in 1:(P+1)) {
  Xp[,i] = Xp[,i]^(i-1)
  Xp[1:length(dt),] = NA
}
Zp = Xp%*%B
plot(dt, xlab = "Period", ylab = "Value", las = 1, type = 'o', 
     pch = 21, bg = "yellow", cex.lab = 1.1, xlim = c(0, length(Zp)))
lines(Z, col = "red", lty = 2, type = "p", pch = 18)
lines(Zp, col = "green3", lty = 2, type = "p", pch = 8)
abline(v = seq(13,length(Zp),12), lty = 2, col = "gray")

legend(x = length(dt)/2, y = min(dt)-21,
       xpd = TRUE,
       xjust = .5,
       yjust = 1,
       horiz = TRUE,
       legend = c("Samples",
                  "Model",
                  "Forecast"),
       lty = c(1, 1, 1),
       col = c("black", "red", "green3"),
       pch = c(21, 18, 8),
       pt.bg = c("yellow", NA, NA),
       bty = "n")

################################################################
##                    MEDIDA DE ACURÁCIA                      ##
################################################################
MAPE = sum(abs((dt-Z)/dt))/length(dt)*100
MAD = sum(abs(dt-Z))/length(dt)
MSD = sum(abs(dt-Z)^2)/length(dt)

Value = c(MAPE,MAD,MSD)  
Modelo <- c("MAPE", "MAD", "MSD")
data.frame(Modelo, Value)
##   Modelo    Value
## 1   MAPE  12.7094
## 2    MAD  12.3885
## 3    MSD 206.5822