library(fpp2)
library(astsa)
library(readxl)
library(tseries)
library(kableExtra)

Se usara un modelo de entrenamiento y un modelo de prueba de h=57 (meses).

mexico<-read_xlsx("Mexico.xlsx")
ts.mexico<-ts(mexico[,c(2)],start = c(1835,1),frequency = 12)

ts.mexico%>%
  window(start=2009)->ts.mexico.test
ts.mexico%>%
  window(end=c(2008,12))->ts.mexico.train
  
ts.mexico%>%
  autoplot(facet=T)

## Warning: Ignoring unknown parameters: facet

decompose(ts.mexico.train)%>%
  autoplot()

Se observa que hay un comportamiento estacional claro, que se repite cada año (por las estaciones del año), asimismo, hay una tendencia en los datos. En los residuos no se observa que sean heterocedasticos.

Ajuste por SARIMA

Ver la cantidad de ajustes de orden 1

ndiffs(ts.mexico.train)

## [1] 1

Ver la cantidad de diferencias estacionales

nsdiffs(ts.mexico.train)

## [1] 1

Se realiza la diferenciacion

ddmexico<-diff(diff(ts.mexico.train),lag = 12)
ddmexico%>%autoplot

Prueba de KPSS para verificar la estacionariedad

kpss<-kpss.test(ddmexico)

## Warning in kpss.test(ddmexico): p-value greater than printed p-value

kpss$p.value

## [1] 0.1

El P Valor es mayor que el nivel de significancia 0.05, por lo que se acepta que la serie diferenciada ya es estacionaria.

Se procede a evaluar los gráficos de autocorrelacion y autocorrelacion parcial

ggAcf(ddmexico,lag.max = 60)

ggPacf(ddmexico,lag.max = 60)

Para la parte estacional, se observa que hay un corte en el ACF en L=12, por lo que se trata de un modelo (0,1,1)[12]. Por otro lado, en la parte no estacional se observa un corte en 1, pero 2 tambien se encuentra cerca del límite, por lo que se usarán (0,1,1), (0,1,2), (1,2,1), (1,2,2) y se elegirá el más adecuado de acuerdo con sus valores AIC

auto.arima(ts.mexico.train)

## Series: ts.mexico.train 
## ARIMA(0,0,5)(1,1,0)[12] with drift 
## 
## Coefficients:
##          ma1     ma2     ma3     ma4     ma5     sar1   drift
##       0.2167  0.1105  0.0427  0.0337  0.0403  -0.4858  0.0009
## s.e.  0.0220  0.0225  0.0223  0.0224  0.0214   0.0192  0.0015
## 
## sigma^2 = 0.6836:  log likelihood = -2548.98
## AIC=5113.96   AICc=5114.03   BIC=5159.07

arima1<-arima(ts.mexico.train, order=c(0,1,1),seasonal = c(0,1,1))
arima2<-arima(ts.mexico.train, order=c(0,1,2),seasonal = c(0,1,1))
arima3<-arima(ts.mexico.train, order=c(1,2,1),seasonal = c(0,1,1))
arima4<-arima(ts.mexico.train, order=c(1,2,2),seasonal = c(0,1,1))
arima5<-arima(ts.mexico.train, order=c(0,0,5),seasonal = c(1,1,0))
aics<-c(arima1$aic,arima2$aic,arima3$aic,arima4$aic,arima5$aic)
which.min(aics)

## [1] 2

El modelo más adecuado es el Sarima (0,1,2)x(0,1,1)[12].

Revision de residuos

sarima.res<-checkresiduals(arima2)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,1,2)(0,1,1)[12]
## Q* = 39.595, df = 21, p-value = 0.008326
## 
## Model df: 3.   Total lags used: 24

sarima.res$p.value

## [1] 0.008326306

El P valor no muestra que el modelo tenga un ajuste adecuado, por los residuos. Sin embargo, se hará una estimación con este método.

