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

# Data processing
# ==============================================================================
import numpy as np
import pandas as pd
from astral.sun import sun
from astral import LocationInfo
from skforecast.datasets import fetch_dataset

# Plots
# ==============================================================================
import matplotlib.pyplot as plt
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from skforecast.plot import plot_residuals, 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})

# Modelling and Forecasting
# ==============================================================================
import xgboost
import lightgbm
import catboost
import sklearn
import shap
from xgboost import XGBRegressor
from lightgbm import LGBMRegressor
from catboost import CatBoostRegressor
from sklearn.ensemble import HistGradientBoostingRegressor
from sklearn.preprocessing import OneHotEncoder, PolynomialFeatures
from skforecast.preprocessing import CalendarFeatures
from feature_engine.timeseries.forecasting import WindowFeatures
from sklearn.feature_selection import RFECV
from sklearn.compose import make_column_transformer, make_column_selector
import skforecast
from skforecast.recursive import ForecasterEquivalentDate, ForecasterRecursive
from skforecast.model_selection import (
    TimeSeriesFold,
    OneStepAheadFold,
    bayesian_search_forecaster,
    backtesting_forecaster,
)
from skforecast.preprocessing import RollingFeatures
from skforecast.feature_selection import select_features
from skforecast.stats import calculate_lag_autocorrelation
from skforecast.metrics import calculate_coverage
from skforecast.exceptions import IgnoredArgumentWarning

# 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__}')
print(f'{color}Versión xgboost: {xgboost.__version__}')
print(f'{color}Versión catboost: {catboost.__version__}')
print(f'{color}Versión pandas: {pd.__version__}')
print(f'{color}Versión numpy: {np.__version__}')

# Descarga de datos
# ==============================================================================
datos = fetch_dataset('bike_sharing', raw=True)

# Preprocesado de datos (estableciendo índice y frecuencia)
# ==============================================================================
datos = datos[
    ['date_time', 'users', 'holiday', 'weather', 'temp', 'atemp', 'hum', 'windspeed']
].copy()
datos['date_time'] = pd.to_datetime(datos['date_time'], format='%Y-%m-%d %H:%M:%S')
datos = datos.set_index('date_time')
datos = datos.asfreq('h')
datos = datos.sort_index()
datos.head()

# Separación de datos en entrenamiento, validación y test
# ==============================================================================
fin_train = '2012-04-30 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)})'
    )

# 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áfico de la serie temporal con zoom
# ==============================================================================
zoom = ('2011-08-01 00:00:00', '2011-08-15 00:00:00')
fig, axs = plt.subplots(2, 1, figsize=(8, 4), gridspec_kw={'height_ratios': [1, 2]})
datos_train['users'].plot(ax=axs[0], label='train', alpha=0.5)
datos_val['users'].plot(ax=axs[0], label='validation', alpha=0.5)
datos_test['users'].plot(ax=axs[0], label='test', alpha=0.5)
axs[0].axvspan(zoom[0], zoom[1], color='blue', alpha=0.7)
axs[0].set_title('Número de usuarios')
axs[0].set_xlabel('')
# Zoom plot
datos.loc[zoom[0] : zoom[1], 'users'].plot(ax=axs[1], color='blue')
axs[1].set_title(f'Zoom: {zoom[0]} to {zoom[1]}', fontsize=10)
plt.tight_layout()
plt.show()

# 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()

# 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()

# Top 10 lags con mayor autocorrelación parcial absoluta
# ==============================================================================
calculate_lag_autocorrelation(
    data    = datos['users'],
    n_lags  = 72,
    sort_by = 'partial_autocorrelation_abs'
).head(10)

# Crear un baseline: valor de la misma hora del día anterior
# ==============================================================================
forecaster = ForecasterEquivalentDate(
                offset    = pd.DateOffset(days=1),
                n_offsets = 1
             )

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

