La
presente entrada intenta explicar como usar la función auto.arima de R, disponible en el paquete forecast,
esta permite identificar de forma automática el modelo arima “que mejor ajuste”
a las características de las series temporales. Es importante que siempre de manera
preliminar es importante estudiar las características de las series temporales,
datos atípicos, componentes estacionales, tendencias y
otros elementos característicos del análisis exploratorio de series.
La función utiliza el algoritmo Hyndman-Khandakar
La función auto.ariama permite utilizarse en su versión más sencilla (aunque particularmente no recomendada) solo requiere utilizar la serie temporal como argumento de la función. Note que en el siguiente ejemplo se utiliza la función para identificar el modelo arima asociado a la serie AirPassengers, disponible en R:
library(forecast)
auto.arima(AirPassengers)
Series: AirPassengers
ARIMA(2,1,1)(0,1,0)[12]
Coefficients:
ar1 ar2
ma1
0.5960 0.2143
-0.9819
s.e. 0.0888 0.0880 0.0292
sigma^2 estimated as 132.3: log
likelihood=-504.92
AIC=1017.85 AICc=1018.17 BIC=1029.35
Note
que la función anterior identifica un modelo ARIMA(ar=2,i=1,ma=1)(sar=0,si=1,mar=0)[12],
autoregresivo de segundo orden (2), integrado de primer orden (1) utilizando el
test de estacionariedad de Kwiatkowski–Phillips–Schmidt–Shin (KPSS)
Ahora, para que la función testee la parte del componente estacional se modifica el argumento seasonal. Mientras que también podemos modificar el criterio de información que se utiliza para determinar el modelo a mejor ajuste.
auto.arima(AirPassengers,
seasonal=T,
ic ="aic")
Series: AirPassengers
ARIMA(0,1,1)(2,1,0)[12]
Coefficients:
ma1 sar1
sar2
-0.3634 -0.1239
0.1911
s.e. 0.0899 0.0934 0.1036
sigma^2 estimated as 133.5: log
likelihood=-505.59
AIC=1019.18 AICc=1019.5 BIC=1030.68
En tercer lugar, se puede también modificar el criterio trace, que, al poner el valor verdadero, imprime el modelo criterio de información asociado a cada uno de los modelos estimados. Esta función particularmente solo la activo cuando necesito generar un reporte en pdf con todos los resultados del script de R para un superior.
auto.arima(AirPassengers,
seasonal=T,
ic ="aic",
trace=T)
ARIMA(2,1,2)(1,1,1)[12] : Inf
ARIMA(0,1,0)(0,1,0)[12] : 1031.508
ARIMA(1,1,0)(1,1,0)[12] : 1020.393
ARIMA(0,1,1)(0,1,1)[12] : 1021.003
ARIMA(1,1,0)(0,1,0)[12] : 1020.394
ARIMA(1,1,0)(2,1,0)[12] : 1019.24
ARIMA(1,1,0)(2,1,1)[12] : Inf
ARIMA(1,1,0)(1,1,1)[12] : Inf
ARIMA(0,1,0)(2,1,0)[12] : 1032.12
ARIMA(2,1,0)(2,1,0)[12] : 1021.12
ARIMA(1,1,1)(2,1,0)[12] : 1021.033
ARIMA(0,1,1)(2,1,0)[12] : 1019.178
ARIMA(0,1,1)(1,1,0)[12] : 1020.425
ARIMA(0,1,1)(2,1,1)[12] : Inf
ARIMA(0,1,1)(1,1,1)[12] : Inf
ARIMA(0,1,2)(2,1,0)[12] : 1021.148
ARIMA(1,1,2)(2,1,0)[12] : 1022.805
Best model: ARIMA(0,1,1)(2,1,0)[12]
Series: AirPassengers
ARIMA(0,1,1)(2,1,0)[12]
Coefficients:
ma1 sar1
sar2
-0.3634 -0.1239
0.1911
s.e. 0.0899 0.0934 0.1036
sigma^2 estimated as 133.5: log
likelihood=-505.59
AIC=1019.18 AICc=1019.5 BIC=1030.68
Otro argumento que es recomendable modificar es stepwise, al colocar falso permite identificar una mayor cantidad de modelo al no hacer la estimación paso a paso. Puntualmente, esto hace ir testeando modelo a modelo según las ganancias obtenidas en los modelos comparados, mientras que al omitir este paso se evalúan todos los modelos y se hace una selección global. Esto porque el autor inteligentemente ha predefinido ciertos valores de los parámetros del modelo, esto para acelerar la búsqueda de los modelos. Observe que aquí se imprime nuevamente la función trece, por lo que, la función busca ahora mayor cantidad de modelos, identificándose en este caso mayor cantidad de parámetros en los modelos.
auto.arima(AirPassengers,
seasonal=T, ic
="aic",
trace=T,
stepwise = FALSE)
ARIMA(0,1,0)(0,1,0)[12] : 1031.508
ARIMA(0,1,0)(0,1,1)[12] : 1030.752
ARIMA(0,1,0)(0,1,2)[12] : 1032.276
ARIMA(0,1,0)(1,1,0)[12] : 1030.408
ARIMA(0,1,0)(1,1,1)[12] : 1032.128
ARIMA(0,1,0)(1,1,2)[12] : 1034.096
ARIMA(0,1,0)(2,1,0)[12] : 1032.12
ARIMA(0,1,0)(2,1,1)[12] : 1034.106
ARIMA(0,1,0)(2,1,2)[12] : Inf
ARIMA(0,1,1)(0,1,0)[12] : 1020.639
ARIMA(0,1,1)(0,1,1)[12] : 1021.003
ARIMA(0,1,1)(0,1,2)[12] : 1019.495
ARIMA(0,1,1)(1,1,0)[12] : 1020.425
ARIMA(0,1,1)(1,1,1)[12] : Inf
ARIMA(0,1,1)(1,1,2)[12] : Inf
ARIMA(0,1,1)(2,1,0)[12] : 1019.178
ARIMA(0,1,1)(2,1,1)[12] : Inf
ARIMA(0,1,1)(2,1,2)[12] : Inf
ARIMA(0,1,2)(0,1,0)[12] : 1022.627
ARIMA(0,1,2)(0,1,1)[12] : 1023.001
ARIMA(0,1,2)(0,1,2)[12] : 1021.454
ARIMA(0,1,2)(1,1,0)[12] : 1022.425
ARIMA(0,1,2)(1,1,1)[12] : Inf
ARIMA(0,1,2)(1,1,2)[12] : Inf
ARIMA(0,1,2)(2,1,0)[12] : 1021.148
ARIMA(0,1,2)(2,1,1)[12] : Inf
ARIMA(0,1,3)(0,1,0)[12] : 1020.204
ARIMA(0,1,3)(0,1,1)[12] : 1020.806
ARIMA(0,1,3)(0,1,2)[12] : 1020.16
ARIMA(0,1,3)(1,1,0)[12] : 1020.362
ARIMA(0,1,3)(1,1,1)[12] : 1020.618
ARIMA(0,1,3)(2,1,0)[12] : 1019.745
ARIMA(0,1,4)(0,1,0)[12] : 1020.696
ARIMA(0,1,4)(0,1,1)[12] : 1020.177
ARIMA(0,1,4)(1,1,0)[12] : 1019.607
ARIMA(0,1,5)(0,1,0)[12] : 1021.399
ARIMA(1,1,0)(0,1,0)[12] : 1020.394
ARIMA(1,1,0)(0,1,1)[12] : 1020.914
ARIMA(1,1,0)(0,1,2)[12] : 1019.493
ARIMA(1,1,0)(1,1,0)[12] : 1020.393
ARIMA(1,1,0)(1,1,1)[12] : Inf
ARIMA(1,1,0)(1,1,2)[12] : Inf
ARIMA(1,1,0)(2,1,0)[12] : 1019.24
ARIMA(1,1,0)(2,1,1)[12] : Inf
ARIMA(1,1,0)(2,1,2)[12] : Inf
ARIMA(1,1,1)(0,1,0)[12] : 1022.394
ARIMA(1,1,1)(0,1,1)[12] : 1022.897
ARIMA(1,1,1)(0,1,2)[12] : 1021.313
ARIMA(1,1,1)(1,1,0)[12] : 1022.356
ARIMA(1,1,1)(1,1,1)[12] : Inf
ARIMA(1,1,1)(1,1,2)[12] : Inf
ARIMA(1,1,1)(2,1,0)[12] : 1021.033
ARIMA(1,1,1)(2,1,1)[12] : Inf
ARIMA(1,1,2)(0,1,0)[12] : 1024.161
ARIMA(1,1,2)(0,1,1)[12] : 1024.718
ARIMA(1,1,2)(0,1,2)[12] : 1023.125
ARIMA(1,1,2)(1,1,0)[12] : Inf
ARIMA(1,1,2)(1,1,1)[12] : Inf
ARIMA(1,1,2)(2,1,0)[12] : 1022.805
ARIMA(1,1,3)(0,1,0)[12] : 1019.253
ARIMA(1,1,3)(0,1,1)[12] : 1019.533
ARIMA(1,1,3)(1,1,0)[12] : 1019.134
ARIMA(1,1,4)(0,1,0)[12] : Inf
ARIMA(2,1,0)(0,1,0)[12] : 1022.394
ARIMA(2,1,0)(0,1,1)[12] : 1022.905
ARIMA(2,1,0)(0,1,2)[12] : 1021.381
ARIMA(2,1,0)(1,1,0)[12] : 1022.373
ARIMA(2,1,0)(1,1,1)[12] : Inf
ARIMA(2,1,0)(1,1,2)[12] : Inf
ARIMA(2,1,0)(2,1,0)[12] : 1021.12
ARIMA(2,1,0)(2,1,1)[12] : Inf
ARIMA(2,1,1)(0,1,0)[12] : 1017.848
ARIMA(2,1,1)(0,1,1)[12] : 1018.36
ARIMA(2,1,1)(0,1,2)[12] : 1017.951
ARIMA(2,1,1)(1,1,0)[12] : 1017.915
ARIMA(2,1,1)(1,1,1)[12] : Inf
ARIMA(2,1,1)(2,1,0)[12] : 1017.658
ARIMA(2,1,2)(0,1,0)[12] : 1019.291
ARIMA(2,1,2)(0,1,1)[12] : 1019.936
ARIMA(2,1,2)(1,1,0)[12] : 1019.547
ARIMA(2,1,3)(0,1,0)[12] : 1019.797
ARIMA(3,1,0)(0,1,0)[12] : 1023.666
ARIMA(3,1,0)(0,1,1)[12] : 1024.441
ARIMA(3,1,0)(0,1,2)[12] : 1023.149
ARIMA(3,1,0)(1,1,0)[12] : 1024.004
ARIMA(3,1,0)(1,1,1)[12] : Inf
ARIMA(3,1,0)(2,1,0)[12] : 1022.819
ARIMA(3,1,1)(0,1,0)[12] : 1019.085
ARIMA(3,1,1)(0,1,1)[12] : 1019.697
ARIMA(3,1,1)(1,1,0)[12] : 1019.328
ARIMA(3,1,2)(0,1,0)[12] : Inf
ARIMA(4,1,0)(0,1,0)[12] : 1021.976
ARIMA(4,1,0)(0,1,1)[12] : 1022.158
ARIMA(4,1,0)(1,1,0)[12] : 1021.53
ARIMA(4,1,1)(0,1,0)[12] : 1023.947
ARIMA(5,1,0)(0,1,0)[12] : 1023.944
Best model: ARIMA(2,1,1)(2,1,0)[12]
Series: AirPassengers
ARIMA(2,1,1)(2,1,0)[12]
Coefficients:
ar1 ar2
ma1 sar1 sar2
0.5741 0.2547
-0.9822 -0.1136 0.1641
s.e. 0.0899 0.0933 0.0286 0.0963 0.1073
sigma^2 estimated as 129.5: log
likelihood=-502.83
AIC=1017.66 AICc=1018.34 BIC=1034.91
Otra recomendación
al momento de utilizar la función es activar el argumento parallel.
Esta permite acelerar el ritmo de búsqueda del modelo, dado que utiliza los diversos núcleos de los modelos, pero tenga pendiente es que la función trace debe apagarse si se desea activar esta función. Además, al hacer un summary sobre el modelo arima podemos verificar una evaluación del modelo dentro del set de datos. Tenga pendiente que estos no son buenos indicadores sobre la capacidad de pronóstico del modelo fuera de muestra o pronósticos.
modeloAR <- auto.arima(AirPassengers,
seasonal=T, ic
="aic",
trace=F,
stepwise = FALSE,
parallel=TRUE)
summary(modeloAR)
Series: AirPassengers
ARIMA(2,1,1)(2,1,0)[12]
Coefficients:
ar1 ar2
ma1 sar1 sar2
0.5741 0.2547
-0.9822 -0.1136 0.1641
s.e. 0.0899 0.0933 0.0286 0.0963 0.1073
sigma^2 estimated as 129.5: log
likelihood=-502.83
AIC=1017.66 AICc=1018.34 BIC=1034.91
Training set error measures:
ME RMSE
MAE MPE MAPE
MASE ACF1
Training set 1.30706 10.64393 7.69726 0.4019023 2.734687 0.2403118 -0.007573175
Finalmente, se puede utilizar el modelo identificado para proyectar los valores de la serie a un determinado horizonte (12 en este caso) fuera de la muestra. Esta función utiliza el modelo, el horizonte a proyección y un nivel de confianza.
pronos <- forecast(modeloAR,h=12,level=95)
autoplot(pronos) + theme_minimal()
auto.arima(AirPassengers,
seasonal=T, ic
="aic",
trace=F,
stepwise = FALSE,
parallel=TRUE,
approximation=FALSE,
allowdrift=TRUE,
allowmean=TRUE,
lambda="auto")
Referencias
Hyndman, R., & Khandakar, Y. (2008). Automatic time series forecasting: The forecast package for R. Journal of Statistical Software. Recuperado el 21 de 12 de 2020, de https://www.jstatsoft.org/article/view/v027i03
Hyndman, R., Athanasopoulos, G., Bergmeir, C., Caceres, G., Chhay, L., O'Hara, M., . . . Wang, E. (2020). Package ‘forecast’. CRAN.
Kwiatkowski, D., Phillips, P., Schmidt, P., Shin,
& Y. (1992). Testing the null hypothesis of stationarity against the
alternative of a unit root. North-Holland: Journal of Econometrics.
Meca, I. y Agulló, O. (2018). Modelos Arima. Disponible en: https://rpubs.com/.
