Forecasting time series with gradient boosting: Skforecast, XGBoost, LightGBM, Scikit-learn and CatBoost

Joaquín Amat Rodrigo, Javier Escobar OrtizFebruary, 2021 (last update November 2024)

More about forecasting

Introduction¶

Gradient boosting models have gained popularity in the machine learning community due to their ability to achieve excellent results in a wide range of use cases, including both regression and classification. Although these models have traditionally been less common in forecasting, they can be highly effective in this domain. Some of the key benefits of using gradient boosting models for forecasting include:

The ease with which exogenous variables can be included in the model, in addition to autoregressive variables.
The ability to capture non-linear relationships between variables.
High scalability, allowing models to handle large volumes of data.
Some implementations allow the inclusion of categorical variables without the need for additional encoding, such as one-hot encoding.

Despite these benefits, the use of machine learning models for forecasting can present several challenges that can make analysts reluctant to use them, the main ones being:

Transforming the data so that it can be used as a regression problem.
Depending on how many future predictions are needed (prediction horizon), an iterative process may be required where each new prediction is based on previous ones.
Model validation requires specific strategies such as backtesting, walk-forward validation or time series cross-validation. Traditional cross-validation cannot be used.

The skforecast library provides automated solutions to these challenges, making it easier to apply and validate machine learning models to forecasting problems. The library supports several advanced gradient boosting models, including XGBoost, LightGBM, Catboost and scikit-learn HistGradientBoostingRegressor. This document shows how to use them to build accurate forecasting models. To ensure a smooth learning experience, an initial exploration of the data is performed. Then, the modeling process is explained step by step, starting with a recursive model utilizing a LightGBM regressor and moving on to a model incorporating exogenous variables and various coding strategies. The document concludes by demonstrating the use of other gradient boosting model implementations, including XGBoost, CatBoost, and the scikit-learn HistGradientBoostingRegressor.

✎ Note

Machine learning models do not always outperform statistical learning models such as AR, ARIMA or Exponential Smoothing. Which one works best depends largely on the characteristics of the use case to which it is applied. Visit ARIMA and SARIMAX models with skforecast to learn more about statistical models.

✎ Note

Additional examples of how using gradient boosting models for forecasting can be found in the documents Forecasting energy demand with machine learning and Global Forecasting Models.

Use case¶

Bicycle sharing is a popular shared transport service that provides bicycles to individuals for short-term use. These systems typically provide bike docks where riders can borrow a bike and return it to any dock belonging to the same system. The docks are equipped with special bike racks that secure the bike and only release it via computer control. One of the major challenges faced by operators of these systems is the need to redistribute bikes to ensure that there are bikes available at all docks, as well as free spaces for returns.

In order to improve the planning and execution of bicycle distribution, it is proposed to create a model capable of forecasting the number of users over the next 36 hours. In this way, at 12:00 every day, the company in charge of managing the system will be able to know the expected users for the rest of the day (12 hours) and the next day (24 hours).

For illustrative purposes, the current example only models a single station, but it can be easily extended to cover multiple stations unsing global multi-series forecasting, thereby improving the management of bike-sharing systems on a larger scale.

Libraries¶

Libraries used in this document.

# 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
from feature_engine.datetime import DatetimeFeatures
from feature_engine.creation import CyclicalFeatures
from feature_engine.timeseries.forecasting import WindowFeatures

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

# 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
from sklearn.preprocessing import OrdinalEncoder
from sklearn.preprocessing import FunctionTransformer
from sklearn.preprocessing import PolynomialFeatures
from sklearn.feature_selection import RFECV
from sklearn.pipeline import make_pipeline
from sklearn.compose import make_column_transformer
from sklearn.compose import make_column_selector

import skforecast
from skforecast.recursive import ForecasterEquivalentDate
from skforecast.recursive import ForecasterRecursive
from skforecast.model_selection import TimeSeriesFold, OneStepAheadFold
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

# Warnings configuration
# ==============================================================================
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 lightgbm: {lightgbm.__version__}")
print(f"{color}Version xgboost: {xgboost.__version__}")
print(f"{color}Version catboost: {catboost.__version__}")
print(f"{color}Version pandas: {pd.__version__}")
print(f"{color}Version numpy: {np.__version__}")

Version skforecast: 0.14.0
Version scikit-learn: 1.5.2
Version lightgbm: 4.4.0
Version xgboost: 2.1.2
Version catboost: 1.2.7
Version pandas: 2.2.3
Version numpy: 1.26.4

⚠ Warning

At the time of writing this document, catboost is only compatible with numpy version lower than 2.0. If you have a higher version, you can downgrade it by running the following command: pip install numpy==1.26.4

Data¶

The data in this document represent the 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, information about weather conditions and holidays is available. The original data was obtained from the UCI Machine Learning Repository and has been previously cleaned (code) applying the following modifications:

Renamed columns with more descriptive names.
Renamed categories of the weather variables. The category of heavy_rain has been combined with that of rain.
Denormalized temperature, humidity and wind variables.
date_time variable created and set as index.
Imputed missing values by forward filling.

# Downloading data
# ==============================================================================
data = fetch_dataset('bike_sharing', raw=True)

bike_sharing
------------
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, information
about weather conditions and holidays is available.
Fanaee-T,Hadi. (2013). Bike Sharing Dataset. UCI Machine Learning Repository.
https://doi.org/10.24432/C5W894.
Shape of the dataset: (17544, 12)

# Preprocessing data (setting index and frequency)
# ==============================================================================
data = data[['date_time', 'users', 'holiday', 'weather', 'temp', 'atemp', 'hum', 'windspeed']]
data['date_time'] = pd.to_datetime(data['date_time'], format='%Y-%m-%d %H:%M:%S')
data = data.set_index('date_time')
if pd.__version__ < '2.2':
    data = data.asfreq('H')
else:
    data = data.asfreq('h')
data = data.sort_index()
data.head()

To facilitate the training of the models, the search for optimal hyperparameters and the evaluation of their predictive accuracy, the data are divided into three separate sets: training, validation and test.

# Split train-validation-test
# ==============================================================================
end_train = '2012-04-30 23:59:00'
end_validation = '2012-08-31 23:59:00'
data_train = data.loc[: end_train, :]
data_val   = data.loc[end_train:end_validation, :]
data_test  = data.loc[end_validation:, :]

print(f"Dates train      : {data_train.index.min()} --- {data_train.index.max()}  (n={len(data_train)})")
print(f"Dates validacion : {data_val.index.min()} --- {data_val.index.max()}  (n={len(data_val)})")
print(f"Dates test       : {data_test.index.min()} --- {data_test.index.max()}  (n={len(data_test)})")

Dates train      : 2011-01-01 00:00:00 --- 2012-04-30 23:00:00  (n=11664)
Dates validacion : 2012-05-01 00:00:00 --- 2012-08-31 23:00:00  (n=2952)
Dates test       : 2012-09-01 00:00:00 --- 2012-12-31 23:00:00  (n=2928)

Data exploration¶

Graphical exploration of time series can be an effective way of identifying trends, patterns, and seasonal variations. This, in turn, helps to guide the selection of the most appropriate forecasting model.

Plot time series¶

Full time series

# Interactive plot of time series
# ==============================================================================
fig = go.Figure()
fig.add_trace(go.Scatter(x=data_train.index, y=data_train['users'], mode='lines', name='Train'))
fig.add_trace(go.Scatter(x=data_val.index, y=data_val['users'], mode='lines', name='Validation'))
fig.add_trace(go.Scatter(x=data_test.index, y=data_test['users'], mode='lines', name='Test'))
fig.update_layout(
    title  = 'Number of users',
    xaxis_title="Time",
    yaxis_title="Users",
    legend_title="Partition:",
    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()

# Static plot of time series with zoom
# ==============================================================================
zoom = ('2011-08-01 00:00:00', '2011-08-15 00:00:00')
fig = plt.figure(figsize=(8, 4))
grid = plt.GridSpec(nrows=8, ncols=1, hspace=0.1, wspace=0)
main_ax = fig.add_subplot(grid[1:3, :])

data_train['users'].plot(ax=main_ax, label='train', alpha=0.5)
data_val['users'].plot(ax=main_ax, label='validation', alpha=0.5)
data_test['users'].plot(ax=main_ax, label='test', alpha=0.5)
min_y = min(data['users'])
max_y = max(data['users'])
main_ax.fill_between(zoom, min_y, max_y, facecolor='blue', alpha=0.5, zorder=0)
main_ax.set_title(f'Number of users: {data.index.min()}, {data.index.max()}', fontsize=10)
main_ax.set_xlabel('')
main_ax.legend(loc='lower center', ncol=3, bbox_to_anchor=(0.5, -0.8))
zoom_ax = fig.add_subplot(grid[5:, :])
data.loc[zoom[0]: zoom[1]]['users'].plot(ax=zoom_ax, color='blue', linewidth=1)
zoom_ax.set_title(f'Number of users: {zoom}', fontsize=10)
zoom_ax.set_xlabel('')
plt.subplots_adjust(hspace=1)
plt.show();

Seasonality plots¶

Seasonal plots are a useful tool for identifying seasonal patterns and trends in a time series. They are created by averaging the values of each season over time and then plotting them against time.

# Annual, weekly and daily seasonality
# ==============================================================================
fig, axs = plt.subplots(2, 2, figsize=(8, 5), sharex=False, sharey=True)
axs = axs.ravel()

# Users distribution by month
data['month'] = data.index.month
data.boxplot(column='users', by='month', ax=axs[0], flierprops={'markersize': 3, 'alpha': 0.3})
data.groupby('month')['users'].median().plot(style='o-', linewidth=0.8, ax=axs[0])
axs[0].set_ylabel('Users')
axs[0].set_title('Users distribution by month', fontsize=9)

# Users distribution by week day
data['week_day'] = data.index.day_of_week + 1
data.boxplot(column='users', by='week_day', ax=axs[1], flierprops={'markersize': 3, 'alpha': 0.3})
data.groupby('week_day')['users'].median().plot(style='o-', linewidth=0.8, ax=axs[1])
axs[1].set_ylabel('Users')
axs[1].set_title('Users distribution by week day', fontsize=9)

# Users distribution by the hour of the day
data['hour_day'] = data.index.hour + 1
data.boxplot(column='users', by='hour_day', ax=axs[2], flierprops={'markersize': 3, 'alpha': 0.3})
data.groupby('hour_day')['users'].median().plot(style='o-', linewidth=0.8, ax=axs[2])
axs[2].set_ylabel('Users')
axs[2].set_title('Users distribution by the hour of the day', fontsize=9)

# Users distribution by week day and hour of the day
mean_day_hour = data.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      = "Day and hour",
    ylabel      = "Number of users"
)
axs[3].title.set_size(9)