# Backtesting
# ==============================================================================
cv = TimeSeriesFold(steps = 36, initial_train_size = len(datos.loc[:fin_validacion]))
metrica_baseline, predicciones = backtesting_forecaster(
                                    forecaster = forecaster,
                                    y          = datos['users'],
                                    cv         = cv,
                                    metric     = 'mean_absolute_error'
                                )
metrica_baseline

# 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
             )

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

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

# 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

# Variables basadas en el calendario, con encoding cíclico
# ==============================================================================
calendar_transformer = CalendarFeatures(
    features = ['month', 'week', 'day_of_week', 'hour'],
    encoding = 'cyclical',
    keep_original_columns = False,
)
variables_calendario = calendar_transformer.fit_transform(datos)
variables_calendario.head(2)

# Variables basadas en la luz solar
# ==============================================================================
location = LocationInfo(
    name      = 'Washington D.C.',
    region    = 'USA',
    latitude  = 38.89,
    longitude = -77.04,
    timezone  = 'America/New_York'
)
# El amanecer y el atardecer se calculan una vez por día y se asignan a cada hora
days = datos.index.normalize().unique()
sun_times = [
    sun(location.observer, date=day, tzinfo=location.timezone) for day in days
]
sunrise_hour = pd.Series([s['sunrise'] for s in sun_times], index=days)
sunset_hour = pd.Series([s['sunset'] for s in sun_times], index=days)
sunrise_hour = sunrise_hour.dt.round('h').dt.hour.reindex(datos.index, method='ffill')
sunset_hour = sunset_hour.dt.round('h').dt.hour.reindex(datos.index, method='ffill')
sunrise_hour_sin = np.sin(2 * np.pi * sunrise_hour / 24)
sunrise_hour_cos = np.cos(2 * np.pi * sunrise_hour / 24)
sunset_hour_sin = np.sin(2 * np.pi * sunset_hour / 24)
sunset_hour_cos = np.cos(2 * np.pi * sunset_hour / 24)

daylight_hours = sunset_hour - sunrise_hour
is_daylight = np.where(
    (datos.index.hour >= sunrise_hour) & (datos.index.hour < sunset_hour), 1, 0,
)

variables_solares = pd.DataFrame({
                        'sunrise_hour_sin': sunrise_hour_sin,
                        'sunrise_hour_cos': sunrise_hour_cos,
                        'sunset_hour_sin': sunset_hour_sin,
                        'sunset_hour_cos': sunset_hour_cos,
                        'daylight_hours': daylight_hours,
                        'is_daylight': is_daylight
                     })
variables_solares.head(2)

# Variables basadas en festivos
# ==============================================================================
variables_festivos = datos[['holiday']].astype(int)
variables_festivos['holiday_previous_day'] = variables_festivos['holiday'].shift(24)
variables_festivos['holiday_next_day'] = variables_festivos['holiday'].shift(-24)
variables_festivos.head(2)

# Variables basadas en temperatura
# ==============================================================================
wf_transformer = WindowFeatures(
    variables = ['temp'],
    window    = ['1D', '7D'],
    functions = ['mean', 'max', 'min'],
    freq      = 'h',
)
variables_temp = wf_transformer.fit_transform(datos[['temp']])
variables_temp.head(2)

# Unión de variables exógenas
# ==============================================================================
assert all(variables_calendario.index == variables_solares.index)
assert all(variables_calendario.index == variables_festivos.index)
assert all(variables_calendario.index == variables_temp.index)
variables_exogenas = pd.concat([
    variables_calendario,
    variables_solares,
    variables_temp,
    variables_festivos
], axis=1)

# Debido a la creación de ventanas móviles, hay valores faltantes al principio
# de la serie. Y debido a holiday_next_day hay valores faltantes al final.
variables_exogenas = variables_exogenas.iloc[7 * 24:, :]
variables_exogenas = variables_exogenas.iloc[:-24, :]
variables_exogenas.head(3)

# Interacción entre variables exógenas
# ==============================================================================
transformer_poly = PolynomialFeatures(
                        degree           = 2,
                        interaction_only = True,
                        include_bias     = False
                    ).set_output(transform='pandas')
