Si te gusta Skforecast , ayúdanos dándonos una estrella en GitHub! ⭐️
Forecasting series temporales con gradient boosting: Skforecast, XGBoost, LightGBM y CatBoost
Más sobre ciencia de datos: cienciadedatos.net
- Skforecast: forecasting series temporales con Python y Scikit Learn
- Predicción (forecasting) de la demanda eléctrica con Python
- Predicción (forecasting) de visitas a página web con machine learning
- Predicción del precio de Bitcoin con Python
- Intervalos de predicción en modelos de forecasting
- Multi-series forecasting
- Predicción de demanda intermitente
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.
✎ Nota
Los modelos de machine learning no siempre superan a los modelos propios del aprendizaje estadístico como AR, ARIMA o Exponential Smoothing. Cuál funciona mejor depende en gran medida de las características del caso de uso al que se apliquen. Consultar Modelos ARIMA y SARIMAX con python para aprender más sobre modelos estadísticos.✎ Nota
Otro ejemplo de cómo utilizar modelos de gradient boosting para forecasting puede encontrarse en el documento Predicción (forecasting) de la demanda energética con machine learning.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:
# Tratamiento de datos
# ==============================================================================
import numpy as np
import pandas as pd
from astral.sun import sun
from astral import LocationInfo
from skforecast.datasets import fetch_dataset
# Gráficos
# ==============================================================================
import matplotlib.pyplot as plt
from statsmodels.graphics.tsaplots import plot_acf
from statsmodels.graphics.tsaplots import plot_pacf
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)
plt.style.use('seaborn-v0_8-darkgrid')
# Modelado y Forecasting
# ==============================================================================
from xgboost import XGBRegressor
from lightgbm import LGBMRegressor
from catboost import CatBoostRegressor
from sklearn.ensemble import HistGradientBoostingRegressor
from sklearn.preprocessing import OneHotEncoder
from sklearn.preprocessing import OrdinalEncoder
from sklearn.preprocessing import FunctionTransformer
from sklearn.preprocessing import PolynomialFeatures
from sklearn.feature_selection import RFECV
from sklearn.pipeline import make_pipeline
from sklearn.compose import make_column_transformer
from sklearn.compose import make_column_selector
from skforecast.ForecasterBaseline import ForecasterEquivalentDate
from skforecast.ForecasterAutoreg import ForecasterAutoreg
from skforecast.model_selection import bayesian_search_forecaster
from skforecast.model_selection import backtesting_forecaster
import shap
shap.initjs()
# 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)
# 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')
datos = datos.asfreq('H')
datos = datos.sort_index()
datos.head()
Con el objetivo de poder entrenar los modelos, hacer búsqueda de los mejores hiperparámetros y evaluar su capacidad predictiva, se reparten los datos en tres conjuntos: entrenamiento, validación y test.
# Separación de datos en entrenamiento, validación y test
# ==============================================================================
fin_train = '2012-03-31 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()} (n={len(datos_train)})")
print(f"Fechas validación : {datos_val.index.min()} --- {datos_val.index.max()} (n={len(datos_val)})")
print(f"Fechas test : {datos_test.index.min()} --- {datos_test.index.max()} (n={len(datos_test)})")
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¶
Serie temporal completa
# Gráfico interactivo de la serie temporal
# ==============================================================================
fig = go.Figure()
fig.add_trace(go.Scatter(x=datos_train.index, y=datos_train['users'], mode='lines', name='Train'))
fig.add_trace(go.Scatter(x=datos_val.index, y=datos_val['users'], mode='lines', name='Validation'))
fig.add_trace(go.Scatter(x=datos_test.index, y=datos_test['users'], mode='lines', name='Test'))
fig.update_layout(
title = 'Número de usuarios de bicicletas',
xaxis_title="Fecha",
yaxis_title="Usiarios",
legend_title="Partición:",
width=800,
height=350,
margin=dict(l=20, r=20, t=35, b=20),
legend=dict(
orientation="h",
yanchor="top",
y=1,
xanchor="left",
x=0.001
)
)
#fig.update_xaxes(rangeslider_visible=True)
fig.show()
# Gráfico de la serie temporal con zoom
# ==============================================================================
zoom = ('2011-08-01 00:00:00','2011-08-15 00:00:00')
fig = plt.figure(figsize=(8, 4))
grid = plt.GridSpec(nrows=8, ncols=1, hspace=0.1, wspace=0)
main_ax = fig.add_subplot(grid[1:3, :])
zoom_ax = fig.add_subplot(grid[5:, :])
datos_train['users'].plot(ax=main_ax, label='train', alpha=0.5)
datos_val['users'].plot(ax=main_ax, label='validation', alpha=0.5)
datos_test['users'].plot(ax=main_ax, label='test', alpha=0.5)
min_y = min(datos['users'])
max_y = max(datos['users'])
main_ax.fill_between(zoom, min_y, max_y, facecolor='blue', alpha=0.5, zorder=0)
main_ax.set_xlabel('')
main_ax.legend(loc='lower center', ncol=3, bbox_to_anchor=(0.5, -0.8))
datos.loc[zoom[0]: zoom[1]]['users'].plot(ax=zoom_ax, color='blue', linewidth=1)
main_ax.set_title(f'Número de usuarios: {datos.index.min()}, {datos.index.max()}', fontsize=10)
zoom_ax.set_title(f'Número de usuarios: {zoom}', fontsize=10)
zoom_ax.set_xlabel('')
plt.subplots_adjust(hspace=1)
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, 5.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])
datos.groupby('month')['users'].median().plot(style='o-', linewidth=0.8, ax=axs[0])
axs[0].set_ylabel('Users')
axs[0].set_title('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])
datos.groupby('week_day')['users'].median().plot(style='o-', linewidth=0.8, ax=axs[1])
axs[1].set_ylabel('Users')
axs[1].set_title('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])
datos.groupby('hour_day')['users'].median().plot(style='o-', linewidth=0.8, ax=axs[2])
axs[2].set_ylabel('Users')
axs[2].set_title('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])
axs[3].set(
title = "Mean users during week",
xticks = [i * 24 for i in range(7)],
xticklabels = ["Mon", "Tue", "Wed", "Thu", "Fri", "Sat", "Sun"],
xlabel = "Día y hora",
ylabel = "Users"
)
fig.suptitle("Gráficos de estacionalidad", fontsize=16)
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)
plt.show()
# Gráfico autocorrelación parcial
# ==============================================================================
fig, ax = plt.subplots(figsize=(5, 2))
plot_pacf(datos['users'], ax=ax, lags=72, method='ywm')
plt.show()
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.
✎ Nota
En las siguientes celdas de código, se entrena un modelo baseline y se evalúa su capacidad predictiva mediante un proceso de backtesting. Si este concepto es nuevo para ti, no te preocupes, se explicará en detalle a lo largo del documento. Por ahora, basta con saber que el proceso de backtesting consiste en entrenar el modelo con una cierta cantidad de datos y evaluar su capacidad predictiva con los datos que el modelo no ha visto. La métrica del error se utilizará como referencia para comparar la capacidad predictiva de los modelos más complejos que se implementarán a lo largo del documento.# 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
# Backtesting
# ==============================================================================
metrica_baseline, predicciones = backtesting_forecaster(
forecaster = forecaster,
y = datos['users'],
steps = 36,
metric = 'mean_absolute_error',
initial_train_size = len(datos.loc[:fin_validacion]),
refit = False,
n_jobs = 'auto',
verbose = False,
show_progress = True
)
print(f"Error de backtest (MAE): {metrica_baseline}")
El modelo baseline alcanza un MAE de 91,7, que se utiliza como referencia para evaluar si merece la pena aplicar los modelos más complejos.
Modelo autoregresivo recursivo 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) 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
# ==============================================================================
forecaster = ForecasterAutoreg(
regressor = LGBMRegressor(random_state=15926, verbose=-1),
lags = 24
)
# Entrenar el forecaster
# ==============================================================================
forecaster.fit(y=datos.loc[:fin_validacion, 'users'])
forecaster
# Predicciones
# ==============================================================================
forecaster.predict(steps=10)
Backtesting¶
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.
# Backtest del modelo con lo datos de test
# ==============================================================================
metrica, predicciones = backtesting_forecaster(
forecaster = forecaster,
y = datos['users'],
steps = 36,
metric = 'mean_absolute_error',
initial_train_size = len(datos[:fin_validacion]),
refit = False,
n_jobs = 'auto',
verbose = True, # Cambiar a False para mostrar menos información
show_progress = True
)
predicciones.head()
# Error de backtest
# ==============================================================================
print(f'Error de backtest (MAE): {metrica}')
El error de predicción del modelo autorregresivo alcanza un MAE inferior al del modelo de referencia.
Optimización de hiperparámetros (tuning)¶
El ForecasterAutoreg 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
# ==============================================================================
# Lags candidatos
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', 400, 1200, step=100),
'max_depth' : trial.suggest_int('max_depth', 3, 10, step=1),
'learning_rate' : trial.suggest_float('learning_rate', 0.01, 0.5),
'reg_alpha' : trial.suggest_float('reg_alpha', 0, 1, step=0.1),
'reg_lambda' : trial.suggest_float('reg_lambda', 0, 1, step=0.1),
}
return search_space
resultados_busqueda, frozen_trial = bayesian_search_forecaster(
forecaster = forecaster,
y = datos.loc[:fin_validacion, 'users'], # Datos test no incluidos
search_space = search_space,
lags_grid = lags_grid,
steps = 36,
refit = False,
metric = 'mean_absolute_error',
initial_train_size = len(datos_train),
fixed_train_size = False,
n_trials = 20, # Aumentar para una búsqueda más exhaustiva
random_state = 123,
return_best = True,
n_jobs = 'auto',
verbose = False,
show_progress = True
)
# Resultados de la búsqueda
# ==============================================================================
resultados_busqueda.head(10)
Al indicar return_best = True, el objeto forecaster se actualiza automáticamente con la mejor configuración encontrada y se entrena con el conjunto de datos completo. Este modelo final puede utilizarse para obtener predicciones sobre nuevos datos.
# Mejor modelo
# ==============================================================================
forecaster
Backtesting con datos de test¶
Una vez identificada la mejor combinación de hiperparámetros utilizando los datos de validación, se evalúa la capacidad predictiva del modelo cuando se aplica al conjunto de test. 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 modelo final con datos de test
# ==============================================================================
metrica, predicciones = backtesting_forecaster(
forecaster = forecaster,
y = datos['users'],
steps = 36,
metric = 'mean_absolute_error',
initial_train_size = len(datos[:fin_validacion]),
refit = False,
n_jobs = 'auto',
verbose = False,
show_progress = True
)
print(f'Error de backtest (MAE): {metrica}')
predicciones.head()
# Gráfico predicciones vs valor real
# ======================================================================================
fig = go.Figure()
trace1 = go.Scatter(x=datos_test.index, y=datos_test['users'], name="test", mode="lines")
trace2 = go.Scatter(x=predicciones.index, y=predicciones['pred'], name="prediction", mode="lines")
fig.add_trace(trace1)
fig.add_trace(trace2)
fig.update_layout(
title="Predicción vs valor real en los datos de test",
xaxis_title="Fecha",
yaxis_title="Usiarios",
width=800,
height=350,
margin=dict(l=20, r=20, t=35, b=20),
legend=dict(
orientation="h",
yanchor="top",
y=1.1,
xanchor="left",
x=0.001
)
)
fig.show()
Tras optimizar los lags y los hiperparámetros, el error de predicción del modelo autorregresivo alcanza un MAE significativamente inferior al del modelo de referencia. A continuación, se evalúa la capacidad predictiva del modelo cuando se añaden variables exógenas.
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.
⚠ Warning
Las variables exógenas deben conocerse en el momento de la predicción. Por ejemplo, si se utiliza la temperatura como variable exógena, el valor de la temperatura para la hora siguiente debe conocerse en el momento de la previsión. Si no se conoce el valor de la temperatura, la predicción no será posible.Las variables meteorológicas deben utilizarse con precaución. Cuando el modelo se pone en producción, las condiciones meteorológicas futuras no se conocen, sino que son predicciones realizadas por los servicios meteorológicos. Al tratarse de predicciones, introducen errores en el modelo de previsión. Como consecuencia, es probable que las predicciones del modelo empeoren. Una forma de anticiparse a este problema, y conocer (no evitar) el rendimiento esperado del modelo, es utilizar las previsiones meteorológicas disponibles en el momento en que se entrena el modelo, en lugar de las condiciones reales registradas.
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
# ==============================================================================
variables_calendario = pd.DataFrame(index=datos.index)
variables_calendario['mes'] = variables_calendario.index.month
variables_calendario['semana_anyo'] = variables_calendario.index.isocalendar().week
variables_calendario['dia_semana'] = variables_calendario.index.day_of_week + 1
variables_calendario['hora_dia'] = variables_calendario.index.hour + 1
# Variables basadas en la luz solar
# ==============================================================================
location = LocationInfo(
name='Washington DC',
region='USA',
timezone='US/Eastern',
latitude=40.516666666666666,
longitude=-77.03333333333333
)
hora_amanecer = [
sun(location.observer, date=date, tzinfo=location.timezone)['sunrise'].hour
for date in datos.index
]
hora_anochecer = [
sun(location.observer, date=date, tzinfo=location.timezone)['sunset'].hour
for date in datos.index
]
variables_solares = pd.DataFrame({
'hora_amanecer': hora_amanecer,
'hora_anochecer': hora_anochecer},
index = datos.index
)
variables_solares['horas_luz_solar'] = (
variables_solares['hora_anochecer'] - variables_solares['hora_amanecer']
)
variables_solares['es_de_dia'] = np.where(
(datos.index.hour >= variables_solares['hora_amanecer']) & \
(datos.index.hour < variables_solares['hora_anochecer']),
1,
0
)
# Variables basadas en festivos
# ==============================================================================
variables_festivos = datos[['holiday']].astype(int)
variables_festivos['holiday_dia_anterior'] = variables_festivos['holiday'].shift(24)
variables_festivos['holiday_dia_siguiente'] = variables_festivos['holiday'].shift(-24)
# Variables basadas en temperatura
# ==============================================================================
variables_temp = datos[['temp']].copy()
variables_temp['temp_roll_mean_1_dia'] = variables_temp['temp'].rolling(24, closed='left').mean()
variables_temp['temp_roll_mean_7_dia'] = variables_temp['temp'].rolling(24*7, closed='left').mean()
variables_temp['temp_roll_max_1_dia'] = variables_temp['temp'].rolling(24, closed='left').max()
variables_temp['temp_roll_min_1_dia'] = variables_temp['temp'].rolling(24, closed='left').min()
variables_temp['temp_roll_max_7_dia'] = variables_temp['temp'].rolling(24*7, closed='left').max()
variables_temp['temp_roll_min_7_dia'] = variables_temp['temp'].rolling(24*7, closed='left').min()
# Merge all exogenous variables
# ==============================================================================
variables_exogenas = pd.concat([
variables_calendario,
variables_solares,
variables_temp,
variables_festivos
], axis=1)
variables_exogenas.head(4)
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.
✎ Nota
Para información más detallada sobre las diferentes formas de tratar las variables cíclicas, consultar el documento Cyclical features in time series forecasting.# Codificación cíclica de las variables de calendario y luz solar
# ==============================================================================
def codificacion_ciclica(datos: pd.Series, longitud_ciclo: int) -> pd.DataFrame:
"""
Codifica una variable cíclica con dos nuevas variables: seno y coseno.
Se asume que el valor mínimo de la variable es 0. El valor máximo de la
variable se pasa como argumento.
Parameters
----------
datos : pd.Series
Serie con la variable a codificar.
longitud_ciclo : int
La longitud del ciclo. Por ejemplo, 12 para meses, 24 para horas, etc.
Este valor se utiliza para calcular el ángulo del seno y coseno.
Returns
-------
resultado : pd.DataFrame
Dataframe con las dos nuevas características: seno y coseno.
"""
seno = np.sin(2 * np.pi * datos/longitud_ciclo)
coseno = np.cos(2 * np.pi * datos/longitud_ciclo)
resultado = pd.DataFrame({
f"{datos.name}_seno": seno,
f"{datos.name}_coseno": coseno
})
return resultado
mes_encoded = codificacion_ciclica(variables_exogenas['mes'], longitud_ciclo=12)
semana_anyo_encoded = codificacion_ciclica(variables_exogenas['semana_anyo'], longitud_ciclo=52)
dia_semana_encoded = codificacion_ciclica(variables_exogenas['dia_semana'], longitud_ciclo=7)
hora_dia_encoded = codificacion_ciclica(variables_exogenas['hora_dia'], longitud_ciclo=24)
hora_amanecer_encoded = codificacion_ciclica(variables_exogenas['hora_amanecer'], longitud_ciclo=24)
hora_anochecer_encoded = codificacion_ciclica(variables_exogenas['hora_anochecer'], longitud_ciclo=24)
variables_ciclicas = pd.concat([
mes_encoded,
semana_anyo_encoded,
dia_semana_encoded,
hora_dia_encoded,
hora_amanecer_encoded,
hora_anochecer_encoded
], axis=1)
variables_exogenas = pd.concat([variables_exogenas, variables_ciclicas], axis=1)
variables_exogenas.head(3)
✎ Nota
Con fines ilustrativos, las variables de calendario y sus transformaciones cíclicas se han creado manualmente. Sin embargo, para un enfoque más eficiente y automatizado, se recomienda utilizar las clases DatetimeFeatures y CyclicalFeatures disponibles en la librería Feature-engine y pasarlos directamente al argumentotransformer_exog del forecaster.
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 = [
'mes_seno',
'mes_coseno',
'semana_anyo_seno',
'semana_anyo_coseno',
'dia_semana_seno',
'dia_semana_coseno',
'hora_dia_seno',
'hora_dia_coseno',
'hora_amanecer_seno',
'hora_amanecer_coseno',
'hora_anochecer_seno',
'hora_anochecer_coseno',
'horas_luz_solar',
'es_de_dia',
'holiday_dia_anterior',
'holiday_dia_siguiente',
'temp_roll_mean_1_dia',
'temp_roll_mean_7_dia',
'temp_roll_max_1_dia',
'temp_roll_min_1_dia',
'temp_roll_max_7_dia',
'temp_roll_min_7_dia',
'temp',
'holiday'
]
variables_poly = transformer_poly.fit_transform(variables_exogenas[poly_cols].dropna())
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(" ", "__")
variables_exogenas = pd.concat([variables_exogenas, variables_poly], axis=1)
variables_exogenas.head(3)
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")
One hot encoding¶
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.
✎ Nota
Es posible aplicar una transformación a todo el conjunto de datos independientemente del forecaster. Sin embargo, es crucial asegurarse de que las transformaciones se aprenden sólo a partir de los datos de entrenamiento para evitar fugas de información (*leakage*). Además, la misma transformación debe aplicarse a los datos de entrada durante la predicción. Por lo tanto, es aconsejable incluir la transformación en el forecaster para que se gestione internamente. Este enfoque garantiza la coherencia en la aplicación de las transformaciones y reduce la probabilidad de errores.# 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 = ForecasterAutoreg(
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)
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.
⚠ Warning
Los cuatro librerías de gradient boosting - LightGBM, HistogramGradientBoosting de scikit-learn, XGBoost y CatBoost - son capaces de manejar directamente características categóricas dentro del modelo. Sin embargo, es importante señalar que cada una tiene sus propias configuraciones, ventajas y potenciales problemas. Para comprender plenamente cómo utilizar cada una, se recomienda consultar esta guía de usuario.Cuando se desplegan modelos en producción, se recomienda evitar la detección automática basada en columnas de tipo
category de pandas. Aunque pandas proporciona una codificación interna para estas columnas, no es consistente entre diferentes conjuntos de datos y puede variar dependiendo de las categorías presentes en cada conjunto de datos. Por lo tanto, es crucial ser consciente de este problema y tomar las medidas adecuadas para garantizar la consistencia en la codificación de las variables categóricas al desplegar modelos en producción. Se puede encontrar más información sobre este problema en Native implementation for categorical features.
# 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_exclude=np.number)
),
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 = ForecasterAutoreg(
regressor = LGBMRegressor(random_state=15926, verbose=-1),
lags = 72,
transformer_exog = transformer_exog,
fit_kwargs = {"categorical_feature": "auto"}
)
Esta es la estrategia que se utilizará en el resto del documento.
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='_seno$|_coseno$').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='left'
)
# Debido a la creación de medias móviles, hay valores NaN al principio de la serie.
# Debido a holiday_dia_siguiente hay valores NaN al final de la serie.
# Las columnas numéricas se convierten a float32.
datos = datos.dropna()
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()
Como ahora se incluyen muchas más variables en el modelo, es posible que los hiperparámetros identificados previamente ya no sean óptimos, por lo que se realiza una nueva búsqueda. Esta vez, sin embargo, la búsqueda se adapta al rango de valores identificados en la búsqueda anterior.
# Búsqueda de hiperparámetros
# ==============================================================================
forecaster = ForecasterAutoreg(
regressor = LGBMRegressor(random_state=15926, verbose=-1),
lags = 72,
transformer_exog = transformer_exog,
fit_kwargs = {"categorical_feature": "auto"}
)
# Lags grid
lags_grid = [[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', 800, 1400, step=100),
'max_depth' : trial.suggest_int('max_depth', 3, 8, step=1),
'learning_rate': trial.suggest_float('learning_rate', 0.01, 0.5),
'reg_alpha' : trial.suggest_float('reg_alpha', 0, 1, step=0.1),
'reg_lambda' : trial.suggest_float('reg_lambda', 0, 1, step=0.1),
}
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],
search_space = search_space,
lags_grid = lags_grid,
steps = 36,
refit = False,
metric = 'mean_absolute_error',
initial_train_size = len(datos_train),
fixed_train_size = False,
n_trials = 20,
random_state = 123,
return_best = True,
n_jobs = 'auto',
verbose = False,
show_progress = True
)
# Backtesting en los datos de test incluyendo las variables exógenas
# ==============================================================================
metrica, predictiones = backtesting_forecaster(
forecaster = forecaster,
y = datos['users'],
exog = datos[exog_cols],
steps = 36,
metric = 'mean_absolute_error',
initial_train_size = len(datos[:fin_validacion]),
refit = False,
n_jobs = 'auto',
verbose = False,
show_progress = True
)
print(f'Error de backtest (MAE): {metrica:.2f}")
# Gráfico predicciones vs valor real
# ======================================================================================
fig = go.Figure()
trace1 = go.Scatter(x=datos_test.index, y=datos_test['users'], name="test", mode="lines")
trace2 = go.Scatter(x=predicciones.index, y=predicciones['pred'], name="prediction", mode="lines")
fig.add_trace(trace1)
fig.add_trace(trace2)
fig.update_layout(
title="Predicción vs valor real en los datos de test",
xaxis_title="Fecha",
yaxis_title="Usiarios",
width=800,
height=350,
margin=dict(l=20, r=20, t=35, b=20),
legend=dict(
orientation="h",
yanchor="top",
y=1.1,
xanchor="left",
x=0.001
)
)
fig.show()
La incorporación de variables exógenas mejora aún más la capacidad predictiva del modelo.
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
# ==============================================================================
params = {
'n_estimators': 900,
'max_depth': 7,
'learning_rate': 0.022655,
'reg_alpha': 0.8,
'reg_lambda': 0.2,
'random_state': 15926,
'verbose': -1
}
forecaster = ForecasterAutoreg(
regressor = LGBMRegressor(**params),
lags = [1, 2, 3, 23, 24, 25, 167, 168, 169],
transformer_exog = transformer_exog,
fit_kwargs = {"categorical_feature": "auto"}
)
forecaster.fit(
y = datos.loc[:fin_validacion, 'users'],
exog = datos.loc[:fin_validacion, exog_cols]
)
Model-specific feature importance¶
# Extraer importancia de los predictores
# ==============================================================================
importancia = forecaster.get_feature_importances()
importancia.sort_values(by='importance', ascending=False).head(10)
⚠ Warning
El métodoget_feature_importances() sólo devuelve valores si el regresor del forecaster tiene el atributo coef_ o feature_importances_, que son los estándares en scikit-learn.
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]
)
display(X_train.head(3))
display(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)
✎ Note
La librería Shap cuenta con varios *Explainers*, cada uno diseñado para un tipo de modelo diferente. Elshap.TreeExplainer explainer se utiliza para modelos basados en árboles, como el LGBMRegressor utilizado en este ejemplo. Para más información, consultar la documentación de SHAP..
# Shap summary plot (top 10)
# ==============================================================================
shap.summary_plot(shap_values, X_train_sample, max_display=10, show=False)
fig, ax = plt.gcf(), plt.gca()
ax.set_title("SHAP Summary plot")
ax.tick_params(labelsize=8)
fig.set_size_inches(8, 4.5)
# Force plot para la primera observación de la muestra
# ==============================================================================
shap.force_plot(explainer.expected_value, shap_values[0,:], X_train_sample.iloc[0,:])
Selecció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.
💡 Tip
La selección de predictores es una herramienta potente para mejorar el rendimiento de los modelos de machine learning. Sin embargo, es computacionalmente costosa y puede requerir mucho tiempo. Dado que el objetivo es encontrar el mejor subconjunto de predictores, no el mejor modelo, no es necesario utilizar todos los datos disponibles ni un modelo muy complejo. En su lugar, se recomienda utilizar un pequeño subconjunto de los datos y un modelo sencillo. Una vez identificados los mejores predictores, el modelo puede entrenarse utilizando todo el conjunto de datos y una configuración más compleja.Por ejemplo, en este caso de uso, el mejor modelo es un LGMBRegressor con 1100 árboles y una profundidad máxima de 10. Sin embargo, para la selección predictores, sólo se utilizan 100 árboles y una profundidad máxima de 5.
# Crear forecaster
# ==============================================================================
forecaster = ForecasterAutoreg(
regressor = LGBMRegressor(random_state=15926, verbose=-1),
lags = [1, 2, 3, 23, 24, 25, 167, 168, 169],
transformer_exog = transformer_exog,
fit_kwargs = {"categorical_feature": "auto"}
)
# Matrices de entrenamiento utilizadas por el forecaster para ajustar el regresor interno
# ==============================================================================
X_train, y_train = forecaster.create_train_X_y(
y = datos_train['users'],
exog = datos_train[exog_cols]
)
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.
# Eliminación recursiva de predictores con validación cruzada
# ==============================================================================
# Muestrear el 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, :]
y_train_sample = y_train.loc[sample]
regressor = LGBMRegressor(
n_estimators = 100,
max_depth = 5,
random_state = 15926,
verbose = -1
)
selector = RFECV(
estimator = regressor,
step = 1,
cv = 3,
min_features_to_select = 25,
n_jobs = -1
)
selector.fit(X_train_sample, y_train_sample)
predictores_seleccionados = selector.get_feature_names_out()
print("Eliminación recursiva de predictores")
print("-------------------------------------")
print(f"Número de predictores disponibles: {X_train.shape[1]}")
print(f"Número de observaciones disponibles: {X_train.shape[0]}")
print(f"Número de observaciones utilizadas para selección de predictores: {X_train_sample.shape[0]}")
print(f"Número de predictores seleccionados: {selector.n_features_}")
print(f"Predictores seleccionados: \n {predictores_seleccionados}")
El forecater se entrena y evalúa de nuevo utilizando el conjunto de predictores seleccionados.
# Crear forecaster con los predictores seleccionados
# ==============================================================================
params = {
'n_estimators': 900,
'max_depth': 7,
'learning_rate': 0.022655,
'reg_alpha': 0.8,
'reg_lambda': 0.2,
'random_state': 15926,
'verbose': -1
}
lags_seleccionados = [
int(predictor.replace("lag_", ""))
for predictor in predictores_seleccionados
if predictor.startswith("lag_")
]
exog_seleccionadas = [
predictor
for predictor in predictores_seleccionados
if not predictor.startswith("lag_")
]
forecaster = ForecasterAutoreg(
regressor = LGBMRegressor(**params),
lags = lags_seleccionados,
transformer_exog = transformer_exog,
fit_kwargs = {"categorical_feature": "auto"}
)
# Backtesting con los predictores seleccionados y los datos de test
# ==============================================================================
metrica_lgbm, predicciones = backtesting_forecaster(
forecaster = forecaster,
y = datos['users'],
exog = datos[exog_seleccionadas],
steps = 24,
metric = 'mean_absolute_error',
initial_train_size = len(datos[:fin_validacion]),
refit = False,
n_jobs = 'auto',
verbose = False,
show_progress = True
)
print(f"Error de backtest (MAE): {metrica_lgbm:.2f}")
El rendimiento del modelo sigue siendo similar al del modelo entrenado con todas las características. 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 características más importantes.
# Actualizar las variables exógenas utilizadas
# ==============================================================================
exog_cols = exog_seleccionadas
XGBoost, CatBoost y 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.
XGBoost: desarrollado por Tianqi Chen.
HistGradientBoostingRegressor: desarrollado por scikit-learn.
CatBoost: desarrollado por Yandex.
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.
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 = ForecasterAutoreg(
regressor = XGBRegressor(
tree_method = 'hist',
enable_categorical = True,
random_state = 123
),
lags = 24,
transformer_exog = transformer_exog
)
# Búsqueda de hiperparámetros
# ==============================================================================
lags_grid = [[1, 2, 3, 23, 24, 25, 167, 168, 169]]
def search_space(trial):
search_space = {
'n_estimators' : trial.suggest_int('n_estimators', 400, 1200, step=100),
'max_depth' : trial.suggest_int('max_depth', 3, 10, step=1),
'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),
}
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],
search_space = search_space,
lags_grid = lags_grid,
steps = 36,
refit = False,
metric = 'mean_absolute_error',
initial_train_size = len(datos_train),
fixed_train_size = False,
n_trials = 20,
random_state = 123,
return_best = True,
n_jobs = 'auto',
verbose = False,
show_progress = True
)
# Backtesting con datos de test
# ==============================================================================
metrica_xgboost, predicciones = backtesting_forecaster(
forecaster = forecaster,
y = datos['users'],
exog = datos[exog_cols],
initial_train_size = len(datos.loc[:fin_validacion]),
fixed_train_size = False,
steps = 36,
refit = False,
metric = 'mean_absolute_error',
n_jobs = 'auto',
verbose = False
)
print(f"Error de backtest (MAE): {metrica_xgboost:.2f}")
HistGradientBoostingRegressor¶
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 = ForecasterAutoreg(
regressor = HistGradientBoostingRegressor(
categorical_features=categorical_features,
random_state=123
),
lags = 24,
transformer_exog = transformer_exog
)
# Busqueda de hiperparámetros
# ==============================================================================
lags_grid = [[1, 2, 3, 23, 24, 25, 167, 168, 169]]
def search_space(trial):
search_space = {
'max_iter' : trial.suggest_int('max_iter', 400, 1200, step=100),
'max_depth' : trial.suggest_int('max_depth', 3, 10, step=1),
'learning_rate' : trial.suggest_float('learning_rate', 0.01, 1),
'min_samples_leaf' : trial.suggest_int('min_samples_leaf', 1, 20, step=1),
'l2_regularization' : trial.suggest_float('l2_regularization', 0, 1)
}
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],
search_space = search_space,
lags_grid = lags_grid,
steps = 36,
refit = False,
metric = 'mean_absolute_error',
initial_train_size = len(datos_train),
fixed_train_size = False,
n_trials = 20,
random_state = 123,
return_best = True,
n_jobs = 'auto',
verbose = False,
show_progress = True
)
# Backtesting con datos de test
# ==============================================================================
metrica_histgb, predicciones = backtesting_forecaster(
forecaster = forecaster,
y = datos['users'],
exog = datos[exog_cols],
initial_train_size = len(datos.loc[:fin_validacion]),
fixed_train_size = False,
steps = 36,
refit = False,
metric = 'mean_absolute_error',
n_jobs = 'auto',
verbose = False
)
print(f"Error de backtest (MAE): {metrica_histgb:.2f}")
CatBoost¶
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 = ForecasterAutoreg(
regressor = CatBoostRegressor(
random_state=123,
silent=True,
allow_writing_files=False,
boosting_type = 'Plain', # Faster training
leaf_estimation_iterations = 3, # Faster training
),
lags = 24,
transformer_exog = one_hot_encoder
)
# Búsqueda de hiperparámetros
# ==============================================================================
lags_grid = [[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, step=1),
'learning_rate' : trial.suggest_float('learning_rate', 0.01, 1),
}
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],
search_space = search_space,
lags_grid = lags_grid,
steps = 36,
refit = False,
metric = 'mean_absolute_error',
initial_train_size = len(datos_train),
fixed_train_size = False,
n_trials = 20,
random_state = 123,
return_best = True,
n_jobs = 'auto',
verbose = False,
show_progress = True
)
# Backtesting con datos de test
# ==============================================================================
metrica_catboost, predicciones = backtesting_forecaster(
forecaster = forecaster,
y = datos['users'],
exog = datos[exog_cols],
initial_train_size = len(datos.loc[:fin_validacion]),
fixed_train_size = False,
steps = 36,
refit = False,
metric = 'mean_absolute_error',
n_jobs = 'auto',
verbose = False
)
print(f"Error de backtest (MAE): {metrica_catboost:.2f}")
Conclusió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, CatBoost y HistGradientBoostingRegressor.
pd.DataFrame({
'Model': ['Baseline', 'LGBMRegressor', 'XGBRegressor',
'HistGradientBoostingRegressor', 'CatBoostRegressor'],
'MAE': [metrica_baseline, metrica_lgbm, metrica_xgboost,
metrica_histgb, metrica_catboost]
}).round(2).sort_values(by='MAE')
🖉 Nota
Para ilustrar el proceso, la búsqueda de hiperparámetros se ha mantenido reducida. Sin embargo, para obtener resultados óptimos, puede ser necesario realizar una búsqueda más amplia para cada modelo.Información de sesión¶
import session_info
session_info.show(html=False)
Bibliografía¶
Hyndman, R.J., & Athanasopoulos, G. (2021) Forecasting: principles and practice, 3rd edition, OTexts: Melbourne, Australia.
Time Series Analysis and Forecasting with ADAM Ivan Svetunkov.
Joseph, M. (2022). Modern time series forecasting with Python: Explore industry-ready time series forecasting using modern machine learning and Deep Learning. Packt Publishing.
Instrucciones para citar¶
¿Cómo citar este documento?
Si utilizas este documento o alguna parte de él, te agradecemos que lo cites. ¡Muchas gracias!
Forecasting series temporales con gradient boosting: Skforecast, XGBoost, LightGBM y Catboost por Joaquín Amat Rodrigo y Javier Escobar Ortiz, disponible bajo licencia Attribution-NonCommercial-ShareAlike 4.0 International en https://www.cienciadedatos.net/py39-forecasting-series-temporales-con-skforecast-xgboost-lightgbm-catboost.html
¿Cómo citar skforecast?
Zenodo:
Amat Rodrigo, Joaquin, & Escobar Ortiz, Javier. (2023). skforecast (v0.10.0). Zenodo. https://doi.org/10.5281/zenodo.8382788
APA:
Amat Rodrigo, J., & Escobar Ortiz, J. (2023). skforecast (Version 0.10.0) [Computer software]. https://doi.org/10.5281/zenodo.8382788
BibTeX:
@software{skforecast, author = {Amat Rodrigo, Joaquin and Escobar Ortiz, Javier}, title = {skforecast}, version = {0.10.0}, month = {9}, year = {2023}, license = {BSD-3-Clause}, url = {https://skforecast.org/}, doi = {10.5281/zenodo.8382788} }
¿Te ha gustado el artículo? Tu ayuda es importante
Mantener un sitio web tiene unos costes elevados, tu contribución me ayudará a seguir generando contenido divulgativo gratuito. ¡Muchísimas gracias! 😊
Este documento creado por Joaquín Amat Rodrigo y Javier Escobar Ortiz tiene licencia Attribution-NonCommercial-ShareAlike 4.0 International.