Herramientas para Pronósticos de Series de Tiempo en Python - Parte 1: SARIMAX

Un clásico arranca un duelo entre diferentes herramientas.

Omar Humberto González Ávila

Última actualización el feb. 28, 2022 20 min de lectura Pronósticos

Fuente: Rustyness en Unsplash

1. Introducción

Ah, el futuro… la última frontera. Cómo economista, hacer estimaciones de lo que pasará en el futuro es pan de cada día, tanto así que hay muchos chistes al respecto. Uno de ellos cuenta, “le pedí su número de teléfono a un economista … y me dio una estimación”. El chiste reside en que el experto en números ni siquiera tiene un número preciso de su propio número telefónico; sin embargo, para mí expresa una situación mucho más simple, y es que por mucho que cualquiera analice los datos, nadie sabe lo que pasará. NADIE¹.

¿Entonces, por qué hacer pronósticos? Porque a través de diferentes herramientas se pueden hacer aproximaciones de lo que pasará en el futuro. En ese sentido en mi entrega anterior probé diferentes herramientas para hacer pronósticos de series de tiempo: a la vieja escuela, nueva escuela y la herramienta Prophet. Cabe señalar que mi objetivo con dicho análisis era comparar los resultados entre diferentes metodologías con base en el pronóstico en sí y no hacer la evaluación con las típicas medidas intra muestra. Sin embargo, no estoy del todo conforme con esos resultados.

En ese sentido, las primeras 2 herramientas que utilicé, en el fondo usan la metodología ARIMA que aprendí de mis profesores y solo Prophet sale de este marco teórico. Por tanto, decidí repetir dicho ejercicio con más y diferentes herramientas, con lo cual, de paso practicaré su implementación en Python y Jupyter.

Arrancaré este segundo round con el modelo SARIMAX, cuya notebook completa puede consultarse en mi repositorio de Github o la pueden ver directamente aquí. Debido a que esta parte quedó muy extensa, me saltaré buena parte del código que esta inspirado en esta publicación de Medium, así que ¡a darle átomos!

2. Recordando el problema

Para esta nueva serie, vuelvo a analizar la Recaudación Federal Participable (RFP), siguiendo la referencia de un artículo desarrollado por mis profesores donde se estiman la elasticidad ingreso de los impuestos federales en México y sus efectos en la misma. La RFP está conformada por distintos ingresos del Gobierno Federal, y es de suma importancia porque a partir ella se determinan las participaciones federales, mismas que son la principal fuente de ingresos de libre disposición que tienen las Entidades Federativas. De manera breve, sólo señalaré que la RFP depende de diversos ingresos del Gobierno Federal, entre los que destacan: Impuesto Sobre la Renta, ISR, que es un impuesto a los ingresos; Impuesto al Valor Agregado, IVA, que es un impuesto al consumo; el Impuesto Especial sobre Producción y Servicios, IEPS, que a su vez se compone de diversos impuestos específicos para disminuir el consumo de ciertos bienes, como son el alcohol, el tabaco y las gasolinas; y finalmente una parte de los ingresos petroleros que recibe el Gobierno Federal vía el Fondo Mexicano del Petróleo para la Estabilización y el Desarrollo.

Al igual que en mi proyecto anterior, los datos para la RFP los tomé de las Estadísticas Oportunas de Finanzas Públicas de donde se pueden consultar las series que componen dicha RFP a partir de 2002. Cabe señalar que existe un desfase de 2 meses entre dichos impuestos (según la RFP) con el período que lo causa, puesto que el ISR de marzo, corresponde a la recaudación de febrero, que a su vez, corresponde con la actividad del mes de enero. Debido a que son muchas las series, solo utilizaré las series más importantes: ISR, IVA, IEPS federal a las gasolinas y diésel, IEPS a tabacos, IEPS a bebidas alcohólicas, IEPS a cervezas, IEPS a las bebidas saborizadas, IEPS a los alimentos con alto contenido calórico, impuesto a la importación e ingresos petroleros. Dado que hay un retraso entre la publicación de las variables y el periodo al que corresponden, no tengo información completa sobre el segundo trimestre de 2021.

