2.1 Modelos ARIMA

# install.packages(c("readxl", "forecast", "tseries", "ggplot2"))

library(readxl)
library(forecast)

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

library(tseries)
library(ggplot2)

data <- read_excel("Water2.xlsx", sheet = "Hoja1")
serie <- ts(data$Data, frequency = 12)  # mensual

autoplot(serie) + 
  ggtitle("Serie de Tiempo - Datos Mensuales") +
  xlab("Tiempo (meses)") + ylab("Valor")

adf.test(serie)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  serie
## Dickey-Fuller = -2.6275, Lag order = 4, p-value = 0.3166
## alternative hypothesis: stationary

serie_diff <- diff(serie)

adf.test(serie_diff)

## Warning in adf.test(serie_diff): p-value smaller than printed p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  serie_diff
## Dickey-Fuller = -5.7771, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationary

par(mfrow = c(1,2))
acf(serie_diff, main="ACF - Serie diferenciada")
pacf(serie_diff, main="PACF - Serie diferenciada")

par(mfrow = c(1,1))

modelo_arima <- Arima(serie, order = c(1, 1, 1))
summary(modelo_arima)

## Series: serie 
## ARIMA(1,1,1) 
## 
## Coefficients:
##          ar1      ma1
##       0.8701  -1.0000
## s.e.  0.0569   0.1454
## 
## sigma^2 = 4.762:  log likelihood = -219.96
## AIC=445.91   AICc=446.16   BIC=453.73
## 
## Training set error measures:
##                      ME     RMSE      MAE      MPE     MAPE      MASE
## Training set -0.1990862 2.149632 1.748978 8.264955 240.9662 0.3765516
##                    ACF1
## Training set 0.06364312

forecast_arima <- forecast(modelo_arima, h = 4)
autoplot(forecast_arima) + ggtitle("Pronóstico ARIMA(1,1,1)")

modelo_sarima <- Arima(serie, order = c(1, 1, 1), seasonal = c(1, 1, 1))
summary(modelo_sarima)

## Series: serie 
## ARIMA(1,1,1)(1,1,1)[12] 
## 
## Coefficients:
##          ar1      ma1     sar1     sma1
##       0.2419  -0.2161  -0.2922  -0.2124
## s.e.  1.1469   1.1491   0.3452   0.4052
## 
## sigma^2 = 6.209:  log likelihood = -204.7
## AIC=419.39   AICc=420.12   BIC=431.78
## 
## Training set error measures:
##                      ME     RMSE      MAE      MPE     MAPE      MASE
## Training set 0.05278957 2.272356 1.698842 -60.3544 300.1729 0.3657574
##                     ACF1
## Training set -0.01185301

forecast_sarima <- forecast(modelo_sarima, h = 4)
autoplot(forecast_sarima) + ggtitle("Pronóstico SARIMA(1,1,1)(1,1,1)[12]")

modelo_alt1 <- Arima(serie, order = c(2, 1, 2)) 
modelo_alt2 <- Arima(serie, order = c(0, 1, 1), seasonal = c(1, 1, 0)) 

summary(modelo_alt1)

## Series: serie 
## ARIMA(2,1,2) 
## 
## Coefficients:
##           ar1      ar2     ma1     ma2
##       -1.2148  -0.6230  1.2736  0.7745
## s.e.   0.2447   0.2623  0.2000  0.2163
## 
## sigma^2 = 4.925:  log likelihood = -219.69
## AIC=449.39   AICc=450.02   BIC=462.41
## 
## Training set error measures:
##                     ME     RMSE      MAE      MPE     MAPE      MASE
## Training set -0.104673 2.163663 1.791323 27.51791 237.6793 0.3856684
##                     ACF1
## Training set -0.01393519

summary(modelo_alt2)

## Series: serie 
## ARIMA(0,1,1)(1,1,0)[12] 
## 
## Coefficients:
##          ma1     sar1
##       0.0196  -0.4580
## s.e.  0.1012   0.0984
## 
## sigma^2 = 6.107:  log likelihood = -204.88
## AIC=415.76   AICc=416.05   BIC=423.19
## 
## Training set error measures:
##                      ME     RMSE      MAE       MPE     MAPE      MASE
## Training set 0.05480725 2.280311 1.701955 -65.77481 292.9239 0.3664276
##                      ACF1
## Training set -0.003247394