More about forecasting in cienciadedatos.net

• ARIMA and SARIMAX models with python
• Time series forecasting with machine learning
• Forecasting time series with gradient boosting: XGBoost, LightGBM and CatBoost
• Forecasting time series with XGBoost
• Probabilistic forecasting
• Forecasting with deep learning
• Forecasting energy demand with machine learning
• Forecasting web traffic with machine learning
• Intermittent demand forecasting
• Modelling time series trend with tree-based models
• Bitcoin price prediction with Python
• Stacking ensemble of machine learning models to improve forecasting
• Interpretable forecasting models
• Mitigating the Impact of Covid on forecasting Models
• Forecasting time series with missing values

Introduction¶

When faced with scenarios involving the prediction of hundreds or thousands of time series, a crucial decision arises: should one develop individual models for each series, or should one use a unified model to handle them all at once?

In Single-Series Modeling (Local Forecasting Model), a separate predictive model is created for each time series. While this method provides a comprehensive understanding of each series, its scalability can be challenged by the need to create and maintain hundreds or thousands of models.

Multi-Series Modeling (Global Forecasting Model) involves building a single predictive model that considers all time series simultaneously. It attempts to capture the core patterns that govern the series, thereby mitigating the potential noise that each series might introduce. This approach is computationally efficient, easy to maintain, and can yield more robust generalizations across time series, albeit potentially at the cost of sacrificing some individual insights.

This document shows how to forecast more than 1000 time series with a single model including exogenous features with some of them having different values per series.

Libraries¶

Libraries used in this document.

# Data management
# ==============================================================================
import numpy as np
import pandas as pd

# Plots
# ==============================================================================
import matplotlib.pyplot as plt
plt.style.use('seaborn-v0_8-darkgrid')

# Forecasting
# ==============================================================================
import skforecast
import lightgbm
from lightgbm import LGBMRegressor
from sklearn.feature_selection import RFECV
from skforecast.recursive import ForecasterRecursiveMultiSeries
from skforecast.model_selection import (
     TimeSeriesFold,
     OneStepAheadFold,
     backtesting_forecaster_multiseries,
     bayesian_search_forecaster_multiseries
)
from skforecast.feature_selection import select_features_multiseries
from skforecast.preprocessing import (
     CalendarFeatures,
     RollingFeatures,
     reshape_series_long_to_dict,
     reshape_exog_long_to_dict
)
from feature_engine.timeseries.forecasting import WindowFeatures
from skforecast.datasets import fetch_dataset

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

color = '\033[1m\033[38;5;208m'
print(f'{color}Version skforecast: {skforecast.__version__}')
print(f'{color}Version lightgbm: {lightgbm.__version__}')

Version skforecast: 0.25.0
Version lightgbm: 4.7.0

Data¶

Data used in this document has been obtained from The Building Data Genome Project 2. The dataset contains information on the energy consumption of more than 1500 buildings. The time range of the time-series data is the two full years (2016 and 2017) and the frequency is hourly measurements of electricity, heating and cooling water, steam, and irrigation meters. Additionally, the dataset includes information on the characteristics of the buildings and the weather conditions. Data has been aggregated to a daily resolution and only the electricity among the different energy sources has been considered.

# Load data
# ==============================================================================
data = fetch_dataset(name='bdg2_daily')
print('Data shape:', data.shape)
data.head(3)

╭─────────────────────────────────── bdg2_daily ───────────────────────────────────╮
│ Description:                                                                     │
│ Daily energy consumption data from the The Building Data Genome Project 2 with   │
│ building metadata and weather data. https://github.com/buds-lab/building-data-   │
│ genome-project-2                                                                 │
│                                                                                  │
│ Source:                                                                          │
│ Miller, C., Kathirgamanathan, A., Picchetti, B. et al. The Building Data Genome  │
│ Project 2, energy meter data from the ASHRAE Great Energy Predictor III          │
│ competition. Sci Data 7, 368 (2020). https://doi.org/10.1038/s41597-020-00712-x  │
│                                                                                  │
│ URL:                                                                             │
│ https://huggingface.co/datasets/skforecast/bdg2_daily/resolve/main/bdg2_daily.pa │
│ rquet                                                                            │
│                                                                                  │
│ Shape: 1153518 rows x 17 columns                                                 │
╰──────────────────────────────────────────────────────────────────────────────────╯

Data shape: (1153518, 17)

