1 jul 2021

Usar la función auto.arima en R

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 (2008), desarrollado en R en el paquete “forecast(Hyndman, et al., 2020)

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) (Kwiatkowski, Phillips, Schmidt, Shin, & Y., 1992), media móvil de primer orden (1), luego, referido a la parte estacional no se identifica efecto sar o sma, mientras que establece que se requiere una diferencia estacional de primer orden. Finalmente, el 12 se refiere a la frecuencia de la serie temporal.

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()


Ojo: finalmente, siempre es recomendable realizar un serio análisis exploratorio para identificar la serie y evaluar las propiedades del residuo del modelo arima estimado. En el siguiente enlace se pueden identificar el resto de argumentos de la función: https://www.rdocumentation.org/packages/forecast/versions/8.15/topics/auto.arima. Como el número máximo de modelos que se pueden identificar. Xreg para agregar regresores exógenos, permitir deriva o allowdrift. Una opción recomienda dentro de los argumentos de la función es (otros detalles pueden verse en https://repositorio.uniandes.edu.co/bitstream/handle/1992/45515/u826935.pdf?sequence=1):

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/.

 

 

 

 

 

Recesión plot en R usando ggplot (recession plot in r)

En el siguiente ejemplo se simular y replica -parcialmente- un ejemplo usado por el FMI para ilustrar la importancia del uso de series de ti...