fig.suptitle("Seasonality plots", fontsize=12)
fig.tight_layout()

There is a clear difference between weekdays and weekends. There is also a clear intra-day pattern, with a different influx of users depending on the time of day.

Autocorrelation plots¶

Auto-correlation plots are a useful tool for identifying the order of an autoregressive model. The autocorrelation function (ACF) is a measure of the correlation between the time series and a lagged version of itself. The partial autocorrelation function (PACF) is a measure of the correlation between the time series and a lagged version of itself, controlling for the values of the time series at all shorter lags. These plots are useful for identifying the lags to be included in the autoregressive model.

# Autocorrelation plot
# ==============================================================================
fig, ax = plt.subplots(figsize=(5, 2))
plot_acf(data['users'], ax=ax, lags=72)
plt.show()

# Partial autocorrelation plot
# ==============================================================================
fig, ax = plt.subplots(figsize=(5, 2))
plot_pacf(data['users'], ax=ax, lags=72, method='ywm')
plt.show()

The results of the autocorrelation study show that there is a significant correlation between the number of users in previous hours, as well as the days before, and the number of users in the future. This means that knowing the number of users during certain periods in the past could be valuable in predicting the number of users in the future.

Baseline¶

When faced with a forecasting problem, it is important to establish a baseline model. This is usually a very simple model that can be used as a reference to assess whether more complex models are worth implementing.

Skforecast allows you to easily create a baseline with its class ForecasterEquivalentDate. This model, also known as Seasonal Naive Forecasting, simply returns the value observed in the same period of the previous season (e.g., the same working day from the previous week, the same hour from the previous day, etc.).

Based on the exploratory analysis performed, the baseline model will be the one that predicts each hour using the value of the same hour on the previous day.

✎ Note

In the following code cells, a baseline forecaster is trained and its predictive ability is evaluated using a backtesting process. If this concept is new to you, do not worry, it will be explained in detail throughout the document. For now, it is sufficient to know that the backtesting process consists of training the model with a certain amount of data and evaluating its predictive ability with the data that the model has not seen. The error metric will be used as a reference to evaluate the predictive ability of the more complex models that will be implemented later.

# Create baseline: value of the same hour of the previous day
# ==============================================================================
forecaster = ForecasterEquivalentDate(
                 offset    = pd.DateOffset(days=1),
                 n_offsets = 1
             )

# Train forecaster
# ==============================================================================
forecaster.fit(y=data.loc[:end_validation, 'users'])
forecaster

======================== 
ForecasterEquivalentDate 
======================== 
Offset: <DateOffset: days=1> 
Number of offsets: 1 
Aggregation function: mean 
Window size: 24 
Training range: [Timestamp('2011-01-01 00:00:00'), Timestamp('2012-08-31 23:00:00')] 
Training index type: DatetimeIndex 
Training index frequency: h 
Creation date: 2024-11-03 20:15:38 
Last fit date: 2024-11-03 20:15:38 
Skforecast version: 0.14.0 
Python version: 3.12.4 
Forecaster id: None

# Backtesting
# ==============================================================================
cv = TimeSeriesFold(
        steps              = 36,
        initial_train_size = len(data.loc[:end_validation]),
        refit              = False
)
metric_baseline, predictions = backtesting_forecaster(
                                   forecaster    = forecaster,
                                   y             = data['users'],
                                   cv            = cv,
                                   metric        = 'mean_absolute_error',
                                   n_jobs        = 'auto',
                                   verbose       = False,
                                   show_progress = True
                               )
metric_baseline

The MAE of the baseline model is used as a reference to assess whether the more complex models are worth implementing.

Recursive multi-step forecasting with LightGBM¶

LightGBM is a highly efficient implementation of the stochastic gradient boosting algorithm, which has become a benchmark in the field of machine learning. The LightGBM library includes its own API as well as the scikit-learn API, making it compatible with skforecast.

First, an ForecasterAutoreg model is trained using past values (lags) and moving average as predictors. Later, exogenous variables are added to the model and the improvement in its performance is assessed. Since Gradient Boosting models have a large number of hyperparameters, a Bayesian Search is performed using the bayesian_search_forecaster() function to find the best combination of hyperparameters and lags. Finally, the predictive ability of the model is evaluated using a backtesting process.

Forecaster¶

# Create forecaster
# ==============================================================================
window_features = RollingFeatures(stats=["mean"], window_sizes=24 * 3)
forecaster = ForecasterRecursive(
                regressor       = LGBMRegressor(random_state=15926, verbose=-1),
                lags            = 24,
                window_features = window_features
             )

# Train forecaster
# ==============================================================================
forecaster.fit(y=data.loc[:end_validation, 'users'])
forecaster

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

2012-09-01 00:00:00    108.331027
2012-09-01 01:00:00     68.562982
2012-09-01 02:00:00     33.499525
2012-09-01 03:00:00     10.027583
2012-09-01 04:00:00      3.037563
2012-09-01 05:00:00     17.162543
2012-09-01 06:00:00     51.059825
2012-09-01 07:00:00    146.940053
2012-09-01 08:00:00    344.320596
2012-09-01 09:00:00    439.738683
Freq: h, Name: pred, dtype: float64

Backtesting¶

In order to have a robust estimate of the predictive ability of the model, a backtesting process is carried out. The backtesting process consists of generating a forecast for each observation in the test set, following the same procedure as it would be done in production, and then comparing the predicted value with the actual value.

# Backtest model on test data
# ==============================================================================
cv = TimeSeriesFold(
        steps              = 36,
        initial_train_size = len(data[:end_validation]),
        refit              = False,
)
metric, predictions = backtesting_forecaster(
                          forecaster    = forecaster,
                          y             = data['users'],
                          cv            = cv,
                          metric        = 'mean_absolute_error',
                          n_jobs        = 'auto',
                          verbose       = False,  # Change to True to see more information
                          show_progress = True
                      )
predictions.head()

# Backtesting error
# ==============================================================================
metric

The prediction error of the autoregressive model achieves a lower MAE than that of the baseline model.

Hyperparameter tuning¶

Hyperparameter tuning is a critical step in developing effective machine learning models. The ForecasterRecursive used in the previous section included the first 24 lags and a LGMBRegressor model with the default hyperparameters. However, there is no reason why these values are the most appropriate. To find the best hyperparameters, a Bayesian Search is performed using the bayesian_search_forecaster() function. The search is carried out using the same backtesting process as before, but each time, the model is trained with different combinations of hyperparameters and lags. It is important to note that the hiperparameter search must be done using the validation set, so the test data is never used.

The search is performed by testing each combination of hyperparameters and lags as follows:

Train the model using only the training set.
The model is evaluated using the validation set via backtesting.
Select the combination of hyperparameters and lags that gives the lowest error.
Train the model again using the best combination found, this time using both the training and validation data.

By following these steps, one can obtain a model with optimized hyperparameters and avoid overfitting.

✎ Note