# Las interacciones se crean únicamente entre las variables cíclicas
poly_cols = [
    'month_sin',
    'month_cos',
    'week_sin',
    'week_cos',
    'day_of_week_sin',
    'day_of_week_cos',
    'hour_sin',
    'hour_cos',
    'sunrise_hour_sin',
    'sunrise_hour_cos',
    'sunset_hour_sin',
    'sunset_hour_cos',
]
variables_poly = transformer_poly.fit_transform(variables_exogenas[poly_cols])
variables_poly = variables_poly.drop(columns=poly_cols)
variables_poly.columns = [f'poly_{col}' for col in variables_poly.columns]
variables_poly.columns = variables_poly.columns.str.replace(' ', '__')
assert all(variables_exogenas.index == variables_poly.index)
variables_exogenas = pd.concat([variables_exogenas, variables_poly], axis=1)
variables_exogenas.head(3)

# Almacenar las variables categóricas como tipo "category"
# ==============================================================================
datos['weather'] = datos['weather'].astype('category')

# Transformación con codificación one-hot
# ==============================================================================
one_hot_encoder = make_column_transformer(
    (
        OneHotEncoder(sparse_output=False, drop='if_binary'),
        make_column_selector(dtype_include=['category', 'object']),
    ),
    remainder='passthrough',
    verbose_feature_names_out=False,
).set_output(transform='pandas')

# Crear un forecaster con un transformer para las variables exógenas
# ==============================================================================
forecaster = ForecasterRecursive(
                estimator        = LGBMRegressor(random_state=15926, verbose=-1),
                lags             = 72,
                window_features  = window_features,
                transformer_exog = one_hot_encoder
             )

# Mostrar matrices de entrenamiento
# ==============================================================================
exog_cols = ['weather']
X_train, y_train = forecaster.create_train_X_y(
                        y    = datos.loc[:fin_validacion, 'users'],
                        exog = datos.loc[:fin_validacion, exog_cols]
                   )
X_train.head(3)

# Crear un forecaster con detección automática de variables categóricas
# ==============================================================================
forecaster = ForecasterRecursive(
                estimator            = LGBMRegressor(random_state=15926, verbose=-1),
                lags                 = 72,
                window_features      = window_features,
                categorical_features = 'auto'
             )

# Selección de variables exógenas a incluir en el modelo
# ==============================================================================
exog_cols = []
# Variables cíclicas: columnas que terminan con _sin o _cos. Sus interacciones
# (columnas que empiezan con poly_) también terminan con _sin o _cos, por lo que
# quedan seleccionadas
exog_cols.extend(variables_exogenas.filter(regex='_sin$|_cos$').columns.tolist())
# Columnas que empiezan con temp_ son seleccionadas
exog_cols.extend(variables_exogenas.filter(regex='^temp_.*').columns.tolist())
# Columnas que empiezan con holiday_ son seleccionadas
exog_cols.extend(variables_exogenas.filter(regex='^holiday_.*').columns.tolist())
exog_cols.extend(['daylight_hours', 'is_daylight', 'temp', 'holiday', 'weather'])

variables_exogenas = variables_exogenas.filter(exog_cols, axis=1)

# Combinar variables exógenas y target en el mismo dataframe
# ==============================================================================
datos = datos[['users', 'weather']].merge(
            variables_exogenas,
            left_index  = True,
            right_index = True,
            how         = 'inner'  # Solo fechas para las que hay datos exógenos
        )
datos = datos.astype({col: np.float32 for col in datos.select_dtypes('number').columns})
datos_train = datos.loc[: fin_train, :].copy()
datos_val   = datos.loc[fin_train:fin_validacion, :].copy()
datos_test  = datos.loc[fin_validacion:, :].copy()

# Backtesting en los datos de test incluyendo las variables exógenas
# ==============================================================================
cv = TimeSeriesFold(steps = 36, initial_train_size = len(datos.loc[:fin_validacion]))
metrica, predicciones = backtesting_forecaster(
                            forecaster = forecaster,
                            y          = datos['users'],
                            exog       = datos[exog_cols],
                            cv         = cv,
                            metric     = 'mean_absolute_error'
                        )
