Forecasting series temporales con LightGBM

Más sobre forecasting en: cienciadedatos.net

• Forecasting series temporales con machine learning
• Modelos ARIMA y SARIMAX
• Forecasting series temporales con gradient boosting: Skforecast, XGBoost, LightGBM y CatBoost
• Forecasting series temporales con XGBoost
• Global Forecasting: Multi-series forecasting
• Forecasting de la demanda eléctrica con machine learning
• Modelos de forecasting globales: Análisis comparativo de modelos de una y múltiples series
• Forecasting con deep learning
• Forecasting de visitas a página web con machine learning
• Forecasting del precio de Bitcoin
• Forecasting probabilístico
• Forecasting de demanda intermitente
• Reducir el impacto del Covid en modelos de forecasting
• Modelar series temporales con tendencia utilizando modelos de árboles

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 proporciona soluciones automatizadas a estos retos, facilitando la aplicación y validación de modelos de machine learning a problemas de forecasting. La librería soporta varios modelos avanzados de gradient boosting, incluyendo LightGBM.

Librerías¶

Librerías utilizadas en este documento.

# Procesado de datos
# ==============================================================================
import numpy as np
import pandas as pd
from skforecast.datasets import fetch_dataset

# Gráficos
# ==============================================================================
import matplotlib.pyplot as plt
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from skforecast.plot import set_dark_theme
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')
plt.rcParams.update({'font.size': 8})

# Modelado y forecasting
# ==============================================================================
import lightgbm
import skforecast
import sklearn
from lightgbm import LGBMRegressor
from sklearn.feature_selection import RFECV
from skforecast.recursive import ForecasterRecursive
from skforecast.model_selection import TimeSeriesFold
from skforecast.model_selection import bayesian_search_forecaster
from skforecast.model_selection import backtesting_forecaster
from skforecast.feature_selection import select_features
from skforecast.preprocessing import RollingFeatures
import shap

# Configuración warnings
# ==============================================================================
import warnings
warnings.filterwarnings('once')

color = '\033[1m\033[38;5;208m'
print(f'{color}Versión skforecast: {skforecast.__version__}')
print(f'{color}Versión scikit-learn: {sklearn.__version__}')
print(f'{color}Versión lightgbm: {lightgbm.__version__}')

Versión skforecast: 0.25.0
Versión scikit-learn: 1.7.2
Versión lightgbm: 4.7.0

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.

# Descarga de datos
# ==============================================================================
datos = fetch_dataset('bike_sharing_extended_features')
datos.head(4)

╭──────────────────────── bike_sharing_extended_features ─────────────────────────╮
│ Description:                                                                    │
│ 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, the dataset   │
│ was enriched by introducing supplementary features. Addition includes calendar- │
│ based variables (day of the week, hour of the day, month, etc.), indicators for │
│ sunlight, incorporation of rolling temperature averages, and the creation of    │
│ polynomial features generated from variable pairs. All cyclic variables are     │
│ encoded using sine and cosine functions to ensure accurate representation.      │
│                                                                                 │
│ Source:                                                                         │
│ Fanaee-T,Hadi. (2013). Bike Sharing Dataset. UCI Machine Learning Repository.   │
│ https://doi.org/10.24432/C5W894.                                                │
│                                                                                 │
│ URL:                                                                            │
│ https://raw.githubusercontent.com/skforecast/skforecast-                        │
│ datasets/main/data/bike_sharing_extended_features.csv                           │
│                                                                                 │
│ Shape: 17352 rows x 90 columns                                                  │
╰─────────────────────────────────────────────────────────────────────────────────╯

Para facilitar el entrenamiento de los modelos, la búsqueda de hiperparámetros óptimos y la evaluación de su precisión predictiva, los datos se dividen en tres conjuntos separados: 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:, :]

particiones = {'train': datos_train, 'validación': datos_val, 'test': datos_test}
for nombre, particion in particiones.items():
    print(
        f'Fechas {nombre:<10} : {particion.index.min()} --- '
        f'{particion.index.max()}  (n={len(particion)})'
    )

