library(dplyr)
library(ggplot2)
library(kableExtra)
library(fpp2)
library(tseries)
library(astsa)
df <- read.csv("https://query.data.world/s/rm3p3z2jzpiccwuu23w6cokypzrhyl", header=TRUE, stringsAsFactors=FALSE)
p<-df$Monthly.car.sales.in.Quebec.1960.1968
cars<-ts(p,start=c(1960,1),frequency = 12)
cars%>%
  autoplot()

cars.test<-window(cars,start=c(1968,1))
cars.train<-window(cars,end=c(1967,12))

Decomponer

decompose(cars.train)%>%
  autoplot()

Sarima

Hay que evaluar los grficos ACF PACF

ndiffs(cars.train)
## [1] 1
nsdiffs(cars.train)
## [1] 1
ddcars<-diff(diff(cars.train,lag=12))
ddcars%>%
  autoplot()

ggAcf(ddcars, lag.max = 60)

ggPacf(ddcars,lag.max = 60)

Estacional (0,1,1) No estacional(1,1,1)(1,1,0)(2,1,1)(3,1,1)

auto.arima(cars)
## Series: cars 
## ARIMA(1,0,1)(0,1,1)[12] with drift 
## 
## Coefficients:
##          ar1      ma1     sma1    drift
##       0.7663  -0.5394  -0.5619  83.3957
## s.e.  0.1685   0.2208   0.1337  12.8709
## 
## sigma^2 = 2328583:  log likelihood = -840.23
## AIC=1690.45   AICc=1691.12   BIC=1703.27
arima1<-Arima(cars.train,order=c(1,1,1),seasonal=c(0,1,1))
arima2<-Arima(cars.train,order=c(2,1,1),seasonal=c(0,1,1))
arima3<-Arima(cars.train,order=c(3,1,1),seasonal=c(0,1,1))
arima4<-Arima(cars.train,order=c(1,0,1),seasonal=c(0,1,1))
arima5<-Arima(cars.train,order=c(1,1,0),seasonal=c(0,1,1))

which.min(c(arima1$aic,arima2$aic,arima3$aic,arima4$aic,arima5$aic))
## [1] 1
forecast(arima1,h=12)->fcast.sarima

autoplot(cars.train)+
  autolayer(fcast.sarima$mean,
            series="Forecast")+
  autolayer(cars.test,
            series="Real")

Calculo del error de pronostico

SSE_SARIMA<-sum((fcast.sarima$mean-cars.test)^2)
SSE_SARIMA
## [1] 42198548

Model Holt Winters Aditivo

holt.winters.add<-hw(cars.train,h=12,seasonal = "additive")
autoplot(cars.train)+
  autolayer(holt.winters.add$mean,
            series="Forecast")+
  autolayer(cars.test,
            series="Real")

Calcular Error

SSE_HW.ADD<-sum((holt.winters.add$mean-cars.test)^2)
SSE_HW.ADD
## [1] 34637591

Model Holt Winters Multiplicativo

holt.winters.mult<-hw(cars.train,h=12,seasonal = "multiplicative")
autoplot(cars.train)+
  autolayer(holt.winters.mult$mean,
            series="Forecast")+
  autolayer(cars.test,
            series="Real")

Calcular Error

SSE_HW.MULT<-sum((holt.winters.mult$mean-cars.test)^2)
SSE_HW.MULT
## [1] 42672484

Naive Estacional

snaive<-snaive(cars.train,h=12)
autoplot(cars.train)+
  autolayer(snaive$mean,
            series="Forecast")+
  autolayer(cars.test,
            series="Real")

## Calcular Error

SSE_NAIVE<-sum((snaive$mean-cars.test)^2)
SSE_NAIVE
## [1] 62974674

Modelo de Regresión Dinámica

xreg.simple<-auto.arima(cars.train,
           xreg = 1:length(cars.train))
xreg.simple
## Series: cars.train 
## Regression with ARIMA(2,0,1)(1,0,0)[12] errors 
## 
## Coefficients:
##           ar1     ar2     ma1    sar1  intercept     xreg
##       -0.2804  0.4025  0.5867  0.8701   9998.568  74.1317
## s.e.   0.1625  0.0941  0.1568  0.0410   1701.219  21.4823
## 
## sigma^2 = 2576450:  log likelihood = -850.39
## AIC=1714.77   AICc=1716.05   BIC=1732.72
fcast.xreg.simple<-forecast(xreg.simple,
         xreg=length(cars.train)+1:12)

autoplot(cars.train)+
  autolayer(fcast.xreg.simple$mean,
            series="Forecast")+
  autolayer(cars.test,
            series="Real")

Calcular Error

SSE_xreg.simple<-sum((fcast.xreg.simple$mean-cars.test)^2)
SSE_xreg.simple
## [1] 34984270

Dinamico de Fourier

fourier<-fourier(cars.train,K=1)
xreg.fourier<-auto.arima(cars.train,
           xreg = cbind(1:length(cars.train),
                        fourier),
           seasonal = F)

Optimizacion de la serie de fourier

  aic<-NULL
for (i in 1:6) {
xreg.fourier<-auto.arima(cars.train,
           xreg = cbind(1:length(cars.train),
                        fourier(cars.train,K=i)),
           seasonal = F)
base<-xreg.fourier$aic
aic<-c(aic,base)    
}
aic
## [1] 1746.028 1690.822 1683.745 1687.327 1674.276 1676.259
which.min(aic)
## [1] 5
fourier<-fourier(cars.train,K=5)
xreg.fourier<-auto.arima(cars.train,
           xreg = cbind(1:length(cars.train),
                        fourier),
           seasonal = F)

## Forecast
hfourier<-fourier(cars.train,K=5,h=12)

fcast.xreg.fourier<-forecast(xreg.fourier,
         xreg=cbind(length(cars.train)+1:12,
                    hfourier))

autoplot(cars.train)+
  autolayer(fcast.xreg.fourier$mean,
            series="Forecast")+
  autolayer(cars.test,
            series="Real")

## Revisar el ajuste del modelo

xreg.fourier$fitted%>%
  autoplot()+
  autolayer(cars.train)

Calcular Error

SSE_xreg.fourier<-sum((fcast.xreg.fourier$mean-cars.test)^2)
SSE_xreg.fourier
## [1] 31658239

Crear el modelo final

fourier<-fourier(cars,K=5)
modelo.final<-auto.arima(cars,
           xreg = cbind(1:length(cars),
                        fourier),
           seasonal = F)

## Forecast
hfourier<-fourier(cars,K=5,h=24)

fcast.xreg.fourier<-forecast(modelo.final,
         xreg=cbind(length(cars.train)+1:24,
                    hfourier))

fcast.xreg.fourier%>%
  autoplot()