Se realiza forecast con el modelo anterior.

sarima.forecast<-forecast(arima2,h=57)

autoplot(
  window(
    ts.mexico.train,start=2000))+
  autolayer(sarima.forecast$mean,
            series="forecast")+
  autolayer(ts.mexico.test,
            series="real")

Se calcula el error de estimacion

SSE_Sarima<-sum((sarima.forecast$mean-ts.mexico.test)^2)
SSE_Sarima

## [1] 28.03144

Ajustando por el metodo NAIVE

naive<-naive(ts.mexico.train,h=57)


autoplot(
  window(
    ts.mexico.train,start=2000))+
  autolayer(naive$mean,
            series="forecast")+
  autolayer(ts.mexico.test,
            series="real")

Se calcula el error de estimacion

SSE_Naive<-sum((naive$mean-ts.mexico.test)^2)
SSE_Naive

## [1] 1073.93

Ajustando por el metodo SNAIVE

snaive<-snaive(ts.mexico.train,h=57)


autoplot(
  window(
    ts.mexico.train,start=2000))+
  autolayer(snaive$mean,
            series="forecast")+
  autolayer(ts.mexico.test,
            series="real")

Se calcula el error de estimacion

SSE_SNaive<-sum((snaive$mean-ts.mexico.test)^2)
SSE_SNaive

## [1] 41.45584

Se observa que este método aumenta sustancialmente el error.

Ajustando por un método naive estacional

ses<-ses(ts.mexico.train,h=57,initial="optimal")
autoplot(
  window(
    ts.mexico.train,start=2000))+
  autolayer(ses$mean,
            series="forecast")+
  autolayer(ts.mexico.test,
            series="real")

Se calcula el error de estimacion

SSE_ses<-sum((ses$mean-ts.mexico.test)^2)
SSE_ses

## [1] 1073.908

Ajustando por un modelo Holt Aditivo

holt.add<-holt(ts.mexico.train,h=57,alpha  = 0.05)
autoplot(
  window(
    ts.mexico.train,start=2000))+
  autolayer(holt.add$mean,
            series="forecast")+
  autolayer(ts.mexico.test,
            series="real")

Calcular el error estandar

SSE_Holt<-sum((holt.add$mean-ts.mexico.test)^2)
SSE_Holt

## [1] 336.0584

Ajustando con un modelo holt winters multiplicativo

hw.mult<-hw(ts.mexico.train,h=57,seasonal = "multiplicative")
autoplot(
  window(
    ts.mexico.train,start=2000))+
  autolayer(hw.mult$mean,
            series="forecast")+
  autolayer(ts.mexico.test,
            series="real")

Se calcula el error de pronostico

SSE_hwmult<-sum((hw.mult$mean-ts.mexico.test)^2)
SSE_hwmult

## [1] 27.33142

Ajustando con un modelo holt winters aditivo

hw.add<-hw(ts.mexico.train,h=57,seasonal = "additive")
autoplot(
  window(
    ts.mexico.train,start=2000))+
  autolayer(hw.add$mean,
            series="forecast")+
  autolayer(ts.mexico.test,
            series="real")

Se calcula el error de pronostico

SSE_hwadd<-sum((hw.add$mean-ts.mexico.test)^2)
SSE_hwadd

## [1] 29.96151

Ajustando con un método ETS

ETS<-forecast(ts.mexico.train, h=57)
autoplot(
  window(
    ts.mexico.train,start=2000))+
  autolayer(ETS$mean,
            series="forecast")+
  autolayer(ts.mexico.test,
            series="real")

Se calcula el error de pronostico

SSE_ETS<-sum((ETS$mean-ts.mexico.test)^2)
SSE_ETS

## [1] 28.96334

Ajustando por un método de regresion dinamico usando fourier sin método determinista

Optimizando K con un ciclo for

 aic<-NULL
