Forecasting series temporales con gradient boosting: Skforecast, XGBoost, LightGBM y CatBoost

Forecasting series temporales con gradient boosting: Skforecast, XGBoost, LightGBM y CatBoost

Joaquín Amat Rodrigo
Diciembre, 2021 (última actualización Enero 2022)

Introducción


Los modelos basados en gradient boosting destacan dentro de la comunidad de machine learning por los buenos resultados que consiguen en multitud de casos de uso, tanto de regresión como de clasificación. Aunque su uso en el forecasting de series temporales ha sido limitado, en los últimos años se ha demostrado que también pueden conseguir resultados muy competitivos. Algunas de las ventajas que presenta el uso de este tipo de modelos son:

  • Facilidad para incluir variables exógenas además de las autorregresivas.

  • Permiten incorporar en el modelo comportamientos no lineales.

  • Elevada escalabilidad, lo que permite aplicarlos cuando se dispone de grandes volúmenes de datos.

A pesar de sus potenciales ventajas, a la hora de aplicar modelos de machine learning a problemas de forecasting aparecen una serie de complicaciones que pueden hacer que el analista sea reticente a su uso, las principales son:

  • Reestructurar los datos de forma que puedan ser utilizados como si de un problema de regresión se tratara.

  • Dependiendo de cuántas predicciones a futuro se quieran obtener (horizonte de predicción), se tienen que programar procesos iterativos en los que, cada nueva predicción, utiliza las predicciones anteriores.

  • La validación de los modelos requiere de estrategias específicas como backtesting, walk-forward validation o time series cross-validation. El uso de la validación cruzada tradicional no aplica en estos casos.

La librería skforecast automatiza muchos de estos procesos, facilitando el uso y validación modelos de machine learning en problemas de forecasting. A lo largo de este documento, se muestra cómo utilizar tres de las implementaciones más avanzadas de modelos gradient boosting: XGBoost, LightGBM y Catboost.

Nota: no siempre los modelos de machine learning superan a modelos modelos propios del aprendizaje estadístico como AR, ARIMA o Exponential Smoothing. Cuál funcione mejor depende en gran medida de las características del caso de uso al que se apliquen.



Caso de uso


Los sistemas de bicicletas compartidas, también conocidos como sistemas de bicicletas públicas, facilitan la disponibilidad automática de bicicletas para que sean utilizadas temporalmente como medio de transporte. La mayoría de estos sistemas permiten recoger una bicicleta y devolverla en un punto diferente (estaciones o dockers), para que el usuario solo necesite tener la bicicleta en su posesión durante el desplazamiento. Uno de los principales retos en la gestión de estos sistemas es la necesidad de redistribuir las bicicletas para intentar que, en todas las estaciones, haya bicicletas disponibles a la vez que espacios libres para devoluciones.

Con el objetivo de mejorar la planificación y ejecución de la distribución de las bicicletas, se plantea crear un modelo capaz de predecir el número de usuarios para las siguientes 36 horas. De esta forma, a las 12h de cada día, la compañía encargada de gestionar las estaciones de alquiler, podrá conocer la demanda prevista el resto del día (12 horas) y el siguiente día (24 horas). Para simplificar el ejemplo mostrado en este documento, se asume que se está modelando una única estación.

Librerías


Las librerías utilizadas en este documento son:

In [1]:
# Tratamiento de datos
# ==============================================================================
import numpy as np
import pandas as pd

# Gráficos
# ==============================================================================
import matplotlib.pyplot as plt
from statsmodels.graphics.tsaplots import plot_acf
from statsmodels.graphics.tsaplots import plot_pacf
import plotly.express as px
plt.style.use('fivethirtyeight')
plt.rcParams['lines.linewidth'] = 1.5
%matplotlib inline

# Modelado y Forecasting
# ==============================================================================
from xgboost import XGBRegressor
from lightgbm import LGBMRegressor
from catboost import CatBoostRegressor

from sklearn.preprocessing import OneHotEncoder
from sklearn.preprocessing import StandardScaler
from sklearn.compose import ColumnTransformer
from sklearn.pipeline import make_pipeline

