Modelos ARIMA y SARIMAX

Introducción

ARIMA (AutoRegressive Integrated Moving Average) y SARIMAX (Seasonal AutoRegressive Integrated Moving Average with eXogenous regressors) son modelos estadísticos ampliamente reconocidos y utilizados para la predicción de series temporales (forecasting). Este modelo consta de tres componentes. El elemento autorregresivo (AR) relaciona el valor actual con valores pasados (lags). El elemento de media móvil (MA) asume que el error de predicción es una combinación lineal de los errores de predicción pasados. Por último, el componente integrado (I) indica que los valores de la serie original han sido reemplazados por la diferencia entre valores consecutivos (y este proceso de diferencia puede haberse realizado más de una vez).

Si bien los modelos ARIMA son ampliamente conocidos, los modelos SARIMAX extienden el marco de ARIMA al incorporar patrones estacionales y variables exógenas.

# Librerías
# ======================================================================================
import numpy as np
import pandas as pd
from io import StringIO
import contextlib
import re
import matplotlib.pyplot as plt
plt.style.use('seaborn-v0_8-darkgrid')

# pmdarima
#!pip install pmdarima
import pmdarima
from pmdarima import ARIMA
from pmdarima import auto_arima

# statsmodels
import statsmodels
from statsmodels.tsa.statespace.sarimax import SARIMAX
from statsmodels.tsa.stattools import adfuller
from statsmodels.tsa.stattools import kpss
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.tsa.seasonal import seasonal_decompose

# skforecast
import skforecast
from skforecast.datasets import fetch_dataset
from skforecast.plot import set_dark_theme
from skforecast.sarimax import Sarimax
from skforecast.recursive import ForecasterSarimax
from skforecast.model_selection import TimeSeriesFold
from skforecast.model_selection import backtesting_sarimax
from skforecast.model_selection import grid_search_sarimax

import warnings
warnings.filterwarnings('once')
C:\Users\WINDOWS 11\Documents\.virtualenvs\r-tensorflow\lib\site-packages\tqdm\auto.py:21: TqdmWarning:

IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html

Datos

El conjunto de datos de este documento es un resumen del consumo mensual de combustible en España.

# Descarga datos
# ======================================================================================
datos = fetch_dataset(name='fuel_consumption', raw=True)
datos = datos[['Fecha', 'Gasolinas']]
datos = datos.rename(columns={'Fecha':'date', 'Gasolinas':'litters'})
datos['date'] = pd.to_datetime(datos['date'], format='%Y-%m-%d')
datos = datos.set_index('date')
datos = datos.loc[:'1990-01-01 00:00:00']
datos = datos.asfreq('MS')
datos = datos['litters']
datos.head(4)
fuel_consumption
----------------
Monthly fuel consumption in Spain from 1969-01-01 to 2022-08-01.
Obtained from Corporación de Reservas Estratégicas de Productos Petrolíferos and
Corporación de Derecho Público tutelada por el Ministerio para la Transición
Ecológica y el Reto Demográfico. https://www.cores.es/es/estadisticas
Shape of the dataset: (644, 6)
date
1969-01-01    166875.2129
1969-02-01    155466.8105
1969-03-01    184983.6699
1969-04-01    202319.8164
Freq: MS, Name: litters, dtype: float64
# Fechas Train-test
# ======================================================================================
set_dark_theme()
fin_train = '1980-01-01 23:59:59'
print(
    f"Fechas train : {datos.index.min()} --- {datos.loc[:fin_train].index.max()}  "
    f"(n={len(datos.loc[:fin_train])})"
)
print(
    f"Fechas test  : {datos.loc[fin_train:].index.min()} --- {datos.loc[:].index.max()}  "
    f"(n={len(datos.loc[fin_train:])})"
)
datos_train = datos.loc[:fin_train]
datos_test  = datos.loc[fin_train:]
datos_train
Fechas train : 1969-01-01 00:00:00 --- 1980-01-01 00:00:00  (n=133)
Fechas test  : 1980-02-01 00:00:00 --- 1990-01-01 00:00:00  (n=120)
date
1969-01-01    166875.2129
1969-02-01    155466.8105
1969-03-01    184983.6699
1969-04-01    202319.8164
1969-05-01    206259.1523
                 ...     
1979-09-01    476677.5163
1979-10-01    487880.0221
1979-11-01    462139.3874
1979-12-01    485646.8776
1980-01-01    413886.2617
Freq: MS, Name: litters, Length: 133, dtype: float64
# Gráfico
# ======================================================================================
fig, ax=plt.subplots(figsize=(8, 4))
datos_train.plot(ax=ax, label='train')
datos_test.plot(ax=ax, label='test')
ax.set_title('Consumo mensual combustible España')
ax.legend();
plt.show()

Análisis Exploratorio

Crear un modelo ARIMA requiere un análisis exploratorio exhaustivo. Este paso crítico sirve de brújula, guiando al analista hacia una comprensión detallada de la dinámica intrínseca de los datos. Antes de entrenar un modelo ARIMA a una serie temporal, es importante realizar un análisis exploratorio para determinar, como mínimo, lo siguiente:

  • Estacionariedad: la estacionariedad significa que las propiedades estadísticas (media, varianza…) permanecen constantes a lo largo del tiempo, por lo que las series temporales con tendencias o estacionalidad no son estacionarias. Dado que ARIMA presupone la estacionariedad de los datos, es esencial someterlos a pruebas rigurosas, como la prueba Dickey-Fuller aumentada, para evaluar que se cumple. Si se constata la no estacionariedad, las series deben diferenciarse hasta alcanzar la estacionariedad. Este análisis ayuda a determinar el valor óptimo del parámetro dd.

  • Análisis de autocorrelación: Graficar las funciones de autocorrelación y autocorrelación parcial (ACF y PACF) para identificar posibles relaciones de rezago (lags) entre los valores de la serie. Este análisis visual ayuda a determinar los términos autorregresivos (AR) y de media móvil (MA) adecuados (pp y qq) para el modelo ARIMA.

  • Descomposición estacional: en los casos donde se sospecha de estacionalidad, descomponer la serie en componentes de tendencia, estacionales y residuales utilizando técnicas como las medias móviles la descomposición estacional de series temporales (STL) puede revelar patrones ocultos y ayudar a identificar la estacionalidad. Este análisis ayuda a determinar los valores óptimos de los parámetros PP, DD, QQ y mm.

Estos análisis exploratorios establecen la base para empezar a construir un modelo ARIMA efectivo que capture los patrones fundamentales y las asociaciones dentro de los datos.

