Si te gusta Skforecast , ayúdanos dándonos una estrella en GitHub! ⭐️
Forecasting series temporales con XGBoost
Más sobre forecasting en: cienciadedatos.net
- Forecasting series temporales con machine learning
- Modelos ARIMA y SARIMAX
- Forecasting series temporales con gradient boosting: 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 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.
# 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.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
# 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__}")
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)
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.
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)})"
)
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
# 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",
width=750,
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.
# 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 = "Usiarios 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)
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.
Modelo autoregresivo recursivo con XGBoost¶
Primero, se entrena un modelo ForecasterRecursive 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¶
# Crear el forecaster
# ==============================================================================
window_features = RollingFeatures(stats=["mean"], window_sizes=24 * 3)
forecaster = ForecasterRecursive(
regressor = XGBRegressor(random_state=15926, enable_categorical=True),
lags = 24,
window_features = window_features
)
# 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
# ==============================================================================
cv = TimeSeriesFold(
steps = 36,
initial_train_size = len(datos[:fin_validacion]),
refit = False,
)
metrica, predicciones = backtesting_forecaster(
forecaster = forecaster,
y = datos['users'],
cv = cv,
metric = 'mean_absolute_error',
n_jobs = 'auto',
verbose = False, # Cambiar a False para mostrar más información
show_progress = True
)
predicciones.head()
# Métrica de backtesting
# ==============================================================================
metrica
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:# 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',
n_jobs = 'auto',
verbose = False,
show_progress = True
)
metrica
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 ForecasterRecursive 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:
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 lags que dé el menor error.
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ón0.12.0 los lags se incluye dentro del espacio de búsqueda (search_space).
# Búsqueda de hiperparámetros
# ==============================================================================
# Crear forecaster
forecaster = ForecasterRecursive(
regressor = XGBRegressor(random_state=15926, enable_categorical=True),
lags = 24,
window_features = window_features
)
# 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
# Folds
cv_searh = TimeSeriesFold(
steps = 36,
initial_train_size = len(datos_train),
refit = False,
)
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,
cv = cv_searh,
metric = 'mean_absolute_error',
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]
# Resultados de la búsqueda de hiperparámetros
# ==============================================================================
resultados_busqueda.head(3)
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
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.
# 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',
n_jobs = 'auto',
verbose = False,
show_progress = True
)
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="Predicción vs valor real en los datos de test",
xaxis_title="Fecha",
yaxis_title="Usiarios",
width=750,
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()