Searching for hyperparameters can take a long time, especially when using a validation strategy based on backtesting (TimeSeriesFold). A faster alternative is to use a validation strategy based on one-step-ahead predictions (OneStepAheadFold). Although this strategy is faster, it may not be as accurate as validation based on backtesting. For a more detailed description of the pros and cons of each strategy, see the section [backtesting vs one-step-ahead](https://skforecast.org/latest/faq/parameters-search-backetesting-vs-one-step-ahead.html).

# Hyperparameters search
# ==============================================================================
# Lags grid
lags_grid = [48, 72, [1, 2, 3, 23, 24, 25, 167, 168, 169]]

# Regressor hyperparameters search space
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),
        'min_data_in_leaf': trial.suggest_int('min_data_in_leaf', 25, 500),
        'learning_rate'   : trial.suggest_float('learning_rate', 0.01, 0.5),
        'feature_fraction': trial.suggest_float('feature_fraction', 0.5, 1, step=0.1),
        'max_bin'         : trial.suggest_int('max_bin', 50, 250, step=25),
        'reg_alpha'       : trial.suggest_float('reg_alpha', 0, 1, step=0.1),
        'reg_lambda'      : trial.suggest_float('reg_lambda', 0, 1, step=0.1),
        'lags'            : trial.suggest_categorical('lags', lags_grid)
    } 
    return search_space

# Folds training and validation
cv = TimeSeriesFold(
        steps              = 36,
        initial_train_size = len(data_train),
        refit              = False,
    )

results_search, frozen_trial = bayesian_search_forecaster(
    forecaster    = forecaster,
    y             = data.loc[:end_validation, 'users'], # Test data not used
    cv            = cv,
    search_space  = search_space,
    metric        = 'mean_absolute_error',
    n_trials      = 20, # Increase this value for a more exhaustive search
    random_state  = 123,
    return_best   = True,
    n_jobs        = 'auto',
    verbose       = False,
    show_progress = True
)

