Más sobre forecasting en: cienciadedatos.net

• Forecasting con machine learning
• ARIMA y SARIMAX
• Forecasting series temporales con gradient boosting
• Modelos de forecasting globales I
• Modelos de forecasting globales II
• Modelos de forecasting globales III
• Modelos de forecasting globales IV
• Modelos de forecasting globales V
• Clustering de series temporales
• Forecasting de la demanda eléctrica
• Predicción de las visitas a una web
• Predicción del precio de Bitcoin
• Forecasting probabilístico I
• Forecasting probabilístico II
• Forecasting de series incompletas
• Forecasting de demanda intermitente
• Forecasting con modelos de deep learning
• Reducir el impacto del covid-19 en modelos de forecasting
• Interpretabilidad en modelos de forecasting
• Forecasting tendencia con modelos de árboles
• Stacking ensemble de modelos de forecasting
• Detección de anomalías en series temporales
• Data leakage en modelos de forecasting preentrenados
• Acelerar modelos de forecasting con GPU

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 se han aplicado diferencias entre valores consecutivos de la serie original (y este proceso de diferenciación 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.

En la notación del modelo ARIMA-SARIMAX, los parámetros $p$, $d$, y $q$ representan las componentes autorregresivas, de diferenciación y de media móvil, respectivamente. $P$, $D$, y $Q$ son las mismas componentes para la parte estacional del modelo y $m$ el número de períodos en cada temporada.

• $p$ es el orden (número de lags temporales) de la parte autorregresiva del modelo.

• $d$ es el grado de diferenciación (el número de veces que se han restado los valores consecutivos de la serie).

• $q$ es el tamaño de la media móvil del modelo.

• $P$ es el orden (número de lags temporales) de la parte estacional del modelo.

• $D$ es el grado de diferenciación de la parte estacional del modelo.

• $Q$ es el tamaño de la media móvil de la parte estacional del modelo.

• $m$ indica al número de períodos en cada temporada.

Cuando los términos $P$, $D$, $Q$, y $m$ son cero y no se incluyen variables exógenas, el modelo SARIMAX es equivalente a un ARIMA.

La librería de python skforecast proporciona dos clases que permiten la creación de modelos de la familia ARIMA:

• skforecast Sarimax: Un wrapper del conocido modelo SARIMAX de la librería statsmodels que sigue la API de scikit-learn. Similar a pmdarima, esta versión ha sido depurada para incluir únicamente los elementos esenciales requeridos por skforecast, lo que se traduce en mejoras significativas en el rendimiento.

• skforecast Arima: Una implementación nativa que también sigue la API de scikit-learn. Esta versión está optimizada para la velocidad mediante la compilación JIT de Numba. Soporta componentes estacionales, variables exógenas y búsqueda automática de parámetros, lo que la convierte en una herramienta versátil para toda la familia ARIMA: ARIMA, SARIMAX y AutoARIMA.

Este documento explora estas librerías y explica cómo construir modelos ARIMA-SARIMAX utilizando cada una de ellas. Además, muestra cómo la clase ForecasterStats amplía las capacidades de los modelos ARIMA-SARIMAX al incorporar funcionalidades como backtesting, ajuste de hiperparámetros, forecasting probabilístico y más. Esta integración permite a los usuarios aprovechar las fortalezas de los modelos ARIMA-SARIMAX mientras se benefician de las características avanzadas proporcionadas por la librería skforecast.

Librerías¶

Librerías utilizadas en este documento.

# Librerías
# ======================================================================================
import numpy as np
import pandas as pd
from io import StringIO
import contextlib
import re
import matplotlib.pyplot as plt

# 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, plot_prediction_intervals
from skforecast.stats import Sarimax, Arima
from skforecast.recursive import ForecasterStats
from skforecast.utils import expand_index
from skforecast.model_selection import TimeSeriesFold, backtesting_stats, grid_search_stats

import warnings
warnings.filterwarnings('once')
warnings.filterwarnings('ignore', category=DeprecationWarning)

color = '\033[1m\033[38;5;208m' 
print(f"{color}Versión skforecast: {skforecast.__version__}")
print(f"{color}Versión statsmodels: {statsmodels.__version__}")
print(f"{color}Versión pandas: {pd.__version__}")
print(f"{color}Versión numpy: {np.__version__}")

Versión skforecast: 0.20.0
Versión statsmodels: 0.14.6
Versión pandas: 2.3.3
Versión numpy: 2.2.6

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']
display(datos.head(4))

╭──────────────────────────────── fuel_consumption ────────────────────────────────╮
│ Description:                                                                     │
│ Monthly fuel consumption in Spain from 1969-01-01 to 2022-08-01.                 │
│                                                                                  │
│ Source:                                                                          │
│ 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            │
│                                                                                  │
│ URL:                                                                             │
│ https://raw.githubusercontent.com/skforecast/skforecast-                         │
│ datasets/main/data/consumos-combustibles-mensual.csv                             │
│                                                                                  │
│ Shape: 644 rows x 6 columns                                                      │
╰──────────────────────────────────────────────────────────────────────────────────╯

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
# ======================================================================================
fin_train = '1983-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:]

# Gráfico
# ======================================================================================
set_dark_theme()
fig, ax=plt.subplots(figsize=(7, 3))
datos_train.plot(ax=ax, label='train')
datos_test.plot(ax=ax, label='test')
ax.set_title('Consumo mensual combustible España')
ax.legend();

Fechas train : 1969-01-01 00:00:00 --- 1983-01-01 00:00:00  (n=169)
Fechas test  : 1983-02-01 00:00:00 --- 1990-01-01 00:00:00  (n=84)

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 si 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 $d$.

• 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 ($p$ y $q$) 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 o 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 $P$, $D$, $Q$ y $m$.

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:

Donde $X_t$ es el valor en el tiempo $t$ y $X_{t-1}$ es el valor en el tiempo $t−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 ($H_O$): La serie tiene una raíz unitaria, no es estacionaria.

• Hipótesis alternativa ($H_A$): 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
# ==============================================================================
datos_diff_1 = datos_train.diff().dropna()
datos_diff_2 = datos_diff_1.diff().dropna()

# Suppress warnings
with warnings.catch_warnings():
    warnings.filterwarnings('ignore')
    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]}')