En cuanto a las variables explicativas, usaré el Producto Interno Bruto (PIB) de México que publica el Instituto Nacional de Estadística y Geografía (INEGI); y de esa misma fuente tomaré el Índice de Precios Implícito (IPI) para convertir las series de ISR e IVA a términos reales (ver siguiente sección). Para el PIB de periodos futuros utilicé las proyecciones del Banco de México de su informe trimestral; de hecho, ya hay una nueva actualización de esta información pero al momento que trabajé en estas proyecciones aún no estaba disponible.

3. Análisis previo de cointegración y quiebre estructural

Este análisis lo hice para mi primer proyecto, pero creo que vale la pena repetirlo, puesto que a través de éste, se confirma que es correcto utilizar el PIB como variable explicativa del ISR y del IVA, pero no del resto de las variables. Adelantando los resultados que mostraré serie por serie, es necesario verificar la existencia de raíz unitaria, porque de lo contrario, podría caer en regresiones espurias (relaciones estadísticas entre variables que no tienen relación alguna).

Para hacer este análisis, apliqué la prueba de contraste aumentado de Dickey Fuller (ADF), con lo que se concluye que todas las series, con excepción del IEPS federal a las gasolinas y diésel y del IEPS a los alimentos, tienen raíz unitaria, pero esta se elimina al tomar las series en primera diferencia (* denota un valor p menor al 0.10, ** denota un valor p menor al 0.05 y *** denota un valor p menor al 0.01):

Serie	Series en niveles		Series en 1a diferencia
	Estadístico ADF-GLS	Decisión	Estadístico ADF	Decisión
PIB Real	-1.14535	Raíz Unitaria	-2.5909***	Estacionaria
ISR	-0.821204	Raíz Unitaria	-5.457999***	Estacionaria
IVA	-1.109062	Raíz Unitaria	-7.997366***	Estacionaria
IEPS Gasolinas	-2.891800*	Estacionaria
IEPS Tabacos	-0.809903	Raíz Unitaria	-6.605603***	Estacionaria
IEPS Bebidas Alcohólicas	-0.640294	Raíz Unitaria	-6.181862***	Estacionaria
IEPS Cervezas	-0.236237	Raíz Unitaria	-6.233054***	Estacionaria
IEPS Bebidas Saborizadas	-0.195569	Raíz Unitaria	-6.610010***	Estacionaria
IEPS Alimentos	-4.737662***	Estacionaria
Impuesto a la Importación	-0.214093	Raíz Unitaria	-5.421195***	Estacionaria
Ingresos Petroleros	-2.313334	Raíz Unitaria	-6.450070***	Estacionaria

3.1 Cointegración de las Series

Una vez que se tienen los resultados de la prueba ADF, se puede verificar si las series estacionarias en efecto están relacionadas y es correcto utilizar a una para explicar a la otra. Para ello, utilizaré la prueba de cointegración de Engle-Granger, donde se verifica primero que las dos series a utilizar tengan raíz unitaria (paso que ya se realizó), luego se verifica si al hacer la regresión por Mínimos Cuadrados Ordinarios (MCO) entre ellas, su residual es estacionario. Los resultados se muestran a continuación:

Variable Dependiente	Variable Explicativa	ADF de res. MCO	Decisión
ISR	PIB Real	-3.6334**	Estacionaria
IVA	PIB Real	-2.77879*	Estacionaria
IEPS Tabacos	PIB Real	-4.66036***	Estacionaria
IEPS Bebidas Alcohólicas	PIB Real	-3.63147**	Estacionaria
IEPS Cervezas	PIB Real	-2.50686	Raíz Unitaria
IEPS Bebidas Saborizadas	PIB Real	-2.78059	Raíz Unitaria
Impuesto a la importación	PIB Real	-1.4614	Raíz Unitaria
Ingresos Petroleros	PIB Real	-1.30229	Raíz Unitaria

Por lo tanto, y dado que los residuales del MCO en el caso del ISR y del IVA con el PIB son estacionarios, concluyó que es correcto utilizar el PIB como variable para explicar el ISR y el IVA, y podría utilizarlo para los IEPS de tabacos y a las bebidas alcohólicas. Sin embargo, solo lo utilizaré para el ISR y el IVA, como es común en la literatura, y dejaré fuera al resto de las series para hacer su análisis mensual.