`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': 3, 'min_data_in_leaf': 236, 'learning_rate': 0.10101592561754799, 'feature_fraction': 0.6, 'max_bin': 225, 'reg_alpha': 0.4, 'reg_lambda': 0.0}
  Backtesting metric: 64.19204873851318

# Search results
# ==============================================================================
results_search.head(3)

Since return_best has been set to True, the forecaster object is automatically updated with the best configuration found and trained on the entire dataset. This final model can then be used for future predictions on new data.

# Best model
# ==============================================================================
forecaster

Backtesting on test data¶

Once the best combination of hyperparameters has been identified using the validation data, the predictive capacity of the model is evaluated when applied to the test set. It is highly recommended to review the documentation for the backtesting_forecaster() function to gain a better understanding of its capabilities. This will help to utilize its full potential to analyze the predictive ability of the model.

# Backtest final model on test data
# ==============================================================================
cv = TimeSeriesFold(
        steps              = 36,
        initial_train_size = len(data[:end_validation]),
        refit              = False,
     )
metric, predictions = backtesting_forecaster(
                          forecaster    = forecaster,
                          y             = data['users'],
                          cv            = cv,
                          metric        = 'mean_absolute_error',
                          n_jobs        = 'auto',
                          verbose       = False,
                          show_progress = True
                      )
display(metric)
predictions.head()

# Plot predictions vs real value
# ======================================================================================
fig = go.Figure()
trace1 = go.Scatter(x=data_test.index, y=data_test['users'], name="test", mode="lines")
trace2 = go.Scatter(x=predictions.index, y=predictions['pred'], name="prediction", mode="lines")
fig.add_trace(trace1)
fig.add_trace(trace2)
fig.update_layout(
    title="Real value vs predicted in test data",
    xaxis_title="Date time",
    yaxis_title="Users",
    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()

After optimizing lags and hyperparameters, the prediction error of the autoregressive model achieves an MAE that is significantly lower than that of the baseline model. Next, the predictive capacity of the model is evaluated when exogenous variables are added.

Exogenous Variables¶

So far, only lagged values of the time series have been used as predictors. However, it is possible to include other variables as predictors. These variables are known as exogenous variables (features) and their use can improve the predictive capacity of the model. A very important point to keep in mind is that the values of the exogenous variables must be known at the time of prediction.

Common examples of exogenous variables are those derived from the calendar, such as the day of the week, month, year, or holidays. Weather variables such as temperature, humidity, and wind also fall into this category, as do economic variables such as inflation and interest rates.

⚠ Warning

Exogenous variables must be known at the time of the forecast. For example, if temperature is used as an exogenous variable, the temperature value for the next hour must be known at the time of the forecast. If the temperature value is not known, the forecast will not be possible.

Weather variables should be used with caution. When the model is deployed into production, future weather conditions are not known, but are predictions made by meteorological services. Because they are predictions, they introduce errors into the forecasting model. As a result, the model's predictions are likely to get worse. One way to anticipate this problem, and to know (not avoid) the expected performance of the model, is to use the weather forecasts available at the time the model is trained, rather than the recorded conditions.

Next, exogenous variables are created based on calendar information, sunrise and sunset times, temperature, and holidays. These new variables are then added to the training, validation and test sets, and used as predictors in the autoregressive model.

Calendar and meteorological features¶

# Calendar features
# ==============================================================================
features_to_extract = [
    'month',
    'week',
    'day_of_week',
    'hour'
]
calendar_transformer = DatetimeFeatures(
    variables='index',
    features_to_extract=features_to_extract,
    drop_original=True,
)
calendar_features = calendar_transformer.fit_transform(data)[features_to_extract]

# Sunlight features
# ==============================================================================
location = LocationInfo(
    name      = 'Washington DC',
    region    = 'USA',
    timezone  = 'US/Eastern',
    latitude  = 40.516666666666666,
    longitude = -77.03333333333333
)
sunrise_hour = [
    sun(location.observer, date=date, tzinfo=location.timezone)['sunrise'].hour
    for date in data.index
]
sunset_hour = [
    sun(location.observer, date=date, tzinfo=location.timezone)['sunset'].hour
    for date in data.index
]
sun_light_features = pd.DataFrame({
                        'sunrise_hour': sunrise_hour,
                        'sunset_hour': sunset_hour}, 
                        index = data.index
                     )
sun_light_features['daylight_hours'] = (
    sun_light_features['sunset_hour'] - sun_light_features['sunrise_hour']
)
sun_light_features["is_daylight"] = np.where(
    (data.index.hour >= sun_light_features["sunrise_hour"])
    & (data.index.hour < sun_light_features["sunset_hour"]),
    1,
    0,
)

# Holiday features
# ==============================================================================
holiday_features = data[['holiday']].astype(int)
holiday_features['holiday_previous_day'] = holiday_features['holiday'].shift(24)
holiday_features['holiday_next_day'] = holiday_features['holiday'].shift(-24)

# Rolling windows of temperature
# ==============================================================================
wf_transformer = WindowFeatures(
    variables   = ["temp"],
    window      = ["1D", "7D"],
    functions   = ["mean", "max", "min"],
    freq        = "h",
)
temp_features = wf_transformer.fit_transform(data[['temp']])


# Merge all exogenous variables
# ==============================================================================
assert all(calendar_features.index == sun_light_features.index)
assert all(calendar_features.index == temp_features.index)
assert all(calendar_features.index == holiday_features.index)
df_exogenous_features = pd.concat([
    calendar_features,
    sun_light_features,
    temp_features,
    holiday_features
], axis=1)

# Due to the creation of moving averages, there are missing values at the beginning
# of the series. And due to holiday_next_day there are missing values at the end.
df_exogenous_features = df_exogenous_features.iloc[7 * 24:, :]
df_exogenous_features = df_exogenous_features.iloc[:-24, :]
df_exogenous_features.head(3)

Features with a cyclical pattern¶

Certain aspects of the calendar, such as hours or days, are cyclical. For example, the hour-day cycle ranges from 0 to 23 hours. These types of features can be handled in several ways, each with its advantages and disadvantages.

One approach is to use the features directly as numeric values without any transformation. This method avoids creating a lot of new features, but may impose an incorrect linear order on the values. For example, hour 23 of one day and hour 00 of the next are very far apart in their linear representation, when in fact there is only one hour difference between them.
Another possibility is to treat cyclical features as categorical variables to avoid imposing a linear order. However, this approach may result in the loss of the cyclical information inherent in the variable.
There is a third way to deal with cyclical features that is often preferred to the other two methods. This involves transforming the features using the sine and cosine of their period. This approach generates only two new features that capture the cyclicality of the data more accurately than the previous two methods because it preserves the natural order of the feature and avoids imposing a linear order.

✎ Note

For a more detailed explanation of the different ways of dealing with cyclical features, see the Cyclical features in time series forecasting document.

# Cliclical encoding of calendar and sunlight features
# ==============================================================================
features_to_encode = [
    "month",
    "week",
    "day_of_week",
    "hour",
    "sunrise_hour",
    "sunset_hour",
]
max_values = {
    "month": 12,
    "week": 52,
    "day_of_week": 6,
    "hour": 23,
    "sunrise_hour": 23,
    "sunset_hour": 23,
}
cyclical_encoder = CyclicalFeatures(
    variables     = features_to_encode,
    max_values    = max_values,
    drop_original = False
)

df_exogenous_features = cyclical_encoder.fit_transform(df_exogenous_features)
df_exogenous_features.head(3)

Feature interaction¶

In many cases, exogenous variables are not isolated. Rather, their effect on the target variable depends on the value of other variables. For example, the effect of temperature on electricity demand depends on the time of day. The interaction between the exogenous variables can be captured by new variables that are obtained by multiplying existing variables together. This interactions are easily obtained with the PolynomialFeatures class from scikit-learn.

# Interaction between exogenous variables
# ==============================================================================
transformer_poly = PolynomialFeatures(
                        degree           = 2,
                        interaction_only = True,
                        include_bias     = False
                    ).set_output(transform="pandas")
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',
    'daylight_hours',
    'is_daylight',
    'holiday_previous_day',
    'holiday_next_day',
    'temp_window_1D_mean',
    'temp_window_1D_min',
    'temp_window_1D_max',
    'temp_window_7D_mean',
    'temp_window_7D_min',
    'temp_window_7D_max',
    'temp',
    'holiday'
]
poly_features = transformer_poly.fit_transform(df_exogenous_features[poly_cols])
poly_features = poly_features.drop(columns=poly_cols)
poly_features.columns = [f"poly_{col}" for col in poly_features.columns]
poly_features.columns = poly_features.columns.str.replace(" ", "__")
assert all(poly_features.index == df_exogenous_features.index)
df_exogenous_features = pd.concat([df_exogenous_features, poly_features], axis=1)
df_exogenous_features.head(3)

Categorical features¶

There are several approaches to incorporating categorical variables into LightGBM (and other gradient boosting frameworks):

One option is to transform the data by converting categorical values to numerical values using methods such as one-hot encoding or ordinal encoding. This approach is applicable to all machine learning models.
Alternatively, LightGBM can handle categorical variables internally without the need for preprocessing. This can be done automatically by setting the parameter categorical_features='auto' and encoding the variables as the category datatype within a Pandas DataFrame. As another option, users can specify the names of the features to be treated as categorical by passing a list of column names to the categorical_features parameter.

There is no one method that is always better than the others. General rules are:

When the cardinality of categorical variables is high (many different values), it is better to use the native support for categorical variables than to use one-hot encoding.
With one-hot encoded data, more split points (i.e., more depth) are needed to recover an equivalent split that could be obtained with a single split point using native handling.
When a categorical variable is converted to multiple dummy variables using one-hot, its importance is diluted, making the analysis of the importance of the features more complex to interpret.

# Store categorical variables as category type
# ==============================================================================
data["weather"] = data["weather"].astype("category")

One hot encoding¶

ColumnTransformers in scikit-learn provide a powerful way to define transformations and apply them to specific features. By encapsulating the transformations in a ColumnTransformer object, it can be passed to a Forecaster using the transformer_exog argument.

✎ Note

It is possible to apply a transformation to the entire dataset independent of the Forecaster. However, it is crucial to ensure that the transformations are learned only from the training data to avoid information leakage. In addition, the same transformation should be applied to the input data during the prediction. Therefore, it is advisable to include the transformation in the Forecaster so that it is handled internally. This approach ensures consistency in the application of transformations and reduces the likelihood of errors.

# One hot encoding transformer
# ==============================================================================
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")

# Create a forecaster with a transformer for exogenous features
# ==============================================================================
forecaster = ForecasterRecursive(
                regressor        = LGBMRegressor(random_state=15926, verbose=-1),
                lags             = 72,
                window_features  = window_features,
                transformer_exog = one_hot_encoder
             )

To examine how the data is transformed, it is possible to use the create_train_X_y method to generate the matrices that the Forecaster uses to train the model. This approach provides insight into the specific data manipulations that occur during the training process.

# View training matrix
# ==============================================================================
exog_features = ['weather']         
X_train, y_train = forecaster.create_train_X_y(
                        y    = data.loc[:end_validation, 'users'],
                        exog = data.loc[:end_validation, exog_features]
                   )
X_train.head(3)

The One Hot Encoder strategy has been shown for didactic purposes. For the rest of the document, however, the native support for categorical variables is used.

Native implementation for categorical features¶

The Gradient Boosting libraries (XGBoost, LightGBM, CatBoost, and HistGradientBoostingRegressor) assume that the input categories are integers starting from 0 up to the number of categories [0, 1, ..., n_categories-1]. In most real cases, categorical variables are not encoded with numbers but with strings, so an intermediate transformation step is necessary. Two options are:

Set columns with categorical variables to the type category. Internally, this data structure consists of an array of categories and an array of integer values (codes) that point to the actual value of the array of categories. That is, internally it is a numeric array with a mapping that relates each value with a category. Models are able to automatically identify the columns of type category and access their internal codes. This approach is applicable to XGBoost, LightGBM and CatBoost.
Preprocess the categorical columns with an OrdinalEncoder to transform their values to integers and then explicitly indicate which features should be treated as categorical.

To use the the automatic detection in skforecast, categorical variables must first be encoded as integers and then stored as the type category again. This is because skforecast internally uses a numeric numpy array to speed up the calculation.

⚠ Warning

The four gradient boosting frameworks – LightGBM, scikit-learn's HistogramGradientBoosting, XGBoost, and CatBoost – are capable of directly handling categorical features within the model. However, it is important to note that each framework has its own configurations, benefits and potential problems. To fully comprehend how to use these frameworks, it is highly recommended to refer to the skforecast user guide for a detailed understanding.

When deploying models in production, it is strongly recommended to avoid using automatic detection based on pandas category type columns. Although pandas provides an internal coding for these columns, it is not consistent across different datasets and may vary depending on the categories present in each dataset. It is therefore crucial to be aware of this issue and to take appropriate measures to ensure consistency in the coding of categorical features when deploying models in production. More details on this issue can be found in Native implementation for categorical features.

# Transformer: Ordinal encoding + cast to category type
# ==============================================================================
pipeline_categorical = make_pipeline(
    OrdinalEncoder(
        dtype=int,
        handle_unknown="use_encoded_value",
        unknown_value=-1,
        encoded_missing_value=-1
    ),
    FunctionTransformer(
        func=lambda x: x.astype('category'),
        feature_names_out= 'one-to-one'
    )
)

transformer_exog = make_column_transformer(
    (
        pipeline_categorical,
        make_column_selector(dtype_include=['category', 'object']),
    ),
    remainder="passthrough",
    verbose_feature_names_out=False,
).set_output(transform="pandas")

When creating a Forecaster with LGBMRegressor, it is necessary to specify how to handle the categorical columns using the fit_kwargs argument. This is because the categorical_feature argument is only specified in the fit method of LGBMRegressor, and not during its initialization.

# Create a forecaster with automatic categorical detection
# ==============================================================================
forecaster = ForecasterRecursive(
                regressor        = LGBMRegressor(random_state=15926, verbose=-1),
                lags             = 72,
                window_features  = window_features,
                transformer_exog = transformer_exog,
                fit_kwargs       = {"categorical_feature": "auto"}
             )

This is the strategy that will be used throughout the rest of this document.

Evaluate model with exogenous features¶

The forecaster is trained again, but this time, the exogenous variables are also included as predictors. For categorical features, the native implementation is used.

# Select exogenous variables to be included in the model
# ==============================================================================
exog_features = []
# Columns that ends with _sin or _cos are selected
exog_features.extend(df_exogenous_features.filter(regex='_sin$|_cos$').columns.tolist())
# columns that start with temp_ are selected
exog_features.extend(df_exogenous_features.filter(regex='^temp_.*').columns.tolist())
# Columns that start with holiday_ are selected
exog_features.extend(df_exogenous_features.filter(regex='^holiday_.*').columns.tolist())
exog_features.extend(['temp', 'holiday', 'weather'])

df_exogenous_features = df_exogenous_features.filter(exog_features, axis=1)

# Merge target and exogenous variables in the same dataframe
# ==============================================================================
data = data[['users', 'weather']].merge(
            df_exogenous_features,
            left_index=True,
            right_index=True,
            how='inner' # To use only dates for which we have all the variables
        )
data = data.astype({col: np.float32 for col in data.select_dtypes("number").columns})
data_train = data.loc[: end_train, :].copy()
data_val   = data.loc[end_train:end_validation, :].copy()
data_test  = data.loc[end_validation:, :].copy()

Because many more features are now included in the model, previously identified hyperparameters may no longer be optimal, so a new hyperparameter search is performed. This time, however, the search is tailored to the range of values identified in the previous search.

# Hyperparameters search
# ==============================================================================
# Create forecaster with a transformer categorical features
forecaster = ForecasterRecursive(
                regressor        = LGBMRegressor(random_state=15926, verbose=-1),
                lags             = 72,
                window_features  = window_features,
                transformer_exog = transformer_exog,
                fit_kwargs       = {"categorical_feature": "auto"}
             )

# Lags grid
lags_grid = [48, 72, [1, 2, 3, 23, 24, 25, 167, 168, 169]]

# Regressor hyperparameters search space
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),
        'min_data_in_leaf': trial.suggest_int('min_data_in_leaf', 25, 500),
        'learning_rate'   : trial.suggest_float('learning_rate', 0.01, 0.5),
        'feature_fraction': trial.suggest_float('feature_fraction', 0.5, 1, step=0.1),
        'max_bin'         : trial.suggest_int('max_bin', 50, 250, step=25),
        'reg_alpha'       : trial.suggest_float('reg_alpha', 0, 1, step=0.1),
        'reg_lambda'      : trial.suggest_float('reg_lambda', 0, 1, step=0.1),
        'lags'            : trial.suggest_categorical('lags', lags_grid)
    } 
    return search_space

# Folds
cv = TimeSeriesFold(
        steps              = 36,
        initial_train_size = len(data_train),
        refit              = False
     )

results_search, frozen_trial = bayesian_search_forecaster(
    forecaster    = forecaster,
    y             = data.loc[:end_validation, 'users'],
    exog          = data.loc[:end_validation, exog_features],
    cv            = cv,
    search_space  = search_space,
    metric        = 'mean_absolute_error',
    n_trials      = 20,
    random_state  = 123,
    return_best   = True,
    n_jobs        = 'auto',
    verbose       = False,
    show_progress = True
)
best_params = results_search['params'].iat[0]
best_params = best_params | {'random_state': 15926, 'verbose': -1}

`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': 8, 'min_data_in_leaf': 207, 'learning_rate': 0.10182210479936714, 'feature_fraction': 1.0, 'max_bin': 250, 'reg_alpha': 0.8, 'reg_lambda': 0.0}
  Backtesting metric: 56.13534870179934