from skforecast.ForecasterAutoreg import ForecasterAutoreg
from skforecast.ForecasterAutoregMultiOutput import ForecasterAutoregMultiOutput
from skforecast.model_selection import grid_search_forecaster
from skforecast.model_selection import backtesting_forecaster

from joblib import dump, load

# Configuración warnings
# ==============================================================================
import warnings
warnings.filterwarnings('ignore')
%config Completer.use_jedi = False

La versión de skforecast utiliza es: skforecast v0.4.2.

Datos


Los datos empleados en este documento representan el uso, a nivel horario, del sistema de alquiler de bicicletas en la ciudad de Washington D.C. durante los años 2011 y 2012. Además del número de usuarios por hora, se dispone de información sobre las condiciones meteorológicas y sobre los días festivos. Los datos originales se han obtenido del UCI Machine Learning Repository y han sido limpiados previamente (código) aplicando las siguientes modificaciones :

  • Columnas renombradas con nombres más descriptivos.

  • Categorías de la variable meteorológica renombradas. La categoría de heavy rain, se ha combinado con la de rain.

  • Variables de temperatura, humedad y viento desnormalizadas.

  • Creada variable date_time y establecida como índice.

  • Imputación de valores missing mediante forward fill.

El data set resultante contiene las columnas:

  • date_time: fecha y hora.

  • month : mes (1 al 12).

  • hour : hora (0 al 23).

  • holiday : si el día es festivo o no (obtenido de from http://dchr.dc.gov/page/holiday-schedule).

  • weekday : día de la semana (Lunes=0, Domingo=6).

  • workingday : si es un día laboral.

  • weathersit : la meteorología del día (clear, mist, rain).

  • temp : temperatura registrada.

  • atemp: sensación térmica.

  • hum: humedad registrada.

  • windspeed: velocidad del viento registrada.

  • users: usuarios totales del servicio de alquiler de bicicletas.

In [2]:
# Descarga de datos
# ==============================================================================
url = ('https://raw.githubusercontent.com/JoaquinAmatRodrigo/Estadistica-machine-'
       'learning-python/master/data/bike_sharing_dataset_clean.csv')
datos = pd.read_csv(url)
datos['date_time'] = pd.to_datetime(datos['date_time'], format='%Y-%m-%d %H:%M:%S')
datos = datos.set_index('date_time')
datos = datos.asfreq('H')
datos = datos.sort_index()
datos
Out[2]:
holiday workingday weather temp atemp hum windspeed users month hour weekday
date_time
2011-01-01 00:00:00 0.0 0.0 clear 9.84 14.395 81.0 0.0000 16.0 1 0 5
2011-01-01 01:00:00 0.0 0.0 clear 9.02 13.635 80.0 0.0000 40.0 1 1 5
2011-01-01 02:00:00 0.0 0.0 clear 9.02 13.635 80.0 0.0000 32.0 1 2 5
2011-01-01 03:00:00 0.0 0.0 clear 9.84 14.395 75.0 0.0000 13.0 1 3 5
2011-01-01 04:00:00 0.0 0.0 clear 9.84 14.395 75.0 0.0000 1.0 1 4 5
... ... ... ... ... ... ... ... ... ... ... ...
2012-12-31 19:00:00 0.0 1.0 mist 10.66 12.880 60.0 11.0014 119.0 12 19 0
2012-12-31 20:00:00 0.0 1.0 mist 10.66 12.880 60.0 11.0014 89.0 12 20 0
2012-12-31 21:00:00 0.0 1.0 clear 10.66 12.880 60.0 11.0014 90.0 12 21 0
2012-12-31 22:00:00 0.0 1.0 clear 10.66 13.635 56.0 8.9981 61.0 12 22 0
2012-12-31 23:00:00 0.0 1.0 clear 10.66 13.635 65.0 8.9981 49.0 12 23 0

17544 rows × 11 columns

Con el objetivo de poder entrenar los modelos, hacer búsqueda de los mejores hiperparámetros y, finalmente, evaluar su capacidad predictiva, se reparten los datos en tres conjuntos: entrenamiento, validación y test.

In [3]:
# Separación datos train-val-test
# ==============================================================================
fin_train = '2012-03-31 23:59:00'
fin_validacion = '2012-08-31 23:59:00'
datos_train = datos.loc[: fin_train, :]
datos_val   = datos.loc[fin_train:fin_validacion, :]
datos_test  = datos.loc[fin_validacion:, :]

print(f"Fechas train      : {datos_train.index.min()} --- {datos_train.index.max()}  (n={len(datos_train)})")
print(f"Fechas validacion : {datos_val.index.min()} --- {datos_val.index.max()}  (n={len(datos_val)})")
print(f"Fechas test       : {datos_test.index.min()} --- {datos_test.index.max()}  (n={len(datos_test)})")
Fechas train      : 2011-01-01 00:00:00 --- 2012-03-31 23:00:00  (n=10944)
Fechas validacion : 2012-04-01 00:00:00 --- 2012-08-31 23:00:00  (n=3672)
Fechas test       : 2012-09-01 00:00:00 --- 2012-12-31 23:00:00  (n=2928)

Exploración gráfica

La exploración gráfica de una serie temporal es muy útil para identificar tendencias, patrones y estacionalidad, lo que resulta muy útil para orientar la selección de lags temporales que pueden ser buenos predictores.

Representación de la serie temporal

In [4]:
# Gráfico serie temporal
# ==============================================================================
fig, ax = plt.subplots(figsize=(11, 4))
datos_train['users'].plot(ax=ax, label='entrenamiento')
datos_val['users'].plot(ax=ax, label='validación')
datos_test['users'].plot(ax=ax, label='test')
ax.set_title('Número de usuarios')
ax.legend();
In [5]:
# Gráfico serie temporal con zoom
# ==============================================================================
zoom = ('2011-08-01 00:00:00','2011-08-15 00:00:00')

fig = plt.figure(figsize=(11, 6))
grid = plt.GridSpec(nrows=8, ncols=1, hspace=0.6, wspace=0)

main_ax = fig.add_subplot(grid[1:3, :])
zoom_ax = fig.add_subplot(grid[5:, :])

datos_train['users'].plot(ax=main_ax, label='entrenamiento', alpha=0.5)
datos_val['users'].plot(ax=main_ax, label='validación', alpha=0.5)
datos_test['users'].plot(ax=main_ax, label='test', alpha=0.5)
min_y = min(datos['users'])
max_y = max(datos['users'])
main_ax.fill_between(zoom, min_y, max_y, facecolor='blue', alpha=0.5, zorder=0)
main_ax.set_xlabel('')
main_ax.legend(loc='lower center', ncol=3, bbox_to_anchor=(0.5, -0.8))

datos.loc[zoom[0]: zoom[1]]['users'].plot(ax=zoom_ax, color='blue', linewidth=2)

main_ax.set_title(f'Número de usuarios: {datos.index.min()}, {datos.index.max()}', fontsize=14)
zoom_ax.set_title(f'Número de usuarios: {zoom}', fontsize=14)
plt.subplots_adjust(hspace=1)
In [6]:
# Gráfico interactivo serie temporal
# ==============================================================================
datos.loc[:fin_train, 'particion'] = 'entrenamiento'
datos.loc[fin_train:fin_validacion, 'particion'] = 'validación'
datos.loc[fin_validacion:, 'particion'] = 'test'

fig = px.line(
    data_frame = datos.reset_index(),
    x      = 'date_time',
    y      = 'users',
    color  = 'particion',
    title  = 'Número de usuarios',
    width  = 900,
    height = 500
)

fig.update_xaxes(
    rangeslider_visible=True,
    rangeselector=dict(
        buttons=list([
            dict(count=1, label="1m", step="month", stepmode="backward"),
            dict(count=6, label="6m", step="month", stepmode="backward"),
            dict(count=1, label="YTD", step="year", stepmode="todate"),
            dict(count=1, label="1y", step="year", stepmode="backward"),
            dict(step="all")
        ])
    )
)

fig.show()

datos=datos.drop(columns='particion')