metrica

# Búsqueda de hiperparámetros
# ==============================================================================
forecaster = ForecasterRecursive(
                estimator            = LGBMRegressor(random_state=15926, verbose=-1),
                lags                 = 72,
                window_features      = window_features,
                categorical_features = 'auto'
             )

# Lags grid
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', 300, 1000, 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),
        'max_bin'         : trial.suggest_int('max_bin', 50, 250),
        '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)
    }

# 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'],
    exog          = datos.loc[:fin_validacion, exog_cols],
    cv            = cv_search,
    search_space  = search_space,
    metric        = 'mean_absolute_error',
    n_trials      = 20,
    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)

# Backtesting en 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'],
                            exog       = datos[exog_cols],
                            cv         = cv,
                            metric     = 'mean_absolute_error'
                        )
metrica

# 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',
    xaxis_title='Date time',
    yaxis_title='Users',
    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()

# 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_cols],
    select_only     = None,
    force_inclusion = None,
    subsample       = 0.5,
    random_state    = 123,
    verbose         = True,
)

# 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
# ==============================================================================
cv = TimeSeriesFold(steps = 36, initial_train_size = len(datos.loc[:fin_validacion]))
metrica_lgbm, predicciones = backtesting_forecaster(
    forecaster = forecaster,
    y          = datos['users'],
    exog       = datos[exog_seleccionadas],
    cv         = cv,
    metric     = 'mean_absolute_error'
)
metrica_lgbm

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

# Crear y entrenar forecaster
# ==============================================================================
forecaster = ForecasterRecursive(
    estimator            = LGBMRegressor(**best_params),
    lags                 = best_lags,
    window_features      = window_features,
    categorical_features = 'auto',
    binner_kwargs        = {'n_bins': 5}
)
forecaster.fit(
    y    = datos.loc[:fin_train, 'users'],
    exog = datos.loc[:fin_train, exog_cols],
    store_in_sample_residuals = True
)

# Predicción de intervalos
# ==============================================================================
# Como el modelo ha sido entrenado con variables exógenas, se tienen que pasar
# para las predicciones.
predicciones = forecaster.predict_interval(
    exog     = datos.loc[fin_train:, exog_cols],
    steps    = 36,
    interval = [0.05, 0.95],
    method   = 'conformal',
)
predicciones.head()

# Backtesting con los datos de validación para obtener residuos out-sample
# ==============================================================================
cv = TimeSeriesFold(steps = 36, initial_train_size = len(datos.loc[:fin_train]))
_, predicciones_val = backtesting_forecaster(
    forecaster = forecaster,
    y          = datos.loc[:fin_validacion, 'users'],
    exog       = datos.loc[:fin_validacion, exog_cols],
    cv         = cv,
    metric     = 'mean_absolute_error'
)

# Distribución de los residuos out-sample
# ==============================================================================
residuals = datos.loc[predicciones_val.index, 'users'] - predicciones_val['pred']
print(pd.Series(np.where(residuals < 0, 'negative', 'positive')).value_counts())
plt.rcParams.update({'font.size': 8})
_ = plot_residuals(
        y_true = datos.loc[predicciones_val.index, 'users'],
        y_pred = predicciones_val['pred'],
        figsize=(7, 4)
    )

# Almacenar residuos out-sample en el forecaster
# ==============================================================================
forecaster.set_out_sample_residuals(
    y_true = datos.loc[predicciones_val.index, 'users'],
    y_pred = predicciones_val['pred']
)

# Backtesting con intervalos de predicción en test usando out-sample residuals
# ==============================================================================
cv = TimeSeriesFold(steps = 36, initial_train_size = len(datos.loc[:fin_validacion]))
metrica, predicciones = backtesting_forecaster(
   forecaster              = forecaster,
   y                       = datos['users'],
   exog                    = datos[exog_cols],
   cv                      = cv,
   metric                  = 'mean_absolute_error',
   interval                = [0.05, 0.95],  # 90% intervalo de predicción
   interval_method         = 'conformal',
   use_in_sample_residuals = False,  # Usar out-sample residuals
   use_binned_residuals    = True,   # Usar residuos condicionados al valor predicho
)
predicciones.head(5)