# Backtesting model with exogenous variables on test data
# ==============================================================================
cv = TimeSeriesFold(
        steps              = 36,
        initial_train_size = len(data[:end_validation]),
        refit              = False,
     )

metric, predictions = backtesting_forecaster(
                        forecaster    = forecaster,
                        y             = data['users'],
                        exog          = data[exog_features],
                        cv            = cv,
                        metric        = 'mean_absolute_error',
                        n_jobs        = 'auto',
                        verbose       = False,
                        show_progress = True
                     )
metric

# Plot predictions vs real value
# ======================================================================================
fig = go.Figure()
trace1 = go.Scatter(x=data_test.index, y=data_test['users'], name="test", mode="lines")
trace2 = go.Scatter(x=predictions.index, y=predictions['pred'], name="prediction", mode="lines")
fig.add_trace(trace1)
fig.add_trace(trace2)
fig.update_layout(
    title="Real value vs predicted in test data",
    xaxis_title="Date time",
    yaxis_title="Users",
    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()

The inclusion of exogenous variables as predictors further improves the predictive capacity of the model.

Feature selection¶

Feature selection is the process of selecting a subset of relevant features for use in model construction. It is an important step in the machine learning process, as it can help to reduce overfitting, improve model accuracy, and reduce training time. Since the underlying regressors of skforecast follow the scikit-learn API, it is possible to apply the feature selection methods available in scikit-learn with the function select_features(). Two of the most popular methods are Recursive Feature Elimination and Sequential Feature Selection.

💡 Tip

Feature selection is a powerful tool for improving the performance of machine learning models. However, it is computationally expensive and can be time-consuming. Since the goal is to find the best subset of features, not the best model, it is not necessary to use the entire data set or a highly complex model. Instead, it is recommended to use a small subset of the data and a simple model. Once the best subset of features has been identified, the model can then be trained using the entire dataset and a more complex configuration.

# Create forecaster
# ==============================================================================
regressor = LGBMRegressor(
    n_estimators = 100,
    max_depth    = 5,
    random_state = 15926,
    verbose      = -1
)
forecaster = ForecasterRecursive(
    regressor        = regressor,
    lags             = [1, 2, 3, 23, 24, 25, 167, 168, 169],
    window_features  = window_features,
    transformer_exog = transformer_exog,
    fit_kwargs       = {"categorical_feature": "auto"}
)

# Recursive feature elimination with cross-validation
# ==============================================================================
selector = RFECV(
    estimator = regressor,
    step      = 1,
    cv        = 3,
    n_jobs    = -1
)
selected_lags, selected_window_features, selected_exog = select_features(
    forecaster      = forecaster,
    selector        = selector,
    y               = data_train['users'],  
    exog            = data_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: 11327
Total number of records used for feature selection: 5663
Number of features available: 99
    Lags            (n=9)
    Window features (n=1)
    Exog            (n=89)
Number of features selected: 22
    Lags            (n=9) : [1, 2, 3, 23, 24, 25, 167, 168, 169]
    Window features (n=1) : ['roll_mean_72']
    Exog            (n=12) : ['weather', 'hour_sin', 'hour_cos', 'poly_month_sin__hour_sin', 'poly_week_sin__hour_sin', 'poly_day_of_week_sin__hour_sin', 'poly_day_of_week_sin__hour_cos', 'poly_hour_sin__hour_cos', 'poly_hour_sin__sunset_hour_cos', 'poly_hour_cos__sunset_hour_sin', 'temp_window_1D_mean', 'temp']

Scikit-learn's RFECV starts by training a model on the initial set of features, and obtaining the importance of each feature (through attributes such as coef_ or feature_importances_). Then, in each round, the least important features are iteratively removed, followed by cross-validation to calculate the performance of the model with the remaining features. This process continues until further feature removal doesn't improve or starts to degrade the performance of the model (based on a chosen metric), or the min_features_to_select is reached.

The final result is an optimal subset of features that ideally balances model simplicity and predictive power, as determined by the cross-validation process.

The forecaster is trained and re-evaluated using the best subset of features.

# Create a forecaster with the selected features
# ==============================================================================
forecaster = ForecasterRecursive(
    regressor        = LGBMRegressor(**best_params),
    lags             = selected_lags,
    window_features  = window_features,
    transformer_exog = transformer_exog,
    fit_kwargs       = {"categorical_feature": "auto"}
)

# Backtesting model with exogenous variables on test data
# ==============================================================================
cv = TimeSeriesFold(
        steps              = 36,
        initial_train_size = len(data[:end_validation]),
        refit              = False,
)
metric_lgbm, predictions = backtesting_forecaster(
    forecaster    = forecaster,
    y             = data['users'],
    exog          = data[selected_exog],
    cv            = cv,
    metric        = 'mean_absolute_error',
    n_jobs        = 'auto',
    verbose       = False,
    show_progress = True
)
metric_lgbm

The performance of the model remains similar to that of the model trained with all the features. However, the model is now much simpler, which will make it faster to train and less prone to overfitting. For the rest of the document, the model will be trained using only the most important features.

# Update exogenous variables
# ==============================================================================
exog_features = selected_exog

Probabilistic forecasting: Prediction intervals¶

A prediction interval defines the interval within which the true value of the target variable can be expected to be found with a given probability. Rob J Hyndman and George Athanasopoulos, in their book Forecasting: Principles and Practice, list multiple ways to estimate prediction intervals, most of which require that the residuals (errors) of the model to be normally distributed. If this cannot be assumed, one can resort to bootstrapping, which requires only that the residuals be uncorrelated. This is one of the methods available in skforecast. A more detailed explanation of prediction intervals can be found in the Probabilistic forecasting: prediction intervals and prediction distribution user guide.

The following code shows how to generate prediction intervals for the autoregressive model. The prediction_interval() function is used to generate the prediction intervals for each predicted step. Then, the backtesting_forecaster() function is used to generate the prediction intervals for the entire test set. The interval argument is used to specify the desired coverage probability of the prediction intervals. In this case, interval is set to [10, 90], which means that the prediction intervals are calculated for the 10th and 90th percentiles, resulting in a theoretical coverage probability of 80%. The n_boot argument is used to specify the number of bootstrap samples to be used in estimating the prediction intervals. The larger the number of samples, the more accurate the prediction intervals will be, but the longer the calculation will take.

# Create and train forecaster
# ==============================================================================
forecaster = ForecasterRecursive(
    regressor        = LGBMRegressor(**best_params),
    lags             = 72,
    window_features  = window_features,
    transformer_exog = transformer_exog,
    fit_kwargs       = {"categorical_feature": "auto"},
    binner_kwargs    = {"n_bins": 20}
)
forecaster.fit(
    y    = data.loc[:end_train, 'users'],
    exog = data.loc[:end_train, exog_features]
)

# Predict intervals
# ==============================================================================
# Since the model has been trained with exogenous variables, they must be provided
# for the prediction.
predictions = forecaster.predict_interval(
    exog     = data.loc[end_train:, exog_features],
    steps    = 24,
    interval = [5, 95], 
)
predictions.head()

By default, intervals are calculated using in-sample residuals (residuals from the training set). However, this can result in intervals that are too narrow (overly optimistic). To avoid this, the set_out_sample_residuals() method is used to specify out-sample residuals computed with a validation set through backtesting.

If the predicted values are passed to set_out_sample_residuals() in addition to the residuals, then the residuals used in the bootstrapping process can be conditioned on the range of values of the predictions. This can help to improve the coverage of the estimated intervals while keeping them as narrow as possible.

# Backtesting on validation data to obtain out-sample residuals
# ==============================================================================
cv = TimeSeriesFold(
        steps              = 24,
        initial_train_size = len(data.loc[:end_train]),
        refit              = False,
)
_, predictions_val = backtesting_forecaster(
    forecaster    = forecaster,
    y             = data.loc[:end_validation, 'users'],
    exog          = data.loc[:end_validation, exog_features],
    cv            = cv,
    metric        = 'mean_absolute_error',
    n_jobs        = 'auto',
    verbose       = False,
    show_progress = True
)

# Out-sample residuals distribution
# ==============================================================================
residuals = data.loc[predictions_val.index, 'users'] - predictions_val['pred']
print(pd.Series(np.where(residuals < 0, 'negative', 'positive')).value_counts())
plt.rcParams.update({'font.size': 8})
_ = plot_residuals(
        y_true = data.loc[predictions_val.index, 'users'],
        y_pred = predictions_val['pred'],
        figsize=(7, 4)
    )

negative    1589
positive    1363
Name: count, dtype: int64

# Store out-sample residuals in the forecaster
# ==============================================================================
forecaster.set_out_sample_residuals(
    y_true = data.loc[predictions_val.index, 'users'],
    y_pred = predictions_val['pred']
)

The backtesting process is then run to estimate the prediction intervals in the test set. The argument use_in_sample_residuals is set to False so that the previously stored out-sample residuals are used and use_binned_residuals so that the residuals used in bootstrapping are selected conditional on the range of prediction values.

⚠ Warning

Estimating prediction intervals using the bootstraping approach is computationally intensive. To reduce the time required for the following example, the test set is reduced.

# Backtesting with prediction intervals in test data using out-sample residuals
# ==============================================================================
end_test = '2012-10-31 23:00:00'
cv = TimeSeriesFold(
        steps              = 24,
        initial_train_size = len(data.loc[:end_validation]),
        refit              = False,
    )
metric, predictions = backtesting_forecaster(
   forecaster              = forecaster,
   y                       = data.loc[:end_test, 'users'],
   exog                    = data.loc[:end_test, exog_features],
   cv                      = cv,
   metric                  = 'mean_absolute_error',
   interval                = [5, 95], # 90% prediction interval
   n_boot                  = 150,
   use_in_sample_residuals = False,  # Use out-sample residuals
   use_binned_residuals    = True,   # Use residuals conditioned on predicted values
   n_jobs                  = 'auto',
   verbose                 = False,
   show_progress           = True
)
predictions.head(5)

# Plot prediction intervals vs real value
# ==============================================================================
fig = go.Figure([
    go.Scatter(name='Prediction', x=predictions.index, y=predictions['pred'], mode='lines'),
    go.Scatter(
        name='Real value', x=data_test.loc[:end_test].index, y=data_test.loc[:end_test, 'users'], mode='lines',
    ),
    go.Scatter(
        name='Upper Bound', x=predictions.index, y=predictions['upper_bound'], mode='lines',
        marker=dict(color="#444"), line=dict(width=0), showlegend=False
    ),
    go.Scatter(
        name='Lower Bound', x=predictions.index, y=predictions['lower_bound'], marker=dict(color="#444"),
        line=dict(width=0), mode='lines', fillcolor='rgba(68, 68, 68, 0.3)', fill='tonexty', showlegend=False
    )
])
fig.update_layout(
    title="Real value vs predicted in test data",
    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),
    # Initial zoom on x axis betwee 1 oct to 10 oct
    xaxis=dict(range=['2012-10-01', '2012-10-10'])
)
fig.show()

# Predicted interval coverage (on test data)
# ==============================================================================
coverage = np.mean(
    np.logical_and(
        data.loc[end_validation:end_test, 'users'] >= predictions["lower_bound"],
        data.loc[end_validation:end_test, 'users'] <= predictions["upper_bound"]
    )
)
area = (predictions['upper_bound'] - predictions['lower_bound']).sum()
print(f"Total area of the interval: {round(area, 2)}")
print(f"Predicted interval coverage: {round(100 * coverage, 2)} %")

Total area of the interval: 579968.65
Predicted interval coverage: 94.4 %

The observed coverage of the intervals is higher than the theoretically expected coverage (90%), the intervals are more conservative than necessary. This means that they contain the true value with a higher probability than expected.

✎ Note

For a more detailed explanation of the probabilistic forecasting features available in skforecast visit: Probabilistic forecasting.

Model explanaibility¶

Due to the complex nature of many modern machine learning models, such as ensemble methods, they often function as black boxes, making it difficult to understand why a particular prediction was made. Explanability techniques aim to demystify these models, providing insight into their inner workings and helping to build trust, improve transparency, and meet regulatory requirements in various domains. Enhancing model explainability not only helps to understand model behavior, but also helps to identify biases, improve model performance, and enable stakeholders to make more informed decisions based on machine learning insights.

Skforecast is compatible with some of the most popular model explainability methods: model-specific feature importances, SHAP values, and partial dependence plots.

# Create and train forecaster
# ==============================================================================
forecaster = ForecasterRecursive(
    regressor        = LGBMRegressor(**best_params),
    lags             = [1, 2, 3, 23, 24, 25, 167, 168, 169],
    window_features  = window_features,
    transformer_exog = transformer_exog,
    fit_kwargs       = {"categorical_feature": "auto"}
)
forecaster.fit(
    y    = data.loc[:end_validation, 'users'],
    exog = data.loc[:end_validation, exog_features]
)

Model-specific feature importance¶

# Extract feature importance
# ==============================================================================
importance = forecaster.get_feature_importances()
importance.head(10)

⚠ Warning

The get_feature_importances() method will only return values if the forecaster's regressor has either the coef_ or feature_importances_ attribute, which is the default in scikit-learn.

Shap values¶

SHAP (SHapley Additive exPlanations) values are a popular method for explaining machine learning models, as they help to understand how variables and values influence predictions visually and quantitatively.

It is possible to generate SHAP-values explanations from skforecast models with just two essential elements:

The internal regressor of the forecaster.
The training matrices created from the time series and used to fit the forecaster.

By leveraging these two components, users can create insightful and interpretable explanations for their skforecast models. These explanations can be used to verify the reliability of the model, identify the most significant factors that contribute to model predictions, and gain a deeper understanding of the underlying relationship between the input variables and the target variable.

# Training matrices used by the forecaster to fit the internal regressor
# ==============================================================================
X_train, y_train = forecaster.create_train_X_y(
                       y    = data.loc[:end_validation, 'users'],
                       exog = data.loc[:end_validation, exog_features]
                   )
display(X_train.head(3))
display(y_train.head(3))

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

# Create SHAP explainer (for three base models)
# ==============================================================================
explainer = shap.TreeExplainer(forecaster.regressor)

# Sample 50% of the data to speed up the calculation
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

Shap library has several explainers, each designed for a different type of model. The shap.TreeExplainer explainer is used for tree-based models, such as the LGBMRegressor used in this example. For more information, see the SHAP documentation.

# 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(10, 4.5)

# Force plot for the first observation
# ==============================================================================
shap.force_plot(explainer.expected_value, shap_values[0, :], X_train_sample.iloc[0, :])

XGBoost, CatBoost and HistGradientBoostingRegressor¶

Since the success of Gradient Boosting as a machine learning algorithm, several implementations have been developed. In addition to LightGBM, three other implementations are very popular: XGBoost, CatBoost and HistGradientBoostingRegressor. All of them are compatible with skforecast.

XGBoost: developed by Tianqi Chen.
HistGradientBoostingRegressor: developed by scikit-learn.
CatBoost: developed by Yandex.

The following sections show how to use these implementations to build forecasting models, with an emphasis on using their native support for categorical features. This time a one-step-ahead validation strategy is used to speed up the search for hyperparameters.

# Folds used for the hyperparameter search and backtesting
# ==============================================================================
cv_search = OneStepAheadFold(initial_train_size = len(data_train))

cv_backtesting = TimeSeriesFold(
                    steps              = 36,
                    initial_train_size = len(data[:end_validation]),
                    refit              = False,
                 )

XGBoost¶

# Transformer: Ordinal encoding + cast to category type
# ==============================================================================
pipeline_categorical = make_pipeline(
    OrdinalEncoder(
        dtype=int,
        handle_unknown="use_encoded_value",
        unknown_value=-1,
        encoded_missing_value=-1
    ),
    FunctionTransformer(
        func=lambda x: x.astype('category'),
        feature_names_out= 'one-to-one'
    )
)

transformer_exog = make_column_transformer(
    (
        pipeline_categorical,
        make_column_selector(dtype_exclude=np.number)
    ),
    remainder="passthrough",
    verbose_feature_names_out=False,
).set_output(transform="pandas")

# Create forecaster
# ==============================================================================
forecaster = ForecasterRecursive(
    regressor = XGBRegressor(tree_method='hist', enable_categorical=True, random_state=123),
    lags = 24,
    window_features  = window_features,
    transformer_exog = transformer_exog
)

# Hyperparameters search
# ==============================================================================
# Lags grid
lags_grid = [48, 72, [1, 2, 3, 23, 24, 25, 167, 168, 169]]

# Regressor hyperparameters search space
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

results_search, frozen_trial = bayesian_search_forecaster(
    forecaster    = forecaster,
    y             = data.loc[:end_validation, 'users'],
    exog          = data.loc[:end_validation, exog_features],
    cv            = cv_search,
    search_space  = search_space,
    metric        = 'mean_absolute_error',
    n_trials      = 20,
    random_state  = 123,
    return_best   = True,
    n_jobs        = 'auto',
    verbose       = False,
    show_progress = True
)


One-step-ahead predictions are used for faster model comparison, but they may not fully represent multi-step prediction performance. It is recommended to backtest the final model for a more accurate multi-step performance estimate.

`Forecaster` refitted using the best-found lags and parameters, and the whole data set: 
  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
 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48] 
  Parameters: {'n_estimators': 900, 'max_depth': 3, 'learning_rate': 0.02321068166899635, 'subsample': 0.6809852212786971, 'colsample_bytree': 0.7638542941294673, 'gamma': 0.3217712671270035, 'reg_alpha': 0.9683999802520469, 'reg_lambda': 0.35358771689006163}
  One-step-ahead metric: 36.8952522277832