# Overall range of available dates
# ==============================================================================
print(
    f'Range of dates available : {data.index.min()} --- {data.index.max()}  '
    f'(n_days={(data.index.max() - data.index.min()).days})'
)

Range of dates available : 2016-01-01 00:00:00 --- 2017-12-31 00:00:00  (n_days=730)

# Ensure index of all series is complete without intermediate gaps
# ==============================================================================
data = (
    data
    .groupby('building_id')
    .apply(lambda group: group.asfreq('D', fill_value=np.nan), include_groups=False)
    .reset_index(level=0)
)
print('Data shape:', data.shape)

Data shape: (1153518, 17)

# Range of available dates per building
# ==============================================================================
available_dates_per_series = (
    data
    .dropna(subset='meter_reading')
    .reset_index()
    .groupby('building_id')
    .agg(
        min_index=('timestamp', 'min'),
        max_index=('timestamp', 'max'),
        n_values=('timestamp', 'nunique')
    )
)
display(available_dates_per_series)
print(f'Unique length of series: {available_dates_per_series.n_values.unique()}')

Unique length of series: [731]

All time series have the same length, starting from January 1, 2016 and ending on December 31, 2017. Some exogenous variables contain missing values. Skforecast does not require that the time series have the same length, and missing values are allowed as long as the underlying estimator can handle them, which is the case for LightGBM, XGBoost, and HistGradientBoostingRegressor.

# Missing values per feature
# ==============================================================================
data.isna().mean().mul(100).round(2)

building_id               0.00
meter_reading             0.00
site_id                   0.00
primaryspaceusage         1.20
sub_primaryspaceusage     1.20
sqm                       0.00
lat                      14.83
lng                      14.83
timezone                  0.00
airTemperature            0.02
cloudCoverage             7.02
dewTemperature            0.03
precipDepth1HR            0.02
precipDepth6HR            0.02
seaLvlPressure            9.56
windDirection             0.02
windSpeed                 0.02
dtype: float64

Exogenous variables¶

Exogenous variables are variables that are external to the time series and can be used as features to improve the forecast. In this case, the exogenous variables used are: the characteristics of the buildings, calendar variables and weather conditions.

Building characteristics¶

One of the key attributes associated with each building is its designated use. This feature may play a crucial role in influencing the energy consumption pattern, as distinct uses can significantly impact both the quantity and timing of energy consumption.

# Number of buildings
# ==============================================================================
print(f"Number of buildings: {data['building_id'].nunique()}")
print(f"Number of building types: {data['primaryspaceusage'].nunique()}")
print(f"Number of building subtypes: {data['sub_primaryspaceusage'].nunique()}")

Number of buildings: 1578
Number of building types: 16
Number of building subtypes: 104

For certain type and subtype categories, there are a limited number of buildings in the dataset. Types with fewer than 100 buildings and subtypes with fewer than 50 buildings are grouped into the "Other" category.

# Group infrequent categories
# ==============================================================================
infrequent_types = (
    data
    .drop_duplicates(subset=['building_id'])['primaryspaceusage']
    .value_counts()
    .loc[lambda x: x < 100]
    .index
    .tolist()
)
infrequent_subtypes = (
    data
    .drop_duplicates(subset=['building_id'])['sub_primaryspaceusage']
    .value_counts()
    .loc[lambda x: x < 50]
    .index
    .tolist()
)

data['primaryspaceusage'] = np.where(
    data['primaryspaceusage'].isin(infrequent_types),
    'Other',
    data['primaryspaceusage']
)
data['sub_primaryspaceusage'] = np.where(
    data['sub_primaryspaceusage'].isin(infrequent_subtypes),
    'Other',
    data['sub_primaryspaceusage']
)

display(data.drop_duplicates(subset=['building_id'])['primaryspaceusage'].value_counts(dropna=False))
display(data.drop_duplicates(subset=['building_id'])['sub_primaryspaceusage'].value_counts(dropna=False))

primaryspaceusage
Education                        604
Office                           296
Entertainment/public assembly    203
Public services                  166
Lodging/residential              149
Other                            141
NaN                               19
Name: count, dtype: int64

sub_primaryspaceusage
Other                          612
Office                         295
College Classroom              131
College Laboratory             116
K-12 School                    109
Dormitory                       91
Primary/Secondary Classroom     84
Education                       67
Library                         54
NaN                             19
Name: count, dtype: int64

Calendar features¶

# Calendar features
# ==============================================================================
features_to_extract = [
    'month',
    'week',
    'day_of_week',
]
calendar_transformer = CalendarFeatures(
                            features = features_to_extract,
                            encoding = 'cyclical',
                            keep_original_columns = True,
                       )