Estacionariedad

Existen varios métodos para evaluar si una serie temporal es estacionaria o no estacionaria:

  • Inspección visual de la serie temporal: inspeccionando visualmente el gráfico de la serie temporal, es posible identificar la presencia de una tendencia o estacionalidad notables. Si se observan estos patrones, es probable que la serie no sea estacionaria.

  • Valores estadísticos: calcular estadísticos como la media y la varianza, de varios segmentos de la serie. Si existen diferencias significativas, la serie no es estacionaria.

  • Pruebas estadísticas: utilizar test estadísticos como la prueba Dickey-Fuller aumentada o la prueba Kwiatkowski-Phillips-Schmidt-Shin (KPSS).

El gráfico generado en el apartado anterior muestra una clara tendencia positiva, lo que indica un aumento constante a lo largo del tiempo. En consecuencia, la media de la serie aumenta con el tiempo, lo que confirma su no estacionariedad.

La diferenciación es una de las técnicas más sencillas para eliminar la tendencia de una serie temporal. Consiste en generar una nueva serie en la que cada valor se calcula como la diferencia entre el valor actual y el valor anterior, es decir, la diferencia entre valores consecutivos. Matemáticamente, la primera diferencia se calcula como:

ΔXt=Xt−Xt−1ΔXt=Xt−Xt−1

Donde XtXt es el valor en el tiempo tt y Xt−1Xt−1 es el valor en el tiempo t−1t−1. Esta es conocida como diferenciación de primer orden. Este proceso se puede repetir si es necesario hasta que se alcance la estacionariedad deseada.

En los siguientes partado, la serie temporal original se somete a una diferenciación de primer y segundo orden y se aplican pruebas estadísticas para determinar si se consigue la estacionariedad.

Prueba de Dickey-Fuller aumentada

La prueba Dickey-Fuller aumentada considera como hipótesis nula que la serie temporal tiene una raíz unitaria, una característica frecuente de las series temporales no estacionarias. Por el contrario, la hipótesis alternativa (bajo la cual se rechaza la hipótesis nula) es que la serie es estacionaria.

  • Hipótesis nula (HOHO): La serie tiene una raíz unitaria, no es estacionaria.

  • Hipótesis alternativa (HAHA): La serie no tiene raíz unitaria, es estacionaria.

Dado que la hipótesis nula supone la presencia de una raíz unitaria, el p-value obtenido debe ser inferior a un nivel de significación determinado, a menudo fijado en 0.05, para rechazar esta hipótesis. Este resultado indica la estacionariedad de la serie. La función adfuller() de la biblioteca Statsmodels permite aplicar la prueba ADF. Su resultado incluye cuatro valores: el p-value, el valor del estadístico, el número de retardos (lags) incluidos en la prueba y los umbrales del valor crítico para tres niveles diferentes de significancia.

Prueba Kwiatkowski-Phillips-Schmidt-Shin (KPSS).

La prueba KPSS comprueba si una serie temporal es estacionaria en torno a una media o una tendencia lineal. En esta prueba, la hipótesis nula es que la serie es estacionaria. Por consiguiente, los p-values pequeños (por ejemplo, inferiores a 0.05) rechazan la hipótesis nula y sugieren que es necesario diferenciar. La librería Statsmodels proporciona una implementación de la prueba KPSS a través de la función kpss().

# Test estacionariedad
# ==============================================================================
warnings.filterwarnings("ignore")

datos_diff_1 = datos_train.diff().dropna()
datos_diff_2 = datos_diff_1.diff().dropna()

print('Test estacionariedad serie original')
print('-------------------------------------')
adfuller_result = adfuller(datos)
kpss_result = kpss(datos)
print(f'ADF Statistic: {adfuller_result[0]}, p-value: {adfuller_result[1]}')
print(f'KPSS Statistic: {kpss_result[0]}, p-value: {kpss_result[1]}')

print('\nTest estacionariedad para serie diferenciada (order=1)')
print('--------------------------------------------------')
adfuller_result = adfuller(datos_diff_1)
kpss_result = kpss(datos.diff().dropna())
print(f'ADF Statistic: {adfuller_result[0]}, p-value: {adfuller_result[1]}')
print(f'KPSS Statistic: {kpss_result[0]}, p-value: {kpss_result[1]}')

print('\nTest estacionariedad para serie diferenciada (order=2)')
print('--------------------------------------------------')
adfuller_result = adfuller(datos_diff_2)
kpss_result = kpss(datos.diff().diff().dropna())
print(f'ADF Statistic: {adfuller_result[0]}, p-value: {adfuller_result[1]}')
print(f'KPSS Statistic: {kpss_result[0]}, p-value: {kpss_result[1]}')

warnings.filterwarnings("default")

# Gráfico series
# ==============================================================================
fig, axs = plt.subplots(nrows=3, ncols=1, figsize=(8, 5), sharex=True)
datos.plot(ax=axs[0], title='Serie original')
datos_diff_1.plot(ax=axs[1], title='Diferenciación orden 1')
datos_diff_2.plot(ax=axs[2], title='Diferenciación orden 2');
plt.show()
Test estacionariedad serie original
-------------------------------------
ADF Statistic: -0.4461298099822783, p-value: 0.9021071923942668
KPSS Statistic: 2.2096370946978383, p-value: 0.01

Test estacionariedad para serie diferenciada (order=1)
--------------------------------------------------
ADF Statistic: -3.641727690032306, p-value: 0.005011605002137535
KPSS Statistic: 0.3132711623572804, p-value: 0.1

Test estacionariedad para serie diferenciada (order=2)
--------------------------------------------------
ADF Statistic: -8.233942641655862, p-value: 5.959599575501017e-13
KPSS Statistic: 0.080656682674821, p-value: 0.1

El p-value obtenido tras la primera diferenciación es estadísticamente significativo acorde al umbral ampliamente reconocido y aceptado de 0.05. Por lo tanto, la selección más adecuada para el parámetro ARIMA dd es 1.

Análisis de autocorrelación

El gráfico de la función de autocorrelación ( Autocorrelation Function ACF) y la función de autocorrelación parcial (Partial Autocorrelation Function (PACF)) de la serie temporal proporciona información útil sobre los posibles valores adecuados de pp y qq. La ACF ayuda a identificar el valor de qq (retardos en la parte de media móvil), mientras que la PACF ayuda a identificar el valor de pp (retardos en la parte autorregresiva).

Función de autocorrelación (ACF)