Lo hago así, porque al agregar el PIB como variable explicativa, necesito agregar el supuesto de cuál será el PIB a futuro y cuál será el IPI. Además evitó perder información por el rezago en los datos, y de hecho, en las series mensuales tengo información hasta el mes de julio. En cambio, al hacer el desfase en las series de ISR e IVA, no tengo información completa para el segundo trimestre de 2021.

3.2 Quiebres Estructurales

Por otro lado, es necesario identificar si las series que se analizan presentan cambios estructurales que pueden deberse a mofificaciones legislativas u otro hechos que sean económicamente relevantes. Para revisar esto, se hace primero una regresión entre las variables que se van a analizar, en este caso, ISR e IVA, y la variable explicativa PIB. Luego se aplican las pruebas CUSUM y CUSUMQ para identificar posibles quiebres. Estos posibles quiebres se verifican con la prueba de Chow:

Variable	Quiebres CUSUM	Quiebres CUSUMQ	Valor p, prueba de Chow
ISR	2008-T1	2012-T1, 2014-T1, 2017-T2	2008-T1*, 2012-T1, 2014-T1, 2017-T2*
IVA	2018-T1	2010-T1	2018-T1, 2010-T1***

Como resultado de las pruebas, la prueba de Chow indica que todos los quiebres son significativos, mientras que en el caso del IVA, solo es significativo el quiebre del primer trimestre de 2010. Con estos resultados tengo que voltear a los cambios institucionales para identificar cuáles cambios coinciden con eventos económicos.

En el caso del IVA, el quiebre significativo coincide con el incremento general de la tasa del IVA del 15% al 16%. Por lo anterior, considero adecuado agregar una variable dummy, que indique la tasa general del impuesto antes de 2010 (15%, indicado con un cero) y después (16%, indicado con un uno). De hecho, antes de 2002, la tasa estuvo cambiando entre 10% y 15%, pero los datos de mi muestra no llegan hasta esas fechas.

Para el caso del ISR, y en cuanto a los quiebres significativos, únicamente identifico la reforma del ISR que entró en vigor a partir de 2014. Como resultado de dicha reforma, se tuvo una Ley del Impuesto Sobre la Renta completamente nueva, por lo que considero adecuado a utilizar una variable dummy que indica el período en el cual entro en vígor la última reforma del ISR (ceros antes de 2014 y a partir de esa fecha, unos).

De hecho, en ambos análisis del ISR y del IVA, también hice el análisis de incluir una variable dummy que señalara que los últimos 3 trimestres de 2020, donde hubo crisis económica (causada por la pandemia del Covid-19), pero resultó no ser significativa, y no la incluí en el modelo final. Esto se debe a que a través de la variable del PIB ya se incluye el efecto económico de la crisis en la estimación.

Para el caso de las otras variables no realicé la prueba CUSUM pues no se tiene variable explícativa. Sin embargo, sí hice un ajuste en el IEPS federal a gasolinas y diésel intencionalmente para dejar fuera la parte de la serie anterior a 2015. Esto lo hago porque justo en ese año hubo una reforma importante en este impuesto que pasó a cobrarse a cuotas por litro, mientras que antes de dicha reforma, el IEPS se utilizaba para controlar el precio de las gasolinas, e incluso podía llegar a ser negativo (subsidio al precio final).

4. Series mensuales

Con el análisis que he presentado hasta el momento, puedo comenzar a utilizar el módelo SARIMAX:

S, de seasonal, porque considera series que pueden tener estacionalidad.
AR, de autoregresivo, por series que dependen de observaciones previas.
I, de integradas, porque los modelos pueden considerar series que se integran, o dependen de si mismas de manera estacional (por ejemplo, anual).
MA, de moving average o media móvil porque considera series que pueden ser medias móviles.
X, por el vector X de variables regresivas que puede incluir.

En general, la metodología que sigue este modelo, es la metodología Box-Jenkins: 1) identificar y seleccionar el modelo con variables son estacionarias; 2) estimar los parámetros que mejor ajusten el modelo ARIMA seleccionado; 3) comprobar el modelo obtenido mediante distintas pruebas y 4) en caso de una incorrecta especificación regresar al paso 1).