# Gráfico intervalos de predicción vs valor real
# ==============================================================================
fig = go.Figure([
    go.Scatter(
        name='Prediction', x=predicciones.index, y=predicciones['pred'], mode='lines'
    ),
    go.Scatter(
        name='Real value', x=datos_test.index, y=datos_test['users'], mode='lines',
    ),
    go.Scatter(
        name='Upper Bound', x=predicciones.index, y=predicciones['upper_bound'],
        mode='lines', marker=dict(color='#444'), line=dict(width=0), showlegend=False
    ),
    go.Scatter(
        name='Lower Bound', x=predicciones.index, y=predicciones['lower_bound'],
        mode='lines', marker=dict(color='#444'), line=dict(width=0),
        fillcolor='rgba(68, 68, 68, 0.3)', fill='tonexty', showlegend=False
    )
])
fig.update_layout(
    title='Valor real vs predicciones en los datos de test',
    xaxis_title='Date time',
    yaxis_title='users',
    width=800,
    height=400,
    margin=dict(l=20, r=20, t=35, b=20),
    hovermode='x',
    legend=dict(orientation='h', yanchor='top', y=1.1, xanchor='left', x=0.001),
    # Zoom inicial en el eje x del 1 al 10 de octubre
    xaxis=dict(range=['2012-10-01', '2012-10-10'])
)
fig.show()

# Cobertura del intervalo en los datos de test
# ==============================================================================
cobertura = calculate_coverage(
              y_true       = datos.loc[fin_validacion:, 'users'],
              lower_bound  = predicciones['lower_bound'],
              upper_bound  = predicciones['upper_bound']
            )
area = (predicciones['upper_bound'] - predicciones['lower_bound']).sum()
print(f'Área total del intervalo: {round(area, 2)}')
print(f'Cobertura del intervalo predicho: {round(100 * cobertura, 2)} %')

# 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_cols]
)

# Extraer importancia de los predictores
# ==============================================================================
importancia = forecaster.get_feature_importances()
importancia.head(10)

# 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_cols]
                    )
display(X_train.head(3))
display(y_train.head(3))

# 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)

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

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_cols].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:],
)

# Particiones utilizadas para la búsqueda de hiperparámetros y backtesting
# ==============================================================================
cv_search = OneStepAheadFold(initial_train_size = len(datos_train))
cv_backtesting = TimeSeriesFold(
    steps = 36, initial_train_size = len(datos.loc[:fin_validacion])
)

# Crear forecaster
# ==============================================================================
forecaster = ForecasterRecursive(
    estimator            = XGBRegressor(tree_method='hist', random_state=123),
    lags                 = 24,
    window_features      = window_features,
    categorical_features = 'auto'
)

# Búsqueda de hiperparámetros
# ==============================================================================
warnings.simplefilter('ignore', category=IgnoredArgumentWarning)

# Se reutiliza el grid de lags (`lags_grid`) definido en la búsqueda de LightGBM
# Espacio de búsqueda de hiperparámetros
def search_space(trial):
    return {
        'n_estimators'    : trial.suggest_int('n_estimators', 300, 1000, step=100),
        'max_depth'       : trial.suggest_int('max_depth', 3, 10),
        '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)
    }

results_search, frozen_trial = bayesian_search_forecaster(
    forecaster   = forecaster,
    y            = datos.loc[:fin_validacion, 'users'],
    exog         = datos.loc[:fin_validacion, exog_cols],
    cv           = cv_search,
    search_space = search_space,
    metric       = 'mean_absolute_error',
    n_trials     = 20
)