La ACF calcula la correlación entre una serie temporal y sus valores retardados (lags). En el contexto de la modelización ARIMA, una caída brusca de la ACF después de unos pocos retardos indica que los datos tienen un orden autorregresivo finito. El retardo en el que cae la ACF proporciona una estimación del valor de qq. Si el ACF muestra un patrón sinusoidal o sinusoidal amortiguado, sugiere la presencia de estacionalidad y requiere la consideración de órdenes estacionales además de órdenes no estacionales.

Función de autocorrelación parcial (PACF)

La PACF mide la correlación entre un valor retardado (lag) y el valor actual de la serie temporal, teniendo en cuenta el efecto de los retardos intermedios. En el contexto de la modelización ARIMA, si la PACF se corta bruscamente después de un determinado retardo, mientras que los valores restantes están dentro del intervalo de confianza, sugiere un modelo AR de ese orden. El desfase en el que se corta el PACF da una idea del valor de pp.

# Grafico de autocorrelación para la serie original y la serie diferenciada
# ==============================================================================
fig, axs = plt.subplots(nrows=2, ncols=1, figsize=(8, 5), sharex=True)
plot_acf(datos, ax=axs[0], lags=50, alpha=0.05)
axs[0].set_title('Autocorrelación serie original')
plot_acf(datos_diff_1, ax=axs[1], lags=50, alpha=0.05)
axs[1].set_title('Autocorrelación serie diferenciada (order=1)');
plt.show()

# Autocorrelación parcial para la serie original y la serie diferenciada
# ==============================================================================
fig, axs = plt.subplots(nrows=2, ncols=1, figsize=(8, 5), sharex=True)
plot_pacf(datos, ax=axs[0], lags=50, alpha=0.05)
axs[0].set_title('Autocorrelación parcial serie original')
plot_pacf(datos_diff_1, ax=axs[1], lags=50, alpha=0.05)
axs[1].set_title('Autoorrelación parcial serie diferenciada (order=1)');
plt.tight_layout();
plt.show()

Acorde a la función de autocorrelación parcial, el valor óptimo para el parámetro pp es 0. Sin embargo, se va a asignar un valor de 1 para proporcionar un componente autorregresivo al modelo. En cuanto al componente qq, la función de autocorrelación sugiere un valor de 1.

Descomposición de series temporales

La descomposición de series temporales consiste en descomponer la serie temporal original en sus componentes fundamentales: la tendencia, la estacionalidad y los residuos. Esta descomposición puede llevarse a cabo de manera aditiva o multiplicativa. Al combinar la descomposición de las series temporales con el análisis de la ACF y la PACF, se obtiene una descripción bastante completa con la que comprender la estructura subyacente de los datos y acotar el valor los parámetros ARIMA más apropiados.

# Descomposición de la serie original y la serie diferenciada
# ==============================================================================
res_decompose = seasonal_decompose(datos, model='additive', extrapolate_trend='freq')
res_descompose_diff_2 = seasonal_decompose(datos_diff_1, model='additive', extrapolate_trend='freq')

fig, axs = plt.subplots(nrows=4, ncols=2, figsize=(8, 6), sharex=True)
res_decompose.observed.plot(ax=axs[0, 0])
axs[0, 0].set_title('Serie original', fontsize=12)
res_decompose.trend.plot(ax=axs[1, 0])
axs[1, 0].set_title('Tendencia', fontsize=12)
res_decompose.seasonal.plot(ax=axs[2, 0])
axs[2, 0].set_title('Estacionalidad', fontsize=12)
res_decompose.resid.plot(ax=axs[3, 0])
axs[3, 0].set_title('Residuos', fontsize=12)
res_descompose_diff_2.observed.plot(ax=axs[0, 1])
axs[0, 1].set_title('Series diferenciadas (order=1)', fontsize=12)
res_descompose_diff_2.trend.plot(ax=axs[1, 1])
axs[1, 1].set_title('Tendencia', fontsize=12)
res_descompose_diff_2.seasonal.plot(ax=axs[2, 1])
axs[2, 1].set_title('Estacionalidad', fontsize=12)
res_descompose_diff_2.resid.plot(ax=axs[3, 1])
axs[3, 1].set_title('Residuos', fontsize=12)
fig.suptitle('Descomposición de la serie original vs serie diferenciada', fontsize=14)
fig.tight_layout();
plt.show()

El patrón recurrente cada 12 meses sugiere una estacionalidad anual, probablemente influenciada por factores vacacionales. El gráfico de ACF respalda aún más la presencia de esta estacionalidad, ya que se observan picos significativos en los lags correspondientes a los intervalos de 12 meses, confirmando la idea de patrones recurrentes.

Conclusiones

Basandose en los resultados del análisis exploratorio, utilizar una combinación de diferenciación de primer orden y diferenciación estacional puede ser el enfoque más apropiado. La diferenciación de primer orden es efectiva para capturar las transiciones entre observaciones y resaltar las fluctuaciones a corto plazo. Al mismo tiempo, la diferenciación estacional, que abarca un período de 12 meses y representa el cambio de un año a otro, captura de manera efectiva los patrones cíclicos inherentes en los datos. Este enfoque nos permite lograr la estacionariedad necesaria para el proceso de modelado ARIMA subsiguiente.

# Diferenciaciación de orden 1 combinada con diferenciación estacional
# ==============================================================================
datos_diff_1_12 = datos_train.diff().diff(12).dropna()

warnings.filterwarnings("ignore")
adfuller_result = adfuller(datos_diff_1_12)
print(f'ADF Statistic: {adfuller_result[0]}, p-value: {adfuller_result[1]}')
kpss_result = kpss(datos_diff_1_12)
print(f'KPSS Statistic: {kpss_result[0]}, p-value: {kpss_result[1]}')
warnings.filterwarnings("default")
ADF Statistic: -4.387457230769962, p-value: 0.0003123773271126866
KPSS Statistic: 0.06291573421251048, p-value: 0.1

Modelo ARIMA-SARIMAX

La siguiente sección muestra cómo entrenar un modelo ARIMA-SARIMAX y PREDECIR valores futuros utilizando cada una de las tres librerías.

Statsmodels

En statsmodels, se diferencia entre el proceso de definir un modelo y entrenarlo. Este enfoque puede resultar familiar para usuarios del lenguaje de programación R, pero puede parecer algo menos convencional para aquellos acostumbrados a librerías como scikit-learn o XGBoost en el ecosistema de Python.

