# Tratamiento de datos
# ===================================================Series de Tiempo y Ciencia de Datos
Introducción https://cienciadedatos.net
Una serie temporal (time series) es una sucesión de datos ordenados cronológicamente, espaciados a intervalos iguales o desiguales. El proceso de forecasting consiste en predecir el valor futuro de una serie temporal, bien modelando la serie únicamente en función de su comportamiento pasado (autorregresivo) o empleando otras variables externas.
Esta guía describe cómo utilizar modelos de regresión de Scikit-learn para realizar forecasting de series temporales. En concreto, se hace uso de Skforecast, una librería que contiene las clases y funciones necesarias para adaptar cualquier modelo de regresión de Scikit-learn a problemas de forecasting.
Skforecast: forecasting series temporales con Python, Machine Learning y Scikit-learn por Joaquín Amat Rodrigo y Javier Escobar Ortiz.
import numpy as np
import pandas as pd
from skforecast.datasets import fetch_dataset
# Gráficos
# ==============================================================================
import matplotlib.pyplot as plt
from skforecast.plot import set_dark_theme
# Modelado y Forecasting
# ==============================================================================
import sklearn
from sklearn.linear_model import Ridge
from sklearn.ensemble import RandomForestRegressor
from sklearn.metrics import mean_squared_error, mean_absolute_error
from sklearn.preprocessing import StandardScaler
import skforecast
from skforecast.recursive import ForecasterRecursive
from skforecast.direct import ForecasterDirect
from skforecast.model_selection import TimeSeriesFold, grid_search_forecaster, backtesting_forecaster
from skforecast.preprocessing import RollingFeatures
from skforecast.utils import save_forecaster, load_forecaster
from skforecast.metrics import calculate_coverage
from skforecast.plot import plot_prediction_intervals
import shap
# Configuración warnings
# ==============================================================================
import warnings
warnings.filterwarnings('once')Datos
Se dispone de una serie temporal con el gasto mensual (millones de dólares) en fármacos con corticoides que tuvo el sistema de salud Australiano entre 1991 y 2008. Se pretende crear un modelo autoregresivo capaz de predecir el futuro gasto mensual. Los datos empleados en los ejemplos de este documento se han obtenido del magnífico libro Forecasting: Principles and Practice by Rob J Hyndman and George Athanasopoulos.
# Descarga de datos
# ==============================================================================
datos = fetch_dataset(name='h2o_exog', raw=True)h2o_exog
--------
Monthly expenditure ($AUD) on corticosteroid drugs that the Australian health
system had between 1991 and 2008. Two additional variables (exog_1, exog_2) are
simulated.
Hyndman R (2023). fpp3: Data for Forecasting: Principles and Practice (3rd
Edition). http://pkg.robjhyndman.com/fpp3package/,
https://github.com/robjhyndman/fpp3package, http://OTexts.com/fpp3.
Shape of the dataset: (195, 4)# Preparación del dato
# ==============================================================================
datos['fecha'] = pd.to_datetime(datos['fecha'], format='%Y-%m-%d')
datos = datos.set_index('fecha')
datos = datos.asfreq('MS')
datos = datos.sort_index()
datos.head() y exog_1 exog_2
fecha
1992-04-01 0.379808 0.958792 1.166029
1992-05-01 0.361801 0.951993 1.117859
1992-06-01 0.410534 0.952955 1.067942
1992-07-01 0.483389 0.958078 1.097376
1992-08-01 0.475463 0.956370 1.122199print(f'Número de filas con missing values: {datos.isnull().any(axis=1).mean()}')Número de filas con missing values: 0.0# Verificar que un índice temporal está completo
# ==============================================================================
fecha_inicio = datos.index.min()
fecha_fin = datos.index.max()
date_range_completo = pd.date_range(start=fecha_inicio, end=fecha_fin, freq=datos.index.freq)
print(f"Índice completo: {(datos.index == date_range_completo).all()}")Índice completo: True# Separación datos train-test
# ==============================================================================
steps = 36
datos_train = datos[:-steps]
datos_test = datos[-steps:]
print(f"Fechas train : {datos_train.index.min()} --- {datos_train.index.max()} (n={len(datos_train)})")Fechas train : 1992-04-01 00:00:00 --- 2005-06-01 00:00:00 (n=159)print(f"Fechas test : {datos_test.index.min()} --- {datos_test.index.max()} (n={len(datos_test)})")Fechas test : 2005-07-01 00:00:00 --- 2008-06-01 00:00:00 (n=36)
set_dark_theme()
fig, ax = plt.subplots(figsize=(8, 5))
datos_train['y'].plot(ax=ax, label='train', fontsize=20)
datos_test['y'].plot(ax=ax, label='test', fontsize=20)
ax.legend();
Forecasting autorregresivo recursivo
Se crea y entrena un modelo ForecasterRecursive a partir de un regresor RandomForestRegressor y una ventana temporal de 6 lags. Esto último significa que, el modelo, utiliza como predictores los 6 meses anteriores.
# Crear y entrenar forecaster
# ==============================================================================
forecaster = ForecasterRecursive(
regressor = RandomForestRegressor(random_state=123),
lags = 6
)
forecaster.fit(y=datos_train['y'])
forecasterForecasterRecursive
General Information
- Regressor: RandomForestRegressor
- Lags: [1 2 3 4 5 6]
- Window features: None
- Window size: 6
- Series name: y
- Exogenous included: False
- Weight function included: False
- Differentiation order: None
- Creation date: 2025-05-27 19:05:25
- Last fit date: 2025-05-27 19:05:25
- Skforecast version: 0.16.0
- Python version: 3.10.11
- Forecaster id: None
Exogenous Variables
-
None
Data Transformations
- Transformer for y: None
- Transformer for exog: None
Training Information
- Training range: [Timestamp('1992-04-01 00:00:00'), Timestamp('2005-06-01 00:00:00')]
- Training index type: DatetimeIndex
- Training index frequency: MS
Regressor Parameters
-
{'bootstrap': True, 'ccp_alpha': 0.0, 'criterion': 'squared_error', 'max_depth': None, 'max_features': 1.0, 'max_leaf_nodes': None, 'max_samples': None, 'min_impurity_decrease': 0.0, 'min_samples_leaf': 1, 'min_samples_split': 2, 'min_weight_fraction_leaf': 0.0, 'monotonic_cst': None, 'n_estimators': 100, 'n_jobs': None, 'oob_score': False, 'random_state': 123, 'verbose': 0, 'warm_start': False}
Fit Kwargs
-
{}
# Predicciones
# ==============================================================================
steps = 36
predicciones = forecaster.predict(steps=steps)
predicciones.head(5)2005-07-01 0.878756
2005-08-01 0.882167
2005-09-01 0.973184
2005-10-01 0.983678
2005-11-01 0.849494
Freq: MS, Name: pred, dtype: float64# Gráfico de predicciones vs valores reales
# ==============================================================================
fig, ax = plt.subplots(figsize=(6, 2.5))
datos_train['y'].plot(ax=ax, label='train')
datos_test['y'].plot(ax=ax, label='test')
predicciones.plot(ax=ax, label='predicciones')
ax.legend();
# Error test
# ==============================================================================
error_mse = mean_squared_error(
y_true = datos_test['y'],
y_pred = predicciones
)
print(f"Error de test (mse): {error_mse}")Error de test (mse): 0.07326833976120374
Ajuste de hiperparámetros (tuning)
El ForecasterRecursive entrenado ha utilizado una ventana temporal de 6 lags y un modelo Random Forest con los hiperparámetros por defecto. Sin embargo, no hay ninguna razón por la que estos valores sean los más adecuados. La librería Skforecast proporciona varias estrategias de búsqueda para encontrar la mejor combinación de hiperparámetros y lags. En este caso, se utiliza la función grid_search_forecaster, que compara los resultados obtenidos con cada combinación de hiperparámetros y lags, e identifica la mejor.
# Búsqueda de hiperparámetros: grid search
# ==============================================================================
forecaster = ForecasterRecursive(
regressor = RandomForestRegressor(random_state=123),
lags = 12 # Este valor será remplazado en el grid search
)
# Particiones de entrenamiento y validación
cv = TimeSeriesFold(
steps = 36,
initial_train_size = int(len(datos_train) * 0.5),
refit = False,
fixed_train_size = False,
)
# Valores candidatos de lags
lags_grid = [10, 20]
# Valores candidatos de hiperparámetros del regresor
param_grid = {
'n_estimators': [100, 250],
'max_depth': [3, 5, 10]
}
resultados_grid = grid_search_forecaster(
forecaster = forecaster,
y = datos_train['y'],
cv = cv,
param_grid = param_grid,
lags_grid = lags_grid,
metric = 'mean_squared_error',
return_best = True,
n_jobs = 'auto',
verbose = False
)
lags grid: 0%| | 0/2 [00:00<?, ?it/s]
params grid: 0%| | 0/6 [00:00<?, ?it/s][A
params grid: 17%|#6 | 1/6 [00:01<00:05, 1.12s/it][A
params grid: 33%|###3 | 2/6 [00:02<00:04, 1.09s/it][A
params grid: 50%|##### | 3/6 [00:03<00:02, 1.03it/s][A
params grid: 67%|######6 | 4/6 [00:04<00:02, 1.02s/it][A
params grid: 83%|########3 | 5/6 [00:04<00:00, 1.05it/s][A
params grid: 100%|##########| 6/6 [00:06<00:00, 1.00it/s][A
[A
lags grid: 50%|##### | 1/2 [00:06<00:06, 6.03s/it]
lags grid: 100%|##########| 2/2 [00:11<00:00, 5.49s/it]
lags grid: 100%|##########| 2/2 [00:11<00:00, 5.57s/it]
`Forecaster` refitted using the best-found lags and parameters, and the whole data set:
Lags: [ 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20]
Parameters: {'max_depth': 3, 'n_estimators': 250}
Backtesting metric: 0.02177319540541341# Resultados de la búsqueda de hiperparámetros
# ==============================================================================
resultados_grid lags ... n_estimators
0 [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14... ... 250
1 [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14... ... 250
2 [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14... ... 250
3 [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14... ... 100
4 [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14... ... 100
5 [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14... ... 100
6 [1, 2, 3, 4, 5, 6, 7, 8, 9, 10] ... 100
7 [1, 2, 3, 4, 5, 6, 7, 8, 9, 10] ... 250
8 [1, 2, 3, 4, 5, 6, 7, 8, 9, 10] ... 100
9 [1, 2, 3, 4, 5, 6, 7, 8, 9, 10] ... 100
10 [1, 2, 3, 4, 5, 6, 7, 8, 9, 10] ... 250
11 [1, 2, 3, 4, 5, 6, 7, 8, 9, 10] ... 250
[12 rows x 6 columns]# Crear y entrenar forecaster con mejores hiperparámetros
# ==============================================================================
regressor = RandomForestRegressor(n_estimators=250, max_depth=3, random_state=123)
forecaster = ForecasterRecursive(
regressor = regressor,
lags = 20
)
forecaster.fit(y=datos_train['y'])
# Predicciones
# ==============================================================================
predicciones = forecaster.predict(steps=steps)
# Gráfico de predicciones vs valores reales
# ==============================================================================
fig, ax = plt.subplots(figsize=(6, 2.5))
datos_train['y'].plot(ax=ax, label='train')
datos_test['y'].plot(ax=ax, label='test')
predicciones.plot(ax=ax, label='predicciones')
ax.legend();
# Error de test
# ==============================================================================
error_mse = mean_squared_error(
y_true = datos_test['y'],
y_pred = predicciones
)
print(f"Error de test (mse) {error_mse}")Error de test (mse) 0.004356831371529945
Backtesting
Para obtener una estimación robusta de la capacidad predictiva del modelo, se lleva a cabo un proceso de backtesting. El backtesting consiste en evaluar el comportamiento de un modelo predictivo al aplicarlo de forma retrospectiva sobre datos históricos. Por lo tanto, es una estrategia de validación que permite cuantificar la capacidad predictiva de un modelo.
# Backtesting partitions using TimeSeriesFold
# ==============================================================================
cv = TimeSeriesFold(
steps = 12 * 3,
initial_train_size = len(datos['y']) - 12 * 9, # Last 9 years are separated for the backtest
window_size = 20,
fixed_train_size = False,
refit = True,
)
cv.split(X=datos['y'], as_pandas=True)Information of folds
--------------------
Number of observations used for initial training: 87
Number of observations used for backtesting: 108
Number of folds: 3
Number skipped folds: 0
Number of steps per fold: 36
Number of steps to exclude between last observed data (last window) and predictions (gap): 0
Fold: 0
Training: 1992-04-01 00:00:00 -- 1999-06-01 00:00:00 (n=87)
Validation: 1999-07-01 00:00:00 -- 2002-06-01 00:00:00 (n=36)
Fold: 1
Training: 1992-04-01 00:00:00 -- 2002-06-01 00:00:00 (n=123)
Validation: 2002-07-01 00:00:00 -- 2005-06-01 00:00:00 (n=36)
Fold: 2
Training: 1992-04-01 00:00:00 -- 2005-06-01 00:00:00 (n=159)
Validation: 2005-07-01 00:00:00 -- 2008-06-01 00:00:00 (n=36)
fold train_start ... test_end_with_gap fit_forecaster
0 0 0 ... 123 True
1 1 0 ... 159 True
2 2 0 ... 195 True
[3 rows x 10 columns]# Backtesting
# ==============================================================================
cv = TimeSeriesFold(
steps = 12 * 3,
initial_train_size = len(datos) - 12 * 9,
fixed_train_size = False,
refit = True,
)
metrica, predicciones_backtest = backtesting_forecaster(
forecaster = forecaster,
y = datos['y'],
cv = cv,
metric = 'mean_squared_error',
verbose = True
)Information of folds
--------------------
Number of observations used for initial training: 87
Number of observations used for backtesting: 108
Number of folds: 3
Number skipped folds: 0
Number of steps per fold: 36
Number of steps to exclude between last observed data (last window) and predictions (gap): 0
Fold: 0
Training: 1992-04-01 00:00:00 -- 1999-06-01 00:00:00 (n=87)
Validation: 1999-07-01 00:00:00 -- 2002-06-01 00:00:00 (n=36)
Fold: 1
Training: 1992-04-01 00:00:00 -- 2002-06-01 00:00:00 (n=123)
Validation: 2002-07-01 00:00:00 -- 2005-06-01 00:00:00 (n=36)
Fold: 2
Training: 1992-04-01 00:00:00 -- 2005-06-01 00:00:00 (n=159)
Validation: 2005-07-01 00:00:00 -- 2008-06-01 00:00:00 (n=36)
0%| | 0/3 [00:00<?, ?it/s]
100%|##########| 3/3 [00:00<00:00, 2999.50it/s]metrica mean_squared_error
0 0.010233
# Gráfico de predicciones de backtest vs valores reales
# ==============================================================================
fig, ax = plt.subplots(figsize=(6, 2.5))
datos.loc[predicciones_backtest.index, 'y'].plot(ax=ax, label='test')
predicciones_backtest.plot(ax=ax, label='predicciones')
ax.legend();
Explicabilidad del modelo
Debido a la naturaleza compleja de muchos de los actuales modelos de machine learning, a menudo funcionan como cajas negras, lo que dificulta entender por qué han hecho una predicción u otra. Las técnicas de explicabilidad pretenden desmitificar estos modelos, proporcionando información sobre su funcionamiento interno y ayudando a generar confianza, mejorar la transparencia y cumplir los requisitos normativos en diversos ámbitos. Mejorar la explicabilidad de los modelos no sólo ayuda a comprender su comportamiento, sino también a identificar sesgos, mejorar su rendimiento y permitir a las partes interesadas tomar decisiones más informadas basadas en los conocimientos del machine learning.
Skforecast es compatible con algunos de los métodos de explicabilidad más populares: model-specific feature importances, SHAP values, and partial dependence plots.
Importancia model-specific
# Importancia predictores
# ==============================================================================
importancia = forecaster.get_feature_importances()
importancia.head(10) feature importance
11 lag_12 0.815564
1 lag_2 0.086286
13 lag_14 0.019047
9 lag_10 0.013819
2 lag_3 0.012943
14 lag_15 0.009637
0 lag_1 0.009141
10 lag_11 0.008130
7 lag_8 0.007377
8 lag_9 0.005268Shap values
Los valores SHAP (SHapley Additive exPlanations) son un método muy utilizado para explicar los modelos de machine learning, ya que ayudan a comprender cómo influyen las variables y los valores en las predicciones de forma visual y cuantitativa.
Se puede obtener un análisis SHAP a partir de modelos skforecast con sólo dos elementos:
El regresor interno del forecaster.
Las matrices de entrenamiento creadas a partir de la serie temporal y variables exógenas, utilizadas para ajustar el pronosticador.
Aprovechando estos dos componentes, los usuarios pueden crear explicaciones interpretables para sus modelos de skforecast. Estas explicaciones pueden utilizarse para verificar la fiabilidad del modelo, identificar los factores más significativos que contribuyen a las predicciones y comprender mejor la relación subyacente entre las variables de entrada y la variable objetivo.
# Matrices de entrenamiento utilizadas por el forecaster para entrenar el regresor
# ==============================================================================
X_train, y_train = forecaster.create_train_X_y(y=datos_train['y'])
# Crear SHAP explainer (para modelos basados en árboles)
# ==============================================================================
explainer = shap.TreeExplainer(forecaster.regressor)
# Se selecciona una muestra del 50% de los datos para acelerar el cálculo
rng = np.random.default_rng(seed=785412)
sample = rng.choice(X_train.index, size=int(len(X_train) * 0.5), replace=False)
X_train_sample = X_train.loc[sample, :]
shap_values = explainer.shap_values(X_train_sample)
# Shap summary plot (top 10)
# ==============================================================================
shap.initjs()<IPython.core.display.HTML object>shap.summary_plot(shap_values, X_train_sample, max_display=15, show=False)
fig, ax = plt.gcf(), plt.gca()
ax.set_title("SHAP Summary plot")
ax.tick_params(labelsize=8)
fig.set_size_inches(8, 7)
# Backtesting indicando que se devuelvan los predictores
# ==============================================================================
cv = TimeSeriesFold(
steps = 12 * 3,
initial_train_size = len(datos) - 12 * 9,
fixed_train_size = False,
refit = True,
)
_, predicciones_backtest = backtesting_forecaster(
forecaster = forecaster,
y = datos['y'],
cv = cv,
return_predictors = True,
metric = 'mean_squared_error',
)
0%| | 0/3 [00:00<?, ?it/s]
100%|##########| 3/3 [00:00<?, ?it/s]predicciones_backtest.head() pred fold lag_1 ... lag_18 lag_19 lag_20
1999-07-01 0.706470 0 0.703872 ... 0.800544 0.994389 0.770522
1999-08-01 0.728661 0 0.706470 ... 0.490557 0.800544 0.994389
1999-09-01 0.792769 0 0.728661 ... 0.524408 0.490557 0.800544
1999-10-01 0.802519 0 0.792769 ... 0.536649 0.524408 0.490557
1999-11-01 0.843620 0 0.802519 ... 0.552091 0.536649 0.524408
[5 rows x 22 columns]# Waterfall para una predicción concreta
# ==============================================================================
iloc_predicted_date = predicciones_backtest.index.get_loc('2005-04-01')
shap_values_single = explainer(predicciones_backtest.iloc[:, 2:])
shap.plots.waterfall(shap_values_single[iloc_predicted_date], show=False)<Axes: >fig, ax = plt.gcf(), plt.gca()
fig.set_size_inches(6, 6)
plt.show()
Forecasting con variables exógenas
En el ejemplo anterior, se han utilizado como predictores únicamente lags de la propia variable predicha. En ciertos escenarios, es posible disponer de información sobre otras variables, cuyo valor a futuro se conoce, y pueden servir como predictoreres adicionales en el modelo.
Siguiendo con el ejemplo anterior, se simula una nueva variable cuyo comportamiento está correlacionado con la serie temporal modelada y que, por lo tanto, se quiere incorporar como predictor. Esto mísmo es aplicable a múltiples variables exógenas.
# Descarga de datos
# ==============================================================================
datos = fetch_dataset(name='h2o_exog', raw=True, verbose=False)
# Preparación del dato
# ==============================================================================
datos['fecha'] = pd.to_datetime(datos['fecha'], format='%Y-%m-%d')
datos = datos.set_index('fecha')
datos = datos.asfreq('MS')
datos = datos.sort_index()
fig, ax = plt.subplots(figsize=(6, 6))
datos['y'].plot(ax=ax, label='y')<Axes: xlabel='fecha'>datos['exog_1'].plot(ax=ax, label='variable exógena')<Axes: xlabel='fecha'>ax.legend();
# Separación datos train-test
# ==============================================================================
steps = 36
datos_train = datos[:-steps]
datos_test = datos[-steps:]
print(
f"Fechas train : {datos_train.index.min()} --- {datos_train.index.max()} (n={len(datos_train)})"
)Fechas train : 1992-04-01 00:00:00 --- 2005-06-01 00:00:00 (n=159)print(
f"Fechas test : {datos_test.index.min()} --- {datos_test.index.max()} (n={len(datos_test)})"
)Fechas test : 2005-07-01 00:00:00 --- 2008-06-01 00:00:00 (n=36)# Crear y entrenar forecaster
# ==============================================================================
forecaster = ForecasterRecursive(
regressor = RandomForestRegressor(random_state=123),
lags = 8
)
forecaster.fit(y=datos_train['y'], exog=datos_train['exog_1'])
# Predicciones
# ==============================================================================
predicciones = forecaster.predict(steps=steps, exog=datos_test['exog_1'])
# Gráfico predicciones vs valores reales
# ==============================================================================
fig, ax = plt.subplots(figsize=(6, 2.5))
datos_train['y'].plot(ax=ax, label='train')<Axes: xlabel='fecha'>datos_test['y'].plot(ax=ax, label='test')<Axes: xlabel='fecha'>predicciones.plot(ax=ax, label='predicciones')<Axes: xlabel='fecha'>ax.legend();
# Error test
# ==============================================================================
error_mse = mean_squared_error(
y_true = datos_test['y'],
y_pred = predicciones
)
print(f"Error de test (mse): {error_mse}")Error de test (mse): 0.03989087922533575
Predictores custom y window features
En determinados escenarios, puede ser interesante incorporar otras características de la serie temporal además de los lags, por ejemplo, la media movil de los últimos n valores puede servir para capturar la tendencia de la serie. El argumento window_features permite incorporar al modelo predictores adicionales, creados a partir de los valores pasados de la serie temporal.
La clase RollingFeatures disponible en skforecast permite crear algunos de los predictores más utilizados:
‘mean’: la media de los n valores anteriores.
‘std’: la desviación estándar de los n valores anteriores.
‘min’: el mínimo de los n valores anteriores.
‘max’: el máximo de los n valores anteriores.
‘sum’: la suma de los n valores anteriores.
‘median’: la mediana de los n valores anteriores.
‘ratio_min_max’: la relación entre el mínimo y el máximo de los n valores anteriores.
‘coef_variation’: el coeficiente de variación de los n valores anteriores.
El usuario puede especificar un tamaño de ventana diferente para cada uno de ellos o el mismo para todos.
# Descarga de datos
# ==============================================================================
datos = fetch_dataset(name='h2o_exog', raw=True, verbose=False)
# Preparación del dato
# ==============================================================================
datos['fecha'] = pd.to_datetime(datos['fecha'], format='%Y-%m-%d')
datos = datos.set_index('fecha')
datos = datos.asfreq('MS')
datos = datos.sort_index()
# Separación datos train-test
# ==============================================================================
steps = 36
datos_train = datos[:-steps]
datos_test = datos[-steps:]
print(f"Fechas train : {datos_train.index.min()} --- {datos_train.index.max()} (n={len(datos_train)})")Fechas train : 1992-04-01 00:00:00 --- 2005-06-01 00:00:00 (n=159)print(f"Fechas test : {datos_test.index.min()} --- {datos_test.index.max()} (n={len(datos_test)})")Fechas test : 2005-07-01 00:00:00 --- 2008-06-01 00:00:00 (n=36)# Window features
# ==============================================================================
window_features = RollingFeatures(
stats = ['mean', 'std', 'min', 'max'],
window_sizes = 20
)
# Crear y entrenar forecaster
# ==============================================================================
forecaster = ForecasterRecursive(
regressor = RandomForestRegressor(random_state=123),
lags = 10,
window_features = window_features,
)
forecaster.fit(y=datos_train['y'])
forecaster===================
ForecasterRecursive
===================
Regressor: RandomForestRegressor
Lags: [ 1 2 3 4 5 6 7 8 9 10]
Window features: ['roll_mean_20', 'roll_std_20', 'roll_min_20', 'roll_max_20']
Window size: 20
Series name: y
Exogenous included: False
Exogenous names: None
Transformer for y: None
Transformer for exog: None
Weight function included: False
Differentiation order: None
Training range: [Timestamp('1992-04-01 00:00:00'), Timestamp('2005-06-01 00:00:00')]
Training index type: DatetimeIndex
Training index frequency: MS
Regressor parameters:
{'bootstrap': True, 'ccp_alpha': 0.0, 'criterion': 'squared_error', 'max_depth':
None, 'max_features': 1.0, 'max_leaf_nodes': None, 'max_samples': None,
'min_impurity_decrease': 0.0, 'min_samples_leaf': 1, 'min_samples_split': 2,
'min_weight_fraction_leaf': 0.0, 'monotonic_cst': None, 'n_estimators': 100,
'n_jobs': None, 'oob_score': False, 'random_state': 123, 'verbose': 0,
'warm_start': False}
fit_kwargs: {}
Creation date: 2025-05-27 19:05:41
Last fit date: 2025-05-27 19:05:41
Skforecast version: 0.16.0
Python version: 3.10.11
Forecaster id: None # Matrices de entrenamiento
# ==============================================================================
X_train, y_train = forecaster.create_train_X_y(y=datos_train['y'])
X_train.head(5) lag_1 lag_2 lag_3 ... roll_std_20 roll_min_20 roll_max_20
fecha ...
1993-12-01 0.699605 0.632947 0.601514 ... 0.122733 0.361801 0.771258
1994-01-01 0.963081 0.699605 0.632947 ... 0.152567 0.361801 0.963081
1994-02-01 0.819325 0.963081 0.699605 ... 0.156751 0.387554 0.963081
1994-03-01 0.437670 0.819325 0.963081 ... 0.155363 0.387554 0.963081
1994-04-01 0.506121 0.437670 0.819325 ... 0.154728 0.387554 0.963081
[5 rows x 14 columns]y_train.head(5)fecha
1993-12-01 0.963081
1994-01-01 0.819325
1994-02-01 0.437670
1994-03-01 0.506121
1994-04-01 0.470491
Freq: MS, Name: y, dtype: float64# Predicciones
# ==============================================================================
steps = 36
predicciones = forecaster.predict(steps=steps)
# Gráfico predicciones vs valores reales
# ==============================================================================
fig, ax = plt.subplots(figsize=(6, 2.5))
datos_train['y'].plot(ax=ax, label='train')<Axes: xlabel='fecha'>datos_test['y'].plot(ax=ax, label='test')<Axes: xlabel='fecha'>predicciones.plot(ax=ax, label='predicciones')<Axes: xlabel='fecha'>ax.legend();
# Error test
# ==============================================================================
error_mse = mean_squared_error(
y_true = datos_test['y'],
y_pred = predicciones
)
print(f"Error de test (mse): {error_mse}")Error de test (mse): 0.04180143590431811
Direct multi-step forecasting
Para conseguir predicciones de varios steps a futuro, los modelos ForecasterRecursive siguen una estrategia de predicción recursiva en la que, cada nueva predicción, se basa en la predicción anterior. Una alternativa es entrenar un modelo para cada uno de los steps que se desea predecir. Esta estrategia, normalmente conocida como direct multi-step forecasting, es computacionalmente más costosa que la recursiva puesto que requiere entrenar varios modelos. Sin embargo, en algunos escenarios, consigue mejores resultados. Este tipo de modelos pueden obtenerse con la clase ForecasterDirect y pueden incluir también window features y variables exógenas.
# Crear forecaster
# ==============================================================================
forecaster = ForecasterDirect(
regressor = Ridge(random_state=123),
transformer_y = StandardScaler(),
steps = 36,
lags = 8 # Este valor será remplazado en el grid search
)
# Búsqueda de hiperparámetros
# ==============================================================================
from skforecast.exceptions import LongTrainingWarning
warnings.simplefilter('ignore', category=LongTrainingWarning)
cv = TimeSeriesFold(
steps = 36,
initial_train_size = int(len(datos_train) * 0.5),
fixed_train_size = False,
refit = False,
)
param_grid = {'alpha': np.logspace(-5, 5, 10)}
lags_grid = [5, 12, 20]
resultados_grid = grid_search_forecaster(
forecaster = forecaster,
y = datos_train['y'],
cv = cv,
param_grid = param_grid,
lags_grid = lags_grid,
metric = 'mean_squared_error',
return_best = True
)
lags grid: 0%| | 0/3 [00:00<?, ?it/s]
params grid: 0%| | 0/10 [00:00<?, ?it/s][A
params grid: 30%|### | 3/10 [00:00<00:00, 28.38it/s][A
params grid: 60%|###### | 6/10 [00:00<00:00, 28.58it/s][A
params grid: 90%|######### | 9/10 [00:00<00:00, 28.55it/s][A
[A
lags grid: 33%|###3 | 1/3 [00:00<00:00, 2.84it/s]
lags grid: 67%|######6 | 2/3 [00:00<00:00, 2.81it/s]
lags grid: 100%|##########| 3/3 [00:01<00:00, 2.78it/s]
lags grid: 100%|##########| 3/3 [00:01<00:00, 2.79it/s]
`Forecaster` refitted using the best-found lags and parameters, and the whole data set:
Lags: [ 1 2 3 4 5 6 7 8 9 10 11 12]
Parameters: {'alpha': 0.2782559402207126}
Backtesting metric: 0.027413948265204574# Resultados de la búsqueda de hiperparámetros
# ==============================================================================
resultados_grid.head() lags ... alpha
0 [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12] ... 0.278256
1 [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12] ... 3.593814
2 [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12] ... 0.021544
3 [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12] ... 0.001668
4 [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12] ... 0.000129
[5 rows x 5 columns]# Predicciones
# ==============================================================================
predicciones = forecaster.predict()
# Gráfico predicciones vs valores reales
# ==============================================================================
fig, ax = plt.subplots(figsize=(6, 2.5))
datos_train['y'].plot(ax=ax, label='train')<Axes: xlabel='fecha'>datos_test['y'].plot(ax=ax, label='test')<Axes: xlabel='fecha'>predicciones.plot(ax=ax, label='predicciones')<Axes: xlabel='fecha'>ax.legend();
# Error test
# ==============================================================================
error_mse = mean_squared_error(y_true = datos_test['y'], y_pred = predicciones)
print(f"Error de test (mse) {error_mse}")Error de test (mse) 0.011792965469623133
Forecasting probabilístico
El forecasting probabilístico es una familia de técnicas que permite predecir la distribución esperada del de la serie en lugar de un único valor. Este tipo de pronóstico proporciona mucha más información, ya que permite la creación de intervalos de predicción.
Un intervalo de predicción define el espacio dentro del cual es de esperar que se encuentre el verdadero valor de yy con una determinada probabilidad. Por ejemplo, es de esperar que el intervalo de predicción (1, 99) contenga el verdadero valor de la predicción con un 98% de probabilidad.
Skforecast implementa varios métodos para el pronóstico probabilístico:
# Descarga de datos
# ==============================================================================
datos = fetch_dataset(name='h2o_exog', raw=True, verbose=False)
# Preparación del dato
# ==============================================================================
datos['fecha'] = pd.to_datetime(datos['fecha'], format='%Y-%m-%d')
datos = datos.set_index('fecha')
datos = datos.asfreq('MS')
datos = datos.sort_index()
# Separación datos train-test
# ==============================================================================
steps = 36
datos_train = datos[:-steps]
datos_test = datos[-steps:]
print(f"Fechas train : {datos_train.index.min()} --- {datos_train.index.max()} (n={len(datos_train)})")Fechas train : 1992-04-01 00:00:00 --- 2005-06-01 00:00:00 (n=159)print(f"Fechas test : {datos_test.index.min()} --- {datos_test.index.max()} (n={len(datos_test)})")Fechas test : 2005-07-01 00:00:00 --- 2008-06-01 00:00:00 (n=36)# Crear y entrenar forecaster
# ==============================================================================
forecaster = ForecasterRecursive(
regressor = Ridge(alpha=0.1, random_state=765),
lags = 15
)
forecaster.fit(y=datos_train['y'], store_in_sample_residuals=True)
# Intervalos de predicción
# ==============================================================================
predicciones = forecaster.predict_interval(
steps = steps,
interval = [5, 95],
method = 'bootstrapping',
n_boot = 150
)
predicciones.head(5) pred lower_bound upper_bound
2005-07-01 0.970598 0.788690 1.059012
2005-08-01 0.990932 0.815214 1.102761
2005-09-01 1.149609 1.064396 1.258934
2005-10-01 1.194584 1.098682 1.313509
2005-11-01 1.231744 1.122570 1.355601# Error de predicción
# ==============================================================================
error_mse = mean_squared_error(
y_true = datos_test['y'],
y_pred = predicciones['pred']
)
print(f"Error de test (mse): {error_mse}")Error de test (mse): 0.010465086161791183# Gráfico
# ==============================================================================
fig, ax = plt.subplots(figsize=(6, 2.5))
plot_prediction_intervals(
predictions = predicciones,
y_true = datos_test,
target_variable = "y",
ax = ax
)
# Cobertura del intervalo predicho
# ==============================================================================
cobertura = calculate_coverage(
y_true = datos.loc[predicciones.index, 'y'],
lower_bound = predicciones['lower_bound'],
upper_bound = predicciones['upper_bound'],
)
print(f"Cobertura del intervalo predicho: {round(100 * cobertura, 2)} %")Cobertura del intervalo predicho: 83.33 %
Métrica custom
En los procesos de backtesting (backtesting_forecaster) y optimización de hiperparámetros (grid_search_forecaster), además de las métricas mean_squared_error, mean_absolute_error y mean_absolute_percentage_error, el usuario puede utilizar cualquier función que desee siempre y cuando cumpla lo siguiente:
Tiene como argumentos:
y_true: verdaderos valores de la serie.y_pred: valores predichos.
Devuelve un valor numérico (
floatoint).El modelo es mejor cuanto menor es la métrica. Esto únicamente es necesario si se quiere que la función
grid_search_forecasterreentrene automáticamente el mejor modelo encontrado.
Gracias a esta flexibilidad, es posible evaluar la capacidad predictiva del modelo con métricas aplicables a escenarios muy diversos. Por ejemplo:
Considerar únicamente determinados meses, días u horas.
Considerar únicamente fechas que sean festivos.
Considerar únicamente el último step del horizonte predicho.
Véase un ejemplo en el que se quiere predecir un horizonte de 12 meses, pero únicamente considerar los últimos 3 meses de cada año para calcular la métrica de interés.
# Métrica custom
# ==============================================================================
def custom_metric(y_true, y_pred):
'''
Calcular el mean_absolute_error utilizando únicamente las predicciones de
los últimos 3 meses del año.
'''
mask = y_true.index.month.isin([10, 11, 12])
metric = mean_absolute_error(y_true[mask], y_pred[mask])
return metric
# Backtesting
# ==============================================================================
metrica, predicciones_backtest = backtesting_forecaster(
forecaster = forecaster,
y = datos['y'],
cv = cv,
metric = custom_metric,
verbose = True
)Information of folds
--------------------
Number of observations used for initial training: 79
Number of observations used for backtesting: 116
Number of folds: 4
Number skipped folds: 0
Number of steps per fold: 36
Number of steps to exclude between last observed data (last window) and predictions (gap): 0
Last fold only includes 8 observations.
Fold: 0
Training: 1992-04-01 00:00:00 -- 1998-10-01 00:00:00 (n=79)
Validation: 1998-11-01 00:00:00 -- 2001-10-01 00:00:00 (n=36)
Fold: 1
Training: No training in this fold
Validation: 2001-11-01 00:00:00 -- 2004-10-01 00:00:00 (n=36)
Fold: 2
Training: No training in this fold
Validation: 2004-11-01 00:00:00 -- 2007-10-01 00:00:00 (n=36)
Fold: 3
Training: No training in this fold
Validation: 2007-11-01 00:00:00 -- 2008-06-01 00:00:00 (n=8)
0%| | 0/4 [00:00<?, ?it/s]
100%|##########| 4/4 [00:00<00:00, 493.71it/s]metrica custom_metric
0 0.171414
Guardar y cargar modelos
Los modelos generados con Skforecast se pueden cargar y guardar usando las librerías Pickle o Joblib. Para facilitar el proceso, dos funciones están disponibles: save_forecaster y load_forecaster. A continuación, se muestra un sencillo ejemplo. Para más información cosultar: skforecast save and load forecaster.
# Crear forecaster
# ==============================================================================
forecaster = ForecasterRecursive(RandomForestRegressor(random_state=123), lags=3)
forecaster.fit(y=datos['y'])
forecaster.predict(steps=3)2008-07-01 0.751967
2008-08-01 0.826505
2008-09-01 0.879444
Freq: MS, Name: pred, dtype: float64# Guardar modelo
# ==============================================================================
save_forecaster(forecaster, file_name='forecaster.joblib', verbose=False)
# Cargar modelo
# ==============================================================================
forecaster_cargado = load_forecaster('forecaster.joblib')===================
ForecasterRecursive
===================
Regressor: RandomForestRegressor
Lags: [1 2 3]
Window features: None
Window size: 3
Series name: y
Exogenous included: False
Exogenous names: None
Transformer for y: None
Transformer for exog: None
Weight function included: False
Differentiation order: None
Training range: [Timestamp('1992-04-01 00:00:00'), Timestamp('2008-06-01 00:00:00')]
Training index type: DatetimeIndex
Training index frequency: MS
Regressor parameters:
{'bootstrap': True, 'ccp_alpha': 0.0, 'criterion': 'squared_error', 'max_depth':
None, 'max_features': 1.0, 'max_leaf_nodes': None, 'max_samples': None,
'min_impurity_decrease': 0.0, 'min_samples_leaf': 1, 'min_samples_split': 2,
'min_weight_fraction_leaf': 0.0, 'monotonic_cst': None, 'n_estimators': 100,
'n_jobs': None, 'oob_score': False, 'random_state': 123, 'verbose': 0,
'warm_start': False}
fit_kwargs: {}
Creation date: 2025-05-27 19:05:45
Last fit date: 2025-05-27 19:05:45
Skforecast version: 0.16.0
Python version: 3.10.11
Forecaster id: None # Predicciones
# ==============================================================================
forecaster_cargado.predict(steps=3)2008-07-01 0.751967
2008-08-01 0.826505
2008-09-01 0.879444
Freq: MS, Name: pred, dtype: float64Uso de modelos en producción
En los proyectos relacionados con forecasting es frecuente que, como resultado de la etapa de experimentación y desarrollo, se genere un modelo. Para que este modelo consiga un impacto real en el negocio, se tiene que poder poner en producción y generar predicciones cada cierto tiempo, con las que tomar decisiones. Esta necesidad ha guiado en gran medida el desarrollo de la librería Skforecast.
Supóngase un caso de uso en el que se han de generar predicciones de forma semanal, por ejemplo, cada lunes el modelo tiene que predecir el resto de la semana. Una forma de conseguir este comportamiento es reentrenando el modelo semanalmente justo antes de que se ejecute la primera predicción y llamar a continuación al método predict del objeto forecaster.
# Crear y entrenar forecaster
# ==============================================================================
forecaster = ForecasterRecursive(
regressor = RandomForestRegressor(random_state=123),
lags = 6
)
forecaster.fit(y=datos_train['y'])
# Predecir con last_window
# ==============================================================================
last_window = datos_test['y'][-6:]
forecaster.predict(last_window=last_window, steps=4)2008-07-01 0.757750
2008-08-01 0.836313
2008-09-01 0.877668
2008-10-01 0.911734
Freq: MS, Name: pred, dtype: float64Casos de uso
Te invitamos a explorar nuestra colección de ejemplos y tutoriales, donde encontrarás múltiples casos de uso que te te brindarán una visión práctica de cómo aplicar esta poderosa biblioteca. Ejemplos y tutoriales.
Información de sesión
#!pip install session_info
import session_info
session_info.show(html=False)-----
matplotlib 3.10.3
numpy 1.26.4
pandas 2.2.3
session_info v1.0.1
shap 0.47.2
skforecast 0.16.0
sklearn 1.6.1
-----
Python 3.10.11 (tags/v3.10.11:7d4cc5a, Apr 5 2023, 00:38:17) [MSC v.1929 64 bit (AMD64)]
Windows-10-10.0.26100-SP0
-----
Session information updated at 2025-05-27 19:05Modelos ARIMA y SARIMAX con Python
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')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)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)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)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])})"
)Fechas train : 1969-01-01 00:00:00 --- 1980-01-01 00:00:00 (n=133)print(
f"Fechas test : {datos.loc[fin_train:].index.min()} --- {datos.loc[:].index.max()} "
f"(n={len(datos.loc[fin_train:])})"
)Fechas test : 1980-02-01 00:00:00 --- 1990-01-01 00:00:00 (n=120)datos_train = datos.loc[:fin_train]
datos_test = datos.loc[fin_train:]
datos_traindate
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=(7, 3))
datos_train.plot(ax=ax, label='train')<Axes: xlabel='date'>datos_test.plot(ax=ax, label='test')<Axes: xlabel='date'>ax.set_title('Consumo mensual combustible España')Text(0.5, 1.0, 'Consumo mensual combustible España')ax.legend();
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')Test estacionariedad serie originalprint('-------------------------------------')-------------------------------------adfuller_result = adfuller(datos)
kpss_result = kpss(datos)
print(f'ADF Statistic: {adfuller_result[0]}, p-value: {adfuller_result[1]}')ADF Statistic: -0.4461298099822807, p-value: 0.9021071923942665print(f'KPSS Statistic: {kpss_result[0]}, p-value: {kpss_result[1]}')KPSS Statistic: 2.2096370946978383, p-value: 0.01print('\nTest estacionariedad para serie diferenciada (order=1)')
Test 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]}')ADF Statistic: -3.641727690032326, p-value: 0.005011605002137185print(f'KPSS Statistic: {kpss_result[0]}, p-value: {kpss_result[1]}')KPSS Statistic: 0.3132711623572804, p-value: 0.1print('\nTest estacionariedad para serie diferenciada (order=2)')
Test 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]}')ADF Statistic: -8.233942641655926, p-value: 5.959599575498749e-13print(f'KPSS Statistic: {kpss_result[0]}, p-value: {kpss_result[1]}')KPSS Statistic: 0.080656682674821, p-value: 0.1warnings.filterwarnings("default")
# Gráfico series
# ==============================================================================
fig, axs = plt.subplots(nrows=3, ncols=1, figsize=(7, 5), sharex=True)
datos.plot(ax=axs[0], title='Serie original')<Axes: title={'center': 'Serie original'}, xlabel='date'>datos_diff_1.plot(ax=axs[1], title='Diferenciación orden 1')<Axes: title={'center': 'Diferenciación orden 1'}, xlabel='date'>datos_diff_2.plot(ax=axs[2], title='Diferenciación orden 2');
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=(6, 4), sharex=True)
plot_acf(datos, ax=axs[0], lags=50, alpha=0.05)<Figure size 600x400 with 2 Axes>axs[0].set_title('Autocorrelación serie original')Text(0.5, 1.0, 'Autocorrelación serie original')plot_acf(datos_diff_1, ax=axs[1], lags=50, alpha=0.05)<Figure size 600x400 with 2 Axes>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)<Figure size 600x300 with 2 Axes>axs[0].set_title('Autocorrelación parcial serie original')Text(0.5, 1.0, 'Autocorrelación parcial serie original')plot_pacf(datos_diff_1, ax=axs[1], lags=50, alpha=0.05)<Figure size 600x300 with 2 Axes>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 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=(9, 6), sharex=True)
res_decompose.observed.plot(ax=axs[0, 0])<Axes: xlabel='date'>axs[0, 0].set_title('Serie original', fontsize=12)Text(0.5, 1.0, 'Serie original')res_decompose.trend.plot(ax=axs[1, 0])<Axes: xlabel='date'>axs[1, 0].set_title('Tendencia', fontsize=12)Text(0.5, 1.0, 'Tendencia')res_decompose.seasonal.plot(ax=axs[2, 0])<Axes: xlabel='date'>axs[2, 0].set_title('Estacionalidad', fontsize=12)Text(0.5, 1.0, 'Estacionalidad')res_decompose.resid.plot(ax=axs[3, 0])<Axes: xlabel='date'>axs[3, 0].set_title('Residuos', fontsize=12)Text(0.5, 1.0, 'Residuos')res_descompose_diff_2.observed.plot(ax=axs[0, 1])<Axes: xlabel='date'>axs[0, 1].set_title('Series diferenciadas (order=1)', fontsize=12)Text(0.5, 1.0, 'Series diferenciadas (order=1)')res_descompose_diff_2.trend.plot(ax=axs[1, 1])<Axes: xlabel='date'>axs[1, 1].set_title('Tendencia', fontsize=12)Text(0.5, 1.0, 'Tendencia')res_descompose_diff_2.seasonal.plot(ax=axs[2, 1])<Axes: xlabel='date'>axs[2, 1].set_title('Estacionalidad', fontsize=12)Text(0.5, 1.0, 'Estacionalidad')res_descompose_diff_2.resid.plot(ax=axs[3, 1])<Axes: xlabel='date'>axs[3, 1].set_title('Residuos', fontsize=12)Text(0.5, 1.0, 'Residuos')fig.suptitle('Descomposición de la serie original vs serie diferenciada', fontsize=14)Text(0.5, 0.98, 'Descomposición de la serie original vs serie diferenciada')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]}')ADF Statistic: -4.387457230769958, p-value: 0.00031237732711269004kpss_result = kpss(datos_diff_1_12)
print(f'KPSS Statistic: {kpss_result[0]}, p-value: {kpss_result[1]}')KPSS Statistic: 0.06291573421251048, p-value: 0.1warnings.filterwarnings("default")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()<class 'statsmodels.iolib.summary.Summary'>
"""
SARIMAX Results
==========================================================================================
Dep. Variable: litters No. Observations: 133
Model: SARIMAX(1, 1, 1)x(1, 1, 1, 12) Log Likelihood -1356.051
Date: ma., 27 may. 2025 AIC 2722.103
Time: 19:05:49 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 3.12e+34. 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.056918
1980-03-01 473997.245805
1980-04-01 489983.091484
1980-05-01 485517.462858
Freq: MS, Name: predicciones_statsmodels, dtype: float64Skforecast
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
ydate
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: float64modelo.fit(y=datos_train)
modelo.summary()<class 'statsmodels.iolib.summary.Summary'>
"""
SARIMAX Results
==========================================================================================
Dep. Variable: litters No. Observations: 133
Model: SARIMAX(1, 1, 1)x(1, 1, 1, 12) Log Likelihood -1356.051
Date: ma., 27 may. 2025 AIC 2722.103
Time: 19:05:49 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 3.12e+34. Standard errors may be unstable.
"""warnings.filterwarnings("default")
# Predictión
# ==============================================================================
predicciones_skforecast = modelo.predict(steps=len(datos_test))
predicciones_skforecast.head(4) pred
1980-02-01 407504.056918
1980-03-01 473997.245805
1980-04-01 489983.091484
1980-05-01 485517.462858pdmarima
# 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)ARIMA(order=(1, 1, 1), seasonal_order=(1, 1, 1, 12))modelo.summary()<class 'statsmodels.iolib.summary.Summary'>
"""
SARIMAX Results
==========================================================================================
Dep. Variable: y No. Observations: 133
Model: SARIMAX(1, 1, 1)x(1, 1, 1, 12) Log Likelihood -1355.749
Date: ma., 27 may. 2025 AIC 2723.498
Time: 19:05:50 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.311385
1980-03-01 472944.444440
1980-04-01 488389.125314
1980-05-01 483432.075663
Freq: MS, Name: predicciones_pdmarima, dtype: float64# Plot predictions
# ==============================================================================
fig, ax = plt.subplots(figsize=(7, 3))
datos_train.plot(ax=ax, label='train')<Axes: xlabel='date'>datos_test.plot(ax=ax, label='test')<Axes: xlabel='date'>predicciones_statsmodels.plot(ax=ax, label='statsmodels')<Axes: xlabel='date'>predicciones_skforecast.columns = ['skforecast']
predicciones_skforecast.plot(ax=ax, label='skforecast')<Axes: xlabel='date'>predicciones_pdmarima.plot(ax=ax, label='pmdarima')<Axes: xlabel='date'>ax.set_title('Predicciones con modelos ARIMA')Text(0.5, 1.0, 'Predicciones con modelos ARIMA')ax.legend();
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.056918
1980-03-01 473997.245805
1980-04-01 489983.091484
1980-05-01 485517.462858
Freq: MS, Name: pred, dtype: float64Backtesting
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 = "auto",
suppress_warnings_fit = True,
verbose = True,
show_progress = True
)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.89it/s]
40%|#### | 4/10 [00:00<00:00, 9.14it/s]
60%|###### | 6/10 [00:00<00:00, 7.86it/s]
70%|####### | 7/10 [00:00<00:00, 7.26it/s]
80%|######## | 8/10 [00:01<00:00, 7.02it/s]
90%|######### | 9/10 [00:01<00:00, 6.59it/s]
100%|##########| 10/10 [00:01<00:00, 6.13it/s]
100%|##########| 10/10 [00:01<00:00, 7.13it/s]metrica mean_absolute_error
0 19611.236349predicciones.head(4) pred
1980-02-01 407504.056918
1980-03-01 473997.245805
1980-04-01 489983.091484
1980-05-01 485517.462858# Gráfico predicciones de backtesting
# ==============================================================================
fig, ax = plt.subplots(figsize=(6, 3))
datos.loc[fin_train:].plot(ax=ax, label='test')<Axes: xlabel='date'>predicciones.plot(ax=ax)<Axes: xlabel='date'>ax.set_title('Predicciones de backtesting con un modelo SARIMAX')Text(0.5, 1.0, 'Predicciones de backtesting con un modelo SARIMAX')ax.legend();
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])})"
)Fechas entrenamiento : 1969-01-01 00:00:00 --- 1976-01-01 00:00:00 (n=85)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])})"
)Fechas validacion : 1976-02-01 00:00:00 --- 1984-01-01 00:00:00 (n=96)print(
f"Fechas test : {datos.loc[fin_val:].index.min()} --- {datos.index.max()} "
f"(n={len(datos.loc[fin_val:])})"
)Fechas test : 1984-02-01 00:00:00 --- 1990-01-01 00:00:00 (n=72)# Gráfico
# ======================================================================================
fig, ax = plt.subplots(figsize=(7, 3))
datos.loc[:fin_train].plot(ax=ax, label='entrenamiento')<Axes: xlabel='date'>datos.loc[fin_train:fin_val].plot(ax=ax, label='validación')<Axes: xlabel='date'>datos.loc[fin_val:].plot(ax=ax, label='test')<Axes: xlabel='date'>ax.set_title('Consumo mensual combustible España')Text(0.5, 1.0, 'Consumo mensual combustible España')ax.legend();
# 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,
)Number of models compared: 45.
params grid: 0%| | 0/45 [00:00<?, ?it/s]
params grid: 9%|8 | 4/45 [00:00<00:01, 31.66it/s]
params grid: 18%|#7 | 8/45 [00:00<00:03, 9.88it/s]
params grid: 22%|##2 | 10/45 [00:01<00:04, 8.10it/s]
params grid: 29%|##8 | 13/45 [00:01<00:03, 10.16it/s]
params grid: 33%|###3 | 15/45 [00:01<00:02, 10.20it/s]
params grid: 38%|###7 | 17/45 [00:02<00:05, 5.39it/s]
params grid: 40%|#### | 18/45 [00:02<00:06, 4.25it/s]
params grid: 47%|####6 | 21/45 [00:02<00:03, 6.50it/s]
params grid: 51%|#####1 | 23/45 [00:03<00:02, 7.65it/s]
params grid: 56%|#####5 | 25/45 [00:03<00:03, 6.24it/s]
params grid: 60%|###### | 27/45 [00:04<00:04, 4.28it/s]
params grid: 64%|######4 | 29/45 [00:04<00:03, 5.19it/s]
params grid: 67%|######6 | 30/45 [00:04<00:02, 5.52it/s]
params grid: 69%|######8 | 31/45 [00:04<00:02, 5.98it/s]
params grid: 71%|#######1 | 32/45 [00:04<00:02, 6.48it/s]
params grid: 73%|#######3 | 33/45 [00:04<00:01, 6.72it/s]
params grid: 76%|#######5 | 34/45 [00:05<00:02, 3.97it/s]
params grid: 78%|#######7 | 35/45 [00:06<00:03, 3.03it/s]
params grid: 80%|######## | 36/45 [00:06<00:03, 2.46it/s]
params grid: 82%|########2 | 37/45 [00:06<00:02, 3.13it/s]
params grid: 84%|########4 | 38/45 [00:06<00:01, 3.88it/s]
params grid: 87%|########6 | 39/45 [00:07<00:01, 4.50it/s]
params grid: 89%|########8 | 40/45 [00:07<00:01, 4.38it/s]
params grid: 91%|#########1| 41/45 [00:07<00:00, 4.32it/s]
params grid: 93%|#########3| 42/45 [00:07<00:00, 4.12it/s]
params grid: 96%|#########5| 43/45 [00:09<00:01, 1.77it/s]
params grid: 98%|#########7| 44/45 [00:10<00:00, 1.26it/s]
params grid: 100%|##########| 45/45 [00:11<00:00, 1.07it/s]
params grid: 100%|##########| 45/45 [00:11<00:00, 3.85it/s]resultados_grid.head(5) params ... trend
0 {'order': (0, 1, 1), 'seasonal_order': (1, 1, ... ... n
1 {'order': (0, 1, 1), 'seasonal_order': (1, 1, ... ... None
2 {'order': (1, 1, 1), 'seasonal_order': (1, 1, ... ... n
3 {'order': (1, 1, 1), 'seasonal_order': (1, 1, ... ... None
4 {'order': (2, 1, 1), 'seasonal_order': (1, 1, ... ... n
[5 rows x 5 columns]# 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.79 sec
ARIMA(0,1,1)(1,1,0)[12] : AIC=3852.231, Time=0.08 sec
ARIMA(0,1,1)(1,1,2)[12] : AIC=3842.593, Time=1.31 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.13 sec
ARIMA(1,1,1)(0,1,2)[12] : AIC=3835.160, Time=0.72 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.07 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.38 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.36 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: 6.207 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.12
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.13
6 (2,1,1) (0,1,1,12) 3836.104 0.38
7 (1,1,2) (0,1,1,12) 3836.107 0.21
8 (1,1,1) (1,1,2,12) 3836.443 1.07
9 (1,1,0) (0,1,1,12) 3836.696 0.08
10 (2,1,2) (0,1,1,12) 3837.435 0.35
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.13
14 (0,1,1) (1,1,2,12) 3842.593 1.04
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.06
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,
)
0%| | 0/6 [00:00<?, ?it/s]
33%|###3 | 2/6 [00:00<00:00, 11.34it/s]
67%|######6 | 4/6 [00:00<00:00, 8.20it/s]
83%|########3 | 5/6 [00:00<00:00, 5.44it/s]
100%|##########| 6/6 [00:00<00:00, 6.21it/s]
100%|##########| 6/6 [00:00<00:00, 6.62it/s]# 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%|###3 | 2/6 [00:00<00:00, 12.67it/s]
67%|######6 | 4/6 [00:00<00:00, 9.25it/s]
83%|########3 | 5/6 [00:00<00:00, 8.65it/s]
100%|##########| 6/6 [00:00<00:00, 7.74it/s]
100%|##########| 6/6 [00:00<00:00, 8.44it/s]# Comparación de métricas
# ==============================================================================
print("Metrica (mean absolute error) del modelo grid search:")Metrica (mean absolute error) del modelo grid search:metrica_m1 mean_absolute_error
0 19803.080694print("Metric (mean_absolute_error) del modelo auto arima:")Metric (mean_absolute_error) del modelo auto arima:metrica_m2 mean_absolute_error
0 20149.352205fig, ax = plt.subplots(figsize=(6, 3))
datos.loc[fin_val:].plot(ax=ax, label='test')<Axes: xlabel='date'>predicciones_m1 = predicciones_m1.rename(columns={'pred': 'grid search'})
predicciones_m2 = predicciones_m2.rename(columns={'pred': 'autoarima'})
predicciones_m1.plot(ax=ax)<Axes: xlabel='date'>predicciones_m2.plot(ax=ax)<Axes: xlabel='date'>ax.set_title('Predicciones de backtesting con un modelo SARIMA')Text(0.5, 1.0, 'Predicciones de backtesting con un modelo SARIMA')ax.legend();
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])})"
)Fechas entrenamiento : 1969-01-01 00:00:00 --- 1980-01-01 00:00:00 (n=133)print(
f"Fechas Last window : {datos.loc[fin_train:].index.min()} --- {datos.index.max()} "
f"(n={len(datos.loc[fin_train:])})"
)Fechas Last window : 1980-02-01 00:00:00 --- 1990-01-01 00:00:00 (n=120)# Gráfico
# ======================================================================================
fig, ax = plt.subplots(figsize=(7, 3))
datos.loc[:fin_train].plot(ax=ax, label='entrenamiento')<Axes: xlabel='date'>datos.loc[fin_train:].plot(ax=ax, label='last window')<Axes: xlabel='date'>ax.set_title('Consumo mensual combustible España')Text(0.5, 1.0, 'Consumo mensual combustible España')ax.legend();
# 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.320778
1990-03-01 693624.449188
1990-04-01 654315.472027
Freq: MS, Name: pred, dtype: float64# Gráfico predicciones
# ======================================================================================
fig, ax = plt.subplots(figsize=(7, 3))
datos.loc[:fin_train].plot(ax=ax, label='entrenamiento')<Axes: xlabel='date'>datos.loc[fin_train:].plot(ax=ax, label='last window')<Axes: xlabel='date'>predicciones.plot(ax=ax, label='predicciones')<Axes: xlabel='date'>ax.set_title('Consumo mensual combustible España')Text(0.5, 1.0, 'Consumo mensual combustible España')ax.legend();
Forecasting series temporales con gradient boosting: Skforecast, XGBoost, LightGBM y CatBoost
Introducción
Los modelos gradient boosting destacan dentro de la comunidad de machine learning debido a su capacidad para lograr excelentes resultados en una amplia variedad de casos de uso, incluyendo tanto la regresión como la clasificación. Aunque su uso en el forecasting de series temporales ha sido limitado, también pueden conseguir resultados muy competitivos en este ámbito. Algunas de las ventajas que presentan los modelos gradient boosting para forecasting son:
La facilidad con que pueden incorporarse al modelo variables exógenas, además de las autorregresivas.
La capacidad de capturar relaciones no lineales entre variables.
Alta escalabilidad, que permite a los modelos manejar grandes volúmenes de datos.
Algunas implementaciones permiten la inclusión de variables categóricas sin necesidad de codificación adicional, como la codificación one-hot.
A pesar de estas ventajas, el uso de modelos de machine learning para forecasting presenta varios retos que pueden hacer que el analista sea reticente a su uso, los principales son:
Reestructurar los datos para poder utilizarlos como si se tratara de un problema de regresión.
Dependiendo de cuántas predicciones futuras se necesiten (horizonte de predicción), puede ser necesario implementar un proceso iterativo en el que cada nueva predicción se base en las anteriores.
La validación de los modelos requiere de estrategias específicas como backtesting, walk-forward validation o time series cross-validation. No puede aplicarse la validación cruzada tradicional.
La librería skforecast ofrece soluciones automatizadas a estos retos, facilitando el uso y la validación de modelos de machine learning en problemas de forecasting. Skforecast es compatible con las implementaciones de gradient boosting más avanzadas, incluyendo XGBoost, LightGBM, Catboost y HistGradientBoostingRegressor. Este documento muestra cómo utilizarlos para construir modelos de forecasting precisos.
Para garantizar una experiencia de aprendizaje fluida, se realiza una exploración inicial de los datos. A continuación, se explica paso a paso el proceso de modelización, empezando por un modelo recursivo que utiliza un regresor LightGBM y pasando por un modelo que incorpora variables exógenas y diversas estrategias de codificación. El documento concluye demostrando el uso de otras implementaciones de modelos de gradient boosting, como XGBoost, CatBoost y el HistGradientBoostingRegressor de scikit-learn.
Caso de uso
Los sistemas de bicicletas compartidas, también conocidos como sistemas de bicicletas públicas, facilitan la disponibilidad automática de bicicletas para que sean utilizadas temporalmente como medio de transporte. La mayoría de estos sistemas permiten recoger una bicicleta y devolverla en un punto diferente (estaciones o dockers), para que el usuario solo necesite tener la bicicleta en su posesión durante el desplazamiento. Uno de los principales retos en la gestión de estos sistemas es la necesidad de redistribuir las bicicletas para intentar que, en todas las estaciones, haya bicicletas disponibles a la vez que espacios libres para devoluciones.
Con el objetivo de mejorar la planificación y ejecución de la distribución de las bicicletas, se plantea crear un modelo capaz de predecir el número de usuarios para las siguientes 36 horas. De esta forma, a las 12h de cada día, la compañía encargada de gestionar las estaciones de alquiler podrá conocer la demanda prevista el resto del día (12 horas) y el siguiente día (24 horas).
A efectos ilustrativos, el ejemplo actual sólo modela una estación, sin embargo, el modelo puede ampliarse para cubrir múltiples estaciones utilizando global multi-series forecasting, mejorando así la gestión de los sistemas de bicicletas compartidas a mayor escala.
Librerías
Las librerías utilizadas en este documento son:
# Data processing
# ==============================================================================
#!pip install astral
#!pip install feature_engine
#!pip install xgboost
#!pip install lightgbm
#!pip install catboost
import numpy as np
import pandas as pd
from astral.sun import sun
from astral import LocationInfo
from skforecast.datasets import fetch_dataset
from feature_engine.datetime import DatetimeFeatures
from feature_engine.creation import CyclicalFeatures
from feature_engine.timeseries.forecasting import WindowFeatures
# Plots
# ==============================================================================
import matplotlib.pyplot as plt
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from skforecast.plot import plot_residuals, calculate_lag_autocorrelation
import plotly.graph_objects as go
import plotly.io as pio
import plotly.offline as poff
pio.templates.default = "seaborn"
poff.init_notebook_mode(connected=True){'text/html': ' <script type="text/javascript">\n window.PlotlyConfig = {MathJaxConfig: \'local\'};\n if (window.MathJax && window.MathJax.Hub && window.MathJax.Hub.Config) {window.MathJax.Hub.Config({SVG: {font: "STIX-Web"}});}\n </script>\n <script type="module">import "https://cdn.plot.ly/plotly-3.0.1.min"</script>\n '}plt.style.use('seaborn-v0_8-darkgrid')
# Modelling and Forecasting
# ==============================================================================
import xgboost
import lightgbm
import catboost
import sklearn
import shap
from xgboost import XGBRegressor
from lightgbm import LGBMRegressor
from catboost import CatBoostRegressor
from sklearn.ensemble import HistGradientBoostingRegressor
from sklearn.preprocessing import (
OneHotEncoder,
OrdinalEncoder,
FunctionTransformer,
PolynomialFeatures,
)
from sklearn.feature_selection import RFECV
from sklearn.pipeline import make_pipeline
from sklearn.compose import make_column_transformer, make_column_selector
import skforecast
from skforecast.recursive import ForecasterEquivalentDate, ForecasterRecursive
from skforecast.model_selection import (
TimeSeriesFold,
OneStepAheadFold,
bayesian_search_forecaster,
backtesting_forecaster,
)
from skforecast.preprocessing import RollingFeatures
from skforecast.feature_selection import select_features
from skforecast.metrics import calculate_coverage
# Configuración warnings
# ==============================================================================
import warnings
warnings.filterwarnings('once')Datos
Los datos empleados en este documento representan el uso, a nivel horario, del sistema de alquiler de bicicletas en la ciudad de Washington D.C. durante los años 2011 y 2012. Además del número de usuarios por hora, se dispone de información sobre las condiciones meteorológicas y sobre los días festivos. Los datos originales se han obtenido del UCI Machine Learning Repository y han sido procesados previamente (código) aplicando las siguientes modificaciones :
Columnas renombradas con nombres más descriptivos.
Categorías de la variable meteorológica renombradas. La categoría de heavy rain, se ha combinado con la de rain.
Variables de temperatura, humedad y viento desnormalizadas.
Creada variable date_time y establecida como índice.
Imputación de valores missing mediante forward fill.
# Descarga de datos
# ==============================================================================
datos = fetch_dataset('bike_sharing', raw=True)bike_sharing
------------
Hourly usage of the bike share system in the city of Washington D.C. during the
years 2011 and 2012. In addition to the number of users per hour, information
about weather conditions and holidays is available.
Fanaee-T,Hadi. (2013). Bike Sharing Dataset. UCI Machine Learning Repository.
https://doi.org/10.24432/C5W894.
Shape of the dataset: (17544, 12)# Preprocesado de datos (estableciendo índice y frecuencia)
# ==============================================================================
datos = datos[['date_time', 'users', 'holiday', 'weather', 'temp', 'atemp', 'hum', 'windspeed']]
datos['date_time'] = pd.to_datetime(datos['date_time'], format='%Y-%m-%d %H:%M:%S')
datos = datos.set_index('date_time')
if pd.__version__ < '2.2':
datos = datos.asfreq('H')
else:
datos = datos.asfreq('h')
datos = datos.sort_index()
datos.head() users holiday weather temp atemp hum windspeed
date_time
2011-01-01 00:00:00 16.0 0.0 clear 9.84 14.395 81.0 0.0
2011-01-01 01:00:00 40.0 0.0 clear 9.02 13.635 80.0 0.0
2011-01-01 02:00:00 32.0 0.0 clear 9.02 13.635 80.0 0.0
2011-01-01 03:00:00 13.0 0.0 clear 9.84 14.395 75.0 0.0
2011-01-01 04:00:00 1.0 0.0 clear 9.84 14.395 75.0 0.0# Separación de datos en entrenamiento, validación y test
# ==============================================================================
fin_train = '2012-04-30 23:59:00'
fin_validacion = '2012-08-31 23:59:00'
datos_train = datos.loc[: fin_train, :]
datos_val = datos.loc[fin_train:fin_validacion, :]
datos_test = datos.loc[fin_validacion:, :]
print(
f"Fechas train : {datos_train.index.min()} --- {datos_train.index.max()} "
f"(n={len(datos_train)})"
)Fechas train : 2011-01-01 00:00:00 --- 2012-04-30 23:00:00 (n=11664)print(
f"Fechas validación : {datos_val.index.min()} --- {datos_val.index.max()} "
f"(n={len(datos_val)})"
)Fechas validación : 2012-05-01 00:00:00 --- 2012-08-31 23:00:00 (n=2952)print(
f"Fechas test : {datos_test.index.min()} --- {datos_test.index.max()} "
f"(n={len(datos_test)})"
)Fechas test : 2012-09-01 00:00:00 --- 2012-12-31 23:00:00 (n=2928)Exploración gráfica
La exploración gráfica de series temporales es una forma eficaz de identificar tendencias, patrones y variaciones estacionales. Esto, a su vez, ayuda a orientar la selección del modelo de forecasting más adecuado.
Representación de la serie temporal
Gráficos de estacionalidad
Los gráficos estacionales son una herramienta útil para identificar patrones y tendencias estacionales en una serie temporal. Se crean promediando los valores de cada estación a lo largo del tiempo y luego trazándolos en función del tiempo.
# Estacionalidad anual, semanal y diaria
# ==============================================================================
fig, axs = plt.subplots(2, 2, figsize=(8, 5), sharex=False, sharey=True)
axs = axs.ravel()
# Distribusión de usuarios por mes
datos['month'] = datos.index.month
datos.boxplot(column='users', by='month', ax=axs[0], flierprops={'markersize': 3, 'alpha': 0.3})<Axes: title={'center': 'users'}, xlabel='month'>datos.groupby('month')['users'].median().plot(style='o-', linewidth=0.8, ax=axs[0])<Axes: title={'center': 'users'}, xlabel='month'>axs[0].set_ylabel('Users')Text(0, 0.5, 'Users')axs[0].set_title('Distribución de usuarios por mes', fontsize=9)Text(0.5, 1.0, 'Distribución de usuarios por mes')# Distribusión de usuarios por día de la semana
datos['week_day'] = datos.index.day_of_week + 1
datos.boxplot(column='users', by='week_day', ax=axs[1], flierprops={'markersize': 3, 'alpha': 0.3})<Axes: title={'center': 'users'}, xlabel='week_day'>datos.groupby('week_day')['users'].median().plot(style='o-', linewidth=0.8, ax=axs[1])<Axes: title={'center': 'users'}, xlabel='week_day'>axs[1].set_ylabel('Users')Text(0, 0.5, 'Users')axs[1].set_title('Distribución de usuarios por día de la semana', fontsize=9)Text(0.5, 1.0, 'Distribución de usuarios por día de la semana')# Distribusión de usuarios por hora del día
datos['hour_day'] = datos.index.hour + 1
datos.boxplot(column='users', by='hour_day', ax=axs[2], flierprops={'markersize': 3, 'alpha': 0.3})<Axes: title={'center': 'users'}, xlabel='hour_day'>datos.groupby('hour_day')['users'].median().plot(style='o-', linewidth=0.8, ax=axs[2])<Axes: title={'center': 'users'}, xlabel='hour_day'>axs[2].set_ylabel('Users')Text(0, 0.5, 'Users')axs[2].set_title('Distribución de usuarios por hora del día', fontsize=9)Text(0.5, 1.0, 'Distribución de usuarios por hora del día')# Distribusión de usuarios por día de la semana y hora del día
mean_day_hour = datos.groupby(["week_day", "hour_day"])["users"].mean()
mean_day_hour.plot(ax=axs[3])<Axes: xlabel='week_day,hour_day'>axs[3].set(
title = "Promedio de usuarios",
xticks = [i * 24 for i in range(7)],
xticklabels = ["Mon", "Tue", "Wed", "Thu", "Fri", "Sat", "Sun"],
xlabel = "Día y hora",
ylabel = "Users"
)[Text(0.5, 1.0, 'Promedio de usuarios'), [<matplotlib.axis.XTick object at 0x0000014DB36EE290>, <matplotlib.axis.XTick object at 0x0000014E28A35480>, <matplotlib.axis.XTick object at 0x0000014E28C12D70>, <matplotlib.axis.XTick object at 0x0000014E28A36B90>, <matplotlib.axis.XTick object at 0x0000014E02DF69B0>, <matplotlib.axis.XTick object at 0x0000014E28C20E50>, <matplotlib.axis.XTick object at 0x0000014E02DF56F0>], [Text(0, 0, 'Mon'), Text(24, 0, 'Tue'), Text(48, 0, 'Wed'), Text(72, 0, 'Thu'), Text(96, 0, 'Fri'), Text(120, 0, 'Sat'), Text(144, 0, 'Sun')], Text(0.5, 0, 'Día y hora'), Text(0, 0.5, 'Users')]axs[3].title.set_size(9)
fig.suptitle("Gráficos de estacionalidad", fontsize=12)Text(0.5, 0.98, 'Gráficos de estacionalidad')fig.tight_layout()
Existe una clara diferencia entre los días entre semana y el fin de semana. También se observa un claro patrón intradiario, con diferente afluencia de usuarios dependiendo de la hora del día.
Gráficos de autocorrelación
Los gráficos de autocorrelación muestran la correlación entre una serie temporal y sus valores pasados. Son una herramienta útil para identificar el orden de un modelo autorregresivo, es decir, los valores pasados (lags) que se deben incluir en el modelo.
La función de autocorrelación (ACF) mide la correlación entre una serie temporal y sus valores pasados. La función de autocorrelación parcial (PACF) mide la correlación entre una serie temporal y sus valores pasados, pero solo después de eliminar las variaciones explicadas por los valores pasados intermedios.
# Gráfico autocorrelación
# ==============================================================================
fig, ax = plt.subplots(figsize=(5, 2))
plot_acf(datos['users'], ax=ax, lags=72)<Figure size 500x200 with 1 Axes>plt.show()
# Gráfico autocorrelación parcial
# ==============================================================================
fig, ax = plt.subplots(figsize=(5, 2))
plot_pacf(datos['users'], ax=ax, lags=72, method='ywm')<Figure size 500x200 with 1 Axes>plt.show()
# Top 10 lags con mayor autocorrelación parcial absoluta
# ==============================================================================
calculate_lag_autocorrelation(
data = datos['users'],
n_lags = 60,
sort_by = "partial_autocorrelation_abs"
).head(10) lag partial_autocorrelation_abs ... autocorrelation_abs autocorrelation
0 1 0.845220 ... 0.845172 0.845172
1 2 0.408349 ... 0.597704 0.597704
2 23 0.355669 ... 0.708470 0.708470
3 22 0.343601 ... 0.520804 0.520804
4 25 0.332366 ... 0.711256 0.711256
5 10 0.272649 ... 0.046483 -0.046483
6 17 0.241984 ... 0.057267 -0.057267
7 19 0.199286 ... 0.159897 0.159897
8 21 0.193404 ... 0.373666 0.373666
9 3 0.182068 ... 0.409680 0.409680
[10 rows x 5 columns]Los resultados del estudio de autocorrelación indican una correlación significativa entre el número de usuarios en las horas anteriores, así como en los días previos. Esto significa que conocer del número de usuarios durante periodos específicos del pasado proporciona información útil para predecir el número de usuarios en el futuro.
Baseline
Al enfrentarse a un problema de forecasting, es recomendable disponer de un modelo de referencia (baseline). Se trata de un modelo muy sencillo que puede utilizarse como referencia para evaluar si merece la pena aplicar modelos más complejos.
Skforecast permite crear fácilmente un modelo de referencia con su clase ForecasterEquivalentDate. Este modelo, también conocido como Seasonal Naive Forecasting, simplemente devuelve el valor observado en el mismo periodo de la temporada anterior (por ejemplo, el mismo día laboral de la semana anterior, la misma hora del día anterior, etc.).
A partir del análisis exploratorio realizado, el modelo de referencia será el que prediga cada hora utilizando el valor de la misma hora del día anterior.
# Crear un baseline: valor de la misma hora del día anterior
# ==============================================================================
forecaster = ForecasterEquivalentDate(
offset = pd.DateOffset(days=1),
n_offsets = 1
)
# Entremaiento del forecaster
# ==============================================================================
forecaster.fit(y=datos.loc[:fin_validacion, 'users'])
forecaster========================
ForecasterEquivalentDate
========================
Offset: <DateOffset: days=1>
Number of offsets: 1
Aggregation function: mean
Window size: 24
Series name: users
Training range: [Timestamp('2011-01-01 00:00:00'), Timestamp('2012-08-31 23:00:00')]
Training index type: DatetimeIndex
Training index frequency: h
Creation date: 2025-05-27 19:06:24
Last fit date: 2025-05-27 19:06:24
Skforecast version: 0.16.0
Python version: 3.10.11
Forecaster id: None # Backtesting
# ==============================================================================
cv = TimeSeriesFold(steps = 36, initial_train_size = len(datos.loc[:fin_validacion]))
metrica_baseline, predicciones = backtesting_forecaster(
forecaster = forecaster,
y = datos['users'],
cv = cv,
metric = 'mean_absolute_error'
)
0%| | 0/82 [00:00<?, ?it/s]
61%|###### | 50/82 [00:00<00:00, 489.25it/s]
100%|##########| 82/82 [00:00<00:00, 497.77it/s]metrica_baseline mean_absolute_error
0 91.668716El MAE del modelo baseline se utiliza como referencia para evaluar si merece la pena aplicar los modelos más complejos.
Forecasting con con LightGBM
LightGBM es una implementación altamente eficiente del algoritmo gradient boosting, que se ha convertido en un referente en el campo del machine learning. La librería LightGBM incluye su propia API, así como la API de scikit-learn, lo que la hace compatible con skforecast.
En primer lugar, se entrena un modelo ForecasterAutoreg utilizando valores pasados de la variable de respuesta (lags) y la media movil como predictores. Posteriormente, se añaden variables exógenas al modelo y se evalúa la mejora de su rendimiento. Dado que los modelos de Gradient Boosting tienen un gran número de hiperparámetros, se realiza una Búsqueda Bayesiana utilizando la función bayesian_search_forecaster() para encontrar la mejor combinación de hiperparámetros y lags. Por último, se evalúa la capacidad predictiva del modelo mediante un proceso de backtesting.
Forecaster
# Crear el forecaster
# ==============================================================================
window_features = RollingFeatures(stats=["mean"], window_sizes=24 * 3)
forecaster = ForecasterRecursive(
regressor = LGBMRegressor(random_state=15926, verbose=-1),
lags = 24,
window_features = window_features
)
# Entrenar el forecaster
# ==============================================================================
forecaster.fit(y=datos.loc[:fin_validacion, 'users'])
forecaster===================
ForecasterRecursive
===================
Regressor: LGBMRegressor
Lags: [ 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24]
Window features: ['roll_mean_72']
Window size: 72
Series name: users
Exogenous included: False
Exogenous names: None
Transformer for y: None
Transformer for exog: None
Weight function included: False
Differentiation order: None
Training range: [Timestamp('2011-01-01 00:00:00'), Timestamp('2012-08-31 23:00:00')]
Training index type: DatetimeIndex
Training index frequency: h
Regressor parameters:
{'boosting_type': 'gbdt', 'class_weight': None, 'colsample_bytree': 1.0,
'importance_type': 'split', 'learning_rate': 0.1, 'max_depth': -1,
'min_child_samples': 20, 'min_child_weight': 0.001, 'min_split_gain': 0.0,
'n_estimators': 100, 'n_jobs': None, 'num_leaves': 31, 'objective': None,
'random_state': 15926, 'reg_alpha': 0.0, 'reg_lambda': 0.0, 'subsample':
1.0, 'subsample_for_bin': 200000, 'subsample_freq': 0, 'verbose': -1}
fit_kwargs: {}
Creation date: 2025-05-27 19:06:24
Last fit date: 2025-05-27 19:06:26
Skforecast version: 0.16.0
Python version: 3.10.11
Forecaster id: None # Predicciones
# ==============================================================================
forecaster.predict(steps=10)2012-09-01 00:00:00 108.331027
2012-09-01 01:00:00 68.562982
2012-09-01 02:00:00 33.499525
2012-09-01 03:00:00 10.027583
2012-09-01 04:00:00 3.037563
2012-09-01 05:00:00 17.162543
2012-09-01 06:00:00 51.059825
2012-09-01 07:00:00 146.940053
2012-09-01 08:00:00 344.320596
2012-09-01 09:00:00 439.738683
Freq: h, Name: pred, dtype: float64Backtesting
Para obtener una estimación robusta de la capacidad predictiva del modelo, se realiza un proceso de backtesting. El proceso de backtesting consiste en generar una predicción para cada observación del conjunto de test, siguiendo el mismo procedimiento que se seguiría si el modelo estuviese en producción, y finalmente comparar el valor predicho con el valor real.
Se recomienda revisar la documentación de la función backtesting_forecaster para comprender mejor sus capacidades. Esto ayudará a utilizar todo su potencial para analizar la capacidad predictiva del modelo.
# Backtest del modelo con lo datos de test
# ==============================================================================
cv = TimeSeriesFold(steps = 36, initial_train_size = len(datos.loc[:fin_validacion]))
metrica, predicciones = backtesting_forecaster(
forecaster = forecaster,
y = datos['users'],
cv = cv,
metric = 'mean_absolute_error'
)
0%| | 0/82 [00:00<?, ?it/s]
7%|7 | 6/82 [00:00<00:01, 54.02it/s]
15%|#4 | 12/82 [00:00<00:01, 56.59it/s]
22%|##1 | 18/82 [00:00<00:01, 57.32it/s]
30%|### | 25/82 [00:00<00:00, 59.65it/s]
39%|###9 | 32/82 [00:00<00:00, 61.75it/s]
48%|####7 | 39/82 [00:00<00:00, 62.39it/s]
56%|#####6 | 46/82 [00:00<00:00, 63.06it/s]
65%|######4 | 53/82 [00:00<00:00, 63.45it/s]
74%|#######4 | 61/82 [00:00<00:00, 66.37it/s]
83%|########2 | 68/82 [00:01<00:00, 65.14it/s]
91%|#########1| 75/82 [00:01<00:00, 66.12it/s]
100%|##########| 82/82 [00:01<00:00, 64.06it/s]
100%|##########| 82/82 [00:01<00:00, 62.90it/s]predicciones.head() pred
2012-09-01 00:00:00 108.331027
2012-09-01 01:00:00 68.562982
2012-09-01 02:00:00 33.499525
2012-09-01 03:00:00 10.027583
2012-09-01 04:00:00 3.037563# Error de backtest
# ==============================================================================
metrica mean_absolute_error
0 76.464247El error de predicción del modelo autorregresivo alcanza un MAE inferior al del modelo de referencia.
Variables exógenas
Hasta ahora, sólo se han utilizado como predictores los valores pasados (lags) de la serie temporal. Sin embargo, es posible incluir otras variables como predictores. Estas variables se conocen como variables exógenas (features) y su uso puede mejorar la capacidad predictiva del modelo. Un punto muy importante que hay que tener en cuenta es que los valores de las variables exógenas deben conocerse en el momento de la predicción.
Ejemplos habituales de variables exógenas son aquellas obtenidas del calendario, como el día de la semana, el mes, el año o los días festivos. Las variables meteorológicas como la temperatura, la humedad y el viento también entran en esta categoría, al igual que las variables económicas como la inflación y los tipos de interés.
Variables de calendario y meteorológicas
A continuación, se crean variables exógenas basadas en información del calendario, las horas de salida y puesta del sol, la temperatura y los días festivos. Estas nuevas variables se añaden a los conjuntos de entrenamiento, validación y test, y se utilizan como predictores en el modelo autorregresivo.
# Variables basadas en el calendario
# ==============================================================================
features_to_extract = [
'month',
'week',
'day_of_week',
'hour'
]
calendar_transformer = DatetimeFeatures(
variables = 'index',
features_to_extract = features_to_extract,
drop_original = True,
)
variables_calendario = calendar_transformer.fit_transform(datos)[features_to_extract]
# Variables basadas en la luz solar
# ==============================================================================
location = LocationInfo(
name = 'Washington DC',
region = 'USA',
timezone = 'US/Eastern',
latitude = 40.516666666666666,
longitude = -77.03333333333333
)
sunrise_hour = [
sun(location.observer, date=date, tzinfo=location.timezone)['sunrise'].hour
for date in datos.index
]
sunset_hour = [
sun(location.observer, date=date, tzinfo=location.timezone)['sunset'].hour
for date in datos.index
]
variables_solares = pd.DataFrame({
'sunrise_hour': sunrise_hour,
'sunset_hour': sunset_hour},
index = datos.index
)
variables_solares['daylight_hours'] = (
variables_solares['sunset_hour'] - variables_solares['sunrise_hour']
)
variables_solares["is_daylight"] = np.where(
(datos.index.hour >= variables_solares["sunrise_hour"])
& (datos.index.hour < variables_solares["sunset_hour"]),
1,
0,
)
# Variables basadas en festivos
# ==============================================================================
variables_festivos = datos[['holiday']].astype(int)
variables_festivos['holiday_previous_day'] = variables_festivos['holiday'].shift(24)
variables_festivos['holiday_next_day'] = variables_festivos['holiday'].shift(-24)
# Variables basadas en temperatura
# ==============================================================================
wf_transformer = WindowFeatures(
variables = ["temp"],
window = ["1D", "7D"],
functions = ["mean", "max", "min"],
freq = "h",
)
variables_temp = wf_transformer.fit_transform(datos[['temp']])
# Unión de variables exógenas
# ==============================================================================
assert all(variables_calendario.index == variables_solares.index)
assert all(variables_calendario.index == variables_festivos.index)
assert all(variables_calendario.index == variables_temp.index)
variables_exogenas = pd.concat([
variables_calendario,
variables_solares,
variables_temp,
variables_festivos
], axis=1)
# Debido a la creación de medias móviles, hay valores faltantes al principio
# de la serie. Y debido a holiday_next_day hay valores faltantes al final.
variables_exogenas = variables_exogenas.iloc[7 * 24:, :]
variables_exogenas = variables_exogenas.iloc[:-24, :]
variables_exogenas.head(3) month week ... holiday_previous_day holiday_next_day
date_time ...
2011-01-08 00:00:00 1 1 ... 0.0 0.0
2011-01-08 01:00:00 1 1 ... 0.0 0.0
2011-01-08 02:00:00 1 1 ... 0.0 0.0
[3 rows x 18 columns]Variables con patrones cíclicos
Algunos aspectos del calendario, como las horas o los días, son cíclicos. Por ejemplo, la hora del día va de 0 a 23 horas. Este tipo de variables pueden tratarse de varias formas, cada una con sus ventajas e inconvenientes.
Un enfoque consiste en utilizar las variables directamente como valores numéricos sin ninguna transformación. Este método evita crear variables nuevas, pero puede imponer un orden lineal incorrecto a los valores. Por ejemplo, la hora 23 de un día y la hora 00 del siguiente están muy alejadas en su representación lineal, cuando en realidad sólo hay una hora de diferencia entre ellas.
Otra posibilidad es tratar las variables cíclicas como variables categóricas para evitar imponer un orden lineal. Sin embargo, este enfoque puede provocar la pérdida de la información cíclica inherente a la variable.
Existe una tercera forma de tratar las variables cíclicas que suele preferirse a los otros dos métodos. Se trata de transformar las variables utilizando el seno y el coseno de su periodo. Este método genera solo dos nuevas variables que captan la ciclicidad de los datos con mayor precisión que los dos métodos anteriores, ya que preserva el orden natural de la variable y evita imponer un orden lineal.
# Codificación cíclica de las variables de calendario y luz solar
# ==============================================================================
features_to_encode = [
"month",
"week",
"day_of_week",
"hour",
"sunrise_hour",
"sunset_hour",
]
max_values = {
"month": 12,
"week": 52,
"day_of_week": 7,
"hour": 24,
"sunrise_hour": 24,
"sunset_hour": 24,
}
cyclical_encoder = CyclicalFeatures(
variables = features_to_encode,
max_values = max_values,
drop_original = False
)
variables_exogenas = cyclical_encoder.fit_transform(variables_exogenas)
variables_exogenas.head(3) month week ... sunset_hour_sin sunset_hour_cos
date_time ...
2011-01-08 00:00:00 1 1 ... -0.866025 -0.5
2011-01-08 01:00:00 1 1 ... -0.866025 -0.5
2011-01-08 02:00:00 1 1 ... -0.866025 -0.5
[3 rows x 30 columns]Interacción entre variables
En muchos casos, las variables exógenas no son independientes. Más bien, su efecto sobre la variable objetivo depende del valor de otras variables. Por ejemplo, el efecto de la temperatura en la demanda de electricidad depende de la hora del día. La interacción entre las variables exógenas puede captarse mediante nuevas variables que se obtienen multiplicando entre sí las variables existentes. Estas interacciones se obtienen fácilmente con la clase PolynomialFeatures de scikit-learn.
# Interacción entre variables exógenas
# ==============================================================================
transformer_poly = PolynomialFeatures(
degree = 2,
interaction_only = True,
include_bias = False
).set_output(transform="pandas")
poly_cols = [
'month_sin',
'month_cos',
'week_sin',
'week_cos',
'day_of_week_sin',
'day_of_week_cos',
'hour_sin',
'hour_cos',
'sunrise_hour_sin',
'sunrise_hour_cos',
'sunset_hour_sin',
'sunset_hour_cos',
'daylight_hours',
'is_daylight',
'holiday_previous_day',
'holiday_next_day',
'temp_window_1D_mean',
'temp_window_1D_min',
'temp_window_1D_max',
'temp_window_7D_mean',
'temp_window_7D_min',
'temp_window_7D_max',
'temp',
'holiday'
]
variables_poly = transformer_poly.fit_transform(variables_exogenas[poly_cols])
variables_poly = variables_poly.drop(columns=poly_cols)
variables_poly.columns = [f"poly_{col}" for col in variables_poly.columns]
variables_poly.columns = variables_poly.columns.str.replace(" ", "__")
assert all(variables_exogenas.index == variables_poly.index)
variables_exogenas = pd.concat([variables_exogenas, variables_poly], axis=1)
variables_exogenas.head(3) month ... poly_temp__holiday
date_time ...
2011-01-08 00:00:00 1 ... 0.0
2011-01-08 01:00:00 1 ... 0.0
2011-01-08 02:00:00 1 ... 0.0
[3 rows x 306 columns]Variables categóricas
Existen varios enfoques para incorporar variables categóricas en LightGBM (y otras implementaciones de gradient boosting):
Una opción es transformar los datos convirtiendo los valores categóricos en valores numéricos utilizando métodos como la codificación one hot la codificación ordinal. Este enfoque es aplicable a todos los modelos de aprendizaje automático.
LightGBM puede manejar variables categóricas internamente sin necesidad de preprocesamiento. Esto puede hacerse automáticamente estableciendo el parámetro
categorical_features='auto'y almacenando las variables con el tipo de datocategorydentro de un Pandas DataFrame. Tambien es posible especificar los nombres de las variables a tratar como categóricas pasando una lista de nombres al parámetrocategorical_features.
No hay un método que sea siempre mejor que los otros. Las reglas generales son:
Cuando la cardinalidad de las variables categóricas es alta (muchos valores diferentes), es mejor utilizar el soporte nativo para variables categóricas que utilizar la codificación one-hot.
Con datos codificados con one hot encoding, se necesitan más puntos de división (es decir, más profundidad) para recuperar una división equivalente a la que podría obtenerse con un solo punto de división utilizando el tratamiento nativo.
Cuando una variable categórica se convierte en múltiples variables dummy utilizando one hot encoding, su importancia se diluye, haciendo que el análisis de la importancia de las características sea más complejo de interpretar.
# Almacenar las variables categoricas como tipo "category"
# ==============================================================================
datos["weather"] = datos["weather"].astype("category")ColumnTransformers en scikit-learn proporcionan una potente forma de definir transformaciones y aplicarlas a variables específicas. Al encapsular las transformaciones en un objeto ColumnTransformer, se puede pasar a un Forecaster utilizando el argumento transformer_exog.
# Transformación con codificación one-hot
# ==============================================================================
one_hot_encoder = make_column_transformer(
(
OneHotEncoder(sparse_output=False, drop='if_binary'),
make_column_selector(dtype_include=['category', 'object']),
),
remainder="passthrough",
verbose_feature_names_out=False,
).set_output(transform="pandas")
# Crear un forecaster con un transformer para las variables exógenas
# ==============================================================================
forecaster = ForecasterRecursive(
regressor = LGBMRegressor(random_state=15926, verbose=-1),
lags = 72,
transformer_exog = one_hot_encoder
)Para examinar cómo se transforman los datos, se puede utilizar el método create_train_X_y y generar las matrices que el forecaster utiliza para entrenar el modelo. Este método permite conocer las manipulaciones específicas de los datos que se producen durante el proceso de entrenamiento.
# Mostrar matrices de entrenamiento
# ==============================================================================
exog_cols = ['weather']
X_train, y_train = forecaster.create_train_X_y(
y = datos.loc[:fin_validacion, 'users'],
exog = datos.loc[:fin_validacion, exog_cols]
)
X_train.head(3) lag_1 lag_2 ... weather_mist weather_rain
date_time ...
2011-01-04 00:00:00 12.0 20.0 ... 0.0 0.0
2011-01-04 01:00:00 5.0 12.0 ... 0.0 0.0
2011-01-04 02:00:00 2.0 5.0 ... 0.0 0.0
[3 rows x 75 columns]La estrategia One Hot Encoder se ha mostrado con fines didácticos. Para el resto del documento, sin embargo, se utiliza el soporte nativo para variables categóricas.
Implementación nativa para variables categóricas
Las librerías de Gradient Boosting (XGBoost, LightGBM, CatBoost y HistGradientBoostingRegressor) asumen que las variables de entrada son números enteros empezando por 0 hasta el número de categorías [0, 1, …, n_categories-1]. En la mayoría de los casos reales, las variables categóricas no se codifican con números sino con cadenas (strings), por lo que es necesario un paso intermedio de transformación. Existen dos opciones:
Asignar a las columnas con variables categóricas el tipo
category. Internamente, esta estructura de datos consiste en una matriz de categorías y una matriz de valores enteros (códigos) que apuntan al valor real de la matriz de categorías. Es decir, internamente es una matriz numérica con un mapeo que relaciona cada valor con una categoría. Los modelos son capaces de identificar automáticamente las columnas de tipocategoryy acceder a sus códigos internos. Este enfoque es aplicable a XGBoost, LightGBM y CatBoost.Preprocesar las columnas categóricas con un
OrdinalEncoderpara transformar sus valores a enteros y luego indicar explícitamente que deben ser tratadas como categóricas.
Para utilizar la detección automática en skforecast, las variables categóricas deben codificarse primero como enteros y luego almacenarse de nuevo como tipo category. Esto se debe a que skforecast utiliza internamente una matriz numérica numpy para acelerar el cálculo.
# Transformación: codificación ordinal + conversión a tipo "category"
# ==============================================================================
pipeline_categorical = make_pipeline(
OrdinalEncoder(
dtype=int,
handle_unknown="use_encoded_value",
unknown_value=-1,
encoded_missing_value=-1
),
FunctionTransformer(
func=lambda x: x.astype('category'),
feature_names_out= 'one-to-one'
)
)
transformer_exog = make_column_transformer(
(
pipeline_categorical,
make_column_selector(dtype_include=['category', 'object']),
),
remainder="passthrough",
verbose_feature_names_out=False,
).set_output(transform="pandas")Cuando se utiliza el regresor LightGBM, se tiene que especificar cómo tratar las variables categóricas utilizando el argumento categorical_feature en el método fit(). Esto es debido a que el argumento categorical_feature sólo se puede especificar en el método fit() y no en la inicialización del regresor.
# Crear un forecaster con detección automática de variables categóricas (LGBMRegressor)
# ==============================================================================
forecaster = ForecasterRecursive(
regressor = LGBMRegressor(random_state=15926, verbose=-1),
lags = 72,
transformer_exog = transformer_exog,
fit_kwargs = {"categorical_feature": "auto"}
)Evaluar el modelo con variables exógenas
Se entrena de nuevo el forecaster, pero esta vez, las variables exógenas también se incluyen como predictores. Para las variables categóricas, se utiliza la implementación nativa.
# Selección de variables exógenas a incluir en el modelo
# ==============================================================================
exog_cols = []
# Columnas que terminan con _seno o _coseno son seleccionadas
exog_cols.extend(variables_exogenas.filter(regex='_sin$|_cos$').columns.tolist())
# Columnas que empiezan con tem_ son seleccionadas
exog_cols.extend(variables_exogenas.filter(regex='^temp_.*').columns.tolist())
# Columnas que empiezan con holiday_ son seleccionadas
exog_cols.extend(variables_exogenas.filter(regex='^holiday_.*').columns.tolist())
exog_cols.extend(['temp', 'holiday', 'weather'])
variables_exogenas = variables_exogenas.filter(exog_cols, axis=1)
# Combinar variables exógenas y target en el mismo dataframe
# ==============================================================================
datos = datos[['users', 'weather']].merge(
variables_exogenas,
left_index = True,
right_index = True,
how = 'inner' # Para utilizar solo fechas para las que hay datos exógenos
)
datos = datos.astype({col: np.float32 for col in datos.select_dtypes("number").columns})
datos_train = datos.loc[: fin_train, :].copy()
datos_val = datos.loc[fin_train:fin_validacion, :].copy()
datos_test = datos.loc[fin_validacion:, :].copy()
# Backtesting en los datos de test incluyendo las variables exógenas
# ==============================================================================
cv = TimeSeriesFold(steps = 36, initial_train_size = len(datos.loc[:fin_validacion]))
metrica, predicciones = backtesting_forecaster(
forecaster = forecaster,
y = datos['users'],
exog = datos[exog_cols],
cv = cv,
metric = 'mean_absolute_error'
)
0%| | 0/81 [00:00<?, ?it/s]
6%|6 | 5/81 [00:00<00:01, 39.81it/s]
11%|#1 | 9/81 [00:00<00:01, 38.27it/s]
16%|#6 | 13/81 [00:00<00:01, 35.19it/s]
21%|## | 17/81 [00:00<00:01, 34.62it/s]
26%|##5 | 21/81 [00:00<00:01, 36.05it/s]
31%|### | 25/81 [00:00<00:01, 36.83it/s]
37%|###7 | 30/81 [00:00<00:01, 37.63it/s]
43%|####3 | 35/81 [00:00<00:01, 37.77it/s]
48%|####8 | 39/81 [00:01<00:01, 37.89it/s]
53%|#####3 | 43/81 [00:01<00:01, 36.35it/s]
59%|#####9 | 48/81 [00:01<00:00, 37.57it/s]
65%|######5 | 53/81 [00:01<00:00, 38.61it/s]
72%|#######1 | 58/81 [00:01<00:00, 38.70it/s]
78%|#######7 | 63/81 [00:01<00:00, 38.84it/s]
83%|########2 | 67/81 [00:01<00:00, 37.20it/s]
88%|########7 | 71/81 [00:01<00:00, 37.86it/s]
93%|#########2| 75/81 [00:01<00:00, 38.07it/s]
98%|#########7| 79/81 [00:02<00:00, 37.01it/s]
100%|##########| 81/81 [00:02<00:00, 37.62it/s]metrica mean_absolute_error
0 48.350039La incorporación de variables exógenas mejora la capacidad predictiva del modelo.
Optimización de hiperparámetros
El ForecasterRecursive entrenado utiliza los primeros 24 lags y un modelo LGMBRegressor con los hiperparámetros por defecto. Sin embargo, no hay ninguna razón por la que estos valores sean los más adecuados. Para encontrar los mejores hiperparámetros, se realiza una Búsqueda Bayesiana con la función bayesian_search_forecaster(). La búsqueda se lleva a cabo utilizando el mismo proceso de backtesting que antes, pero cada vez, el modelo se entrena con diferentes combinaciones de hiperparámetros y lags. Es importante señalar que la búsqueda de hiperparámetros debe realizarse utilizando el conjunto de validación, nunca con los datos de test.
La búsqueda se realiza probando cada combinación de hiperparámetros y retardos del siguiente modo:
Entrenar el modelo utilizando sólo el conjunto de entrenamiento.
El modelo se evalúa utilizando el conjunto de validación mediante backtesting.
Seleccionar la combinación de hiperparámetros y retardos que proporcione el menor error.
- Volver a entrenar el modelo con la mejor combinación encontrada, esta vez utilizando tanto los datos de entrenamiento como los de validación.
Siguiendo estos pasos, se puede obtener un modelo con hiperparámetros optimizados y evitar el sobreajuste.
# Búsqueda de hiperparámetros
# ==============================================================================
forecaster = ForecasterRecursive(
regressor = LGBMRegressor(random_state=15926, verbose=-1),
lags = 72,
window_features = window_features,
transformer_exog = transformer_exog,
fit_kwargs = {"categorical_feature": "auto"}
)
# Lags grid
lags_grid = [48, 72, (1, 2, 3, 23, 24, 25, 167, 168, 169)]
# Espacio de búsqueda de hiperparámetros
def search_space(trial):
search_space = {
'n_estimators' : trial.suggest_int('n_estimators', 300, 1000, step=100),
'max_depth' : trial.suggest_int('max_depth', 3, 10),
'min_data_in_leaf': trial.suggest_int('min_data_in_leaf', 25, 500),
'learning_rate' : trial.suggest_float('learning_rate', 0.01, 0.5),
'feature_fraction': trial.suggest_float('feature_fraction', 0.5, 1),
'max_bin' : trial.suggest_int('max_bin', 50, 250),
'reg_alpha' : trial.suggest_float('reg_alpha', 0, 1),
'reg_lambda' : trial.suggest_float('reg_lambda', 0, 1),
'lags' : trial.suggest_categorical('lags', lags_grid)
}
return search_space
# Particiones de entrenamiento y validación
cv_search = TimeSeriesFold(steps = 36, initial_train_size = len(datos_train))
results_search, frozen_trial = bayesian_search_forecaster(
forecaster = forecaster,
y = datos.loc[:fin_validacion, 'users'],
exog = datos.loc[:fin_validacion, exog_cols],
cv = cv_search,
search_space = search_space,
metric = 'mean_absolute_error',
n_trials = 20,
return_best = True
)
0%| | 0/20 [00:00<?, ?it/s]
Best trial: 0. Best value: 69.3725: 0%| | 0/20 [00:02<?, ?it/s]
Best trial: 0. Best value: 69.3725: 5%|5 | 1/20 [00:02<00:54, 2.86s/it]
Best trial: 1. Best value: 65.189: 5%|5 | 1/20 [00:05<00:54, 2.86s/it]
Best trial: 1. Best value: 65.189: 10%|# | 2/20 [00:05<00:52, 2.90s/it]
Best trial: 2. Best value: 63.5914: 10%|# | 2/20 [00:08<00:52, 2.90s/it]
Best trial: 2. Best value: 63.5914: 15%|#5 | 3/20 [00:08<00:48, 2.88s/it]
Best trial: 3. Best value: 58.8332: 15%|#5 | 3/20 [00:11<00:48, 2.88s/it]
Best trial: 3. Best value: 58.8332: 20%|## | 4/20 [00:11<00:44, 2.80s/it]
Best trial: 3. Best value: 58.8332: 20%|## | 4/20 [00:14<00:44, 2.80s/it]
Best trial: 3. Best value: 58.8332: 25%|##5 | 5/20 [00:14<00:42, 2.85s/it]
Best trial: 3. Best value: 58.8332: 25%|##5 | 5/20 [00:17<00:42, 2.85s/it]
Best trial: 3. Best value: 58.8332: 30%|### | 6/20 [00:17<00:39, 2.82s/it]
Best trial: 3. Best value: 58.8332: 30%|### | 6/20 [00:19<00:39, 2.82s/it]
Best trial: 3. Best value: 58.8332: 35%|###5 | 7/20 [00:19<00:36, 2.80s/it]
Best trial: 3. Best value: 58.8332: 35%|###5 | 7/20 [00:22<00:36, 2.80s/it]
Best trial: 3. Best value: 58.8332: 40%|#### | 8/20 [00:22<00:32, 2.74s/it]
Best trial: 3. Best value: 58.8332: 40%|#### | 8/20 [00:24<00:32, 2.74s/it]
Best trial: 3. Best value: 58.8332: 45%|####5 | 9/20 [00:24<00:29, 2.69s/it]
Best trial: 3. Best value: 58.8332: 45%|####5 | 9/20 [00:27<00:29, 2.69s/it]
Best trial: 3. Best value: 58.8332: 50%|##### | 10/20 [00:27<00:26, 2.70s/it]
Best trial: 3. Best value: 58.8332: 50%|##### | 10/20 [00:30<00:26, 2.70s/it]
Best trial: 3. Best value: 58.8332: 55%|#####5 | 11/20 [00:30<00:23, 2.60s/it]
Best trial: 3. Best value: 58.8332: 55%|#####5 | 11/20 [00:32<00:23, 2.60s/it]
Best trial: 3. Best value: 58.8332: 60%|###### | 12/20 [00:32<00:20, 2.55s/it]
Best trial: 3. Best value: 58.8332: 60%|###### | 12/20 [00:35<00:20, 2.55s/it]
Best trial: 3. Best value: 58.8332: 65%|######5 | 13/20 [00:35<00:18, 2.63s/it]
Best trial: 3. Best value: 58.8332: 65%|######5 | 13/20 [00:37<00:18, 2.63s/it]
Best trial: 3. Best value: 58.8332: 70%|####### | 14/20 [00:37<00:15, 2.54s/it]
Best trial: 3. Best value: 58.8332: 70%|####### | 14/20 [00:40<00:15, 2.54s/it]
Best trial: 3. Best value: 58.8332: 75%|#######5 | 15/20 [00:40<00:12, 2.57s/it]
Best trial: 3. Best value: 58.8332: 75%|#######5 | 15/20 [00:42<00:12, 2.57s/it]
Best trial: 3. Best value: 58.8332: 80%|######## | 16/20 [00:42<00:10, 2.54s/it]
Best trial: 16. Best value: 57.4977: 80%|######## | 16/20 [00:45<00:10, 2.54s/it]
Best trial: 16. Best value: 57.4977: 85%|########5 | 17/20 [00:45<00:07, 2.59s/it]
Best trial: 16. Best value: 57.4977: 85%|########5 | 17/20 [00:48<00:07, 2.59s/it]
Best trial: 16. Best value: 57.4977: 90%|######### | 18/20 [00:48<00:05, 2.60s/it]
Best trial: 16. Best value: 57.4977: 90%|######### | 18/20 [00:50<00:05, 2.60s/it]
Best trial: 16. Best value: 57.4977: 95%|#########5| 19/20 [00:50<00:02, 2.58s/it]
Best trial: 19. Best value: 55.0165: 95%|#########5| 19/20 [00:53<00:02, 2.58s/it]
Best trial: 19. Best value: 55.0165: 100%|##########| 20/20 [00:53<00:00, 2.56s/it]
Best trial: 19. Best value: 55.0165: 100%|##########| 20/20 [00:53<00:00, 2.66s/it]
`Forecaster` refitted using the best-found lags and parameters, and the whole data set:
Lags: [ 1 2 3 23 24 25 167 168 169]
Parameters: {'n_estimators': 300, 'max_depth': 8, 'min_data_in_leaf': 84, 'learning_rate': 0.03208367008069202, 'feature_fraction': 0.6966856646946507, 'max_bin': 141, 'reg_alpha': 0.22914980380651362, 'reg_lambda': 0.15917712164349596}
Backtesting metric: 55.01654818651542best_params = results_search['params'].iat[0]
best_params = best_params | {'random_state': 15926, 'verbose': -1}
best_lags = results_search['lags'].iat[0]
# Search results
# ==============================================================================
results_search.head(3) lags ... reg_lambda
0 [1, 2, 3, 23, 24, 25, 167, 168, 169] ... 0.159177
1 [1, 2, 3, 23, 24, 25, 167, 168, 169] ... 0.014016
2 [1, 2, 3, 23, 24, 25, 167, 168, 169] ... 0.222861
[3 rows x 11 columns]# Backtesting en los datos
# ==============================================================================
cv = TimeSeriesFold(steps = 36, initial_train_size = len(datos[:fin_validacion]))
metrica, predicciones = backtesting_forecaster(
forecaster = forecaster,
y = datos['users'],
exog = datos[exog_cols],
cv = cv,
metric = 'mean_absolute_error'
)
0%| | 0/81 [00:00<?, ?it/s]
5%|4 | 4/81 [00:00<00:01, 38.92it/s]
11%|#1 | 9/81 [00:00<00:01, 39.15it/s]
17%|#7 | 14/81 [00:00<00:02, 28.31it/s]
23%|##3 | 19/81 [00:00<00:01, 33.06it/s]
28%|##8 | 23/81 [00:00<00:01, 34.41it/s]
33%|###3 | 27/81 [00:00<00:01, 34.78it/s]
38%|###8 | 31/81 [00:00<00:01, 35.22it/s]
43%|####3 | 35/81 [00:01<00:01, 35.69it/s]
49%|####9 | 40/81 [00:01<00:01, 37.50it/s]
54%|#####4 | 44/81 [00:01<00:00, 37.86it/s]
59%|#####9 | 48/81 [00:01<00:00, 37.29it/s]
64%|######4 | 52/81 [00:01<00:00, 37.77it/s]
69%|######9 | 56/81 [00:01<00:00, 37.97it/s]
74%|#######4 | 60/81 [00:01<00:00, 37.14it/s]
79%|#######9 | 64/81 [00:01<00:00, 36.46it/s]
84%|########3 | 68/81 [00:01<00:00, 35.57it/s]
89%|########8 | 72/81 [00:02<00:00, 35.10it/s]
94%|#########3| 76/81 [00:02<00:00, 35.30it/s]
100%|##########| 81/81 [00:02<00:00, 36.73it/s]
100%|##########| 81/81 [00:02<00:00, 35.92it/s]metrica mean_absolute_error
0 46.937169Selección de predictores
La selección de predictores (feature selection) es el proceso de identificar un subconjunto de predictores relevantes para su uso en la creación del modelo. Es un paso importante en el proceso de machine learning, ya que puede ayudar a reducir el sobreajuste, mejorar la precisión del modelo y reducir el tiempo de entrenamiento. Dado que los regresores subyacentes de skforecast siguen la API de scikit-learn, es posible utilizar los métodos de selección de predictores disponibles en scikit-learn. Dos de los métodos más populares son Recursive Feature Elimination y Sequential Feature Selection.
# Crear forecaster
# ==============================================================================
regressor = LGBMRegressor(
n_estimators = 100,
max_depth = 5,
random_state = 15926,
verbose = -1
)
forecaster = ForecasterRecursive(
regressor = regressor,
lags = best_lags,
window_features = window_features,
transformer_exog = transformer_exog,
fit_kwargs = {"categorical_feature": "auto"}
)
# Eliminación recursiva de predictores con validación cruzada
# ==============================================================================
warnings.filterwarnings("ignore", message="X does not have valid feature names.*")
selector = RFECV(
estimator = regressor,
step = 1,
cv = 3,
)
lags_seleccionados, window_features_seleccionadas, exog_seleccionadas = select_features(
forecaster = forecaster,
selector = selector,
y = datos_train['users'],
exog = datos_train[exog_cols],
select_only = None,
force_inclusion = None,
subsample = 0.5,
random_state = 123,
verbose = True,
)Recursive feature elimination (RFECV)
-------------------------------------
Total number of records available: 11327
Total number of records used for feature selection: 5663
Number of features available: 99
Lags (n=9)
Window features (n=1)
Exog (n=89)
Number of features selected: 60
Lags (n=9) : [1, 2, 3, 23, 24, 25, 167, 168, 169]
Window features (n=1) : ['roll_mean_72']
Exog (n=50) : ['weather', 'month_sin', 'week_sin', 'week_cos', 'day_of_week_sin', 'hour_sin', 'hour_cos', 'poly_month_sin__week_cos', 'poly_month_sin__day_of_week_sin', 'poly_month_sin__day_of_week_cos', 'poly_month_sin__hour_sin', 'poly_month_sin__hour_cos', 'poly_month_cos__week_sin', 'poly_month_cos__day_of_week_sin', 'poly_month_cos__day_of_week_cos', 'poly_month_cos__hour_sin', 'poly_month_cos__hour_cos', 'poly_week_sin__day_of_week_sin', 'poly_week_sin__day_of_week_cos', 'poly_week_sin__hour_sin', 'poly_week_sin__hour_cos', 'poly_week_sin__sunrise_hour_cos', 'poly_week_cos__day_of_week_sin', 'poly_week_cos__day_of_week_cos', 'poly_week_cos__hour_sin', 'poly_week_cos__hour_cos', 'poly_week_cos__sunrise_hour_cos', 'poly_day_of_week_sin__day_of_week_cos', 'poly_day_of_week_sin__hour_sin', 'poly_day_of_week_sin__hour_cos', 'poly_day_of_week_sin__sunrise_hour_sin', 'poly_day_of_week_sin__sunset_hour_sin', 'poly_day_of_week_cos__hour_sin', 'poly_day_of_week_cos__hour_cos', 'poly_day_of_week_cos__sunset_hour_sin', 'poly_hour_sin__hour_cos', 'poly_hour_sin__sunrise_hour_sin', 'poly_hour_sin__sunrise_hour_cos', 'poly_hour_sin__sunset_hour_sin', 'poly_hour_sin__sunset_hour_cos', 'poly_hour_cos__sunrise_hour_sin', 'poly_hour_cos__sunrise_hour_cos', 'poly_hour_cos__sunset_hour_sin', 'temp_window_1D_mean', 'temp_window_1D_max', 'temp_window_1D_min', 'temp_window_7D_mean', 'temp_window_7D_min', 'temp', 'holiday']El RFECV de scikit-learn empieza entrenando un modelo con todos los predictores disponibles y calculando la importancia de cada uno en base a los atributos como coef_ o feature_importances_. A continuación, se elimina el predictor menos importante y se realiza una validación cruzada para calcular el rendimiento del modelo con los predictores restantes. Este proceso se repite hasta que la eliminación de predictores adicionales no mejora la metrica de rendimiento elegida o se alcanza el min_features_to_select.
El resultado final es un subconjunto de predictores que idealmente equilibra la simplicidad del modelo y su capacidad predictiva, determinada por el proceso de validación cruzada.
El forecater se entrena y evalúa de nuevo utilizando el conjunto de predictores seleccionados.
# Crear forecaster con los predictores seleccionados
# ==============================================================================
forecaster = ForecasterRecursive(
regressor = LGBMRegressor(**best_params),
lags = lags_seleccionados,
window_features = window_features,
transformer_exog = transformer_exog,
fit_kwargs = {"categorical_feature": "auto"}
)
# Backtesting con los predictores seleccionados y los datos de test
# ==============================================================================
cv = TimeSeriesFold(steps = 36, initial_train_size = len(datos[:fin_validacion]))
metrica_lgbm, predicciones = backtesting_forecaster(
forecaster = forecaster,
y = datos['users'],
exog = datos[exog_seleccionadas],
cv = cv,
metric = 'mean_absolute_error'
)
0%| | 0/81 [00:00<?, ?it/s]
5%|4 | 4/81 [00:00<00:01, 39.43it/s]
11%|#1 | 9/81 [00:00<00:01, 40.47it/s]
17%|#7 | 14/81 [00:00<00:01, 38.58it/s]
23%|##3 | 19/81 [00:00<00:01, 40.03it/s]
30%|##9 | 24/81 [00:00<00:01, 39.73it/s]
35%|###4 | 28/81 [00:00<00:01, 38.98it/s]
40%|###9 | 32/81 [00:00<00:01, 38.96it/s]
44%|####4 | 36/81 [00:00<00:01, 38.18it/s]
49%|####9 | 40/81 [00:01<00:01, 38.37it/s]
54%|#####4 | 44/81 [00:01<00:00, 38.00it/s]
59%|#####9 | 48/81 [00:01<00:00, 38.23it/s]
64%|######4 | 52/81 [00:01<00:00, 38.53it/s]
70%|####### | 57/81 [00:01<00:00, 41.26it/s]
77%|#######6 | 62/81 [00:01<00:00, 41.72it/s]
83%|########2 | 67/81 [00:01<00:00, 42.08it/s]
89%|########8 | 72/81 [00:01<00:00, 40.89it/s]
95%|#########5| 77/81 [00:01<00:00, 40.32it/s]
100%|##########| 81/81 [00:02<00:00, 39.75it/s]metrica_lgbm mean_absolute_error
0 46.751648El rendimiento del modelo sigue siendo similar al del modelo entrenado con todas las variables. Sin embargo, el modelo es ahora mucho más simple, lo que hará que sea más rápido de entrenar y menos propenso al sobreajuste. Para el resto del documento, el modelo se entrenará utilizando sólo las variables seleccionadas.
# Actualizar las variables exógenas utilizadas
# ==============================================================================
exog_cols = exog_seleccionadasForecasting probabilístico: intervalos de predicción
Un intervalo de predicción define el intervalo dentro del cual es de esperar que se encuentre el verdadero valor de la variable respuesta con una determinada probabilidad. Skforecast implementa varios métodos para el forecasting probabilístico:
El siguiente código muestra cómo generar intervalos de predicción utilizando el método coformal prediction.
Con la función
prediction_interval()se predicen los intervalos para los siguientes n step.Con la función
backtesting_forecaster()se predicen los intervalos de todo el conjunto de test.
El argumento interval se utiliza para especificar la probabilidad de cobertura deseada de los intervalos. En este caso, interval se establece en [5, 95], lo que significa una cobertura teórica del 90%.
# Crear y entrenar forecaster
# ==============================================================================
forecaster = ForecasterRecursive(
regressor = LGBMRegressor(**best_params),
lags = [1, 2, 3, 23, 24, 25, 167, 168, 169],
window_features = window_features,
transformer_exog = transformer_exog,
fit_kwargs = {"categorical_feature": "auto"},
binner_kwargs = {"n_bins": 5}
)
forecaster.fit(
y = datos.loc[:fin_train, 'users'],
exog = datos.loc[:fin_train, exog_cols],
store_in_sample_residuals = True
)
# Predicción de intervalos
# ==============================================================================
# Como el modelo ha sido entrenado con variables exógenas, se tienen que pasar
# para las predicciones.
predicciones = forecaster.predict_interval(
exog = datos.loc[fin_train:, exog_cols],
steps = 36,
interval = [5, 95],
method = 'conformal',
)
predicciones.head() pred lower_bound upper_bound
2012-05-01 00:00:00 23.714525 0.062383 47.366666
2012-05-01 01:00:00 8.618651 0.297713 16.939588
2012-05-01 02:00:00 5.222926 -3.098012 13.543864
2012-05-01 03:00:00 5.053190 -3.267748 13.374128
2012-05-01 04:00:00 6.758192 -1.562746 15.079130Por defecto, los intervalos se calculan utilizando los residuos in-sample (residuos del conjunto de entrenamiento). Sin embargo, esto puede dar lugar a intervalos demasiado estrechos (demasiado optimistas). Para evitarlo, se utiliza el método set_out_sample_residuals() para almacenar residuos out-sample calculados mediante backtesting con un conjunto de validación.
Si además de los residuos, se le pasan las correspondientes predicciones al método set_out_sample_residuals(), entonces los residuos utilizados en el proceso de bootstrapping se seleccionan condicionados al rango de valores de las predicciones. Esto puede ayudar a conseguir intervalos con mayor covertura a la vez que se mantienen lo más estrechos posibles.
# Backtesting con los datos de validación para obtener residuos out-sample
# ==============================================================================
cv = TimeSeriesFold(steps = 36, initial_train_size = len(datos.loc[:fin_train]))
_, predicciones_val = backtesting_forecaster(
forecaster = forecaster,
y = datos.loc[:fin_validacion, 'users'],
exog = datos.loc[:fin_validacion, exog_cols],
cv = cv,
metric = 'mean_absolute_error',
)
0%| | 0/82 [00:00<?, ?it/s]
6%|6 | 5/82 [00:00<00:01, 40.56it/s]
12%|#2 | 10/82 [00:00<00:01, 40.05it/s]
18%|#8 | 15/82 [00:00<00:01, 41.60it/s]
24%|##4 | 20/82 [00:00<00:01, 40.78it/s]
30%|### | 25/82 [00:00<00:01, 43.12it/s]
37%|###6 | 30/82 [00:00<00:01, 42.53it/s]
43%|####2 | 35/82 [00:00<00:01, 41.63it/s]
49%|####8 | 40/82 [00:00<00:01, 40.75it/s]
55%|#####4 | 45/82 [00:01<00:00, 39.39it/s]
61%|###### | 50/82 [00:01<00:00, 40.16it/s]
67%|######7 | 55/82 [00:01<00:00, 41.41it/s]
73%|#######3 | 60/82 [00:01<00:00, 40.69it/s]
79%|#######9 | 65/82 [00:01<00:00, 40.32it/s]
85%|########5 | 70/82 [00:01<00:00, 40.69it/s]
91%|#########1| 75/82 [00:01<00:00, 40.79it/s]
98%|#########7| 80/82 [00:01<00:00, 40.89it/s]
100%|##########| 82/82 [00:02<00:00, 40.92it/s]# Distribución de los residuos out-sample
# ==============================================================================
residuals = datos.loc[predicciones_val.index, 'users'] - predicciones_val['pred']
print(pd.Series(np.where(residuals < 0, 'negative', 'positive')).value_counts())positive 1721
negative 1231
Name: count, dtype: int64plt.rcParams.update({'font.size': 8})
_ = plot_residuals(
y_true = datos.loc[predicciones_val.index, 'users'],
y_pred = predicciones_val['pred'],
figsize=(7, 4)
)
# Almacenar residuos out-sample en el forecaster
# ==============================================================================
forecaster.set_out_sample_residuals(
y_true = datos.loc[predicciones_val.index, 'users'],
y_pred = predicciones_val['pred']
)
A continuación, se ejecuta el proceso de backtesting para estimar los intervalos de predicción en el conjunto de test. Se indica el argumento use_in_sample_residuals en False para que se utilicen los residuos out-sample almacenados previamente y use_binned_residuals para que el los residuos utilizados en se seleccionen condicionados al rango de valores de las predicciones.
# Backtesting con intervalos de predicción en test usando out-sample residuals
# ==============================================================================
cv = TimeSeriesFold(steps = 36, initial_train_size = len(datos.loc[:fin_validacion]))
metrica, predicciones = backtesting_forecaster(
forecaster = forecaster,
y = datos['users'],
exog = datos[exog_cols],
cv = cv,
metric = 'mean_absolute_error',
interval = [5, 95], # 90% intervalo de predicción
interval_method = 'conformal',
use_in_sample_residuals = False, # Usar out-sample residuals
use_binned_residuals = True, # Usar residuos condicionados al valor predicho
)
predicciones.head(5)
# Cobertura del intervalo en los datos de test
# ==============================================================================Explicabilidad del modelo
Debido a la naturaleza compleja de muchos de los actuales modelos de machine learning, a menudo funcionan como cajas negras, lo que dificulta entender por qué han hecho una predicción u otra. Las técnicas de explicabilidad pretenden desmitificar estos modelos, proporcionando información sobre su funcionamiento interno y ayudando a generar confianza, mejorar la transparencia y cumplir los requisitos normativos en diversos ámbitos. Mejorar la explicabilidad de los modelos no sólo ayuda a comprender su comportamiento, sino también a identificar sesgos, mejorar su rendimiento y permitir a las partes interesadas tomar decisiones más informadas basadas en los conocimientos del machine learning
Skforecast es compatible con algunos de los métodos de explicabilidad más populares: model-specific feature importances, SHAP values, and partial dependence plots.
# Crear y entrenar el forecaster
# ==============================================================================
forecaster = ForecasterRecursive(
regressor = LGBMRegressor(**best_params),
lags = [1, 2, 3, 23, 24, 25, 167, 168, 169],
window_features = window_features,
transformer_exog = transformer_exog,
fit_kwargs = {"categorical_feature": "auto"}
)
forecaster.fit(
y = datos.loc[:fin_validacion, 'users'],
exog = datos.loc[:fin_validacion, exog_cols]
)
# Extraer importancia de los predictores
# ==============================================================================
importancia = forecaster.get_feature_importances()
importancia.head(10)Shap values
Los valores SHAP (SHapley Additive exPlanations) son un método muy utilizado para explicar los modelos de machine learning, ya que ayudan a comprender cómo influyen las variables y los valores en las predicciones de forma visual y cuantitativa.
Se puede obtener un análisis SHAP a partir de modelos skforecast con sólo dos elementos:
El regresor interno del forecaster.
Las matrices de entrenamiento creadas a partir de la serie temporal y variables exógenas, utilizadas para ajustar el pronosticador.
Aprovechando estos dos componentes, los usuarios pueden crear explicaciones interpretables para sus modelos de skforecast. Estas explicaciones pueden utilizarse para verificar la fiabilidad del modelo, identificar los factores más significativos que contribuyen a las predicciones y comprender mejor la relación subyacente entre las variables de entrada y la variable objetivo.
# Matrices de entrenamiento utilizadas por el forecaster para entrenar el regresor
# ==============================================================================
X_train, y_train = forecaster.create_train_X_y(
y = datos.loc[:fin_validacion, 'users'],
exog = datos.loc[:fin_validacion, exog_cols]
)
X_train.head(3)
y_train.head(3)
# Crear SHAP explainer (para modelos basados en árboles)
# ==============================================================================
explainer = shap.TreeExplainer(forecaster.regressor)
# Se selecciona una muestra del 50% de los datos para acelerar el cálculo
rng = np.random.default_rng(seed=785412)
sample = rng.choice(X_train.index, size=int(len(X_train)*0.5), replace=False)
X_train_sample = X_train.loc[sample, :]
shap_values = explainer.shap_values(X_train_sample)
# Backtesting indicando que se devuelvan los predictores
# ==============================================================================
cv = TimeSeriesFold(steps = 36, initial_train_size = len(datos.loc[:fin_validacion]))
metrica, predicciones = backtesting_forecaster(
forecaster = forecaster,
y = datos['users'],
exog = datos[exog_cols],
cv = cv,
metric = 'mean_absolute_error',
return_predictors = True,
)
predicciones.head(3)
# Waterfall para una predicción concreta
# ==============================================================================
predicciones = predicciones.astype(datos[exog_cols].dtypes) # Asegurar tipos
iloc_predicted_date = predicciones.index.get_loc('2012-10-06 12:00:00')
shap_values_single = explainer(predicciones.iloc[:, 2:])
shap.plots.waterfall(shap_values_single[iloc_predicted_date], show=False)
fig, ax = plt.gcf(), plt.gca()
fig.set_size_inches(8, 3.5)
ax_list = fig.axes
ax = ax_list[0]
ax.tick_params(labelsize=8)
ax.set
plt.show()
# Force plot para una predicción concreta
# ==============================================================================
shap.force_plot(
base_value = shap_values_single.base_values[iloc_predicted_date],
shap_values = shap_values_single.values[iloc_predicted_date],
features = predicciones.iloc[iloc_predicted_date, 2:],
)XGBoost, CatBoost, HistGradientBoostingRegressor
Desde el éxito del Gradient Boosting como algoritmo de machine learning, se han desarrollado varias implementaciones. Además de LightGBM, otras tres muy populares son: XGBoost, CatBoost e HistGradientBoostingRegressor. Todas ellas son compatibles con skforecast.
Las siguientes secciones muestran cómo utilizar estas implementaciones para crear modelos de forecasting, haciendo hincapié en el uso de su soporte nativo para características categóricas. Esta vez se utiliza una estrategia de validación one-step-ahead para acelerar la búsqueda de hiperparámetros.
# Particiones utilizadas para la búsqueda de hiperparámetros y backtesting
# ==============================================================================
cv_search = OneStepAheadFold(initial_train_size = len(datos_train))
cv_backtesting = TimeSeriesFold(steps = 36, initial_train_size = len(datos[:fin_validacion]))XGBoost
# Encoding ordinal + conversión a tipo category
# ==============================================================================
pipeline_categorical = make_pipeline(
OrdinalEncoder(
dtype=int,
handle_unknown="use_encoded_value",
unknown_value=-1,
encoded_missing_value=-1
),
FunctionTransformer(
func=lambda x: x.astype('category'),
feature_names_out= 'one-to-one'
)
)
transformer_exog = make_column_transformer(
(
pipeline_categorical,
make_column_selector(dtype_exclude=np.number)
),
remainder="passthrough",
verbose_feature_names_out=False,
).set_output(transform="pandas")# Crear forecaster
# ==============================================================================
forecaster = ForecasterRecursive(
regressor = XGBRegressor(tree_method='hist', enable_categorical=True, random_state=123),
lags = 24,
window_features = window_features,
transformer_exog = transformer_exog
)
# Búsqueda de hiperparámetros
# ==============================================================================
lags_grid = [48, 72, [1, 2, 3, 23, 24, 25, 167, 168, 169]]
def search_space(trial):
search_space = {
'n_estimators' : trial.suggest_int('n_estimators', 300, 1000, step=100),
'max_depth' : trial.suggest_int('max_depth', 3, 10),
'learning_rate' : trial.suggest_float('learning_rate', 0.01, 1),
'subsample' : trial.suggest_float('subsample', 0.1, 1),
'colsample_bytree': trial.suggest_float('colsample_bytree', 0.1, 1),
'gamma' : trial.suggest_float('gamma', 0, 1),
'reg_alpha' : trial.suggest_float('reg_alpha', 0, 1),
'reg_lambda' : trial.suggest_float('reg_lambda', 0, 1),
'lags' : trial.suggest_categorical('lags', lags_grid)
}
return search_space
results_search, frozen_trial = bayesian_search_forecaster(
forecaster = forecaster,
y = datos.loc[:fin_validacion, 'users'],
exog = datos.loc[:fin_validacion, exog_cols],
cv = cv_search,
search_space = search_space,
metric = 'mean_absolute_error',
n_trials = 20
)
# Backtesting con datos de test
# ==============================================================================
metrica_xgboost, predicciones = backtesting_forecaster(
forecaster = forecaster,
y = datos['users'],
exog = datos[exog_cols],
cv = cv_backtesting,
metric = 'mean_absolute_error'
)
metrica_xgboostHistGradientBoostingRegressor
Al crear un forecaster utilizando HistogramGradientBoosting, los nombres de las columnas categóricas deben especificarse durante la instanciación del regresor pasándolos como una lista al argumento categorical_feature.
# Codificación ordinal
# ==============================================================================
# Se utiliza un ColumnTransformer para transformar variables categóricas
# (no numéricas) utilizando codificación ordinal. Las variables numéricas
# no se modifican. Los valores missing se codifican como -1. Si se encuentra una
# nueva categoría en el conjunto de test, se codifica como -1.
categorical_features = ['weather']
transformer_exog = make_column_transformer(
(
OrdinalEncoder(
dtype=int,
handle_unknown="use_encoded_value",
unknown_value=-1,
encoded_missing_value=-1
),
categorical_features
),
remainder="passthrough",
verbose_feature_names_out=False,
).set_output(transform="pandas")
# Crear forecaster
# ==============================================================================
forecaster = ForecasterRecursive(
regressor = HistGradientBoostingRegressor(
categorical_features=categorical_features,
random_state=123
),
lags = 24,
window_features = window_features,
transformer_exog = transformer_exog
)
# Busqueda de hiperparámetros
# ==============================================================================
lags_grid = [48, 72, [1, 2, 3, 23, 24, 25, 167, 168, 169]]
def search_space(trial):
search_space = {
'max_iter' : trial.suggest_int('max_iter', 300, 1000, step=100),
'max_depth' : trial.suggest_int('max_depth', 3, 10),
'learning_rate' : trial.suggest_float('learning_rate', 0.01, 1),
'min_samples_leaf' : trial.suggest_int('min_samples_leaf', 1, 20),
'l2_regularization' : trial.suggest_float('l2_regularization', 0, 1),
'lags' : trial.suggest_categorical('lags', lags_grid)
}
return search_space
results_search, frozen_trial = bayesian_search_forecaster(
forecaster = forecaster,
y = datos.loc[:fin_validacion, 'users'],
exog = datos.loc[:fin_validacion, exog_cols],
cv = cv_search,
search_space = search_space,
metric = 'mean_absolute_error',
n_trials = 20
)
# Backtesting con datos de test
# ==============================================================================
metrica_histgb, predicciones = backtesting_forecaster(
forecaster = forecaster,
y = datos['users'],
exog = datos[exog_cols],
cv = cv_backtesting,
metric = 'mean_absolute_error'
)
metrica_histgbCatBoost
Desafortunadamente, la versión actual de skforecast no es compatible con el manejo nativo de varaibles categóricas de CatBoost. El problema surge porque CatBoost sólo acepta variables categóricas codificadas como enteros, mientras que skforecast convierte los datos de entrada a float para utilizar matrices numpy y así acelerar el proceso. Para evitar esta limitación, es necesario aplicar one-hot encoding o label encoding a las variables categóricas antes de utilizarlas con CatBoost.
# One hot encoding
# ==============================================================================
one_hot_encoder = make_column_transformer(
(
OneHotEncoder(sparse_output=False, drop='if_binary'),
make_column_selector(dtype_exclude=np.number),
),
remainder="passthrough",
verbose_feature_names_out=False,
).set_output(transform="pandas")
# Crear forecaster
# ==============================================================================
forecaster = ForecasterRecursive(
regressor = CatBoostRegressor(
random_state=123,
silent=True,
allow_writing_files=False,
boosting_type = 'Plain', # Faster training
leaf_estimation_iterations = 3, # Faster training
),
lags = 24,
window_features = window_features,
transformer_exog = one_hot_encoder
)
# Búsqueda de hiperparámetros
# ==============================================================================
lags_grid = [48, 72, [1, 2, 3, 23, 24, 25, 167, 168, 169]]
def search_space(trial):
search_space = {
'n_estimators' : trial.suggest_int('n_estimators', 100, 1000, step=100),
'max_depth' : trial.suggest_int('max_depth', 3, 10),
'learning_rate' : trial.suggest_float('learning_rate', 0.01, 1),
'lags' : trial.suggest_categorical('lags', lags_grid)
}
return search_space
results_search, frozen_trial = bayesian_search_forecaster(
forecaster = forecaster,
y = datos.loc[:fin_validacion, 'users'],
exog = datos.loc[:fin_validacion, exog_cols],
cv = cv_search,
search_space = search_space,
metric = 'mean_absolute_error',
n_trials = 20
)
# Backtesting con datos de test
# ==============================================================================
metrica_catboost, predicciones = backtesting_forecaster(
forecaster = forecaster,
y = datos['users'],
exog = datos[exog_cols],
cv = cv_backtesting,
metric = 'mean_absolute_error'
)
metrica_catboostConclusión
Utilizar modelos Gradient Boosting en problemas de forecasting es muy sencillo gracias a las funcionalidades ofrecidas por skforecast.
Como se ha mostrado en este documento, la incorporación de variables exógenas como predictores puede mejorar en gran medida la capacidad predictiva del modelo.
Las variables categóricas pueden incluirse facilmente como variables exógenas sin necesidad de preprocesamiento. Esto es posible gracias al soporte nativo para variables categóricas ofrecido por LightGBM, XGBoost y HistGradientBoostingRegressor.
metricas = pd.concat(
[metrica_baseline, metrica_lgbm, metrica_xgboost, metrica_histgb, metrica_catboost],
axis=0,
)
metricas.index = [
"Baseline",
"LGBMRegressor",
"XGBRegressor",
"HistGradientBoostingRegressor",
"CatBoostRegressor",
]
metricas.round(2).sort_values(by="mean_absolute_error")