data = calendar_transformer.fit_transform(data)
print(f'Data shape: {data.shape}')
data.head(3)

Data shape: (1153518, 23)

Meteorological features¶

Meteorological variables are recorded at the site level, meaning that weather data varies by building location, even at the same timestamp. In other words, while the exogenous variables are consistent across all series, their values differ by location.

# Values of meteorological features for a given date in each site
# ==============================================================================
data.loc['2016-01-01'].groupby('site_id', observed=True).agg(
    {
        'airTemperature': 'first',
        'cloudCoverage': 'first',
        'dewTemperature': 'first',
        'precipDepth1HR': 'first',
        'precipDepth6HR': 'first',
        'seaLvlPressure': 'first',
        'windDirection': 'first',
        'windSpeed': 'first',
    }
)

Skforecast allows you to include different exogenous variables and/or different values for each series in the dataset (more details in the next section).

Rolling features for exogenous variables¶

# Rolling features of meteorological features
# ==============================================================================
wf_transformer = WindowFeatures(
    variables      = ['airTemperature', 'windSpeed'],
    window         = ['7D', '14D'],
    functions      = ['mean'],
    freq           = 'D',
    missing_values = 'ignore',
    drop_na        = False,
)
data = data.groupby('building_id').apply(wf_transformer.fit_transform, include_groups=False).reset_index(level=0)
print(f'Data shape: {data.shape}')
data.head(3)

Data shape: (1153518, 27)

Categorical features¶

Since version 0.22.0, skforecast provides a built-in categorical_features parameter that automatically handles encoding and natively configures gradient boosting estimators (XGBoost, LightGBM, CatBoost and HistGradientBoosting), requiring no manual encoder pipelines or estimator-specific parameters.

Modelling and forecasting¶

ForecasterRecursiveMultiSeries enables modeling of time series with varying lengths and different exogenous variables.

When the series have different lengths, the data must be transformed into either a dictionary or a DataFrame with a MultiIndex. If using a dictionary, the keys should be the series names and the values the corresponding series data. If using a MultiIndex DataFrame, the first level of the index represents the series name, while the second level corresponds to the DateTime temporal index.

The same applies to exogenous variables when they vary across series. These variables should be transformed into a dictionary or a MultiIndex DataFrame, where the dictionary keys are the series names and the values are the respective exogenous variables.

Skforecast simplifies this data transformation process with the functions reshape_series_long_to_dict and reshape_exog_long_to_dict.

If all series have the same length and share the same exogenous variables, dictionaries are not required. In this case, the series can be provided as a single DataFrame with each series as a separate column, and the exogenous variables can be passed as a DataFrame matching the series length.

# Exogenous features selected for modeling
# ==============================================================================
exog_features = [
    'primaryspaceusage',
    'sub_primaryspaceusage',
    'timezone',
    'sqm',
    'airTemperature',
    'cloudCoverage',
    'dewTemperature',
    'precipDepth1HR',
    'precipDepth6HR',
    'seaLvlPressure',
    'windDirection',
    'windSpeed',
    'day_of_week_sin',
    'day_of_week_cos',
    'week_sin',
    'week_cos',
    'month_sin',
    'month_cos',
    'airTemperature_window_7D_mean',
    'windSpeed_window_7D_mean',
    'airTemperature_window_14D_mean',
    'windSpeed_window_14D_mean',
]

# Transform series and exog to dictionaries
# ==============================================================================
series_dict = reshape_series_long_to_dict(
    data      = data.reset_index(),
    series_id = 'building_id',
    index     = 'timestamp',
    values    = 'meter_reading',
    freq      = 'D'
)

exog_dict = reshape_exog_long_to_dict(
    data      = data[exog_features + ['building_id']].reset_index(),
    series_id = 'building_id',
    index     = 'timestamp',
    freq      = 'D'
)

To train the models, search for optimal hyperparameters, and evaluate their predictive performance, the data is divided into three separate sets: training, validation, and test.

# Partition data in train and test
# ==============================================================================
data = data.sort_index()
end_train = '2017-08-31 23:59:00'
end_validation = '2017-10-31 23:59:00'
# Only the partitions used later in the document are stored
series_dict_train = {k: v.loc[:end_train] for k, v in series_dict.items()}
exog_dict_train = {k: v.loc[:end_train] for k, v in exog_dict.items()}
exog_dict_valid = {k: v.loc[end_train:end_validation] for k, v in exog_dict.items()}

