Forecasting series temporales con XGBoost

Si te gusta  Skforecast , ayúdanos dándonos una estrella en  GitHub! ⭐️

Forecasting series temporales con XGBoost

Joaquín Amat Rodrigo, Javier Escobar Ortiz
Junio, 2024 (última actualización Agosto 2024)

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 librearí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 XGBoost.

✎ Nota

Este documento es un resumen de una guía más completa sobre el uso de modelos de *gradient boosting* para el forecasting de series temporales. La guía completa está disponible en Forecasting series temporales con gradient boosting: Skforecast, XGBoost, LightGBM y CatBoost.

Librearías

Librerías utilizadas en este documento.

In [1]:
# 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
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
# ==============================================================================
import xgboost
import skforecast
import sklearn
from xgboost import XGBRegressor
from sklearn.feature_selection import RFECV
from skforecast.ForecasterAutoreg import ForecasterAutoreg
from skforecast.model_selection import bayesian_search_forecaster
from skforecast.model_selection import backtesting_forecaster
from skforecast.model_selection import select_features
import shap

# Warnings
# ==============================================================================
import warnings
warnings.filterwarnings('once')

color = '\033[1m\033[38;5;208m'
print(f"{color}Version skforecast: {skforecast.__version__}")
print(f"{color}Version scikit-learn: {sklearn.__version__}")
print(f"{color}Version xgboost: {xgboost.__version__}")
print(f"{color}Version pandas: {pd.__version__}")
print(f"{color}Version numpy: {np.__version__}")
Version skforecast: 0.13.0
Version scikit-learn: 1.5.1
Version xgboost: 2.1.0
Version pandas: 2.2.2
Version numpy: 2.0.1

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.

In [2]:
# Descarga de datos
# ==============================================================================
datos = fetch_dataset('bike_sharing_extended_features')
datos.head(4)
bike_sharing_extended_features
------------------------------
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.
Fanaee-T,Hadi. (2013). Bike Sharing Dataset. UCI Machine Learning Repository.
https://doi.org/10.24432/C5W894.
Shape of the dataset: (17352, 90)
Out[2]:
users weather 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 holiday_previous_day holiday_next_day temp holiday
date_time
2011-01-08 00:00:00 25.0 mist 0.5 0.866025 0.120537 0.992709 -0.781832 0.62349 0.258819 0.965926 ... 8.063334 10.127976 9.02 6.56 18.86 4.92 0.0 0.0 7.38 0.0
2011-01-08 01:00:00 16.0 mist 0.5 0.866025 0.120537 0.992709 -0.781832 0.62349 0.500000 0.866025 ... 8.029166 10.113334 9.02 6.56 18.86 4.92 0.0 0.0 7.38 0.0
2011-01-08 02:00:00 16.0 mist 0.5 0.866025 0.120537 0.992709 -0.781832 0.62349 0.707107 0.707107 ... 7.995000 10.103572 9.02 6.56 18.86 4.92 0.0 0.0 7.38 0.0
2011-01-08 03:00:00 7.0 rain 0.5 0.866025 0.120537 0.992709 -0.781832 0.62349 0.866025 0.500000 ... 7.960833 10.093809 9.02 6.56 18.86 4.92 0.0 0.0 7.38 0.0

4 rows × 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.

In [3]:
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()}  "
    f"(n={len(datos_train)})"
)
print(
    f"Fechas validación : {datos_val.index.min()} --- {datos_val.index.max()}  "
    f"(n={len(datos_val)})"
)
print(
    f"Fechas test       : {datos_test.index.min()} --- {datos_test.index.max()}  "
    f"(n={len(datos_test)})"
)
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)

Exploracion 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

In [4]:
# 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á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.

In [5]:
# Estacionalidad anual, semanal y diaria
# ==============================================================================
fig, axs = plt.subplots(2, 2, figsize=(8, 5), sharex=False, sharey=True)
axs = axs.ravel()

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

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

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

# 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"
)
axs[3].title.set_size(10)

fig.suptitle("Gráficos de estacionalidad", fontsize=14)
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.

