Se usará la base datos llamada H02
data("h02")
plot(h02)
autoplot(h02) + ylab("Retail index") + xlab("Year")
#Nuestra gráfica muestra que si es estacionaria
plot(decompose(h02))
adf.test(h02)
## 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
# El valor del P value es de 0.01, y debe ser <.05, por lo que si es estacionaria
#Como la gráfica mkuestar una tendencia ascendente, se aplicara la diferenciación
h02 %>% diff(lag=12) %>% ggtsdisplay()
#Vamos a probar con el modelo arima
modelo1 = arima(h02, order=c(3,1,3), seasonal = c(3,1,2))
modelo2 = arima(h02, order=c(4,1,4), seasonal = c(3,1,2))
modelo3 = arima(h02, order=c(4,1,3), seasonal = c(4,1,2))
modelo4 = arima(h02, order=c(4,1,5), seasonal = c(4,1,2))
modelo1$aic
## [1] -562.7761
modelo2$aic
## [1] -559.7355
modelo3$aic
## [1] -563.6873
modelo4$aic
## [1] -559.8206
#De todos nuestros modelos , el mejor por ser mas bajo es -559.8206. El modelo 4
modelo1 = arima(h02, order=c(3,1,3), seasonal = c(3,1,2))
modelo2 = arima(h02, order=c(4,1,4), seasonal = c(3,1,2))
modelo3 = arima(h02, order=c(4,1,3), seasonal = c(4,1,2))
modelo4 = arima(h02, order=c(4,1,5), seasonal = c(4,1,2))
modelo1$aic
## [1] -562.7761
modelo2$aic
## [1] -559.7355
modelo3$aic
## [1] -563.6873
modelo4$aic
## [1] -559.8206
#De todos nuestros modelos , el mejor por ser mas bajo es -559.8206. El modelo 4
modelo4
##
## Call:
## arima(x = h02, order = c(4, 1, 5), seasonal = c(4, 1, 2))
##
## Coefficients:
## Warning in sqrt(diag(x$var.coef)): NaNs produced
## ar1 ar2 ar3 ar4 ma1 ma2 ma3 ma4 ma5
## 0.1307 0.3885 0.3538 -0.1812 -0.885 -0.069 -0.1088 0.053 0.0617
## s.e. 0.0034 NaN 0.0100 NaN NaN NaN 0.0433 NaN NaN
## sar1 sar2 sar3 sar4 sma1 sma2
## 0.3574 -0.8284 -0.1875 -0.3274 -0.8867 0.8277
## s.e. 0.0386 NaN 0.0126 0.0040 0.0678 0.0262
##
## sigma^2 estimated as 0.002425: log likelihood = 295.91, aic = -559.82
euretail %>%
arima(order=c(4,1,5), seasonal=c(4,1,2)) %>%
residuals() %>% ggtsdisplay()
checkresiduals(modelo4)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(4,1,5)(4,1,2)[12]
## Q* = 18.884, df = 9, p-value = 0.02619
##
## Model df: 15. Total lags used: 24
forecast(modelo4, h=8)
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jul 2008 1.0123085 0.9489702 1.0756467 0.9154409 1.1091760
## Aug 2008 0.9866853 0.9214662 1.0519044 0.8869412 1.0864293
## Sep 2008 1.1295408 1.0579455 1.2011361 1.0200452 1.2390364
## Oct 2008 1.1458582 1.0688670 1.2228493 1.0281104 1.2636059
## Nov 2008 1.1229446 1.0454764 1.2004128 1.0044672 1.2414220
## Dec 2008 1.2280962 1.1472736 1.3089188 1.1044887 1.3517037
## Jan 2009 1.2133509 1.1312732 1.2954285 1.0878239 1.3388778
## Feb 2009 0.6950121 0.6120592 0.7779651 0.5681466 0.8218777
autoplot(forecast(modelo4, h=8))