print(
    f'Range of dates available : {data.index.min()} --- {data.index.max()} '
    f'(n_days={(data.index.max() - data.index.min()).days})'
)
print(
    f'  Dates for training    : {data.loc[: end_train, :].index.min()} --- {data.loc[: end_train, :].index.max()} '
    f'(n_days={(data.loc[: end_train, :].index.max() - data.loc[: end_train, :].index.min()).days})'
)
print(
    f'  Dates for validation  : {data.loc[end_train:end_validation, :].index.min()} --- {data.loc[end_train:end_validation, :].index.max()} '
    f'(n_days={(data.loc[end_train:end_validation, :].index.max() - data.loc[end_train:end_validation, :].index.min()).days})'
)
print(
    f'  Dates for test        : {data.loc[end_validation:, :].index.min()} --- {data.loc[end_validation:, :].index.max()} '
    f'(n_days={(data.loc[end_validation:, :].index.max() - data.loc[end_validation:, :].index.min()).days})'
)

Range of dates available : 2016-01-01 00:00:00 --- 2017-12-31 00:00:00 (n_days=730)
  Dates for training    : 2016-01-01 00:00:00 --- 2017-08-31 00:00:00 (n_days=608)
  Dates for validation  : 2017-09-01 00:00:00 --- 2017-10-31 00:00:00 (n_days=60)
  Dates for test        : 2017-11-01 00:00:00 --- 2017-12-31 00:00:00 (n_days=60)

Hyperparameter tuning¶

Hyperparameter and lag tuning involves systematically testing different values or combinations of hyperparameters (and/or lags) to find the optimal configuration that gives the best performance. The skforecast library provides two different methods to evaluate each candidate configuration, implemented through the TimeSeriesFold and OneStepAheadFold classes:

• Backtesting (TimeSeriesFold): In this method, the model predicts several steps ahead in each iteration, using the same forecast horizon and retraining frequency strategy that would be used if the model were deployed (including the option of no retraining at all, i.e. refit=False, used in the examples below for speed). This simulates a real forecasting scenario where the model is retrained and updated over time.

• One-Step Ahead (OneStepAheadFold): Evaluates the model using only one-step-ahead predictions. This method is faster because it requires fewer iterations, but it only tests the model's performance in the immediate next time step (t+1).

Each method uses a different evaluation strategy, so they may produce different results. However, in the long run, both methods are expected to converge to similar selections of optimal hyperparameters. It is recommended to backtest the final model for a more accurate multi-step performance estimate.

In this example, the search is performed with the bayesian_search_forecaster_multiseries function, which relies on Optuna to efficiently explore the hyperparameter space. To keep the runtime of this example short, only 10 trials are evaluated (n_trials=10); in a real application, a much larger number of trials is recommended to properly explore the search space.

# Create forecaster
# ==============================================================================
window_features = RollingFeatures(stats=['mean', 'min', 'max'], window_sizes=7)
forecaster = ForecasterRecursiveMultiSeries(
                estimator            = LGBMRegressor(random_state=8520, verbose=-1),
                lags                 = 14,
                window_features      = window_features,
                categorical_features = 'auto',
                encoding             = 'ordinal'
            )

# Bayesian search with OneStepAheadFold
# ==============================================================================
def search_space(trial):
    search_space  = {
        'lags'             : trial.suggest_categorical('lags', [31, 62]),
        'n_estimators'     : trial.suggest_int('n_estimators', 200, 800, step=100),
        'learning_rate'    : trial.suggest_float('learning_rate', 0.01, 0.3, log=True),
        'num_leaves'       : trial.suggest_int('num_leaves', 15, 255, log=True),
        'max_depth'        : trial.suggest_int('max_depth', 3, 8),
        'min_child_samples': trial.suggest_int('min_child_samples', 5, 200, log=True),
        'colsample_bytree' : trial.suggest_float('colsample_bytree', 0.6, 1.0),
        'subsample'        : trial.suggest_float('subsample', 0.6, 1.0),
        'subsample_freq'   : trial.suggest_int('subsample_freq', 0, 5),
        'reg_alpha'        : trial.suggest_float('reg_alpha', 1e-3, 10.0, log=True),
        'reg_lambda'       : trial.suggest_float('reg_lambda', 1e-3, 10.0, log=True),
    }

    return search_space

cv = OneStepAheadFold(initial_train_size='2017-08-31')  # last date included in the training set