Fechas train      : 2011-01-08 00:00:00 --- 2012-03-31 23:00:00  (n=10776)
Fechas validación : 2012-04-01 00:00:00 --- 2012-08-31 23:00:00  (n=3672)
Fechas test       : 2012-09-01 00:00:00 --- 2012-12-30 23:00:00  (n=2904)

Exploración de datos¶

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()
particiones = {'Train': datos_train, 'Validation': datos_val, 'Test': datos_test}
for nombre, particion in particiones.items():
    fig.add_trace(
        go.Scatter(
            x=particion.index, y=particion['users'], mode='lines', name=nombre
        )
    )
fig.update_layout(
    title  = 'Número de usuarios',
    xaxis_title='Fecha',
    yaxis_title='Usuarios',
    legend_title='Partición:',
    width=800,
    height=400,
    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áficos de estacionalidad¶

Los gráficos estacionales son una herramienta útil para identificar patrones y tendencias estacionales en una serie temporal. Se crean agrupando las observaciones por estación (por ejemplo, mes, día de la semana u hora del día) y comparando la distribución de los valores dentro de cada grupo. En este caso, se utilizan diagramas de caja (boxplots) junto con la mediana de cada grupo.

# Estacionalidad anual, semanal y diaria
# ==============================================================================
set_dark_theme()
fig, axs = plt.subplots(2, 2, figsize=(8, 5), sharex=False, sharey=True)
axs = axs.ravel()
flierprops = dict(
    marker='o', markerfacecolor='white', markeredgecolor='black', markersize=4
)

# Distribución de usuarios por mes
datos['month'] = datos.index.month
datos.boxplot(column='users', by='month', ax=axs[0], flierprops=flierprops)
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', fontsize=10)

# Distribució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=flierprops)
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', fontsize=10)

# Distribució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=flierprops)
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', fontsize=10)

# Distribució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       = 'Usuarios promedio',
    xticks      = [i * 24 for i in range(7)],
    xticklabels = ['Mon', 'Tue', 'Wed', 'Thu', 'Fri', 'Sat', 'Sun'],
    xlabel      = 'Día y hora',
    ylabel      = 'Users'
)
axs[3].title.set_size(10)

fig.suptitle('Gráficos de estacionalidad', fontsize=12)
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, fft=True)
plt.show()

# Gráfico autocorrelación parcial
# ==============================================================================
fig, ax = plt.subplots(figsize=(5, 2))
plot_pacf(datos['users'], ax=ax, lags=72, method='burg')
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 el número de usuarios durante periodos específicos del pasado proporciona información útil para predecir el número de usuarios en el futuro.

Modelo recursivo con LightGBM¶

En primer lugar, se entrena un modelo ForecasterRecursive utilizando valores pasados de la variable de respuesta (lags) y la media móvil 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(
    estimator            = LGBMRegressor(random_state=15926, verbose=-1),
    lags                 = 72,
    window_features      = window_features,
    categorical_features = 'auto'
)

# Entrenar el forecaster
# ==============================================================================
forecaster.fit(y=datos.loc[:fin_validacion, 'users'])
forecaster

# Predicciones
# ==============================================================================
forecaster.predict(steps=10)

2012-09-01 00:00:00    116.543722
2012-09-01 01:00:00     81.892246
2012-09-01 02:00:00     43.328288
2012-09-01 03:00:00     14.088216
2012-09-01 04:00:00      7.404299
2012-09-01 05:00:00     13.773338
2012-09-01 06:00:00     47.292544
2012-09-01 07:00:00    136.366345
2012-09-01 08:00:00    276.658280
2012-09-01 09:00:00    377.346129
Freq: h, Name: pred, dtype: float64

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.

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 los 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'
                        )
predicciones.head()

# Error de backtest
# ==============================================================================
metrica