El proceso comienza con la definición del modelo, que incluye los parámetros configurables y el conjunto de datos de entrenamiento. Cuando se invoca al método de ajuste (fit). En lugar de modificar el objeto modelo, como es típico en las librerías de Python, statsmodels crea un nuevo objeto SARIMAXResults. Este objeto no solo encapsula detalles esenciales como los residuos y los parámetros aprendidos, sino que también proporciona las herramientas necesarias para generar predicciones.

# Modelo SARIMAX con statsmodels.Sarimax
# ==============================================================================
warnings.filterwarnings("ignore", category=UserWarning, message='Non-invertible|Non-stationary')
modelo = SARIMAX(endog = datos_train, order = (1, 1, 1), seasonal_order = (1, 1, 1, 12))
modelo_res = modelo.fit(disp=0)
warnings.filterwarnings("default")
modelo_res.summary()
SARIMAX Results
Dep. Variable: litters No. Observations: 133
Model: SARIMAX(1, 1, 1)x(1, 1, 1, 12) Log Likelihood -1356.051
Date: Thu, 29 May 2025 AIC 2722.103
Time: 11:52:18 BIC 2736.040
Sample: 01-01-1969 HQIC 2727.763
- 01-01-1980
Covariance Type: opg
coef std err z P>|z| [0.025 0.975]
ar.L1 -0.4972 0.134 -3.707 0.000 -0.760 -0.234
ma.L1 -0.0096 0.146 -0.066 0.947 -0.295 0.276
ar.S.L12 0.0465 0.162 0.288 0.774 -0.270 0.364
ma.S.L12 -0.3740 0.203 -1.847 0.065 -0.771 0.023
sigma2 3.291e+08 1.06e-09 3.1e+17 0.000 3.29e+08 3.29e+08
Ljung-Box (L1) (Q): 5.13 Jarque-Bera (JB): 18.12
Prob(Q): 0.02 Prob(JB): 0.00
Heteroskedasticity (H): 1.26 Skew: -0.42
Prob(H) (two-sided): 0.46 Kurtosis: 4.71


Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
[2] Covariance matrix is singular or near-singular, with condition number 1.33e+33. Standard errors may be unstable.

El resumen del modelo muestra mucha información sobre el proceso de ajuste:

  • Estadísticas de Ajuste del Modelo: Esta parte incluye varias estadísticas que ayudan a evaluar qué tan bien el modelo se ajusta a los datos observados:

    • Log-Likelihood (Logaritmo de la Verosimilitud): Una medida de qué tan bien el modelo explica los datos observados, donde valores más negativos indican un ajuste deficiente a los datos y valores más cercanos a cero indican un mejor ajuste.

    • AIC (Criterio de Información de Akaike): Una métrica de bondad de ajuste que equilibra el ajuste del modelo con su complejidad. Cuanto menor el valor de AIC mejor es el modelo.

    • BIC (Criterio de Información Bayesiano): Similar al AIC, pero penaliza más la complejidad del modelo. Al igual que con el AIC, valores más bajos de BIC indican un mejor ajuste.

    • HQIC (Criterio de Información de Hannan-Quinn): Otro criterio de selección de modelo, similar al AIC y al BIC.

  • Coeficientes: Esta tabla lista los coeficientes estimados para los parámetros del modelo. Incluye tanto los parámetros autoregresivos (AR) como los parámetros de media móvil (MA), así como cualquier variable exógena si se incluyen en el modelo. También incluye los errores estándar asociados con los coeficientes estimados para indicar la incertidumbre de dichas estimaciones, sus p-values, que se utilizan para evaluar la significancia de cada coeficiente, y el intervalo de confianza del 95%.

  • Diagnósticos del modelo: Esta sección proporciona información sobre los residuos. Las diferencias entre los valores observados (valores de entrenamiento) y los valores predichos por el modelo.

    • Prueba Ljung-Box: Una prueba de autocorrelación en los residuos.

    • Prueba de Jarque-Bera: Una prueba de normalidad de los residuos.

    • Asimetría y curtosis: Medidas de la forma de la distribución de los residuos.

# Predicción
# ==============================================================================
predicciones_statsmodels = modelo_res.get_forecast(steps=len(datos_test)).predicted_mean
predicciones_statsmodels.name = 'predicciones_statsmodels'
predicciones_statsmodels.head(4)
1980-02-01    407504.056933
1980-03-01    473997.245805
1980-04-01    489983.091489
1980-05-01    485517.462861
Freq: MS, Name: predicciones_statsmodels, dtype: float64

Skforecast

Skforecast adapta el modelo statsmodels.SARIMAX a la API de scikit-learn.

# Modelo SARIMAX con skforecast.Sarimax
# ==============================================================================
warnings.filterwarnings("ignore", category=UserWarning, message='Non-invertible|Non-stationary')
modelo = Sarimax(order=(1, 1, 1), seasonal_order=(1, 1, 1, 12))
y=datos_train
y
modelo.fit(y=datos_train)
modelo.summary()
warnings.filterwarnings("default")
# Predictión
# ==============================================================================
predicciones_skforecast = modelo.predict(steps=len(datos_test))
predicciones_skforecast.head(4)
pred
1980-02-01 407504.056933
1980-03-01 473997.245805
1980-04-01 489983.091489
1980-05-01 485517.462861

pdmarima

# Modelo SARIMAX con pdmarima.Sarimax
# ==============================================================================
warnings.filterwarnings("ignore", message=".*force_all_finite.*", category=FutureWarning)
modelo = ARIMA(order=(1, 1, 1), seasonal_order=(1, 1, 1, 12))
modelo.fit(y=datos_train)
modelo.summary()
SARIMAX Results
Dep. Variable: y No. Observations: 133
Model: SARIMAX(1, 1, 1)x(1, 1, 1, 12) Log Likelihood -1355.749
Date: Thu, 29 May 2025 AIC 2723.498
Time: 11:52:19 BIC 2740.223
Sample: 01-01-1969 HQIC 2730.290
- 01-01-1980
Covariance Type: opg
coef std err z P>|z| [0.025 0.975]
intercept -474.5820 1101.722 -0.431 0.667 -2633.917 1684.753
ar.L1 -0.4896 0.138 -3.554 0.000 -0.760 -0.220
ma.L1 -0.0211 0.151 -0.139 0.889 -0.317 0.275
ar.S.L12 0.0545 0.164 0.331 0.740 -0.268 0.377
ma.S.L12 -0.3841 0.204 -1.884 0.060 -0.784 0.015
sigma2 3.289e+08 0.002 1.84e+11 0.000 3.29e+08 3.29e+08
Ljung-Box (L1) (Q): 4.90 Jarque-Bera (JB): 18.55
Prob(Q): 0.03 Prob(JB): 0.00
Heteroskedasticity (H): 1.27 Skew: -0.43
Prob(H) (two-sided): 0.46 Kurtosis: 4.73


Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
[2] Covariance matrix is singular or near-singular, with condition number 1.73e+27. Standard errors may be unstable. 
# Prediction
# ==============================================================================
warnings.filterwarnings("ignore", message=".*force_all_finite.*", category=FutureWarning)
predicciones_pdmarima = modelo.predict(len(datos_test))
predicciones_pdmarima.name = 'predicciones_pdmarima'
predicciones_pdmarima.head(4)
1980-02-01    406998.311391
1980-03-01    472944.444482
1980-04-01    488389.125304
1980-05-01    483432.075696
Freq: MS, Name: predicciones_pdmarima, dtype: float64
# Plot predictions
# ==============================================================================
fig, ax = plt.subplots(figsize=(8, 4))
datos_train.plot(ax=ax, label='train')
datos_test.plot(ax=ax, label='test')
predicciones_statsmodels.plot(ax=ax, label='statsmodels')
predicciones_skforecast.columns = ['skforecast']
predicciones_skforecast.plot(ax=ax, label='skforecast')
predicciones_pdmarima.plot(ax=ax, label='pmdarima')
ax.set_title('Predicciones con modelos ARIMA')
ax.legend();
plt.show()