# Backtesting model with exogenous variables on test data
# ==============================================================================
metric_xgboost, predictions = backtesting_forecaster(
    forecaster    = forecaster,
    y             = data['users'],
    exog          = data[exog_features],
    cv            = cv_backtesting,
    metric        = 'mean_absolute_error',
    n_jobs        = 'auto',
    verbose       = False,
)
metric_xgboost

HistGradientBoostingRegressor¶

When creating a Forecaster using HistogramGradientBoosting, the names of the categorical columns should be specified during the regressor instantiation by passing them as a list to the categorical_feature argument.

# Transformer: ordinal encoding
# ==============================================================================
# A ColumnTransformer is used to transform categorical features (no numerical)
# using ordinal encoding. Numerical features are left untouched. Missing values
# are encoded as -1. If a new category is found in the test set, it is encoded
# as -1.
categorical_features = ['weather']
transformer_exog = make_column_transformer(
    (
        OrdinalEncoder(
            dtype=int,
            handle_unknown="use_encoded_value",
            unknown_value=-1,
            encoded_missing_value=-1
        ),
        categorical_features
    ),
    remainder="passthrough",
    verbose_feature_names_out=False,
).set_output(transform="pandas")

# Create forecaster
# ==============================================================================
forecaster = ForecasterRecursive(
    regressor = HistGradientBoostingRegressor(
                    categorical_features=categorical_features,
                    random_state=123
                ),
    lags = 24,
    window_features = window_features,
    transformer_exog = transformer_exog
)