El modelo autorregresivo, que utiliza los primeros 72 lags y la media móvil de las últimas 72 horas como predictores, alcanza un MAE en test de 71.0.

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 exógenas incluidas en el modelo
# ==============================================================================
exog_features = [
    'month_sin',
    'month_cos',
    'week_of_year_sin',
    'week_of_year_cos',
    'week_day_sin',
    'week_day_cos',
    'hour_day_sin',
    'hour_day_cos',
    'sunrise_hour_sin',
    'sunrise_hour_cos',
    'sunset_hour_sin',
    'sunset_hour_cos',
    'holiday_previous_day',
    'holiday_next_day',
    'temp_roll_mean_1_day',
    'temp_roll_mean_7_day',
    'temp_roll_max_1_day',
    'temp_roll_min_1_day',
    'temp_roll_max_7_day',
    'temp_roll_min_7_day',
    'temp',
    'holiday'
]

# Backtesting en los datos de test incluyendo las variables exógenas
# ==============================================================================
metrica, predicciones = backtesting_forecaster(
                            forecaster = forecaster,
                            y          = datos['users'],
                            exog       = datos[exog_features],
                            cv         = cv,
                            metric     = 'mean_absolute_error'
                        )
metrica

La incorporación de variables exógenas como predictores mejora la capacidad predictiva del modelo: el MAE se reduce de 71.0 a 59.4. Dado que ambos forecasters utilizan los mismos lags, window features y estimador, la mejora puede atribuirse a las variables exógenas.

Optimización de hiperparámetros¶

La optimización de hiperparámetros es un paso crítico en el desarrollo de modelos de machine learning efectivos. El ForecasterRecursive utilizado en las secciones anteriores incluía los primeros 72 lags y un modelo LGBMRegressor 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:

1. Entrenar el modelo utilizando sólo el conjunto de entrenamiento.

2. El modelo se evalúa utilizando el conjunto de validación mediante backtesting.

3. Seleccionar la combinación de hiperparámetros y retardos que proporcione el menor error.

4. 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(
    estimator            = LGBMRegressor(random_state=15926, verbose=-1),
    lags                 = 72,
    window_features      = window_features,
    categorical_features = 'auto'
)

# 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):
    return {
        'n_estimators'    : trial.suggest_int('n_estimators', 400, 1200, 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, step=0.1),
        'max_bin'         : trial.suggest_int('max_bin', 50, 250, step=25),
        'reg_alpha'       : trial.suggest_float('reg_alpha', 0, 1, step=0.1),
        'reg_lambda'      : trial.suggest_float('reg_lambda', 0, 1, step=0.1),
        'lags'            : trial.suggest_categorical('lags', lags_grid)
    }

# 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'], # Datos de test no utilizados
    exog          = datos.loc[:fin_validacion, exog_features],
    cv            = cv_search,
    search_space  = search_space,
    metric        = 'mean_absolute_error',
    n_trials      = 20, # Aumentar para una búsqueda más exhaustiva
    return_best   = True
)

best_params = results_search['params'].iat[0]
best_params = best_params | {'random_state': 15926, 'verbose': -1}
best_lags   = results_search['lags'].iat[0]

# Resultados de la búsqueda
# ==============================================================================
results_search.head(3)

Al indicar return_best = True, el objeto forecaster se actualiza automáticamente con la mejor configuración encontrada y se reentrena con todos los datos proporcionados a la búsqueda (conjuntos de entrenamiento y validación). El conjunto de test no se utiliza en ningún momento, por lo que puede emplearse para obtener una estimación no sesgada del rendimiento del modelo final.

# Mejor modelo
# ==============================================================================
forecaster

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.