results_search, best_trial = bayesian_search_forecaster_multiseries(
    forecaster        = forecaster,
    series            = {k: v.loc[:end_validation] for k, v in series_dict.items()},
    exog              = {k: v.loc[:end_validation, exog_features] for k, v in exog_dict.items()},
    cv                = cv,
    search_space      = search_space,
    n_trials          = 10,
    random_state      = 123,
    metric            = 'mean_absolute_error',
    suppress_warnings = True
)

best_params = results_search.at[0, 'params']
best_lags = results_search.at[0, 'lags']
print(f'Best lags: {best_lags}')
print(f'Best params: {best_params}')
results_search.head(3)

Best 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
 49 50 51 52 53 54 55 56 57 58 59 60 61 62]
Best params: {'n_estimators': 300, 'learning_rate': 0.032160750764041963, 'num_leaves': 63, 'max_depth': 6, 'min_child_samples': 7, 'colsample_bytree': 0.6523579802656323, 'subsample': 0.7287922425873232, 'subsample_freq': 3, 'reg_alpha': 2.4323434680180993, 'reg_lambda': 0.16331624244561055}

Backtesting on test data¶

After selecting the best combination of lags and hyperparameters, the model's predictive performance is evaluated on the test set (November and December 2017) using the backtesting_forecaster_multiseries function. In each fold, the forecaster predicts 7 days at a time without retraining (refit=False). Since the search was run with the default return_best=True, the forecaster is automatically updated with the best configuration found, so it can be passed directly to the backtesting function.

# Backtesting on test set
# ==============================================================================
cv = TimeSeriesFold(
    initial_train_size = '2017-10-31',  # last date of training + validation
    steps              = 7,
    refit              = False
)
metrics, predictions = backtesting_forecaster_multiseries(
                          forecaster        = forecaster,
                          series            = series_dict,
                          exog              = exog_dict,
                          cv                = cv,
                          metric            = 'mean_absolute_error',
                          suppress_warnings = True
                      )

display(predictions.head())
display(metrics)

# Aggregated metric for all buildings
# ==============================================================================
average_metric_all_buildings = metrics.query("levels == 'average'")['mean_absolute_error'].item()
errors_all_buildings = pd.merge(
    left     = data[['building_id', 'meter_reading']].reset_index(),
    right    = predictions.rename_axis('timestamp').reset_index(),
    left_on  = ['timestamp', 'building_id'],
    right_on = ['timestamp', 'level'],
    how      = 'inner',
    validate = '1:1'
).assign(error=lambda df: df['meter_reading'] - df['pred'])

sum_abs_errors_all_buildings = errors_all_buildings['error'].abs().sum()
sum_bias_all_buildings = errors_all_buildings['error'].sum()
print(f'Average mean absolute error for all buildings: {average_metric_all_buildings:.0f}')
print(f'Sum of absolute errors for all buildings (x 10,000): {sum_abs_errors_all_buildings/10000:.0f}')
print(f'Bias (x 10,000): {sum_bias_all_buildings/10000:.0f}')

Average mean absolute error for all buildings: 311
Sum of absolute errors for all buildings (x 10,000): 2990
Bias (x 10,000): -316

# Plot predictions vs real value for 2 random buildings
# ==============================================================================
rng = np.random.default_rng(14793)
n_buildings = 2
selected_buildings = rng.choice(data['building_id'].unique(), size=n_buildings, replace=False)

fig, axs = plt.subplots(n_buildings, 1, figsize=(7, 4.5), sharex=True)
axs = axs.flatten()

for i, building in enumerate(selected_buildings):
    data.query('building_id == @building').loc[predictions.index.unique(), 'meter_reading'].plot(ax=axs[i], label='test')
    predictions.query('level == @building')['pred'].plot(ax=axs[i], label='predictions')
    axs[i].set_title(f'Building {building}', fontsize=10)
    axs[i].set_xlabel('')
    axs[i].legend()

fig.tight_layout()
plt.show();

Feature selection¶

Feature selection is the process of selecting a subset of relevant features (variables, predictors) for use in model construction. Feature selection techniques are used for several reasons: to simplify models to make them easier to interpret, to reduce training time, to avoid the curse of dimensionality, to improve generalization by reducing overfitting (formally, variance reduction), and others.

Skforecast is compatible with the feature selection methods implemented in the scikit-learn library. There are several methods for feature selection, but the most common are:

• Recursive feature elimination (RFE)

• Sequential Feature Selection (SFS)