In [6]:
# Gráfico autocorrelación
# ==============================================================================
fig, ax = plt.subplots(figsize=(5, 2))
plot_acf(datos['users'], ax=ax, lags=72)
plt.show()
In [7]:
# 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.

Modelo autoregresivo recursivo con XGBoost

Primero, se entrena un modelo ForecasterAutoreg utilizando valores pasados (lags) de la variable respuesta como predictores. Posteriormente, se añaden variables exógenas al modelo y se evalúa la mejora en su rendimiento. Dado que los modelos de Gradient Boosting tienen un elevado 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. Finalmente, se evalúa la capacidad predictiva del modelo utilizando un proceso de backtesting.

Forecaster

In [8]:
# Crear el forecaster
# ==============================================================================
forecaster = ForecasterAutoreg(
    regressor = XGBRegressor(random_state=15926, enable_categorical=True),
    lags      = 24
)

# Entrenar el forecaster
# ==============================================================================
forecaster.fit(y=datos.loc[:fin_validacion, 'users'])
forecaster
Out[8]:
================= 
ForecasterAutoreg 
================= 
Regressor: XGBRegressor(base_score=None, booster=None, callbacks=None,
             colsample_bylevel=None, colsample_bynode=None,
             colsample_bytree=None, device=None, early_stopping_rounds=None,
             enable_categorical=True, eval_metric=None, feature_types=None,
             gamma=None, grow_policy=None, importance_type=None,
             interaction_constraints=None, learning_rate=None, max_bin=None,
             max_cat_threshold=None, max_cat_to_onehot=None,
             max_delta_step=None, max_depth=None, max_leaves=None,
             min_child_weight=None, missing=nan, monotone_constraints=None,
             multi_strategy=None, n_estimators=None, n_jobs=None,
             num_parallel_tree=None, random_state=15926, ...) 
Lags: [ 1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24] 
Transformer for y: None 
Transformer for exog: None 
Window size: 24 
Weight function included: False 
Differentiation order: None 
Exogenous included: False 
Exogenous variables names: None 
Training range: [Timestamp('2011-01-08 00:00:00'), Timestamp('2012-08-31 23:00:00')] 
Training index type: DatetimeIndex 
Training index frequency: h 
Regressor parameters: {'objective': 'reg:squarederror', 'base_score': None, 'booster': None, 'callbacks': None, 'colsample_bylevel': None, 'colsample_bynode': None, 'colsample_bytree': None, 'device': None, 'early_stopping_rounds': None, 'enable_categorical': True, 'eval_metric': None, 'feature_types': None, 'gamma': None, 'grow_policy': None, 'importance_type': None, 'interaction_constraints': None, 'learning_rate': None, 'max_bin': None, 'max_cat_threshold': None, 'max_cat_to_onehot': None, 'max_delta_step': None, 'max_depth': None, 'max_leaves': None, 'min_child_weight': None, 'missing': nan, 'monotone_constraints': None, 'multi_strategy': None, 'n_estimators': None, 'n_jobs': None, 'num_parallel_tree': None, 'random_state': 15926, 'reg_alpha': None, 'reg_lambda': None, 'sampling_method': None, 'scale_pos_weight': None, 'subsample': None, 'tree_method': None, 'validate_parameters': None, 'verbosity': None} 
fit_kwargs: {} 
Creation date: 2024-08-10 18:46:42 
Last fit date: 2024-08-10 18:46:42 
Skforecast version: 0.13.0 
Python version: 3.12.4 
Forecaster id: None 
In [9]:
# Predicciones
# ==============================================================================
forecaster.predict(steps=10)
Out[9]:
2012-09-01 00:00:00    123.423668
2012-09-01 01:00:00     75.687805
2012-09-01 02:00:00     43.349522
2012-09-01 03:00:00     15.084871
2012-09-01 04:00:00      3.778617
2012-09-01 05:00:00     13.915284
2012-09-01 06:00:00     44.951176
2012-09-01 07:00:00    132.655685
2012-09-01 08:00:00    298.680359
2012-09-01 09:00:00    435.629669
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.