Sin embargo, en la sección anterior, concluí que sólo utilizaré el PIB como variable explicativa para mis series de ISR e IVA. Por tanto, analizaré estas series desde su perspectiva mensual. Solo mostraré el análisis completo para la primera serie, y en el resto, dejaré solo las gráficas con los pronósticos.

4.1 IEPS Cervezas

La serie de IEPS a la venta de cervezas, y su descomposición en sus componentes de tendencia, estacional y residuales, se ve como sigue:

Para analizar si esta serie es estacionaria, se puede analizar el correlograma y el autocorrelograma, donde se ve que la serie tiene fuerte autocorrelación y autocorrelación parcial:

Sin embargo, para poder concluirlo adecuadamente, apliqué la prueba ADF que falla en rechazar que la serie tenga raíz unitaria. Con ello concluyó que la serie tiene raíz unitaria.

Results of Dickey-Fuller Test:
Test Statistic                  -0.236237
p-value                          0.934062
#Lags Used                      12.000000
Number of Observations Used    222.000000
Critical Value (1%)             -3.460154
Critical Value (5%)             -2.874649
Critical Value (10%)            -2.573757
dtype: float64

Por tanto, es necesario generar la serie en primeras diferencias, donde vuelvo a realizar la prueba ADF rechazó la hipótesis de que la serie tiene raíz unitaria, en favor de la hipótesis que la serie es estacionaria. En otras palabras, puedo concluir que la serie en primeras diferencias sí es estacionaria (como adelante desde la sección 3).

Results of Dickey-Fuller Test:
Test Statistic                -6.233054e+00
p-value                        4.900097e-08
#Lags Used                     1.500000e+01
Number of Observations Used    2.060000e+02
Critical Value (1%)           -3.462499e+00
Critical Value (5%)           -2.875675e+00
Critical Value (10%)          -2.574304e+00
dtype: float64

Luego, con este comando del paquete pmdarima, automáticamente se evalúan diferentes propuestas de modelos ARIMA, seleccionando el mejor modelo con base en el criterio de Akaike:

model = pm.auto_arima(ieps_cervezas, d=1, D=1,
                      seasonal=True, m=12, trend='c', 
                      start_p=0, start_q=0, max_order=9, test='adf', stepwise=True, trace=True)

Performing stepwise search to minimize aic
 ARIMA(0,1,0)(1,1,1)[12]             : AIC=3292.063, Time=2.28 sec
 ARIMA(0,1,0)(0,1,0)[12]             : AIC=3360.414, Time=0.13 sec
 ARIMA(1,1,0)(1,1,0)[12]             : AIC=3291.291, Time=0.87 sec
 ARIMA(0,1,1)(0,1,1)[12]             : AIC=3217.454, Time=0.69 sec
 ARIMA(0,1,1)(0,1,0)[12]             : AIC=3296.440, Time=0.34 sec
 ARIMA(0,1,1)(1,1,1)[12]             : AIC=3219.324, Time=1.14 sec
 ARIMA(0,1,1)(0,1,2)[12]             : AIC=3219.155, Time=5.85 sec
 ARIMA(0,1,1)(1,1,0)[12]             : AIC=3270.511, Time=2.53 sec
 ARIMA(0,1,1)(1,1,2)[12]             : AIC=inf, Time=7.40 sec
 ARIMA(0,1,0)(0,1,1)[12]             : AIC=3290.610, Time=1.59 sec
 ARIMA(1,1,1)(0,1,1)[12]             : AIC=3213.099, Time=2.84 sec
 ARIMA(1,1,1)(0,1,0)[12]             : AIC=3298.430, Time=0.52 sec
 ARIMA(1,1,1)(1,1,1)[12]             : AIC=3214.994, Time=4.20 sec
 ARIMA(1,1,1)(0,1,2)[12]             : AIC=3214.897, Time=10.21 sec
 ARIMA(1,1,1)(1,1,0)[12]             : AIC=inf, Time=1.10 sec
 ARIMA(1,1,1)(1,1,2)[12]             : AIC=3213.140, Time=9.82 sec
 ARIMA(1,1,0)(0,1,1)[12]             : AIC=3245.988, Time=0.76 sec
 ARIMA(2,1,1)(0,1,1)[12]             : AIC=inf, Time=2.96 sec
 ARIMA(1,1,2)(0,1,1)[12]             : AIC=inf, Time=1.92 sec
 ARIMA(0,1,2)(0,1,1)[12]             : AIC=3215.092, Time=1.21 sec
 ARIMA(2,1,0)(0,1,1)[12]             : AIC=3229.670, Time=3.08 sec
 ARIMA(2,1,2)(0,1,1)[12]             : AIC=inf, Time=1.68 sec
 ARIMA(1,1,1)(0,1,1)[12] intercept   : AIC=3213.099, Time=3.52 sec