• Feature selection based on threshold (SelectFromModel)

Skforecast provides the function select_features_multiseries to apply any of these scikit-learn selectors directly on the lags, window features, exogenous variables and calendar features of a ForecasterRecursiveMultiSeries.

# Feature selection (autoregressive and exog) with scikit-learn RFECV
# ==============================================================================
warnings.filterwarnings('ignore', message='X does not have valid feature names.*')
estimator = LGBMRegressor(n_estimators=50, num_leaves=8, max_bin=63, random_state=15926, verbose=-1, n_jobs=-1)
selector = RFECV(estimator=estimator, step=3, cv=3, n_jobs=1)
selected_lags, selected_window_features, selected_exog, _ = select_features_multiseries(
    forecaster      = forecaster,
    selector        = selector,
    series          = {k: v.loc[:end_validation] for k, v in series_dict.items()},
    exog            = {k: v.loc[:end_validation, exog_features] for k, v in exog_dict.items()},
    select_only     = None,
    force_inclusion = None,
    subsample       = 0.2,
    random_state    = 123,
    verbose         = True,
)

╭──────────────────────────────── MissingValuesWarning ────────────────────────────────╮
│ NaNs detected in `X_train`. Some estimators do not allow NaN values during training. │
│ If you want to drop them, set `forecaster.dropna_from_series = True`.                │
│                                                                                      │
│ Category : skforecast.exceptions.MissingValuesWarning                                │
│ Location :                                                                           │
│ C:\Users\Joaquin\miniconda3\envs\skforecast_24_py13\Lib\site-packages\skforecast\rec │
│ ursive\_forecaster_recursive_multiseries.py:1456                                     │
│ Suppress : warnings.simplefilter('ignore', category=MissingValuesWarning)            │
╰──────────────────────────────────────────────────────────────────────────────────────╯

Recursive feature elimination (RFECV)
-------------------------------------
Total number of records available: 959424
Total number of records used for feature selection: 191884
Number of features available: 87
    Lags            (n=62)
    Window features (n=3)
    Exog            (n=22)
    Calendar        (n=0)
Number of features selected: 45
    Lags            (n=28) : [1, 5, 7, 8, 9, 10, 13, 14, 15, 20, 21, 23, 26, 28, 34, 35, 37, 42, 48, 49, 54, 56, 57, 58, 59, 60, 61, 62]
    Window features (n=3) : ['roll_mean_7', 'roll_min_7', 'roll_max_7']
    Exog            (n=14) : ['primaryspaceusage', 'sub_primaryspaceusage', 'timezone', 'sqm', 'airTemperature', 'cloudCoverage', 'dewTemperature', 'windDirection', 'day_of_week_sin', 'day_of_week_cos', 'airTemperature_window_7D_mean', 'windSpeed_window_7D_mean', 'airTemperature_window_14D_mean', 'windSpeed_window_14D_mean']
    Calendar        (n=0) : []

# Backtesting forecaster with selected features
# ==============================================================================
selected_stats = [feature.split('_')[1] for feature in selected_window_features]
window_features = RollingFeatures(stats=selected_stats, window_sizes=7)
forecaster = ForecasterRecursiveMultiSeries(
                estimator            = LGBMRegressor(**best_params, random_state=8520, verbose=-1),
                lags                 = selected_lags,
                window_features      = window_features,
                categorical_features = 'auto',
                encoding             = 'ordinal'
            )
cv = TimeSeriesFold(
    initial_train_size = '2017-10-31',  # last date of training + validation
    steps              = 7,
    refit              = False
)
metrics, predictions = backtesting_forecaster_multiseries(
                          forecaster        = forecaster,
                          series            = series_dict,
                          exog              = {k: v[selected_exog] for k, v in exog_dict.items()},
                          cv                = cv,
                          metric            = 'mean_absolute_error',
                          suppress_warnings = True
                      )

display(predictions.head())
display(metrics)

The number of lags, window features, and exogenous features included in the model has been reduced with almost no impact on the model performance.

Note that the calendar features shown as Calendar (n=0) in the feature selection output are not evaluated as their own group here because they were passed as regular exogenous variables rather than through the forecaster's calendar_features parameter; their cyclical columns are evaluated as part of the exogenous set instead.

Predict new series (unknown series)¶

ForecasterRecursiveMultiSeries allows predicting unknown series (series not seen during the training process). Two scenarios may occur:

• There is historical data for the unknown series: The user must provide a DataFrame with the historical data required to create the features (lags and window features) used during the prediction process. This DataFrame must contain one column per series to be predicted. It is possible to include both series seen during training and new (unseen) series.

• There is no historical data for the unknown series: This scenario is similar, but the last_window for the new series is composed entirely of NaN values. This approach is only possible when using estimators that support missing values, since all lagged features will be NaN. If no exogenous variables are available either, all model inputs are NaN, so the model returns the same prediction for every step and series, as shown in the second example below.

In both cases, exogenous variables may be available or not. If they are not available, they will be automatically set to NaN.

# Predict new series with available historical data
# ==============================================================================
# Create the last window with the last observations of the new series
forecaster.fit(
    series = series_dict_train,
    exog   = {k: v[selected_exog] for k, v in exog_dict_train.items()},
)
rng = np.random.default_rng(67584)
last_window_new = pd.DataFrame(
    {
        'series_new_1': rng.normal(loc=1000, scale=500, size=62),
        'series_new_2': rng.normal(loc=1000, scale=500, size=62),
    },
    index=pd.date_range('2017-07-01', periods=62, freq='D'),
)

# If exogenous variables are not available for the new series, they are automatically set to NaN
forecaster.predict(
    steps       = 5,
    last_window = last_window_new,
    exog        = {k: v[selected_exog] for k, v in exog_dict_valid.items()}
)

╭──────────────────────────────── MissingValuesWarning ────────────────────────────────╮
│ NaNs detected in `X_train`. Some estimators do not allow NaN values during training. │
│ If you want to drop them, set `forecaster.dropna_from_series = True`.                │
│                                                                                      │
│ Category : skforecast.exceptions.MissingValuesWarning                                │
│ Location :                                                                           │
│ C:\Users\Joaquin\miniconda3\envs\skforecast_24_py13\Lib\site-packages\skforecast\rec │
│ ursive\_forecaster_recursive_multiseries.py:1456                                     │
│ Suppress : warnings.simplefilter('ignore', category=MissingValuesWarning)            │
╰──────────────────────────────────────────────────────────────────────────────────────╯

╭──────────────────────────────── UnknownLevelWarning ─────────────────────────────────╮
│ `levels` {'series_new_2', 'series_new_1'} were not included in training. Unknown     │
│ levels are encoded as NaN, which may cause the prediction to fail if the estimator   │
│ does not accept NaN values.                                                          │
│                                                                                      │
│ Category : skforecast.exceptions.UnknownLevelWarning                                 │
│ Location :                                                                           │
│ C:\Users\Joaquin\miniconda3\envs\skforecast_24_py13\Lib\site-packages\skforecast\uti │
│ ls\utils.py:1258                                                                     │
│ Suppress : warnings.simplefilter('ignore', category=UnknownLevelWarning)             │
╰──────────────────────────────────────────────────────────────────────────────────────╯

╭───────────────────────────────── MissingExogWarning ─────────────────────────────────╮
│ `exog` does not contain keys for levels {'series_new_2', 'series_new_1'}. Missing    │
│ levels are filled with NaN. Most of machine learning models do not allow missing     │
│ values. Prediction method may fail.                                                  │
│                                                                                      │
│ Category : skforecast.exceptions.MissingExogWarning                                  │
│ Location :                                                                           │
│ C:\Users\Joaquin\miniconda3\envs\skforecast_24_py13\Lib\site-packages\skforecast\uti │
│ ls\utils.py:1384                                                                     │
│ Suppress : warnings.simplefilter('ignore', category=MissingExogWarning)              │
╰──────────────────────────────────────────────────────────────────────────────────────╯

# Predict new series without available historical data
# ==============================================================================
# Create the last window of the new series with NaN values
last_window_new = pd.DataFrame(
    {
        'series_new_1': np.nan,
        'series_new_2': np.nan,
    },
    index=pd.date_range('2017-07-01', periods=62, freq='D'),
)

# If exogenous variables are not available for the new series, they are automatically set to NaN
forecaster.predict(
    steps       = 5,
    last_window = last_window_new,
    exog        = {k: v[selected_exog] for k, v in exog_dict_valid.items()}
)