# Hyperparameters search
# ==============================================================================
# Lags grid
lags_grid = [48, 72, [1, 2, 3, 23, 24, 25, 167, 168, 169]]

# Regressor hyperparameters search space
def search_space(trial):
    search_space  = {
        'max_iter'          : trial.suggest_int('max_iter', 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),
        'min_samples_leaf'  : trial.suggest_int('min_samples_leaf', 1, 20, step=1),
        'l2_regularization' : trial.suggest_float('l2_regularization', 0, 1),
        'lags'              : trial.suggest_categorical('lags', lags_grid)
    } 
    return search_space

results_search, frozen_trial = bayesian_search_forecaster(
    forecaster    = forecaster,
    y             = data.loc[:end_validation, 'users'],
    exog          = data.loc[:end_validation, exog_features],
    cv            = cv_search,
    search_space  = search_space,
    metric        = 'mean_absolute_error',
    n_trials      = 20,
    random_state  = 123,
    return_best   = True,
    n_jobs        = 'auto',
    verbose       = False,
    show_progress = True
)


One-step-ahead predictions are used for faster model comparison, but they may not fully represent multi-step prediction performance. It is recommended to backtest the final model for a more accurate multi-step performance estimate.

`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: {'max_iter': 1200, 'max_depth': 9, 'learning_rate': 0.038248393152946335, 'min_samples_leaf': 18, 'l2_regularization': 0.2878398874455151}
  One-step-ahead metric: 38.43500088649492

# Backtesting model with exogenous variables on test data
# ==============================================================================
metric_histgb, predictions = backtesting_forecaster(
    forecaster    = forecaster,
    y             = data['users'],
    exog          = data[exog_features],
    cv            = cv_backtesting,
    metric        = 'mean_absolute_error',
    n_jobs        = 'auto',
    verbose       = False,
    show_progress = True
)
metric_histgb

CatBoost¶

Unfortunately, the current version of skforecast is not compatible with CatBoost's built-in handling of categorical features. The issue arises because CatBoost only accepts categorical features as integers, while skforecast converts input data to floats for faster computation using numpy arrays in the internal prediction process. To work around this limitation, it is necessary to apply either the one-hot encoding or label encoding strategy to the categorical features before using them with CatBoost.

# One hot encoding
# ==============================================================================
one_hot_encoder = make_column_transformer(
    (
        OneHotEncoder(sparse_output=False, drop='if_binary'),
        make_column_selector(dtype_exclude=np.number),
    ),
    remainder="passthrough",
    verbose_feature_names_out=False,
).set_output(transform="pandas")

# Create forecaster
# ==============================================================================
forecaster = ForecasterRecursive(
    regressor = CatBoostRegressor(
                    random_state=123, 
                    silent=True, 
                    allow_writing_files=False,
                    boosting_type = 'Plain',         # Faster training
                    leaf_estimation_iterations = 3,  # Faster training
                ),
    lags = 24,
    window_features = window_features,
    transformer_exog = one_hot_encoder
)

# Hyperparameters search
# ==============================================================================
# Lags grid
lags_grid = [48, 72, [1, 2, 3, 23, 24, 25, 167, 168, 169]]

# Regressor hyperparameters search space
def search_space(trial):
    search_space  = {
        'n_estimators'  : trial.suggest_int('n_estimators', 100, 1000, step=100),
        'max_depth'     : trial.suggest_int('max_depth', 3, 10, step=1),
        'learning_rate' : trial.suggest_float('learning_rate', 0.01, 1),
        'lags'          : trial.suggest_categorical('lags', lags_grid)
    } 
    return search_space

results_search, frozen_trial = bayesian_search_forecaster(
    forecaster    = forecaster,
    y             = data.loc[:end_validation, 'users'],
    exog          = data.loc[:end_validation, exog_features],
    cv            = cv_search, 
    search_space  = search_space,
    metric        = 'mean_absolute_error',
    n_trials      = 20,
    random_state  = 123,
    return_best   = True,
    n_jobs        = 'auto',
    verbose       = False,
    show_progress = True
)


One-step-ahead predictions are used for faster model comparison, but they may not fully represent multi-step prediction performance. It is recommended to backtest the final model for a more accurate multi-step performance estimate.

`Forecaster` refitted using the best-found lags and parameters, and the whole data set: 
  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
 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48] 
  Parameters: {'n_estimators': 100, 'max_depth': 3, 'learning_rate': 0.703242177562586}
  One-step-ahead metric: 38.182225356181206

# Backtesting on test data
# ==============================================================================
metric_catboost, predictions = backtesting_forecaster(
    forecaster = forecaster,
    y          = data['users'],
    exog       = data[exog_features],
    cv         = cv_backtesting,
    metric     = 'mean_absolute_error',
    n_jobs     = 'auto',
    verbose    = False
)
metric_catboost

Conclusion¶

The use of Gradient Boosting models in forecasting problems is very easy thanks to the functionalities offered by skforecast.
As shown in this document, the inclusion of exogenous variables as predictors can significantly improve forecasting performance.
Categorical features can be included as exogenous variables without the need for preprocessing. This is possible thanks to the native support for categorical features offered by LightGBM, XGBoost and HistGradientBoostingRegressor.

metrics = pd.concat(
    [metric_baseline, metric_lgbm, metric_xgboost, metric_histgb, metric_catboost]
)
metrics.index = [
    "Baseline",
    "LGBMRegressor",
    "XGBRegressor",
    "HistGradientBoostingRegressor",
    "CatBoostRegressor",
]
metrics.round(2).sort_values(by="mean_absolute_error")

✎ Note

To illustrate the process, the hyperparameter search has been kept small. However, to achieve optimal results, it may be necessary to perform a more extensive search for each model.

Session information¶

import session_info
session_info.show(html=False)

-----
astral              3.2
catboost            1.2.7
feature_engine      1.8.0
lightgbm            4.4.0
matplotlib          3.9.0
numpy               1.26.4
optuna              4.0.0
pandas              2.2.3
plotly              5.23.0
session_info        1.0.0
shap                0.46.0
skforecast          0.14.0
sklearn             1.5.2
statsmodels         0.14.2
xgboost             2.1.2
-----
IPython             8.25.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:03:56) [MSC v.1929 64 bit (AMD64)]
Windows-11-10.0.26100-SP0
-----
Session information updated at 2024-11-03 20:31

Bibliography¶

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

Time Series Analysis and Forecasting with ADAM Ivan Svetunkov.

Joseph, M. (2022). Modern time series forecasting with Python: Explore industry-ready time series forecasting using modern machine learning and Deep Learning. Packt Publishing.

Citation¶

How to cite this document

If you use this document or any part of it, please acknowledge the source, thank you!

Forecasting time series with gradient boosting: Skforecast, XGBoost, LightGBM and CatBoost by Joaquín Amat Rodrigo and Javier Escobar Ortiz, available under a Attribution-NonCommercial-ShareAlike 4.0 International at https://www.cienciadedatos.net/documentos/py39-forecasting-time-series-with-skforecast-xgboost-lightgbm-catboost.html

How to cite skforecast

If you use skforecast for a scientific publication, we would appreciate it if you cite the published software.

Zenodo:

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

APA:

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

BibTeX:

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

Did you like the article? Your support is important

Website maintenance has high cost, your contribution will help me to continue generating free educational content. Many thanks! 😊

This work by Joaquín Amat Rodrigo and Javier Escobar Ortiz is licensed under a Attribution-NonCommercial-ShareAlike 4.0 International.

Allowed:

Share: copy and redistribute the material in any medium or format.
Adapt: remix, transform, and build upon the material.

Under the following terms:

Attribution: You must give appropriate credit, provide a link to the license, and indicate if changes were made. You may do so in any reasonable manner, but not in any way that suggests the licensor endorses you or your use.
NonCommercial: You may not use the material for commercial purposes.
ShareAlike: If you remix, transform, or build upon the material, you must distribute your contributions under the same license as the original.

	month	week	day_of_week	hour	sunrise_hour	sunset_hour	daylight_hours	is_daylight	temp	temp_window_1D_mean	...	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
date_time
2011-01-08 00:00:00	1	1	5	0	7	16	9	0	7.38	8.063333	...	0.120537	0.992709	-0.866025	0.5	0.000000	1.000000	0.942261	-0.33488	-0.942261	-0.33488
2011-01-08 01:00:00	1	1	5	1	7	16	9	0	7.38	8.029167	...	0.120537	0.992709	-0.866025	0.5	0.269797	0.962917	0.942261	-0.33488	-0.942261	-0.33488
2011-01-08 02:00:00	1	1	5	2	7	16	9	0	7.38	7.995000	...	0.120537	0.992709	-0.866025	0.5	0.519584	0.854419	0.942261	-0.33488	-0.942261	-0.33488

	month	week	day_of_week	hour	sunrise_hour	sunset_hour	daylight_hours	is_daylight	temp	temp_window_1D_mean	...	poly_temp_window_7D_mean__temp_window_7D_min	poly_temp_window_7D_mean__temp_window_7D_max	poly_temp_window_7D_mean__temp	poly_temp_window_7D_mean__holiday	poly_temp_window_7D_min__temp_window_7D_max	poly_temp_window_7D_min__temp	poly_temp_window_7D_min__holiday	poly_temp_window_7D_max__temp	poly_temp_window_7D_max__holiday	poly_temp__holiday
date_time
2011-01-08 00:00:00	1	1	5	0	7	16	9	0	7.38	8.063333	...	49.829643	191.013631	74.744464	0.0	92.7912	36.3096	0.0	139.1868	0.0	0.0
2011-01-08 01:00:00	1	1	5	1	7	16	9	0	7.38	8.029167	...	49.757600	190.737467	74.636400	0.0	92.7912	36.3096	0.0	139.1868	0.0	0.0
2011-01-08 02:00:00	1	1	5	2	7	16	9	0	7.38	7.995000	...	49.709571	190.553357	74.564357	0.0	92.7912	36.3096	0.0	139.1868	0.0	0.0

	lag_1	lag_2	lag_3	lag_4	lag_5	lag_6	lag_7	lag_8	lag_9	lag_10	...	lag_67	lag_68	lag_69	lag_70	lag_71	lag_72	roll_mean_72	weather_clear	weather_mist	weather_rain
date_time
2011-01-04 00:00:00	12.0	20.0	52.0	52.0	110.0	157.0	157.0	76.0	72.0	77.0	...	1.0	1.0	13.0	32.0	40.0	16.0	43.638889	1.0	0.0	0.0
2011-01-04 01:00:00	5.0	12.0	20.0	52.0	52.0	110.0	157.0	157.0	76.0	72.0	...	2.0	1.0	1.0	13.0	32.0	40.0	43.486111	1.0	0.0	0.0
2011-01-04 02:00:00	2.0	5.0	12.0	20.0	52.0	52.0	110.0	157.0	157.0	76.0	...	3.0	2.0	1.0	1.0	13.0	32.0	42.958333	1.0	0.0	0.0

	pred	lower_bound	upper_bound
2012-05-01 00:00:00	23.011572	-7.434268	50.121755
2012-05-01 01:00:00	10.131827	-17.232515	49.026826
2012-05-01 02:00:00	1.680696	-24.857192	36.325549
2012-05-01 03:00:00	-5.929731	-27.910258	34.723111
2012-05-01 04:00:00	6.919658	-16.680601	47.093409

	pred	lower_bound	upper_bound
2012-09-01 00:00:00	137.327934	28.001861	244.224522
2012-09-01 01:00:00	82.382643	22.327932	282.897515
2012-09-01 02:00:00	47.575894	8.702706	300.032463
2012-09-01 03:00:00	16.475415	1.863276	277.280476
2012-09-01 04:00:00	4.755445	2.936639	340.840768

	users	holiday	weather	temp	atemp	hum	windspeed
date_time
2011-01-01 00:00:00	16.0	0.0	clear	9.84	14.395	81.0	0.0
2011-01-01 01:00:00	40.0	0.0	clear	9.02	13.635	80.0	0.0
2011-01-01 02:00:00	32.0	0.0	clear	9.02	13.635	80.0	0.0
2011-01-01 03:00:00	13.0	0.0	clear	9.84	14.395	75.0	0.0
2011-01-01 04:00:00	1.0	0.0	clear	9.84	14.395	75.0	0.0

	pred
2012-09-01 00:00:00	108.331027
2012-09-01 01:00:00	68.562982
2012-09-01 02:00:00	33.499525
2012-09-01 03:00:00	10.027583
2012-09-01 04:00:00	3.037563

	lags	params	mean_absolute_error	n_estimators	max_depth	min_data_in_leaf	learning_rate	feature_fraction	max_bin	reg_alpha
0	[1, 2, 3, 23, 24, 25, 167, 168, 169]	{'n_estimators': 400, 'max_depth': 3, 'min_dat...	64.192049	400.0	3.0	236.0	0.101016	0.6	225.0	0.4
1	[1, 2, 3, 23, 24, 25, 167, 168, 169]	{'n_estimators': 400, 'max_depth': 3, 'min_dat...	64.498000	400.0	3.0	222.0	0.083926	1.0	250.0	0.4
2	[1, 2, 3, 23, 24, 25, 167, 168, 169]	{'n_estimators': 400, 'max_depth': 3, 'min_dat...	65.844639	400.0	3.0	261.0	0.100016	0.6	200.0	0.8

	pred
2012-09-01 00:00:00	115.505608
2012-09-01 01:00:00	91.534104
2012-09-01 02:00:00	60.349164
2012-09-01 03:00:00	29.668854
2012-09-01 04:00:00	10.544213

	feature	importance
0	lag_1	966
7	lag_168	727
4	lag_24	656
8	lag_169	452
1	lag_2	375
6	lag_167	375
5	lag_25	361
2	lag_3	358
3	lag_23	357
9	roll_mean_72	311

	lag_1	lag_2	lag_3	lag_23	lag_24	lag_25	lag_167	lag_168	lag_169	roll_mean_72	...	hour_cos	poly_month_sin__hour_sin	poly_week_sin__hour_sin	poly_day_of_week_sin__hour_sin	poly_day_of_week_sin__hour_cos	poly_hour_sin__hour_cos	poly_hour_sin__sunset_hour_cos	poly_hour_cos__sunset_hour_sin	temp_window_1D_mean	temp
date_time
2011-01-15 01:00:00	28.0	27.0	36.0	1.0	5.0	14.0	16.0	16.0	25.0	55.736111	...	0.962917	0.134898	0.064567	-0.233651	-0.833911	0.259792	-0.018412	-0.960672	6.594167	6.56
2011-01-15 02:00:00	20.0	28.0	27.0	1.0	1.0	5.0	7.0	16.0	16.0	55.930556	...	0.854419	0.259792	0.124345	-0.449973	-0.739949	0.443943	-0.035458	-0.852428	6.696667	6.56
2011-01-15 03:00:00	12.0	20.0	28.0	1.0	1.0	1.0	1.0	7.0	16.0	56.083333	...	0.682553	0.365418	0.174900	-0.632923	-0.591108	0.498834	-0.049874	-0.680962	6.799167	6.56

	mean_absolute_error
LGBMRegressor	49.22
HistGradientBoostingRegressor	49.81
XGBRegressor	61.65
CatBoostRegressor	62.41
Baseline	91.67