for (i in 1:6) {
xreg<-auto.arima(ts.mexico.train,xreg=fourier(ts.mexico.train,K=i),seasonal = F)
aicb<-xreg$aic
aic<-c(aic,aicb) 
}
K<-which.min(aic)
K

## [1] 5

El valor óptimo es K=5

fourier<-fourier(ts.mexico.train,K=K)
modfourier<-auto.arima(ts.mexico.train,xreg=fourier,seasonal = F)
hfourier<-fourier(ts.mexico.train,K=K,h=57)
xregfourier<-forecast(modfourier,xreg = hfourier)

autoplot(
  window(
    ts.mexico.train,start=2000))+
  autolayer(xregfourier$mean,
            series="forecast")+
  autolayer(ts.mexico.test,
            series="real")

Se calcula el error de pronostico

SSE_Fourier<-sum((xregfourier$mean-ts.mexico.test)^2)
SSE_Fourier

## [1] 27.58734

Ajustando con un Modelo de regresion dinamico (variable adcional y parte modelada por fourier)

Usar variables regresoras

ts.uncrt<-ts(mexico[,c(3)],start = c(1835,1),frequency = 12)

ts.uncrt%>%
  window(start=2009)->ts.uncrt.test
ts.uncrt%>%
  window(end=c(2008,12))->ts.uncrt.train

Modelar la regresion dinamica

uncrt.train<-as.data.frame(ts.uncrt.train)$averagetemperatureuncertainty
dinamic<-auto.arima(ts.mexico.train,
           xreg=cbind(
             uncrt.train,
             fourier),
           seasonal = F)

## Se procede a ejecutar el pronostico
uncrt.test<-as.data.frame(ts.uncrt.test)$averagetemperatureuncertainty

dinamicfcast<-forecast(dinamic,xreg=cbind(
  uncrt.test,
  hfourier
))

autoplot(
  window(
    ts.mexico.train,start=2000))+
  autolayer(dinamicfcast$mean,
            series="forecast")+
  autolayer(ts.mexico.test,
            series="real")

Calcular el Error

SSE_Dinamico<-sum((dinamicfcast$mean-ts.mexico.test)^2)
SSE_Dinamico

## [1] 27.77237

Comparativa de Métodos

Metodo	SSE
ETS	28.96334
Fourier	27.58734
Holt	336.05843
Holt Winters Aditivo	29.96151
Holt Wnters Multiplicativo	27.33142
Naive	1073.93015
Naive Estacional	41.45584
SARIMA	28.03144
SS	1073.90822
Dinamico	27.77237

Modelos de Series de Tiempo

Pedro Pablo Morales Ortiz

2022-10-28

Ajuste por SARIMA

Ver la cantidad de ajustes de orden 1

Ver la cantidad de diferencias estacionales

Se realiza la diferenciacion

Prueba de KPSS para verificar la estacionariedad

Se procede a evaluar los gráficos de autocorrelacion y autocorrelacion parcial

Revision de residuos

Se realiza forecast con el modelo anterior.

Se calcula el error de estimacion

Ajustando por el metodo NAIVE

Se calcula el error de estimacion

Ajustando por el metodo SNAIVE

Se calcula el error de estimacion

Ajustando por un método naive estacional

Se calcula el error de estimacion

Ajustando por un modelo Holt Aditivo

Calcular el error estandar

Ajustando con un modelo holt winters multiplicativo

Se calcula el error de pronostico

Ajustando con un modelo holt winters aditivo

Se calcula el error de pronostico

Ajustando con un método ETS

Se calcula el error de pronostico

Ajustando por un método de regresion dinamico usando fourier sin método determinista

Optimizando K con un ciclo for

Se calcula el error de pronostico

Ajustando con un Modelo de regresion dinamico (variable adcional y parte modelada por fourier)

Usar variables regresoras

Modelar la regresion dinamica

Calcular el Error

Comparativa de Métodos