ForecasterSarimax

La clase ForecasterSarimax permite entrenar y validar modelos ARIMA y SARIMAX utilizando la API de skforecast. Dado que ForecasterSarimax sigue la misma API que los otros Forecasters disponibles en la librería, es muy fácil hacer una comparación robusta del rendimiento de modelos ARIMA-SARIMAX frente a otros modelos de machine learning como Random Forest or Gradient Boosting.

Entrenamiento - Predicción

El proceso de entrenamiento y predicción sigue una API similar a la de scikit-learn. Puedes encontrar más detalles en la guía del usuario de ForecasterSarimax.

# Modelo ARIMA con ForecasterSarimax y skforecast Sarimax
# ==============================================================================
forecaster = ForecasterSarimax(
                 regressor=Sarimax(order=(1, 1, 1), seasonal_order=(1, 1, 1, 12))
             )
forecaster.fit(y=datos_train, suppress_warnings=True)

# Predicción
predicciones = forecaster.predict(steps=len(datos_test))
predicciones.head(4)
1980-02-01    407504.056933
1980-03-01    473997.245805
1980-04-01    489983.091489
1980-05-01    485517.462861
Freq: MS, Name: pred, dtype: float64

Backtesting

El siguiente ejemplo muestra el uso de backtesting para evaluar el rendimiento del modelo SARIMAX al generar predicciones para los 12 meses siguientes en un plan anual. En este contexto, se genera una previsión al final de cada mes de diciembre, prediciendo valores para los 12 meses siguientes.

# Backtest forecaster
# ==============================================================================
forecaster = ForecasterSarimax(
                 regressor=Sarimax(
                                order=(1, 1, 1),
                                seasonal_order=(1, 1, 1, 12),
                                maxiter=200
                            )
             )
cv = TimeSeriesFold(
        steps              = 12,
        initial_train_size = len(datos_train),
        refit              = True,
        fixed_train_size   = False,
)
metrica, predicciones = backtesting_sarimax(
                            forecaster            = forecaster,
                            y                     = datos,
                            cv                    = cv,
                            metric                = 'mean_absolute_error',
                            n_jobs                = 1,
                            suppress_warnings_fit = True,
                            verbose               = True,
                            show_progress         = True
                        )
metrica
predicciones.head(4)
Information of folds
--------------------
Number of observations used for initial training: 133
Number of observations used for backtesting: 120
    Number of folds: 10
    Number skipped folds: 0 
    Number of steps per fold: 12
    Number of steps to exclude between last observed data (last window) and predictions (gap): 0

Fold: 0
    Training:   1969-01-01 00:00:00 -- 1980-01-01 00:00:00  (n=133)
    Validation: 1980-02-01 00:00:00 -- 1981-01-01 00:00:00  (n=12)
Fold: 1
    Training:   1969-01-01 00:00:00 -- 1981-01-01 00:00:00  (n=145)
    Validation: 1981-02-01 00:00:00 -- 1982-01-01 00:00:00  (n=12)
Fold: 2
    Training:   1969-01-01 00:00:00 -- 1982-01-01 00:00:00  (n=157)
    Validation: 1982-02-01 00:00:00 -- 1983-01-01 00:00:00  (n=12)
Fold: 3
    Training:   1969-01-01 00:00:00 -- 1983-01-01 00:00:00  (n=169)
    Validation: 1983-02-01 00:00:00 -- 1984-01-01 00:00:00  (n=12)
Fold: 4
    Training:   1969-01-01 00:00:00 -- 1984-01-01 00:00:00  (n=181)
    Validation: 1984-02-01 00:00:00 -- 1985-01-01 00:00:00  (n=12)
Fold: 5
    Training:   1969-01-01 00:00:00 -- 1985-01-01 00:00:00  (n=193)
    Validation: 1985-02-01 00:00:00 -- 1986-01-01 00:00:00  (n=12)
Fold: 6
    Training:   1969-01-01 00:00:00 -- 1986-01-01 00:00:00  (n=205)
    Validation: 1986-02-01 00:00:00 -- 1987-01-01 00:00:00  (n=12)
Fold: 7
    Training:   1969-01-01 00:00:00 -- 1987-01-01 00:00:00  (n=217)
    Validation: 1987-02-01 00:00:00 -- 1988-01-01 00:00:00  (n=12)
Fold: 8
    Training:   1969-01-01 00:00:00 -- 1988-01-01 00:00:00  (n=229)
    Validation: 1988-02-01 00:00:00 -- 1989-01-01 00:00:00  (n=12)
Fold: 9
    Training:   1969-01-01 00:00:00 -- 1989-01-01 00:00:00  (n=241)
    Validation: 1989-02-01 00:00:00 -- 1990-01-01 00:00:00  (n=12)

  0%|          | 0/10 [00:00<?, ?it/s] 20%|██        | 2/10 [00:00<00:00, 15.59it/s] 40%|████      | 4/10 [00:00<00:00,  8.77it/s] 60%|██████    | 6/10 [00:00<00:00,  7.58it/s] 70%|███████   | 7/10 [00:00<00:00,  7.04it/s] 80%|████████  | 8/10 [00:01<00:00,  6.82it/s] 90%|█████████ | 9/10 [00:01<00:00,  6.45it/s]100%|██████████| 10/10 [00:01<00:00,  6.02it/s]100%|██████████| 10/10 [00:01<00:00,  6.95it/s]