# Gráfico series
# ==============================================================================
fig, axs = plt.subplots(nrows=3, ncols=1, figsize=(7, 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');

Test estacionariedad serie original
-------------------------------------
ADF Statistic: -0.44612980998227997, p-value: 0.9021071923942665
KPSS Statistic: 2.2096370946978383, p-value: 0.01

Test estacionariedad para serie diferenciada (order=1)
--------------------------------------------------
ADF Statistic: -2.7436320757907513, p-value: 0.06678751673391206
KPSS Statistic: 0.313271162357279, p-value: 0.1

Test estacionariedad para serie diferenciada (order=2)
--------------------------------------------------
ADF Statistic: -10.490170523414776, p-value: 1.1524989419523892e-18
KPSS Statistic: 0.08065668267482215, 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 $d$ 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 $p$ y $q$. La ACF ayuda a identificar el valor de $q$ (retardos en la parte de media móvil), mientras que la PACF ayuda a identificar el valor de $p$ (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 $q$. 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 $p$.

# Grafico de autocorrelación para la serie original y la serie diferenciada
# ==============================================================================
fig, axs = plt.subplots(nrows=2, ncols=1, figsize=(6, 4), 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)');

# Autocorrelación parcial para la serie original y la serie diferenciada
# ==============================================================================
fig, axs = plt.subplots(nrows=2, ncols=1, figsize=(6, 3), 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();

Acorde a la función de autocorrelación parcial, el valor óptimo para el parámetro $p$ es 0. Sin embargo, se va a asignar un valor de 1 para proporcionar un componente autorregresivo al modelo. En cuanto al componente $q$, 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=(9, 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();

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.86234654148435, p-value: 4.126322028795106e-05
KPSS Statistic: 0.29553731751033135, 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()

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'
display(predicciones_statsmodels.head(4))

1983-02-01    408707.295485
1983-03-01    465933.670647
1983-04-01    506906.544811
1983-05-01    463579.002354
Freq: MS, Name: predicciones_statsmodels, dtype: float64

Skforecast Sarimax¶

La clase Sarimax de 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))
modelo.fit(y=datos_train)
modelo.summary()

# Predictión
# ==============================================================================
predicciones_skforecast_sarimax = modelo.predict(steps=len(datos_test))
predicciones_skforecast_sarimax = pd.DataFrame(predicciones_skforecast_sarimax, index=datos_test.index)
predicciones_skforecast_sarimax.columns = ['skforecast']
display(predicciones_skforecast_sarimax.head(4))

Skforecast Arima¶

Desde la versión 0.20.0 de skforecast, la implementación nativa de modelos ARIMA está disponible a través de la clase Arima. Esta implementación está optimizada para la velocidad mediante la compilación JIT de Numba. Soporta componentes estacionales, variables exógenas y búsqueda automática de parámetros, lo que la convierte en una herramienta versátil para toda la familia ARIMA: ARIMA, SARIMAX y AutoARIMA.

# Modelo ARIMA con skforecast.Arima
# ==============================================================================
modelo = Arima(order=(1, 1, 1), seasonal_order=(1, 1, 1), m=12)
modelo.fit(y=datos_train, suppress_warnings=True)
modelo.summary()

ARIMA Model Summary
============================================================
Model     : Arima(1,1,1)(1,1,1)[12]
Method    : ARIMA(1,1,1)(1,1,1)[12]
Converged : False

Coefficients:
------------------------------------------------------------
  ar1            :    -0.3083  (SE:   0.1137, t:    -2.71)
  ma1            :    -0.5583  (SE:   0.1093, t:    -5.11)
  sar1           :    -0.1390  (SE:   0.1160, t:    -1.20)
  sma1           :    -0.5581  (SE:   0.1010, t:    -5.53)

Model fit statistics:
  sigma^2:             249194557.351875
  Log-likelihood:      -1733.55
  AIC:                 3477.11
  BIC:                 N/A

Residual statistics:
  Mean:                -1167.159495
  Std Dev:             15166.561117
  MAE:                 10860.580394
  RMSE:                15166.599660

Time Series Summary Statistics:
Number of observations: 169
  Mean:                 384743.1773
  Std Dev:              108126.6689
  Min:                  155466.8105
  25%:                  303667.7591
  Median:               397278.0241
  75%:                  466194.3073
  Max:                  605073.0143

# Predictcion
# ==============================================================================
predicciones_skforecast_arima = modelo.predict(steps=len(datos_test))
pred_index = expand_index(index=datos_train.index, steps=len(datos_test))
predicciones_skforecast_arima = pd.Series(predicciones_skforecast_arima, index=pred_index)
predicciones_skforecast_arima.head(4)

1983-02-01    414750.833100
1983-03-01    465097.100031
1983-04-01    504456.783383
1983-05-01    472048.144760
Freq: MS, dtype: float64

Comparación de predicciones¶

# Gráfico de las predicciones
# ==============================================================================
fig, ax = plt.subplots(figsize=(7, 3))
datos_train.plot(ax=ax, label='train')
datos_test.plot(ax=ax, label='test')
predicciones_statsmodels.plot(ax=ax, label='statsmodels')
predicciones_skforecast_sarimax.plot(ax=ax, label='skforecast sarimax')
predicciones_skforecast_arima.plot(ax=ax, label='skforecast arima')
ax.set_title('Predictions with ARIMA models')
ax.legend();

ForecasterStats¶

La clase ForecasterStats amplía las capacidades de los modelos ARIMA-SARIMAX al incorporar funcionalidades como backtesting, ajuste de hiperparámetros, forecasting probabilístico y más. Esta integración permite a los usuarios aprovechar las fortalezas de los modelos ARIMA-SARIMAX mientras se benefician de las características avanzadas proporcionadas por la librería skforecast.

Entrenamiento¶

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

# ForecasterStats con Arima
# ==============================================================================
forecaster = ForecasterStats(
                 estimator=Arima(order=(1, 1, 1), seasonal_order=(1, 1, 1), m=12)
             )
forecaster.fit(y=datos_train, suppress_warnings=True)
forecaster

Predicción¶

Una vez que el modelo está ajustado, puede utilizarse para predecir observaciones futuras. Es importante tener en cuenta que este tipo de modelos requiere que las predicciones sigan inmediatamente después de los datos de entrenamiento; por lo tanto, el pronóstico comienza justo después del último valor observado.

# Predicciones
# ==============================================================================
steps = len(datos_test)
predicciones = forecaster.predict(steps=steps)
predicciones

1983-02-01    414750.833100
1983-03-01    465097.100031
1983-04-01    504456.783383
1983-05-01    472048.144760
1983-06-01    494874.801161
                  ...      
1989-09-01    533179.537771
1989-10-01    520348.872652
1989-11-01    472951.496372
1989-12-01    510776.528620
1990-01-01    465935.769505
Freq: MS, Name: pred, Length: 84, dtype: float64

Intervalos de predicción¶

El método predict_interval permite calcular intervalos de predicción. Los usuarios pueden especificar el nivel de confianza del intervalo estimado utilizando el argumento alpha o interval.

# Intervalos de predicción
# ==============================================================================
predicciones_intervalo = forecaster.predict_interval(steps=steps, alpha=0.05)
predicciones_intervalo.head(3)

# Gráfico de los intervalos de predicción
# ==============================================================================
fig, ax = plt.subplots(figsize=(6, 3))
plot_prediction_intervals(
    predictions         = predicciones_intervalo,
    y_true              = datos_test,
    target_variable     = "litters",
    title               = "Prediction intervals in test data",
    kwargs_fill_between = {'color': 'white', 'alpha': 0.3, 'zorder': 1},
    ax                  = ax
)

Importancia de predictores¶

El método get_feature_importances devuelve los coeficientes del modelo, que pueden utilizarse como medida de la importancia de los predictores en el modelo SARIMAX o ARIMA.

# Importancia de predictores
# ==============================================================================
modelo.get_feature_importances()

Predicciones sobre los datos de entrenamiento (Predicciones In-sample)¶

Predicciones sobre los datos de entrenamiento son cruciales para evaluar la precisión y efectividad del modelo. Al comparar los valores predichos con los valores observados reales en el conjunto de datos de entrenamiento, se puede evaluar qué tan bien el modelo ha aprendido los patrones y tendencias subyacentes en los datos. Esta comparación ayuda a comprender el rendimiento del modelo e identificar áreas donde puede necesitar mejoras o ajustes. En esencia, actúan como un espejo, reflejando cómo el modelo interpreta y reconstruye los datos históricos sobre los que fue entrenado.

Las predicciones de las observaciones utilizadas para ajustar el modelo se almacenan en el atributo fitted_values_ del objeto Arima.

# Predicciones In-sample
# ==============================================================================
# Mostrar solo los primeros 5 valores
modelo.fitted_values_[:5]

array([166778.86748562, 155441.28198886, 184949.38956579, 202282.72293861,
       206228.03030826])

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.

# Crear forecaster
# ==============================================================================
forecaster = ForecasterStats(
                 estimator=Arima(order=(1, 1, 1), seasonal_order=(1, 1, 1), m=12)
             )

# Backtest forecaster
# ==============================================================================
cv = TimeSeriesFold(
        steps              = 12,
        initial_train_size = len(datos_train),
        refit              = True,
        fixed_train_size   = False,
)
metricas, predicciones = backtesting_stats(
                            forecaster        = forecaster,
                            y                 = datos,
                            cv                = cv,
                            metric            = ['mean_absolute_error', 'mean_absolute_percentage_error'],
                            suppress_warnings = True,
                            verbose           = True,
                            show_progress     = True
                        )

Information of folds
--------------------
Number of observations used for initial training: 169
Number of observations used for backtesting: 84
    Number of folds: 7
    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 -- 1983-01-01 00:00:00  (n=169)
    Validation: 1983-02-01 00:00:00 -- 1984-01-01 00:00:00  (n=12)
Fold: 1
    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: 2
    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: 3
    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: 4
    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: 5
    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: 6
    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/7 [00:00<?, ?it/s]

# Predicciones de backtesting
# ==============================================================================
predicciones.head(4)

# Métricas de backtesting
# ==============================================================================
metricas

# Gráfico predicciones de backtesting
# ==============================================================================
fig, ax = plt.subplots(figsize=(6, 3))
datos.loc[fin_train:].plot(ax=ax, label='test')
predicciones['pred'].plot(ax=ax)
ax.set_title('Predicciones de backtesting (particiones de 12 meses)')
ax.legend();

Búsqueda 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. Entre estos métodos, dos enfoques ampliamente utilizados son:

• Criterios estadísticos: Las métricas de criterios de información, como el Criterio de Información de Akaike (AIC) o el Criterio de Información Bayesiano (BIC), utilizan diferentes penalizaciones sobre la estimación de máxima verosimilitud del modelo como medida de la bondad de ajuste. La ventaja de utilizar estas métricas es que se calculan únicamente con los datos de entrenamiento, lo que elimina la necesidad de realizar predicciones sobre nuevos datos. Como resultado, el proceso de optimización se acelera considerablemente. El conocido algoritmo Auto Arima utiliza este enfoque.

• Técnicas de validación: El uso de técnicas de validación, especialmente el backtesting, es otra estrategia efectiva. El backtesting consiste en evaluar el rendimiento del modelo utilizando datos históricos para simular las condiciones del mundo real. Esto ayuda a validar la eficacia de los hiperparámetros en diferentes escenarios, proporcionando una evaluación práctica de su viabilidad.

En el primer enfoque, los cálculos se basan únicamente en los datos de entrenamiento, lo que elimina la necesidad separar un partición de datos de validación. Esto hace que el proceso de optimización sea muy rápido. Sin embargo, es importante señalar que las métricas basadas en los criterios de información sólo miden la calidad relativa de los modelos. Esto significa que los modelos probados podrían ser deficientes. Por lo tanto, el modelo final seleccionado debe someterse a una fase de backtesting. En esta fase se calcula una métrica (como MAE, MSE, MAPE, etc.) que valida su rendimiento en una escala con sentido para el caso de uso.

Por otro lado, el segundo enfoque -las técnicas de validación- suele requerir más tiempo, ya que el modelo debe entrenarse y luego evaluarse con nuevos datos. Sin embargo, los resultados generados suelen ser más robustos y las métricas derivadas pueden proporcionar información más profunda.

AutoARIMA¶

AutoArima es un algoritmo diseñado para automatizar la selección de los hiperparámetros óptimos de un modelo ARIMA. El algoritmo evalúa de forma sistemática distintas combinaciones de parámetros no estacionales ((p, d, q)), parámetros estacionales ((P, D, Q)) y el período estacional ((m)) para identificar la configuración que mejor se ajusta a los datos según un criterio específico, normalmente el Criterio de Información de Akaike (AIC) o el Criterio de Información Bayesiano (BIC).

Con esta estrategia, los cálculos se basan únicamente en los datos de entrenamiento, lo que elimina la necesidad de una partición adicional de los datos, requerida por otras estrategias como grid search. Esto hace que el proceso de optimización sea extremadamente eficiente. No obstante, es importante tener en cuenta que los criterios de información solo miden la calidad relativa de los modelos dentro del espacio de búsqueda definido. Un modelo con el AIC más bajo puede seguir siendo un mal ajuste en términos absolutos. Por ello, el modelo seleccionado debería idealmente someterse a una fase de backtesting. En esta fase se calculan métricas de error de predicción (como MAE, MSE o MAPE) para validar el rendimiento en una escala significativa e interpretable.

La clase Arima de Skforecast activa la funcionalidad AutoArima cuando los argumentos order o seasonal_order se establecen en None. Una vez ajustado el modelo, los parámetros óptimos pueden consultarse a través del atributo best_params_. Para todas las predicciones posteriores, el modelo utiliza automáticamente la configuración óptima identificada durante el proceso de ajuste.

# Skforecast Auto Arima
# ==============================================================================
auto_arima = Arima(
    order          = None,  # Must be None to use AutoArima
    seasonal_order = None,  # Must be None to use AutoArima
    start_p        = 0,
    start_q        = 0,
    max_p          = 5,
    max_q          = 5,
    max_P          = 5,
    max_Q          = 2,
    max_order      = 5,
    max_d          = 2,
    max_D          = 1,
    ic             = "aic",
    m              = 12,
    stepwise       = True,  # True for faster results
    trace          = True,  # True for detailed information of the process
)
auto_arima.fit(y=datos_train, suppress_warnings=True)
auto_arima

Fitting models using approximations...

 ARIMA(p,d,q)(P,D,Q)[m]                     : aic
 ------------------------------------------------
 ARIMA(0,1,0)(1,1,1)[12]                    : 3564.3887
 ARIMA(0,1,0)(0,1,0)[12]                    : 3637.1363
 ARIMA(1,1,0)(1,1,0)[12]                    : 3506.9117
 ARIMA(0,1,1)(0,1,1)[12]                    : 3481.3432
 ARIMA(0,1,1)(0,1,0)[12]                    : 3542.1014
 ARIMA(0,1,1)(1,1,1)[12]                    : 3482.0043
 ARIMA(0,1,1)(0,1,2)[12]                    : 3481.2735
 ARIMA(0,1,1)(1,1,2)[12]                    : 3482.5762
 ARIMA(0,1,0)(0,1,2)[12]                    : 3563.7959
 ARIMA(1,1,1)(0,1,2)[12]                    : 3476.4900
 ARIMA(1,1,1)(0,1,1)[12]                    : 3476.4887
 ARIMA(1,1,1)(0,1,0)[12]                    : 3531.8215
 ARIMA(1,1,1)(1,1,1)[12]                    : 3477.1082
 ARIMA(1,1,1)(1,1,0)[12]                    : 3493.3277
 ARIMA(1,1,1)(1,1,2)[12]                    : 3477.9947
 ARIMA(1,1,0)(0,1,1)[12]                    : 3491.3027
 ARIMA(2,1,1)(0,1,1)[12]                    : 3477.5383
 ARIMA(1,1,2)(0,1,1)[12]                    : 3477.8528
 ARIMA(0,1,0)(0,1,1)[12]                    : 3563.8825
 ARIMA(0,1,2)(0,1,1)[12]                    : 3478.2811
 ARIMA(2,1,0)(0,1,1)[12]                    : 3482.5359
 ARIMA(2,1,2)(0,1,1)[12]                    : 3479.4579

Now re-fitting the best model(s) without approximations...
 ARIMA(1,1,1)(0,1,1)[12]                    : 3476.4887

Best model found: ARIMA(1,1,1)(0,1,1)[12] with aic: 3476.488679996082

AutoArima(1,1,1)(0,1,1)[12]

# Parametros óptimos AutoARIMA
# ==============================================================================
print(f"Best parameters: {auto_arima.best_params_}")

Best parameters: {'order': (1, 1, 1), 'seasonal_order': (0, 1, 1), 'm': 12}

Puede ser de interés capturar la traza generada por la función auto_arima para permitir una exploración más exhaustiva de los resultados. La implementación actual imprime los resultados, pero es posible capturarlos y almacenarlos en un marco de datos estructurado de Pandas.

# Capture auto_arima trace in a pandas dataframe
# ==============================================================================
buffer = StringIO()
with contextlib.redirect_stdout(buffer):
    auto_arima = Arima(
        order          = None,  # Must be None to use AutoArima
        seasonal_order = None,  # Must be None to use AutoArima
        start_p        = 0,
        start_q        = 0,
        max_p          = 5,
        max_q          = 5,
        max_P          = 5,
        max_Q          = 2,
        max_order      = 5,
        max_d          = 2,
        max_D          = 1,
        ic             = "aic",
        m              = 12,
        stepwise       = True,  # True for faster results
        trace          = True,  # True for detailed information of the process
    )
    auto_arima.fit(y=datos_train, suppress_warnings=True)

trace_autoarima = buffer.getvalue()
# Pattern matches: ARIMA(p,d,q)(P,D,Q)[m] : aic_value
pattern = r"ARIMA\((\d+),(\d+),(\d+)\)\((\d+),(\d+),(\d+)\)\[(\d+)\]\s+:\s+([\d\.]+)"
matches = re.findall(pattern, trace_autoarima)

results = pd.DataFrame(
    matches, columns=["p", "d", "q", "P", "D", "Q", "m", "AIC"]
)
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", "AIC"]]
results["AIC"] = results["AIC"].astype(float)
results = results.sort_values(by="AIC").drop_duplicates().reset_index(drop=True)
results.head()

Grid search con backtesting¶

Aunque es menos común en el contexto de modelos estadísticos, también es posible utilizar la búsqueda en cuadrícula (grid search) y la búsqueda aleatoria (random search) para encontrar hiperparámetros óptimos. Es crucial realizar la búsqueda utilizando un conjunto de datos de validación, en lugar del conjunto de datos de prueba, para asegurar una evaluación precisa e imparcial del rendimiento del modelo.

La partición de entrenamiento se divide nuevamente para separar una partición de validación, que se utiliza para la búsqueda de hiperparámetros. La partición de prueba permanece intacta y se reserva exclusivamente para la evaluación final del rendimiento del modelo.

# Split Train-validation-test data
# ======================================================================================
end_train = '1979-01-01 23:59:59'
end_val   = '1983-01-01 23:59:59'

print(
    f"Fechas train      : {datos.index.min()} --- {datos.loc[:end_train].index.max()}  "
    f"(n={len(datos.loc[:end_train])})"
)
print(
    f"Fechas validación : {datos.loc[end_train:].index.min()} --- {datos.loc[:end_val].index.max()}  "
    f"(n={len(datos.loc[end_train:end_val])})"
)
print(
    f"Fechas test      : {datos.loc[end_val:].index.min()} --- {datos.index.max()}  "
    f"(n={len(datos.loc[end_val:])})"
)

# Plot
# ==============================================================================
fig, ax = plt.subplots(figsize=(7, 3))
datos.loc[:end_train].plot(ax=ax, label='train')
datos.loc[end_train:end_val].plot(ax=ax, label='validación')
datos.loc[end_val:].plot(ax=ax, label='test')
ax.set_title('Consumo mensual de combustible en España')
ax.legend();

Fechas train      : 1969-01-01 00:00:00 --- 1979-01-01 00:00:00  (n=121)
Fechas validación : 1979-02-01 00:00:00 --- 1983-01-01 00:00:00  (n=48)
Fechas test      : 1983-02-01 00:00:00 --- 1990-01-01 00:00:00  (n=84)

# Grid search basado en backtesting
# ==============================================================================
forecaster = ForecasterStats(
                 estimator=Arima(order=(1, 0, 0))  # Este valor se reemplaza en el grid search
             )

param_grid = {
    'order': [(0, 1, 0), (0, 1, 1), (1, 1, 0), (1, 1, 1), (2, 1, 1)],
    'seasonal_order': [(0, 0, 0,), (0, 1, 0), (1, 1, 1)],
    'm': [1, 12]
}

cv = TimeSeriesFold(steps=12, initial_train_size=len(datos.loc[:end_train]))

resultados_grid = grid_search_stats(
                   forecaster        = forecaster,
                   y                 = datos.loc[:end_val],
                   param_grid        = param_grid,
                   cv                = cv,
                   metric            = 'mean_absolute_error',
                   return_best       = True,
                   n_jobs            = 'auto',
                   suppress_warnings = True,
                   verbose           = False,
                   show_progress     = True
                 )
resultados_grid.head(5)

Number of models compared: 30.

params grid:   0%|          | 0/30 [00:00<?, ?it/s]

`Forecaster` refitted using the best-found parameters, and the whole data set: 
  Parameters: {'m': 12, 'order': (0, 1, 0), 'seasonal_order': (1, 1, 1)}
  Backtesting metric: 17746.910354365504

Se comparan los dos modelos candidatos: el seleccionado por grid search basado en backtesting con un error absoluto medio, y el seleccionado por auto arima basado en el AIC, al realizar predicciones para los próximos tres años en intervalos de 12 meses.

# Predicciones de backtesting con el mejor modelo según el grid search
# ==============================================================================
forecaster = ForecasterStats(
                 estimator=Arima(order=(0, 1, 0), seasonal_order=(1, 1, 1), m=12),
             )
cv = TimeSeriesFold(
        steps              = 12,
        initial_train_size = len(datos.loc[:end_val]),
        refit              = True,
)
metrica_m1, predicciones_m1 = backtesting_stats(
                                forecaster        = forecaster,
                                y                 = datos,
                                cv                = cv,
                                metric            = 'mean_absolute_error',
                                suppress_warnings = True,
                              )

# Predicciones de backtesting con el mejor modelo según auto arima
# ==============================================================================
forecaster = ForecasterStats(
                 estimator=Arima(order=(1, 1, 1), seasonal_order=(0, 1, 1), m=12),
             )
metrica_m2, predicciones_m2 = backtesting_stats(
                                forecaster        = forecaster,
                                y                 = datos,
                                cv                = cv,
                                metric            = 'mean_absolute_error',
                                suppress_warnings = True
                              )

  0%|          | 0/7 [00:00<?, ?it/s]

  0%|          | 0/7 [00:00<?, ?it/s]

# Comparación de métricas
# ==============================================================================
print("Metrica (mean absolute error) del modelo grid search:")
display(metrica_m1)
print("Metric (mean_absolute_error) del modelo auto arima:")
display(metrica_m2)

fig, ax = plt.subplots(figsize=(6, 3))
datos.loc[end_val:].plot(ax=ax, label='test')
predicciones_m1 = predicciones_m1.rename(columns={'pred': 'grid search'})
predicciones_m2 = predicciones_m2.rename(columns={'pred': 'auto arima'})
predicciones_m1['grid search'].plot(ax=ax)
predicciones_m2['auto arima'].plot(ax=ax)
ax.set_title('Predicciones de backtesting')
ax.legend();

Metrica (mean absolute error) del modelo grid search:

Metric (mean_absolute_error) del modelo auto arima:

Para este conjunto de datos, la configuración identificada por la técnica Auto ARIMA logra mejores resultados.

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 $t$, se necesita el error de la predicción en $t-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.

A pesar de los considerables esfuerzos y avances para acelerar el proceso de entrenamiento de estos modelos, no siempre es factible volver a entrenar el modelo antes de cada etapa de pedicción, ya sea por limitaciones de tiempo o recursos computacionales insuficientes para acceder repetidamente a los datos históricos. Un enfoque intermedio es alimentar al modelo con datos desde la última observación de entrenamiento hasta el inicio de la fase de predicción. Esto permite estimar las predicciones intermedias y, como resultado, los errores necesarios. Por ejemplo, supóngase que un modelo se entrenó hace 20 días con datos diarios de los últimos tres años. Al generar nuevas predicciones, sólo se necesitarían los 20 valores más recientes, en lugar del conjunto de datos históricos completo (365 * 3 + 20).

Integrar nuevos datos en el modelo puede ser complejo, la clase ForecasterStats automatiza el proceso mediante el argumento last_window de su método predict.

# Split data Train - Last window - Test
# ==============================================================================
end_train = '1983-01-01 23:59:59'
end_last_window = '1988-01-01 23:59:59'

print(
    f"Fechas train       : {datos.index.min()} --- {datos.loc[:end_train].index.max()}  "
    f"(n={len(datos.loc[:end_train])})"
)
print(
    f"Fechas last window : {datos.loc[end_train:].index.min()} --- {datos.loc[:end_last_window].index.max()}  "
    f"(n={len(datos.loc[end_train:end_last_window])})"
)
print(
    f"Fechas test        : {datos.loc[end_last_window:].index.min()} --- {datos.index.max()}  "
    f"(n={len(datos.loc[end_last_window:])})"
)
datos_train       = datos.loc[:end_train]
datos_last_window = datos.loc[end_train:end_last_window]
datos_test        = datos.loc[end_last_window:]

# Plot
# ======================================================================================
fig, ax = plt.subplots(figsize=(7, 3))
datos_train.plot(ax=ax, label='train')
datos_last_window.plot(ax=ax, label='last window')
datos_test.plot(ax=ax, label='test')
ax.legend();

Fechas train       : 1969-01-01 00:00:00 --- 1983-01-01 00:00:00  (n=169)
Fechas last window : 1983-02-01 00:00:00 --- 1988-01-01 00:00:00  (n=60)
Fechas test        : 1988-02-01 00:00:00 --- 1990-01-01 00:00:00  (n=24)

El Forecaster se entrena utilizando datos hasta '1983-01-01 23:59:59' y luego utilizará la información restante como la última ventana de observaciones para generar nuevas predicciones. Esto significa que las predicciones no comenzarán inmediatamente después de la última observación de entrenamiento, sino después de la última observación en la last_window proporcionada.

# Entrenar modelo con datos desde 1969-01-01 hasta 1983-01-01 23:59:59
# ==============================================================================
forecaster = ForecasterStats(
                 estimator=Sarimax(order=(1, 1, 1), seasonal_order=(0, 1, 1, 12)),
             )
forecaster.fit(y= datos_train, suppress_warnings=True)

A continuación, se predicien los siguientes 12 valores de la serie.

# Predicción utilizando last window
# ==============================================================================
predicciones = forecaster.predict(
                  steps       = 12,
                  last_window = datos.loc[end_train:]
              )
predicciones.head(3)

1990-02-01    580151.145790
1990-03-01    692490.243788
1990-04-01    655918.187923
Freq: MS, Name: pred, dtype: float64

# Gráfico predicciones
# ======================================================================================
fig, ax = plt.subplots(figsize=(7, 3))
datos.loc[:end_train].plot(ax=ax, label='entrenamiento')
datos.loc[end_train:].plot(ax=ax, label='last window')
predicciones.plot(ax=ax, label='predicciones')
ax.set_title('Consumo mensual combustible España')
ax.legend();

Optimización de memoria¶

Cuando se trabaja con modelos ARIMA-SARIMAX en entornos de producción donde es necesario almacenar muchos modelos ajustados pero solo se requieren capacidades de predicción (no diagnósticos), se puede reducir significativamente el uso de memoria con el método reduce_memory(). Esto es especialmente útil al trabajar con conjuntos de datos grandes o al desplegar modelos en entornos con recursos limitados. Este método elimina los valores ajustados y los residuos dentro de la muestra, que solo son necesarios para fines de diagnóstico pero no para generar pronósticos.

# Comparar tamaño antes y después de reduce_memory()
# ==============================================================================
import sys

def total_model_size(model):
    size = sys.getsizeof(model)
    for attr_name in dir(model):
        if attr_name.startswith('_'):
            continue
        try:
            attr = getattr(model, attr_name)
            size += sys.getsizeof(attr)
        except Exception:
            pass
    return size

model=Arima(order=(1, 1, 1), seasonal_order=(0, 1, 1), m=12)
model.fit(y=datos_train, suppress_warnings=True)
model_size_before = total_model_size(model)
print(f"Memory before reduce_memory(): {model_size_before / 1024:.3f} KB")

# Reduce memory
model.reduce_memory()
model_size_after = total_model_size(model)
print(f"Memory after reduce_memory(): {model_size_after / 1024:.3f} KB")
print(f"Memory reduction: {(1 - model_size_after / model_size_before) * 100:.1f}%")

Memory before reduce_memory(): 6.812 KB
Memory after reduce_memory(): 3.843 KB
Memory reduction: 43.6%

Información de sesión¶

import session_info
session_info.show(html=False)

-----
matplotlib          3.10.8
numpy               2.2.6
pandas              2.3.3
session_info        v1.0.1
skforecast          0.20.0
statsmodels         0.14.6
-----
IPython             9.8.0
jupyter_client      7.4.9
jupyter_core        5.9.1
notebook            6.5.7
-----
Python 3.13.11 | packaged by Anaconda, Inc. | (main, Dec 10 2025, 21:28:48) [GCC 14.3.0]
Linux-6.14.0-37-generic-x86_64-with-glibc2.39
-----
Session information updated at 2026-02-01 20:18

Bibliografía¶

Hyndman, R.J., & Athanasopoulos, G. (2021) Forecasting: principles and practice, 3rd edition, OTexts: Melbourne, Australia.

Time Series Analysis and Forecasting with ADAM Ivan Svetunkov

Python for Finance: Mastering Data-Driven Finance

Forecasting: theory and practice

Instrucciones para citar¶

¿Cómo citar este documento?

Si utilizas este documento o alguna parte de él, te agradecemos que lo cites. ¡Muchas gracias!

Modelos ARIMA y SARIMAX con Python por Joaquín Amat Rodrigo y Javier Escobar Ortiz, disponible con licencia Attribution-NonCommercial-ShareAlike 4.0 International (CC BY-NC-SA 4.0 DEED) en https://www.cienciadedatos.net/documentos/py51-modelos-arima-sarimax-python.html

¿Cómo citar skforecast?

Si utilizas skforecast, te agradeceríamos mucho que lo cites. ¡Muchas gracias!

Zenodo:

Amat Rodrigo, Joaquin, & Escobar Ortiz, Javier. (2024). skforecast (v0.20.0). Zenodo. https://doi.org/10.5281/zenodo.8382788

APA:

Amat Rodrigo, J., & Escobar Ortiz, J. (2024). skforecast (Version 0.20.0) [Computer software]. https://doi.org/10.5281/zenodo.8382788

BibTeX:

@software{skforecast, author = {Amat Rodrigo, Joaquin and Escobar Ortiz, Javier}, title = {skforecast}, version = {0.20.0}, month = {01}, year = {2026}, license = {BSD-3-Clause}, url = {https://skforecast.org/}, doi = {10.5281/zenodo.8382788} }

¿Te ha gustado el artículo? Tu ayuda es importante

Tu contribución me ayudará a seguir generando contenido divulgativo gratuito. ¡Muchísimas gracias! 😊

Este documento creado por Joaquín Amat Rodrigo y Javier Escobar Ortiz tiene licencia Attribution-NonCommercial-ShareAlike 4.0 International.

Se permite:

• Compartir: copiar y redistribuir el material en cualquier medio o formato.

• Adaptar: remezclar, transformar y crear a partir del material.

Bajo los siguientes términos:

• Atribución: Debes otorgar el crédito adecuado, proporcionar un enlace a la licencia e indicar si se realizaron cambios. Puedes hacerlo de cualquier manera razonable, pero no de una forma que sugiera que el licenciante te respalda o respalda tu uso.

• No-Comercial: No puedes utilizar el material para fines comerciales.

• Compartir-Igual: Si remezclas, transformas o creas a partir del material, debes distribuir tus contribuciones bajo la misma licencia que el original.