In [10]:
# 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            = False, # Cambiar a False para mostrar más información
    show_progress      = True
)
predicciones.head()
Out[10]:
pred
2012-09-01 00:00:00 123.423668
2012-09-01 01:00:00 75.687805
2012-09-01 02:00:00 43.349522
2012-09-01 03:00:00 15.084871
2012-09-01 04:00:00 3.778617
In [11]:
# Métrica de backtesting
# ==============================================================================
metrica
Out[11]:
mean_absolute_error
0 73.466462

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.

✎ Nota

For a more detailed explanation of how to use exogenous variables, such as categorical features, visit: Para una explicación más detallada de cómo utilizar variables exógenas, visitar:
In [12]:
# 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'
]
In [13]:
# Backtesting en los datos de test incluyendo las variables exógenas
# ==============================================================================
metrica, predicciones = backtesting_forecaster(
    forecaster         = forecaster,
    y                  = datos['users'],
    exog               = datos[exog_features],
    steps              = 36,
    metric             = 'mean_absolute_error',
    initial_train_size = len(datos[:fin_validacion]),
    refit              = False,
    n_jobs             = 'auto',
    verbose            = False,
    show_progress      = True
)
metrica
Out[13]:
mean_absolute_error
0 62.119768

Adding exogenous variables to the model improves its predictive capacity.

Optimización de hiperparámetros (tuning)

La optimización de hiperparámetros es un paso crítico en el desarrollo de modelos de machine learning efectivos. El ForecasterAutoreg utilizado en las secciones anteriores incluía un modelo XGBoost con los hiperparámetros por defecto. Sin embargo, no hay razón para pensar que estos valores son los más adecuados. Para encontrar los mejores hiperparámetros, se realiza una Búsqueda Bayesiana utilizando la función bayesian_search_forecaster(). La búsqueda se lleva a cabo utilizando el mismo proceso de backtesting que antes, pero cada vez se entrena el modelo con diferentes combinaciones de hiperparámetros y lags. Es importante tener en cuenta que la búsqueda de hiperparámetros debe hacerse utilizando el conjunto de validación, sin incluir el conjunto de test.

La búsqueda se realiza probando cada combinación de hiperparámetros y lags de la siguiente manera:

  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 lags que dé el menor error.

  4. Entrenar el modelo de nuevo 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.

✎ Nota

Desde la versión 0.12.0 los lags se incluye dentro del espacio de búsqueda (search_space).
In [14]:
# Búsqueda de hiperparámetros
# ==============================================================================

# Crear forecaster
forecaster = ForecasterAutoreg(
    regressor = XGBRegressor(random_state=15926, enable_categorical=True),
    lags      = 72
)

# 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, 1),
        'subsample'       : trial.suggest_float('subsample', 0.1, 1),
        'colsample_bytree': trial.suggest_float('colsample_bytree', 0.1, 1),
        'gamma'           : trial.suggest_float('gamma', 0, 1),
        'reg_alpha'       : trial.suggest_float('reg_alpha', 0, 1),
        'reg_lambda'      : trial.suggest_float('reg_lambda', 0, 1),
        'lags'            : trial.suggest_categorical('lags', lags_grid)
    } 
    return search_space

resultados_busqueda, frozen_trial = bayesian_search_forecaster(
    forecaster         = forecaster,
    y                  = datos.loc[:fin_validacion, 'users'], # Datos test no incluidos
    exog               = datos.loc[:fin_validacion, exog_features],
    search_space       = search_space,
    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
)
best_params = resultados_busqueda['params'].iat[0]
best_lags = resultados_busqueda['lags'].iat[0]
`Forecaster` refitted using the best-found lags and parameters, and the whole data set: 
  Lags: [  1   2   3  23  24  25 167 168 169] 
  Parameters: {'n_estimators': 400, 'max_depth': 10, 'learning_rate': 0.1987029596852645, 'subsample': 0.8452581056448839, 'colsample_bytree': 0.6490537930286489, 'gamma': 0.42743952450393, 'reg_alpha': 0.7787952214581544, 'reg_lambda': 0.17127802671184533}
  Backtesting metric: 64.3772965454393