╭──────────────────────────────── UnknownLevelWarning ─────────────────────────────────╮
│ `levels` {'series_new_2', 'series_new_1'} were not included in training. Unknown     │
│ levels are encoded as NaN, which may cause the prediction to fail if the estimator   │
│ does not accept NaN values.                                                          │
│                                                                                      │
│ Category : skforecast.exceptions.UnknownLevelWarning                                 │
│ Location :                                                                           │
│ C:\Users\Joaquin\miniconda3\envs\skforecast_24_py13\Lib\site-packages\skforecast\uti │
│ ls\utils.py:1258                                                                     │
│ Suppress : warnings.simplefilter('ignore', category=UnknownLevelWarning)             │
╰──────────────────────────────────────────────────────────────────────────────────────╯

╭──────────────────────────────── MissingValuesWarning ────────────────────────────────╮
│ `last_window` has missing values. Most of machine learning models do not allow       │
│ missing values. Prediction method may either raise an error or return NaN            │
│ predictions.                                                                         │
│                                                                                      │
│ Category : skforecast.exceptions.MissingValuesWarning                                │
│ Location :                                                                           │
│ C:\Users\Joaquin\miniconda3\envs\skforecast_24_py13\Lib\site-packages\skforecast\uti │
│ ls\utils.py:1338                                                                     │
│ Suppress : warnings.simplefilter('ignore', category=MissingValuesWarning)            │
╰──────────────────────────────────────────────────────────────────────────────────────╯

╭───────────────────────────────── MissingExogWarning ─────────────────────────────────╮
│ `exog` does not contain keys for levels {'series_new_2', 'series_new_1'}. Missing    │
│ levels are filled with NaN. Most of machine learning models do not allow missing     │
│ values. Prediction method may fail.                                                  │
│                                                                                      │
│ Category : skforecast.exceptions.MissingExogWarning                                  │
│ Location :                                                                           │
│ C:\Users\Joaquin\miniconda3\envs\skforecast_24_py13\Lib\site-packages\skforecast\uti │
│ ls\utils.py:1384                                                                     │
│ Suppress : warnings.simplefilter('ignore', category=MissingExogWarning)              │
╰──────────────────────────────────────────────────────────────────────────────────────╯

Clustering time series¶

The idea behind modeling multiple series at the same time is to be able to capture the main patterns that govern the series, thereby reducing the impact of the potential noise that each series may have. This means that series that behave similarly may benefit from being modeled together. One way to identify potential groups of series is to perform a clustering study prior to modeling. If clear groups are identified as a result of clustering, it is appropriate to model each of them separately.

Clustering is an unsupervised analysis technique that groups a set of observations into clusters that contain observations that are considered homogeneous, while observations in different clusters are considered heterogeneous. Algorithms that cluster time series can be divided into two groups: those that use a transformation to create features prior to clustering (feature-driven time series clustering), and those that work directly on the time series (elastic distance measures).

• Feature-driven time series clustering: Features describing structural characteristics are extracted from each time series, then these features are fed into arbitrary clustering algorithms. These features are obtained by applying statistical operations that best capture the underlying characteristics: trend, seasonality, periodicity, serial correlation, skewness, kurtosis, chaos, nonlinearity, and self-similarity.

• Elastic distance measures: This approach works directly on the time series, adjusting or "realigning" the series in comparison to each other. The best known of this family of measures is Dynamic Time Warping (DTW).

For a detailed example of how time series clustering can improve the forecasting models, see Clustering Time Series to Improve Forecasting Models.

Session information¶

import session_info
session_info.show(html=False)

-----
feature_engine      1.9.4
lightgbm            4.7.0
matplotlib          3.10.9
numpy               2.4.6
optuna              4.9.0
pandas              2.3.3
session_info        v1.0.1
skforecast          0.25.0
sklearn             1.7.2
-----
IPython             9.15.0
jupyter_client      8.9.1
jupyter_core        5.9.1
-----
Python 3.13.14 | packaged by conda-forge | (main, Jun 12 2026, 09:44:26) [MSC v.1944 64 bit (AMD64)]
Windows-11-10.0.26200-SP0
-----
Session information updated at 2026-09-24 15:54

Citation¶

How to cite this document

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

Forecasting at scale: modeling thousands of time series with a single global model by Joaquín Amat Rodrigo and Javier Escobar Ortiz, available under a CC BY-NC-SA 4.0 at https://www.cienciadedatos.net/documentos/py59-scalable-forecasting-models.html

How to cite skforecast

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

Zenodo:

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

APA:

Amat Rodrigo, J., & Escobar Ortiz, J. (2026). skforecast (Version 0.25.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.25.0}, month = {09}, year = {2026}, license = {BSD-3-Clause}, url = {https://skforecast.org/}, doi = {10.5281/zenodo.8382788} }

Did you like the article? Your support is important

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.