pred
1980-02-01 407504.056933
1980-03-01 473997.245805
1980-04-01 489983.091489
1980-05-01 485517.462861 
# Gráfico predicciones de backtesting
# ==============================================================================
fig, ax = plt.subplots(figsize=(8, 4))
datos.loc[fin_train:].plot(ax=ax, label='test')
predicciones.plot(ax=ax)
ax.set_title('Predicciones de backtesting con un modelo SARIMAX')
ax.legend();
plt.show()

Busqueda de hiperparámetros p, d, q

El análisis exploratorio ha reducido el espacio de búsqueda para los hiperparámetros óptimos del modelo. Sin embargo, para determinar definitivamente los valores más apropiados, es esencial utilizar métodos de búsqueda estratégicos.

# Train-validation-test
# ======================================================================================
fin_train = '1976-01-01 23:59:59'
fin_val = '1984-01-01 23:59:59'
print(
    f"Fechas entrenamiento : {datos.index.min()} --- {datos.loc[:fin_train].index.max()}  "
    f"(n={len(datos.loc[:fin_train])})"
)
print(
    f"Fechas validacion    : {datos.loc[fin_train:].index.min()} --- {datos.loc[:fin_val].index.max()}  "
    f"(n={len(datos.loc[fin_train:fin_val])})"
)
print(
    f"Fechas test          : {datos.loc[fin_val:].index.min()} --- {datos.index.max()}  "
    f"(n={len(datos.loc[fin_val:])})"
)

# Gráfico
# ======================================================================================
fig, ax = plt.subplots(figsize=(8, 4))
datos.loc[:fin_train].plot(ax=ax, label='entrenamiento')
datos.loc[fin_train:fin_val].plot(ax=ax, label='validación')
datos.loc[fin_val:].plot(ax=ax, label='test')
ax.set_title('Consumo mensual combustible España')
ax.legend();
plt.show()
Fechas entrenamiento : 1969-01-01 00:00:00 --- 1976-01-01 00:00:00  (n=85)
Fechas validacion    : 1976-02-01 00:00:00 --- 1984-01-01 00:00:00  (n=96)
Fechas test          : 1984-02-01 00:00:00 --- 1990-01-01 00:00:00  (n=72)

# Grid search basado en backtesting
# ==============================================================================
forecaster = ForecasterSarimax(
                 regressor=Sarimax(
                                order=(1, 1, 1), # Placeholder replaced in the grid search
                                maxiter=500
                            )
             )

param_grid = {
    'order': [(0, 1, 0), (0, 1, 1), (1, 1, 0), (1, 1, 1), (2, 1, 1)],
    'seasonal_order': [(0, 0, 0, 0), (0, 1, 0, 12), (1, 1, 1, 12)],
    'trend': [None, 'n', 'c']
}

cv = TimeSeriesFold(
        steps              = 12,
        initial_train_size = len(datos_train),
        refit              = True,
        fixed_train_size   = False,
    )

resultados_grid = grid_search_sarimax(
                   forecaster            = forecaster,
                   y                     = datos.loc[:fin_val],
                   cv                    = cv,
                   param_grid            = param_grid,
                   metric                = 'mean_absolute_error',
                   return_best           = False,
                   suppress_warnings_fit = True,
               )
resultados_grid.head(5)
Number of models compared: 45.
params grid:   0%|          | 0/45 [00:00<?, ?it/s]params grid:   9%|▉         | 4/45 [00:00<00:01, 31.12it/s]params grid:  18%|█▊        | 8/45 [00:00<00:04,  9.01it/s]params grid:  22%|██▏       | 10/45 [00:01<00:04,  7.56it/s]params grid:  29%|██▉       | 13/45 [00:01<00:03,  9.66it/s]params grid:  33%|███▎      | 15/45 [00:01<00:03,  9.70it/s]params grid:  38%|███▊      | 17/45 [00:02<00:05,  5.17it/s]params grid:  40%|████      | 18/45 [00:02<00:06,  4.08it/s]params grid:  47%|████▋     | 21/45 [00:03<00:03,  6.26it/s]params grid:  51%|█████     | 23/45 [00:03<00:02,  7.40it/s]params grid:  56%|█████▌    | 25/45 [00:03<00:03,  6.07it/s]params grid:  60%|██████    | 27/45 [00:04<00:04,  4.17it/s]params grid:  62%|██████▏   | 28/45 [00:04<00:03,  4.63it/s]params grid:  67%|██████▋   | 30/45 [00:04<00:02,  5.43it/s]params grid:  69%|██████▉   | 31/45 [00:04<00:02,  5.83it/s]params grid:  71%|███████   | 32/45 [00:05<00:02,  6.25it/s]params grid:  73%|███████▎  | 33/45 [00:05<00:01,  6.44it/s]params grid:  76%|███████▌  | 34/45 [00:05<00:02,  3.90it/s]params grid:  78%|███████▊  | 35/45 [00:06<00:03,  2.98it/s]params grid:  80%|████████  | 36/45 [00:06<00:03,  2.42it/s]params grid:  82%|████████▏ | 37/45 [00:07<00:02,  3.06it/s]params grid:  84%|████████▍ | 38/45 [00:07<00:01,  3.80it/s]params grid:  87%|████████▋ | 39/45 [00:07<00:01,  4.40it/s]params grid:  89%|████████▉ | 40/45 [00:07<00:01,  4.25it/s]params grid:  91%|█████████ | 41/45 [00:07<00:00,  4.21it/s]params grid:  93%|█████████▎| 42/45 [00:08<00:00,  3.90it/s]params grid:  96%|█████████▌| 43/45 [00:09<00:01,  1.72it/s]params grid:  98%|█████████▊| 44/45 [00:10<00:00,  1.24it/s]params grid: 100%|██████████| 45/45 [00:12<00:00,  1.06it/s]params grid: 100%|██████████| 45/45 [00:12<00:00,  3.74it/s]
params mean_absolute_error order seasonal_order trend
0 {'order': (0, 1, 1), 'seasonal_order': (1, 1, ... 18789.897510 (0, 1, 1) (1, 1, 1, 12) n
1 {'order': (0, 1, 1), 'seasonal_order': (1, 1, ... 18789.897510 (0, 1, 1) (1, 1, 1, 12) None
2 {'order': (1, 1, 1), 'seasonal_order': (1, 1, ... 19897.376863 (1, 1, 1) (1, 1, 1, 12) n
3 {'order': (1, 1, 1), 'seasonal_order': (1, 1, ... 19897.376863 (1, 1, 1) (1, 1, 1, 12) None
4 {'order': (2, 1, 1), 'seasonal_order': (1, 1, ... 20176.735382 (2, 1, 1) (1, 1, 1, 12) n 
# Auto arima: seleccion basada en AIC
# ==============================================================================
warnings.filterwarnings("ignore", message=".*force_all_finite.*", category=FutureWarning)
modelo = auto_arima(
            y                 = datos.loc[:fin_val],
            start_p           = 0,
            start_q           = 0,
            max_p             = 3,
            max_q             = 3,
            seasonal          = True,
            test              = 'adf',
            m                 = 12,   # periodicidad de la estacionalidad
            d                 = None, # El algoritmo determina 'd'
            D                 = None, # El algoritmo determina 'D'
            trace             = True,
            error_action      = 'ignore',
            suppress_warnings = True,
            stepwise          = True
        )
Performing stepwise search to minimize aic
 ARIMA(0,1,0)(1,1,1)[12]             : AIC=3903.204, Time=0.07 sec
 ARIMA(0,1,0)(0,1,0)[12]             : AIC=3942.897, Time=0.01 sec
 ARIMA(1,1,0)(1,1,0)[12]             : AIC=3846.786, Time=0.06 sec
 ARIMA(0,1,1)(0,1,1)[12]             : AIC=3840.318, Time=0.09 sec
 ARIMA(0,1,1)(0,1,0)[12]             : AIC=3873.797, Time=0.03 sec
 ARIMA(0,1,1)(1,1,1)[12]             : AIC=3841.882, Time=0.13 sec
 ARIMA(0,1,1)(0,1,2)[12]             : AIC=3841.572, Time=0.69 sec
 ARIMA(0,1,1)(1,1,0)[12]             : AIC=3852.231, Time=0.09 sec
 ARIMA(0,1,1)(1,1,2)[12]             : AIC=3842.593, Time=1.04 sec
 ARIMA(0,1,0)(0,1,1)[12]             : AIC=3904.615, Time=0.05 sec
 ARIMA(1,1,1)(0,1,1)[12]             : AIC=3834.135, Time=0.14 sec
 ARIMA(1,1,1)(0,1,0)[12]             : AIC=3866.187, Time=0.03 sec
 ARIMA(1,1,1)(1,1,1)[12]             : AIC=3835.564, Time=0.14 sec
 ARIMA(1,1,1)(0,1,2)[12]             : AIC=3835.160, Time=0.71 sec
 ARIMA(1,1,1)(1,1,0)[12]             : AIC=3844.410, Time=0.09 sec
 ARIMA(1,1,1)(1,1,2)[12]             : AIC=3836.443, Time=1.10 sec
 ARIMA(1,1,0)(0,1,1)[12]             : AIC=3836.696, Time=0.08 sec
 ARIMA(2,1,1)(0,1,1)[12]             : AIC=3836.104, Time=0.30 sec
 ARIMA(1,1,2)(0,1,1)[12]             : AIC=3836.107, Time=0.21 sec
 ARIMA(0,1,2)(0,1,1)[12]             : AIC=3834.320, Time=0.12 sec
 ARIMA(2,1,0)(0,1,1)[12]             : AIC=3834.277, Time=0.11 sec
 ARIMA(2,1,2)(0,1,1)[12]             : AIC=3837.435, Time=0.35 sec
 ARIMA(1,1,1)(0,1,1)[12] intercept   : AIC=3835.455, Time=0.15 sec

Best model:  ARIMA(1,1,1)(0,1,1)[12]          
Total fit time: 5.804 seconds
# Captura de los resultados de auto_arima en un DataFrame
# ==============================================================================
warnings.filterwarnings("ignore", message=".*force_all_finite.*", category=FutureWarning)
buffer = StringIO()
with contextlib.redirect_stdout(buffer):
    auto_arima(
            y                 = datos.loc[:fin_val],
            start_p           = 0,
            start_q           = 0,
            max_p             = 3,
            max_q             = 3,
            seasonal          = True,
            test              = 'adf',
            m                 = 12,  # periodicidad de la estacionalidad
            d                 = None, # el algoritmo determina 'd'
            D                 = None, # el algoritmo determina 'D'
            trace             = True,
            error_action      = 'ignore',
            suppress_warnings = True,
            stepwise          = True
        )
trace_autoarima = buffer.getvalue()
pattern = r"ARIMA\((\d+),(\d+),(\d+)\)\((\d+),(\d+),(\d+)\)\[(\d+)\]\s+(intercept)?\s+:\s+AIC=([\d\.]+), Time=([\d\.]+) sec"
matches = re.findall(pattern, trace_autoarima)
results = pd.DataFrame(
    matches, columns=["p", "d", "q", "P", "D", "Q", "m", "intercept", "AIC", "Time"]
)
results["order"] = results[["p", "d", "q"]].apply(
    lambda x: f"({x.iloc[0]},{x.iloc[1]},{x.iloc[2]})", axis=1
)
results["seasonal_order"] = results[["P", "D", "Q", "m"]].apply(
    lambda x: f"({x.iloc[0]},{x.iloc[1]},{x.iloc[2]},{x.iloc[3]})", axis=1
)
results = results[["order", "seasonal_order", "intercept", "AIC", "Time"]]
results.sort_values(by="AIC").reset_index(drop=True)
order seasonal_order intercept AIC Time
0 (1,1,1) (0,1,1,12) 3834.135 0.14
1 (2,1,0) (0,1,1,12) 3834.277 0.11
2 (0,1,2) (0,1,1,12) 3834.320 0.11
3 (1,1,1) (0,1,2,12) 3835.160 0.72
4 (1,1,1) (0,1,1,12) intercept 3835.455 0.15
5 (1,1,1) (1,1,1,12) 3835.564 0.14
6 (2,1,1) (0,1,1,12) 3836.104 0.30
7 (1,1,2) (0,1,1,12) 3836.107 0.21
8 (1,1,1) (1,1,2,12) 3836.443 1.08
9 (1,1,0) (0,1,1,12) 3836.696 0.08
10 (2,1,2) (0,1,1,12) 3837.435 0.36
11 (0,1,1) (0,1,1,12) 3840.318 0.09
12 (0,1,1) (0,1,2,12) 3841.572 0.69
13 (0,1,1) (1,1,1,12) 3841.882 0.14
14 (0,1,1) (1,1,2,12) 3842.593 1.03
15 (1,1,1) (1,1,0,12) 3844.410 0.09
16 (1,1,0) (1,1,0,12) 3846.786 0.06
17 (0,1,1) (1,1,0,12) 3852.231 0.09
18 (1,1,1) (0,1,0,12) 3866.187 0.03
19 (0,1,1) (0,1,0,12) 3873.797 0.03
20 (0,1,0) (1,1,1,12) 3903.204 0.07
21 (0,1,0) (0,1,1,12) 3904.615 0.05
22 (0,1,0) (0,1,0,12) 3942.897 0.01 
# Predicciones de backtesting con el mejor modelo según el grid search
# ==============================================================================
forecaster = ForecasterSarimax(
                 regressor=Sarimax(order=(0, 1, 1), seasonal_order=(1, 1, 1, 12), maxiter=500),
             )
cv = TimeSeriesFold(
        steps              = 12,
        initial_train_size = len(datos.loc[:fin_val]),
        refit              = True,
)
metrica_m1, predicciones_m1 = backtesting_sarimax(
                                forecaster            = forecaster,
                                y                     = datos,
                                cv                    = cv,
                                metric                = 'mean_absolute_error',
                                suppress_warnings_fit = True,
                              )

# Predicciones de backtesting con el mejor modelo según auto arima
# ==============================================================================
forecaster = ForecasterSarimax(
                 regressor=Sarimax(order=(1, 1, 1), seasonal_order=(0, 1, 1, 12), maxiter=500),
             )
metrica_m2, predicciones_m2 = backtesting_sarimax(
                                forecaster            = forecaster,
                                y                     = datos,
                                cv                    = cv,
                                metric                = 'mean_absolute_error',
                                suppress_warnings_fit = True
                              )
  0%|          | 0/6 [00:00<?, ?it/s] 33%|███▎      | 2/6 [00:00<00:00, 11.15it/s] 67%|██████▋   | 4/6 [00:00<00:00,  7.61it/s] 83%|████████▎ | 5/6 [00:00<00:00,  5.18it/s]100%|██████████| 6/6 [00:00<00:00,  5.89it/s]100%|██████████| 6/6 [00:00<00:00,  6.28it/s]
  0%|          | 0/6 [00:00<?, ?it/s] 33%|███▎      | 2/6 [00:00<00:00, 12.08it/s] 67%|██████▋   | 4/6 [00:00<00:00,  8.93it/s] 83%|████████▎ | 5/6 [00:00<00:00,  8.29it/s]100%|██████████| 6/6 [00:00<00:00,  7.43it/s]100%|██████████| 6/6 [00:00<00:00,  8.10it/s]
# Comparación de métricas
# ==============================================================================
print("Metrica (mean absolute error) del modelo grid search:")
metrica_m1
print("Metric (mean_absolute_error) del modelo auto arima:")
metrica_m2

fig, ax = plt.subplots(figsize=(8, 4))
datos.loc[fin_val:].plot(ax=ax, label='test')
predicciones_m1 = predicciones_m1.rename(columns={'pred': 'grid search'})
predicciones_m2 = predicciones_m2.rename(columns={'pred': 'autoarima'})
predicciones_m1.plot(ax=ax)
predicciones_m2.plot(ax=ax)
ax.set_title('Predicciones de backtesting con un modelo SARIMA')
ax.legend();
plt.show()
Metrica (mean absolute error) del modelo grid search:
Metric (mean_absolute_error) del modelo auto arima:

Variables exógenas

La implementación de ARIMA-SARIMAX ofrece statsmodels tiene una característica valiosa: la capacidad de integrar variables exógenas junto con la serie temporal principal que se está considerando. El único requisito para incluir una variable exógena es conocer el valor de la variable también durante el período de predicción. La adición de variables exógenas se realiza utilizando el argumento exog.

Utilizar un ARIMA-SARIMAX ya entrenado

Realizar predicciones con un modelo ARIMA se complica cuando los datos del horizonte de previsión no siguen inmediatamente al último valor observado durante la fase de entrenamiento. Esta complejidad se debe al componente de media móvil (MA), que utiliza como predictores los en errores de las prediciones anteriores. Por lo tanto, para predecir en el tiempo tt, se necesita el error de la predicción en t−1t−1. Si esta predicción no está disponible, el error correspondiente permanece inaccesible. Por esta razón, en la mayoría de los casos, los modelos ARIMA se vuelven a entrenar cada vez que se necesitan hacer predicciones.

# División datos Train - Last window
# ==============================================================================
fin_train = '1980-01-01 23:59:59'
                      
print(
    f"Fechas entrenamiento : {datos.index.min()} --- {datos.loc[:fin_train].index.max()}  "
    f"(n={len(datos.loc[:fin_train])})"
)
print(
    f"Fechas Last window  : {datos.loc[fin_train:].index.min()} --- {datos.index.max()}  "
    f"(n={len(datos.loc[fin_train:])})"
)

# Gráfico
# ======================================================================================
fig, ax = plt.subplots(figsize=(8, 4))
datos.loc[:fin_train].plot(ax=ax, label='entrenamiento')
datos.loc[fin_train:].plot(ax=ax, label='last window')
ax.set_title('Consumo mensual combustible España')
ax.legend();
plt.show()
Fechas entrenamiento : 1969-01-01 00:00:00 --- 1980-01-01 00:00:00  (n=133)
Fechas Last window  : 1980-02-01 00:00:00 --- 1990-01-01 00:00:00  (n=120)

# Entrenar modelo con datos desde 1969-01-01 hasta 1980-01-01
# ==============================================================================
forecaster = ForecasterSarimax(
                regressor = Sarimax(
                    order          = (0, 1, 1),
                    seasonal_order = (1, 1, 1, 12),
                    maxiter        = 500
                )
)
forecaster.fit(y=datos.loc[:fin_train])
# Predicción utilizando last window
# ==============================================================================
predicciones = forecaster.predict(
                  steps       = 12,
                  last_window = datos.loc[fin_train:]
              )
predicciones.head(3)
1990-02-01    580893.320803
1990-03-01    693624.449212
1990-04-01    654315.472036
Freq: MS, Name: pred, dtype: float64
# Gráfico predicciones
# ======================================================================================
fig, ax = plt.subplots(figsize=(8, 4))
datos.loc[:fin_train].plot(ax=ax, label='entrenamiento')
datos.loc[fin_train:].plot(ax=ax, label='last window')
predicciones.plot(ax=ax, label='predicciones')
ax.set_title('Consumo mensual combustible España')
ax.legend();
plt.show()