# Backtest modelo final con datos de test
# ==============================================================================
metrica, predicciones = backtesting_forecaster(
    forecaster = forecaster,
    y          = datos['users'],
    exog       = datos[exog_features],
    cv         = cv,
    metric     = 'mean_absolute_error'
)
display(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='Valor real vs predicciones en los datos de test',
    xaxis_title='Fecha',
    yaxis_title='Usuarios',
    width=800,
    height=400,
    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()

Con los hiperparámetros y lags seleccionados, el MAE en test se mantiene prácticamente igual, pasando de 59.4 a 59.6. La optimización de hiperparámetros no garantiza una mejora en datos nuevos: la búsqueda se ha mantenido muy reducida (20 iteraciones) y la configuración se elige según su rendimiento en el periodo de validación, que puede diferir del periodo de test. En un proyecto real, se recomienda realizar una búsqueda más exhaustiva.

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 estimadores 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. Skforecast facilita su uso con forecasters mediante la función select_features() (guía de usuario).

# Crear forecaster
# ==============================================================================
estimator = LGBMRegressor(
    n_estimators = 100,
    max_depth    = 5,
    random_state = 15926,
    verbose      = -1
)
forecaster = ForecasterRecursive(
    estimator            = estimator,
    lags                 = best_lags,
    window_features      = window_features,
    categorical_features = 'auto'
)

# Eliminación recursiva de predictores con validación cruzada
# ==============================================================================
warnings.filterwarnings('ignore', message='X does not have valid feature names.*')
selector = RFECV(
    estimator = estimator,
    step      = 1,
    cv        = 3,
)
lags_seleccionados, wf_seleccionadas, exog_seleccionadas, _ = select_features(
    forecaster      = forecaster,
    selector        = selector,
    y               = datos_train['users'],
    exog            = datos_train[exog_features],
    select_only     = None,
    force_inclusion = None,
    subsample       = 0.5,
    random_state    = 123,
    verbose         = True,
)

Recursive feature elimination (RFECV)
-------------------------------------
Total number of records available: 10607
Total number of records used for feature selection: 5303
Number of features available: 32
    Lags            (n=9)
    Window features (n=1)
    Exog            (n=22)
    Calendar        (n=0)
Number of features selected: 25
    Lags            (n=9) : [1, 2, 3, 23, 24, 25, 167, 168, 169]
    Window features (n=1) : ['roll_mean_72']
    Exog            (n=15) : ['month_sin', 'month_cos', 'week_of_year_sin', 'week_of_year_cos', 'week_day_sin', 'week_day_cos', 'hour_day_sin', 'hour_day_cos', 'temp_roll_mean_1_day', 'temp_roll_mean_7_day', 'temp_roll_max_1_day', 'temp_roll_min_1_day', 'temp_roll_max_7_day', 'temp_roll_min_7_day', 'temp']
    Calendar        (n=0) : []

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 métrica 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 forecaster se entrena y evalúa de nuevo utilizando el conjunto de predictores seleccionados.

# Crear forecaster con los predictores seleccionados
# ==============================================================================
# Las window features se incluyen solo si han sido seleccionadas
wf_finales = window_features if wf_seleccionadas else None
forecaster = ForecasterRecursive(
    estimator            = LGBMRegressor(**best_params),
    lags                 = lags_seleccionados,
    window_features      = wf_finales,
    categorical_features = 'auto'
)

# Backtesting con los predictores seleccionados y los datos de test
# ==============================================================================
metrica, predicciones = backtesting_forecaster(
    forecaster = forecaster,
    y          = datos['users'],
    exog       = datos[exog_seleccionadas],
    cv         = cv,
    metric     = 'mean_absolute_error'
)
metrica

El número de predictores se ha reducido de 32 a 25 y el MAE en test ha mejorado ligeramente, de 59.6 a 59.2. El modelo es ahora mucho más simple, lo que hace que sea más rápido de entrenar y menos propenso al sobreajuste, sin pérdida de capacidad predictiva. Para el resto del documento, el modelo se entrenará utilizando sólo las variables exógenas seleccionadas.

# Actualizar las variables exógenas utilizadas
# ==============================================================================
exog_features = exog_seleccionadas

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(
    estimator            = LGBMRegressor(**best_params),
    lags                 = best_lags,
    window_features      = window_features,
    categorical_features = 'auto'
)
forecaster.fit(
    y    = datos.loc[:fin_validacion, 'users'],
    exog = datos.loc[:fin_validacion, exog_features]
)

Model-specific feature importance¶

# 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 estimador interno del forecaster.

• Las matrices de entrenamiento creadas a partir de la serie temporal y variables exógenas, utilizadas para ajustar el forecaster.

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 estimador
# ==============================================================================
X_train, y_train = forecaster.create_train_X_y(
                        y    = datos.loc[:fin_validacion, 'users'],
                        exog = datos.loc[:fin_validacion, exog_features]
                    )
display(X_train.head(3))
display(y_train.head(3))

date_time
2011-01-15 01:00:00    20.0
2011-01-15 02:00:00    12.0
2011-01-15 03:00:00     8.0
Freq: h, Name: y, dtype: float64

# Crear SHAP explainer (para modelos basados en árboles)
# ==============================================================================
explainer = shap.TreeExplainer(forecaster.estimator)

# 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()
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, colors='white')
fig.set_size_inches(8, 4.5)