Best model:  ARIMA(1,1,1)(0,1,1)[12]          
Total fit time: 66.918 seconds

model.summary()

SARIMAX Results
Dep. Variable:	y	No. Observations:	235
Model:	SARIMAX(1, 1, 1)x(0, 1, 1, 12)	Log Likelihood	-1601.550
Date:	Thu, 12 Aug 2021	AIC	3213.099
Time:	14:11:01	BIC	3230.113
Sample:	0	HQIC	3219.968
	- 235
Covariance Type:	opg

	coef	std err	z	P>\|z\|	[0.025	0.975]
intercept	0.1059	0.660	0.161	0.872	-1.187	1.399
ar.L1	0.2516	0.049	5.177	0.000	0.156	0.347
ma.L1	-0.9086	0.024	-37.853	0.000	-0.956	-0.862
ma.S.L12	-0.7614	0.044	-17.459	0.000	-0.847	-0.676
sigma2	1.031e+05	5163.368	19.973	0.000	9.3e+04	1.13e+05

Ljung-Box (L1) (Q):	0.18	Jarque-Bera (JB):	825.59
Prob(Q):	0.67	Prob(JB):	0.00
Heteroskedasticity (H):	9.40	Skew:	-0.10
Prob(H) (two-sided):	0.00	Kurtosis:	12.45

Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).

A diferencia de la ocasión anterior, ahora trabajaré con la serie completa (sin partirla en parte para entrenar y parte para hacer pruebas), para que el siguiente comando me de los mejores coeficientes para explicar la serie, con base en el modelo que obtuve en el paso anterior:

model = SARIMAX(ieps_cervezas,order=(1,1,1),seasonal_order=(0,1,1,12))
results = model.fit()
results.summary()

SARIMAX Results
Dep. Variable:	ieps_cervezas	No. Observations:	235
Model:	SARIMAX(1, 1, 1)x(0, 1, 1, 12)	Log Likelihood	-1601.564
Date:	Thu, 12 Aug 2021	AIC	3211.128
Time:	14:11:03	BIC	3224.738
Sample:	01-01-2002	HQIC	3216.623
	- 07-01-2021
Covariance Type:	opg

	coef	std err	z	P>\|z\|	[0.025	0.975]
ar.L1	0.2509	0.044	5.750	0.000	0.165	0.336
ma.L1	-0.9078	0.024	-37.752	0.000	-0.955	-0.861
ma.S.L12	-0.7610	0.041	-18.506	0.000	-0.842	-0.680
sigma2	1.032e+05	5026.424	20.523	0.000	9.33e+04	1.13e+05

Ljung-Box (L1) (Q):	0.18	Jarque-Bera (JB):	825.72
Prob(Q):	0.67	Prob(JB):	0.00
Heteroskedasticity (H):	9.37	Skew:	-0.09
Prob(H) (two-sided):	0.00	Kurtosis:	12.45

Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).

Los estadísticos para diagnosticar este modelo se ven bastante prometedores:

Los pronósticos intra muestra de este modelo parecen ajustar adecuadamente:

Mientras que los pronósticos hacia adelante también se ven bastante prometedores. Cabe señalar que esta es la parte que más me interesa para cada serie.

4.2 IEPS Bebidas Alcohólicas

Luego de repetir esta metodología para el caso del IEPS a las bebidas alcohólicas, obtuve el siguiente pronóstico para esta serie: