data("h02")
plot(h02)

autoplot(h02) + ylab("Retail index") + xlab("Year")

plot(decompose(h02))

adf.test(h02) # p-value<.05 para que sea estacionaria
## Warning in adf.test(h02): p-value smaller than printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  h02
## Dickey-Fuller = -9.5147, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary

Se observa un P valor menor a 0.05 por lo que la serie ya es estacionaria y no necesita que se le aplique ninguna diferencia.

h02 %>% ggtsdisplay()

Posteriormente ejecutamos las tres opciones de los posibles modelos, se sabe por el comportamiento de las gráficas anteriores que es necesario usar un modelo AR.

modelo1 = arima(h02, order=c(4,0,0), seasonal = c(0,0,2))
modelo2 = arima(h02, order=c(5,0,0), seasonal = c(0,0,2))
modelo3 = arima(h02, order=c(6,0,0), seasonal = c(0,0,2))
modelo4 = arima(h02, order=c(7,0,0), seasonal = c(0,0,2))

modelo1$aic
## [1] -405.3088
modelo2$aic
## [1] -404.5826
modelo3$aic
## [1] -402.6676
modelo4$aic
## [1] -403.8973
modelo1
## 
## Call:
## arima(x = h02, order = c(4, 0, 0), seasonal = c(0, 0, 2))
## 
## Coefficients:
##          ar1     ar2     ar3      ar4    sma1    sma2  intercept
##       0.5382  0.2199  0.0727  -0.1362  1.0034  0.5145     0.7575
## s.e.  0.0724  0.0803  0.0813   0.0721  0.0782  0.0615     0.0454
## 
## sigma^2 estimated as 0.006896:  log likelihood = 210.65,  aic = -405.31

Se observa que el mejor modelo a utilizar es el 1. Posteriormente ejecutamos los residuales de la serie de tiempo:

h02 %>%
  arima(order=c(4,0,0), seasonal=c(0,0,2)) %>%
  residuals() %>% ggtsdisplay()

checkresiduals(h02)
## Warning in modeldf.default(object): Could not find appropriate degrees of
## freedom for this model.

forecast(modelo1, h=8)
##          Point Forecast     Lo 80     Hi 80     Lo 95    Hi 95
## Jul 2008      0.9387970 0.8323716 1.0452224 0.7760334 1.101561
## Aug 2008      0.9424444 0.8215827 1.0633060 0.7576024 1.127286
## Sep 2008      1.0086171 0.8761455 1.1410887 0.8060193 1.211215
## Oct 2008      0.9310920 0.7896672 1.0725168 0.7148014 1.147383
## Nov 2008      1.0466699 0.9024509 1.1908889 0.8261060 1.267234
## Dec 2008      1.0349198 0.8889965 1.1808431 0.8117494 1.258090
## Jan 2009      1.0223869 0.8757552 1.1690186 0.7981331 1.246641
## Feb 2009      0.8052106 0.6583636 0.9520576 0.5806275 1.029794
autoplot(forecast(modelo1, h=8))