In [15]:
# Resultados de la búsqueda de hiperparámetros
# ==============================================================================
resultados_busqueda.head(5)
Out[15]:
lags params mean_absolute_error n_estimators max_depth learning_rate subsample colsample_bytree gamma reg_alpha reg_lambda
17 [1, 2, 3, 23, 24, 25, 167, 168, 169] {'n_estimators': 400, 'max_depth': 10, 'learni... 64.377297 400.0 10.0 0.198703 0.845258 0.649054 0.427440 0.778795 0.171278
19 [1, 2, 3, 23, 24, 25, 167, 168, 169] {'n_estimators': 500, 'max_depth': 9, 'learnin... 67.255268 500.0 9.0 0.158141 0.874604 0.644145 0.456653 0.748756 0.200001
15 [1, 2, 3, 23, 24, 25, 167, 168, 169] {'n_estimators': 400, 'max_depth': 9, 'learnin... 68.130566 400.0 9.0 0.016108 0.724424 0.663652 0.133031 0.808592 0.151677
16 [1, 2, 3, 23, 24, 25, 167, 168, 169] {'n_estimators': 400, 'max_depth': 10, 'learni... 68.204207 400.0 10.0 0.163782 0.764628 0.634598 0.402922 0.788758 0.014016
13 [1, 2, 3, 23, 24, 25, 167, 168, 169] {'n_estimators': 900, 'max_depth': 8, 'learnin... 69.459926 900.0 8.0 0.012168 0.748844 0.775209 0.085690 0.899431 0.012284

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.

In [16]:
# Mejor modelo
# ==============================================================================
forecaster
Out[16]:
================= 
ForecasterAutoreg 
================= 
Regressor: XGBRegressor(base_score=None, booster=None, callbacks=None,
             colsample_bylevel=None, colsample_bynode=None,
             colsample_bytree=0.6490537930286489, device=None,
             early_stopping_rounds=None, enable_categorical=True,
             eval_metric=None, feature_types=None, gamma=0.42743952450393,
             grow_policy=None, importance_type=None,
             interaction_constraints=None, learning_rate=0.1987029596852645,
             max_bin=None, max_cat_threshold=None, max_cat_to_onehot=None,
             max_delta_step=None, max_depth=10, max_leaves=None,
             min_child_weight=None, missing=nan, monotone_constraints=None,
             multi_strategy=None, n_estimators=400, n_jobs=None,
             num_parallel_tree=None, random_state=15926, ...) 
Lags: [  1   2   3  23  24  25 167 168 169] 
Transformer for y: None 
Transformer for exog: None 
Window size: 169 
Weight function included: False 
Differentiation order: None 
Exogenous included: True 
Exogenous variables names: ['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'] 
Training range: [Timestamp('2011-01-08 00:00:00'), Timestamp('2012-08-31 23:00:00')] 
Training index type: DatetimeIndex 
Training index frequency: h 
Regressor parameters: {'objective': 'reg:squarederror', 'base_score': None, 'booster': None, 'callbacks': None, 'colsample_bylevel': None, 'colsample_bynode': None, 'colsample_bytree': 0.6490537930286489, 'device': None, 'early_stopping_rounds': None, 'enable_categorical': True, 'eval_metric': None, 'feature_types': None, 'gamma': 0.42743952450393, 'grow_policy': None, 'importance_type': None, 'interaction_constraints': None, 'learning_rate': 0.1987029596852645, 'max_bin': None, 'max_cat_threshold': None, 'max_cat_to_onehot': None, 'max_delta_step': None, 'max_depth': 10, 'max_leaves': None, 'min_child_weight': None, 'missing': nan, 'monotone_constraints': None, 'multi_strategy': None, 'n_estimators': 400, 'n_jobs': None, 'num_parallel_tree': None, 'random_state': 15926, 'reg_alpha': 0.7787952214581544, 'reg_lambda': 0.17127802671184533, 'sampling_method': None, 'scale_pos_weight': None, 'subsample': 0.8452581056448839, 'tree_method': None, 'validate_parameters': None, 'verbosity': None} 
fit_kwargs: {} 
Creation date: 2024-08-10 18:46:49 
Last fit date: 2024-08-10 18:48:45 
Skforecast version: 0.13.0 
Python version: 3.12.4 
Forecaster id: None 

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 configuraciones. Esto ayudará a utilizar todo su potencial para analizar la capacidad predictiva del modelo.

In [17]:
# Backtest modelo final con datos de test
# ==============================================================================
metric, predicciones = backtesting_forecaster(
    forecaster         = forecaster,
    y                  = datos['users'],
    exog               = datos[exog_features],
    steps              = 36,
    metric             = 'mean_absolute_error',
    initial_train_size = len(datos[:fin_validacion]),
    refit              = False,
    n_jobs             = 'auto',
    verbose            = False,
    show_progress      = True
)
display(metrica)
predicciones.head()
mean_absolute_error
0 62.119768
Out[17]:
pred
2012-09-01 00:00:00 144.362183
2012-09-01 01:00:00 111.565620
2012-09-01 02:00:00 91.786720
2012-09-01 03:00:00 39.881084
2012-09-01 04:00:00 12.092457
In [18]:
# 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 se reduce.

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.
In [19]:
# Crear forecaster
# ==============================================================================
regressor = XGBRegressor(
    n_estimators       = 100,
    max_depth          = 5,
    random_state       = 15926,
    enable_categorical = True
)
forecaster = ForecasterAutoreg(
    regressor = regressor,
    lags      = [1, 2, 3, 23, 24, 25, 167, 168, 169],
)

# Eliminación recursiva de predictores con validación cruzada
# ==============================================================================
selector = RFECV(
    estimator              = regressor,
    step                   = 1,
    cv                     = 3,
    min_features_to_select = 10,
    n_jobs                 = -1
)
lags_seleccionados, 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: 31
    Autoreg (n=9)
    Exog    (n=22)
Number of features selected: 14
    Autoreg (n=8) : [1, 2, 3, 23, 24, 167, 168, 169]
    Exog    (n=6) : ['week_day_sin', 'week_day_cos', 'hour_day_sin', 'hour_day_cos', 'temp', 'holiday']

El RFECV de scikit-learn empieza entrenando un modelo con todos los predictores disponibles y calculando la importancia de cada uno en base a los atributos como coef_ o feature_importances_. A continuación, se elimina el predictor menos importante y se realiza una validación cruzada para calcular el rendimiento del modelo con los predictores restantes. Este proceso se repite hasta que la eliminación de predictores adicionales no mejora la metrica de rendimiento elegida o se alcanza el min_features_to_select.

El resultado final es un subconjunto de predictores que idealmente equilibra la simplicidad del modelo y su capacidad predictiva, determinada por el proceso de validación cruzada.

El forecater se entrena y evalúa de nuevo utilizando el conjunto de predictores seleccionados.

In [20]:
# Crear forecaster con los predictores seleccionados
# ==============================================================================
forecaster = ForecasterAutoreg(
    regressor = XGBRegressor(**best_params, random_state=15926, enable_categorical=True),
    lags      = lags_seleccionados
)

# Backtesting con los predictores seleccionados y los datos de test
# ==============================================================================
metrica, predicciones = backtesting_forecaster(
    forecaster         = forecaster,
    y                  = datos['users'],
    exog               = datos[exog_seleccionadas],
    steps              = 36,
    metric             = 'mean_absolute_error',
    initial_train_size = len(datos[:fin_validacion]),
    refit              = False,
    n_jobs             = 'auto',
    verbose            = False,
    show_progress      = True
)
metrica
Out[20]:
mean_absolute_error
0 56.845887

El rendimiento del modelo sigue siendo similar al del modelo entrenado con todas las variables. Sin embargo, el modelo es ahora mucho más simple, lo que hará que sea más rápido de entrenar y menos propenso al sobreajuste. Para el resto del documento, el modelo se entrenará utilizando sólo las características más importantes.

In [21]:
# 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 más populares para la interpretabilidad de modelos: importancia específica del modelo, valores SHAP y gráficos de dependencia parcial.

In [22]:
# Create and train forecaster
# ==============================================================================
forecaster = ForecasterAutoreg(
    regressor = XGBRegressor(**best_params, random_state=15926, enable_categorical=True),
    lags      = lags_seleccionados
)
forecaster.fit(
    y    = datos.loc[:fin_validacion, 'users'],
    exog = datos.loc[:fin_validacion, exog_features]
)

Model-specific feature importance

In [23]:
# Extraer importancia de los predictores
# ==============================================================================
importancia = forecaster.get_feature_importances()
importancia.head(10)
Out[23]:
feature importance
0 lag_1 0.366732
6 lag_168 0.221177
4 lag_24 0.111368
10 hour_day_sin 0.076504
11 hour_day_cos 0.072937
5 lag_167 0.028321
7 lag_169 0.027888
13 holiday 0.022521
2 lag_3 0.015814
1 lag_2 0.014428

Warning

El método get_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 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.

In [24]:
# 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_features]
                    )
display(X_train.head(3))
display(y_train.head(3))
lag_1 lag_2 lag_3 lag_23 lag_24 lag_167 lag_168 lag_169 week_day_sin week_day_cos hour_day_sin hour_day_cos temp holiday
date_time
2011-01-15 01:00:00 28.0 27.0 36.0 1.0 5.0 16.0 16.0 25.0 -0.781832 0.62349 0.500000 0.866025 6.56 0.0
2011-01-15 02:00:00 20.0 28.0 27.0 1.0 1.0 7.0 16.0 16.0 -0.781832 0.62349 0.707107 0.707107 6.56 0.0
2011-01-15 03:00:00 12.0 20.0 28.0 1.0 1.0 1.0 7.0 16.0 -0.781832 0.62349 0.866025 0.500000 6.56 0.0
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
In [25]:
# 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. El shap.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..
In [26]:
# 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)
fig.set_size_inches(7, 3)
In [27]:
# Force plot para la primera observación de la muestra
# ==============================================================================
shap.force_plot(explainer.expected_value, shap_values[0,:], X_train_sample.iloc[0,:])
Out[27]:
Visualization omitted, Javascript library not loaded!
Have you run `initjs()` in this notebook? If this notebook was from another user you must also trust this notebook (File -> Trust notebook). If you are viewing this notebook on github the Javascript has been stripped for security. If you are using JupyterLab this error is because a JupyterLab extension has not yet been written.

Información de sesión

In [28]:
import session_info
session_info.show(html=False)
-----
matplotlib          3.9.1
numpy               2.0.1
optuna              3.6.1
pandas              2.2.2
plotly              5.23.0
session_info        1.0.0
shap                0.46.0
skforecast          0.13.0
sklearn             1.5.1
statsmodels         0.14.2
xgboost             2.1.0
-----
IPython             8.26.0
jupyter_client      8.6.2
jupyter_core        5.7.2
notebook            6.4.12
-----
Python 3.12.4 | packaged by Anaconda, Inc. | (main, Jun 18 2024, 15:12:24) [GCC 11.2.0]
Linux-5.15.0-1066-aws-x86_64-with-glibc2.31
-----
Session information updated at 2024-08-10 18:49

Bibliografía

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

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

¿Cómo citar skforecast?

Zenodo:

Amat Rodrigo, Joaquin, & Escobar Ortiz, Javier. (2023). skforecast (v0.13.0). Zenodo. https://doi.org/10.5281/zenodo.8382787

APA:

Amat Rodrigo, J., & Escobar Ortiz, J. (2023). skforecast (Version 0.13.0) [Computer software]. https://doi.org/10.5281/zenodo.8382787

BibTeX:

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


¿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! 😊


Creative Commons Licence
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.

  • NoComercial: No puedes utilizar el material para fines comerciales.

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