Los valores SHAP no solo permiten interpretar el comportamiento general del modelo, sino que también son una herramienta poderosa para analizar predicciones individuales. Esto resulta especialmente útil cuando se quiere entender cómo se ha generado una predicción específica y qué variables han contribuido a ella.

Para llevar a cabo este análisis, es necesario acceder a los valores de los predictores (lags, window features y variables exógenas) en el momento de la predicción. Esto puede lograrse utilizando el método create_predict_X() o bien activando el argumento return_predictors=True en la función backtesting_forecaster().

Supóngase que se quiere entender la predicción obtenida durante el backtesting para la fecha 2012-10-06 12:00:00.

# 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_features],
                        cv                = cv,
                        metric            = 'mean_absolute_error',
                        return_predictors = True,
                     )

Al indicar return_predictors=True, se obtiene un DataFrame con el valor predicho ('pred'), la partición en la que se encuentra ('fold') y el valor de los predictores (lags, window features y variables exógenas) utilizados para realizar cada predicción.

predicciones.head(3)

# Waterfall para una predicción concreta
# ==============================================================================
# Asegurar que los tipos son los mismos que en los datos de entrenamiento
predicciones = predicciones.astype(datos[exog_features].dtypes)
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 = plt.gcf()
fig.set_size_inches(8, 3.5)
fig.axes[0].tick_params(labelsize=8)
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:],
)

Información de sesión¶

import session_info
session_info.show(html=False)

-----
lightgbm            4.7.0
matplotlib          3.10.9
numpy               2.4.6
optuna              4.9.0
pandas              2.3.3
plotly              6.9.0
session_info        v1.0.1
shap                0.52.0
skforecast          0.25.0
sklearn             1.7.2
statsmodels         0.14.6
-----
IPython             9.15.0
jupyter_client      8.9.1
jupyter_core        5.9.1
-----
Python 3.13.14 | packaged by conda-forge | (main, Jun 12 2026, 09:44:26) [MSC v.1944 64 bit (AMD64)]
Windows-11-10.0.26200-SP0
-----
Session information updated at 2026-09-24 12:06

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 LightGBM por Joaquín Amat Rodrigo y Javier Escobar Ortiz, disponible bajo licencia Attribution-NonCommercial-ShareAlike 4.0 International en https://www.cienciadedatos.net/documentos/py58-forecasting-series-temporales-con-lightgbm.html

¿Cómo citar skforecast?

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

Zenodo:

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

APA:

Amat Rodrigo, J., & Escobar Ortiz, J. (2026). skforecast (Version 0.25.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.25.0}, month = {09}, year = {2026}, license = {BSD-3-Clause}, url = {https://skforecast.org/}, doi = {10.5281/zenodo.8382788} }

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

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

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

Se permite:

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

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

Bajo los siguientes términos:

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

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

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