# Backtesting con datos de test
# ==============================================================================
metrica_xgboost, predicciones = backtesting_forecaster(
    forecaster = forecaster,
    y          = datos['users'],
    exog       = datos[exog_cols],
    cv         = cv_backtesting,
    metric     = 'mean_absolute_error'
)
metrica_xgboost

# Crear forecaster
# ==============================================================================
forecaster = ForecasterRecursive(
    estimator            = HistGradientBoostingRegressor(random_state=123),
    lags                 = 24,
    window_features      = window_features,
    categorical_features = 'auto'
)

# Búsqueda de hiperparámetros
# ==============================================================================
# Se reutiliza el grid de lags (`lags_grid`) definido en la búsqueda de LightGBM
# Espacio de búsqueda de hiperparámetros
def search_space(trial):
    return {
        'max_iter'          : trial.suggest_int('max_iter', 300, 1000, step=100),
        'max_depth'         : trial.suggest_int('max_depth', 3, 10),
        'learning_rate'     : trial.suggest_float('learning_rate', 0.01, 1),
        'min_samples_leaf'  : trial.suggest_int('min_samples_leaf', 1, 20),
        'l2_regularization' : trial.suggest_float('l2_regularization', 0, 1),
        'lags'              : trial.suggest_categorical('lags', lags_grid)
    }

results_search, frozen_trial = bayesian_search_forecaster(
    forecaster    = forecaster,
    y             = datos.loc[:fin_validacion, 'users'],
    exog          = datos.loc[:fin_validacion, exog_cols],
    cv            = cv_search,
    search_space  = search_space,
    metric        = 'mean_absolute_error',
    n_trials      = 20
)

# Backtesting con datos de test
# ==============================================================================
metrica_histgb, predicciones = backtesting_forecaster(
    forecaster = forecaster,
    y          = datos['users'],
    exog       = datos[exog_cols],
    cv         = cv_backtesting,
    metric     = 'mean_absolute_error'
)
metrica_histgb

# Crear forecaster
# ==============================================================================
forecaster = ForecasterRecursive(
    estimator            = CatBoostRegressor(
                                random_state=123,
                                silent=True,
                                allow_writing_files=False,
                                boosting_type = 'Plain',         # Faster training
                                leaf_estimation_iterations = 3,  # Faster training
                                one_hot_max_size = 10,           # Faster training
                            ),
    lags                 = 24,
    window_features      = window_features,
    categorical_features = 'auto'
)

# Búsqueda de hiperparámetros
# ==============================================================================
# Se reutiliza el grid de lags (`lags_grid`) definido en la búsqueda de LightGBM
# Espacio de búsqueda de hiperparámetros
def search_space(trial):
    return {
        'n_estimators'  : trial.suggest_int('n_estimators', 100, 1000, step=100),
        'max_depth'     : trial.suggest_int('max_depth', 3, 10),
        'learning_rate' : trial.suggest_float('learning_rate', 0.01, 1),
        'lags'          : trial.suggest_categorical('lags', lags_grid)
    }

results_search, frozen_trial = bayesian_search_forecaster(
    forecaster    = forecaster,
    y             = datos.loc[:fin_validacion, 'users'],
    exog          = datos.loc[:fin_validacion, exog_cols],
    cv            = cv_search,
    search_space  = search_space,
    metric        = 'mean_absolute_error',
    n_trials      = 20
)

# Backtesting con datos de test
# ==============================================================================
metrica_catboost, predicciones = backtesting_forecaster(
    forecaster = forecaster,
    y          = datos['users'],
    exog       = datos[exog_cols],
    cv         = cv_backtesting,
    metric     = 'mean_absolute_error'
)
metrica_catboost

# Comparación del MAE en test
# ==============================================================================
metricas = pd.concat(
    [metrica_baseline, metrica_lgbm, metrica_xgboost, metrica_histgb, metrica_catboost]
)
metricas.index = [
    'Baseline',
    'LGBMRegressor',
    'XGBRegressor',
    'HistGradientBoostingRegressor',
    'CatBoostRegressor',
]
metricas.round(2).sort_values(by='mean_absolute_error')

import session_info
session_info.show(html=False)