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 deep learning model including exogenous features.

Libraries

Libraries used in this document.

import os
os.environ["KERAS_BACKEND"] = "torch"  # Options: "tensorflow", "jax", or "torch"
import keras

# Data management
# ==============================================================================
import numpy as np
import pandas as pd
from skforecast.datasets import fetch_dataset

# Plots
# ==============================================================================
import matplotlib.pyplot as plt
from skforecast.plot import set_dark_theme

# Forecasting
# ==============================================================================
import skforecast
from skforecast.deep_learning import ForecasterRnn
from skforecast.deep_learning.utils import create_and_compile_model
from skforecast.model_selection import TimeSeriesFold
from skforecast.model_selection import backtesting_forecaster_multiseries
from feature_engine.datetime import DatetimeFeatures
from feature_engine.creation import CyclicalFeatures
from sklearn.preprocessing import MinMaxScaler
import keras
from keras.optimizers import Adam
from keras.losses import MeanSquaredError
from keras.callbacks import EarlyStopping
import torch

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

print('Versi0n skforecast:', skforecast.__version__)
print('Version keras:', keras.__version__)
print('Version torch:', torch.__version__)

Versi0n skforecast: 0.19.0
Version keras: 3.12.0
Version torch: 2.8.0

Data

Data used in this document has been obtained from the The Building Data Genome Project 2 https://github.com/buds-lab/building-data-genome-project-2. The dataset contains information on the energy consumption of more than 1500 buildings. The time range of the times-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 of 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://drive.google.com/file/d/1KHYopzclKvS1F6Gt6GoJWKnxiuZ2aqen/view?usp=shari │
│ ng                                                                               │
│                                                                                  │
│ Shape: 1153518 rows x 17 columns                                                 │
╰──────────────────────────────────────────────────────────────────────────────────╯

Data shape: (1153518, 17)

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

Rage 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, drop=False)
)
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 lenght of series: {available_dates_per_series.n_values.unique()}")

Unique lenght of series: [731]

All time series have the same length, starting from January 1, 2016 and ending on December 31, 2017.

ForecasterRNN requires the input data to be a pandas DataFrame where each column represents a different time series and the index is a pandas DatetimeIndex.

# Reshape data to wide format where each column is a time series
# ==============================================================================
series = data.reset_index().pivot(
    index='timestamp',
    columns='building_id',
    values='meter_reading'
).sort_index()
series = series.asfreq('D', fill_value=np.nan)
series.head(3)

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.

# Calendar features
# ==============================================================================
calendar_features = series.index.to_frame()
features_to_extract = [
    'month',
    'week',
    'day_of_week',
]
calendar_transformer = DatetimeFeatures(
                            variables           = 'index',
                            features_to_extract = features_to_extract,
                            drop_original       = False,
                       )
calendar_features = calendar_transformer.fit_transform(calendar_features)
calendar_features = calendar_features.drop(columns=['timestamp'])
print(f"Data shape: {calendar_features.shape}")
calendar_features.head(3)

Data shape: (731, 3)

# Cyclical encoding of calendar features
# ==============================================================================
features_to_encode = [
    "month",
    "week",
    "day_of_week",
]
max_values = {
    "month": 12,
    "week": 52,
    "day_of_week": 6,
}
cyclical_encoder = CyclicalFeatures(
                        variables     = features_to_encode,
                        max_values    = max_values,
                        drop_original = True
                   )

calendar_features = cyclical_encoder.fit_transform(calendar_features)
print(f"Data shape: {calendar_features.shape}")
calendar_features.head(3)

Data shape: (731, 6)

# Values of meteorological features for a guiven date in each site
# ==============================================================================
meteo_features = (
    data
    .reset_index()
    .drop_duplicates(subset=["timestamp", "building_id"])
    .groupby(by="timestamp", observed=True).agg(
        {
            "airTemperature": "median",
            "cloudCoverage":  "median",
            "dewTemperature": "median",
            "precipDepth1HR": "median",
            "precipDepth6HR": "median",
            "seaLvlPressure": "median",
            "windDirection":  "median",
            "windSpeed":      "median",
        }
    )
)
meteo_features.head(3)

Modelling and forecasting

# Exogenous features selected for modeling
# ==============================================================================
exog_features = pd.merge(
    left=calendar_features,
    right=meteo_features,
    left_index=True,
    right_index=True,
    validate="one_to_one"
)
exog_features = exog_features.asfreq('D', fill_value=np.nan)
exog_features.head(3)

# Check that series and exogenous features are aligned
# ==============================================================================
series.index.equals(exog_features.index)

True

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
# ==============================================================================
end_train = '2017-08-31 23:59:00'
end_validation = '2017-10-31 23:59:00'
series_train = series.loc[: end_train,]
series_val = series.loc[end_train: end_validation,]
series_test = series.loc[end_validation:,]
exog_train = exog_features.loc[: end_train,]
exog_val = exog_features.loc[end_train: end_validation,]
exog_test = exog_features.loc[end_validation:,]

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

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

Deep learning model and forecaster

# Create deep learning model
# ==============================================================================
levels = series_train.columns.tolist()
lags = 48

model = create_and_compile_model(
    series                  = series[levels], # DataFrame with all series (predictors)
    exog                    = exog_features, # DataFrame with exogenous features
    levels                  = levels,        # List of series to be predicted by the model
    lags                    = lags, 
    steps                   = 7,
    recurrent_layer         = "LSTM",
    recurrent_units         = [256, 128, 64, 32],
    recurrent_layers_kwargs = {"activation": "tanh"},
    dense_units             = [32],
    compile_kwargs          = {"optimizer": Adam(learning_rate=0.01), "loss": MeanSquaredError()},
    model_name              = "global_model_lstm"
)

model.summary()    

keras version: 3.12.0
Using backend: torch
torch version: 2.8.0

Model: "global_model_lstm"

┏━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━┓
┃ Layer (type)        ┃ Output Shape      ┃    Param # ┃ Connected to      ┃
┡━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━┩
│ series_input        │ (None, 48, 1578)  │          0 │ -                 │
│ (InputLayer)        │                   │            │                   │
├─────────────────────┼───────────────────┼────────────┼───────────────────┤
│ lstm_1 (LSTM)       │ (None, 48, 256)   │  1,879,040 │ series_input[0][… │
├─────────────────────┼───────────────────┼────────────┼───────────────────┤
│ lstm_2 (LSTM)       │ (None, 48, 128)   │    197,120 │ lstm_1[0][0]      │
├─────────────────────┼───────────────────┼────────────┼───────────────────┤
│ lstm_3 (LSTM)       │ (None, 48, 64)    │     49,408 │ lstm_2[0][0]      │
├─────────────────────┼───────────────────┼────────────┼───────────────────┤
│ lstm_4 (LSTM)       │ (None, 32)        │     12,416 │ lstm_3[0][0]      │
├─────────────────────┼───────────────────┼────────────┼───────────────────┤
│ repeat_vector       │ (None, 7, 32)     │          0 │ lstm_4[0][0]      │
│ (RepeatVector)      │                   │            │                   │
├─────────────────────┼───────────────────┼────────────┼───────────────────┤
│ exog_input          │ (None, 7, 14)     │          0 │ -                 │
│ (InputLayer)        │                   │            │                   │
├─────────────────────┼───────────────────┼────────────┼───────────────────┤
│ concat_exog         │ (None, 7, 46)     │          0 │ repeat_vector[0]… │
│ (Concatenate)       │                   │            │ exog_input[0][0]  │
├─────────────────────┼───────────────────┼────────────┼───────────────────┤
│ dense_td_1          │ (None, 7, 32)     │      1,504 │ concat_exog[0][0] │
│ (TimeDistributed)   │                   │            │                   │
├─────────────────────┼───────────────────┼────────────┼───────────────────┤
│ output_dense_td_la… │ (None, 7, 1578)   │     52,074 │ dense_td_1[0][0]  │
│ (TimeDistributed)   │                   │            │                   │
└─────────────────────┴───────────────────┴────────────┴───────────────────┘

 Total params: 2,191,562 (8.36 MB)

 Trainable params: 2,191,562 (8.36 MB)

 Non-trainable params: 0 (0.00 B)

# Forecaster Creation
# ==============================================================================
forecaster = ForecasterRnn(
    estimator=model,
    levels=levels,
    lags=lags,
    transformer_series=MinMaxScaler(),
    transformer_exog=MinMaxScaler(),
    fit_kwargs={
        "epochs": 50, 
        "batch_size": 128, 
        "callbacks": [
            EarlyStopping(monitor="val_loss", patience=10, restore_best_weights=True)
        ],  # Callback to stop training when it is no longer learning.
        "series_val": series_val[levels],  # Validation data for model training.
        "exog_val": exog_val
    }
)

# Fit forecaster
# ==============================================================================
forecaster.fit(series = series_train[levels], exog = exog_train)
forecaster

Using 'torch' backend with device: cuda

Epoch 1/50
5/5 ━━━━━━━━━━━━━━━━━━━━ 5s 942ms/step - loss: 0.2168 - val_loss: 19.0227
Epoch 2/50

5/5 ━━━━━━━━━━━━━━━━━━━━ 4s 756ms/step - loss: 0.0608 - val_loss: 18.9965
Epoch 3/50
5/5 ━━━━━━━━━━━━━━━━━━━━ 3s 687ms/step - loss: 0.0577 - val_loss: 18.9876
Epoch 4/50
5/5 ━━━━━━━━━━━━━━━━━━━━ 3s 673ms/step - loss: 0.0466 - val_loss: 18.9904
Epoch 5/50
5/5 ━━━━━━━━━━━━━━━━━━━━ 4s 775ms/step - loss: 0.0453 - val_loss: 18.9890
Epoch 6/50
5/5 ━━━━━━━━━━━━━━━━━━━━ 4s 688ms/step - loss: 0.0424 - val_loss: 18.9826
Epoch 7/50
5/5 ━━━━━━━━━━━━━━━━━━━━ 3s 672ms/step - loss: 0.0406 - val_loss: 18.9868
Epoch 8/50
5/5 ━━━━━━━━━━━━━━━━━━━━ 3s 670ms/step - loss: 0.0394 - val_loss: 18.9831
Epoch 9/50
5/5 ━━━━━━━━━━━━━━━━━━━━ 4s 766ms/step - loss: 0.0388 - val_loss: 18.9875
Epoch 10/50
5/5 ━━━━━━━━━━━━━━━━━━━━ 3s 678ms/step - loss: 0.0382 - val_loss: 18.9872
Epoch 11/50
5/5 ━━━━━━━━━━━━━━━━━━━━ 3s 678ms/step - loss: 0.0377 - val_loss: 18.9869
Epoch 12/50
5/5 ━━━━━━━━━━━━━━━━━━━━ 3s 684ms/step - loss: 0.0371 - val_loss: 18.9875
Epoch 13/50
5/5 ━━━━━━━━━━━━━━━━━━━━ 3s 706ms/step - loss: 0.0368 - val_loss: 18.9927
Epoch 14/50
5/5 ━━━━━━━━━━━━━━━━━━━━ 4s 736ms/step - loss: 0.0363 - val_loss: 18.9932
Epoch 15/50
5/5 ━━━━━━━━━━━━━━━━━━━━ 5s 1s/step - loss: 0.0361 - val_loss: 18.9926
Epoch 16/50
5/5 ━━━━━━━━━━━━━━━━━━━━ 6s 1s/step - loss: 0.0358 - val_loss: 18.9969

When the magnitude of the loss is very different between training and validation, it is recommended to plot the training and validation loss using different y-axes. This approach allows for a clearer visualization of the loss trends, making it easier to identify overfitting or underfitting issues in the model.

# Track training and overfitting (different axes for loss and val_loss)
# ==============================================================================
set_dark_theme()
fig, ax = plt.subplots(figsize=(7, 3))
epochs = np.arange(len(forecaster.history_['loss']))
p1, = ax.plot(epochs, forecaster.history_['loss'], color='tab:blue', label='loss')
ax.set_ylabel('loss', color='tab:blue')
ax.tick_params(axis='y', labelcolor='tab:blue')

ax2 = ax.twinx()
p2, = ax2.plot(epochs, forecaster.history_['val_loss'], color='tab:orange', label='val_loss')
ax2.set_ylabel('val_loss', color='tab:orange')
ax2.tick_params(axis='y', labelcolor='tab:orange')

ax.legend(handles=[p1, p2], loc='upper right')
ax.set_title('Training and Validation Loss over Epochs')
plt.show()

Backtesting on test data

# Backtesting on test set
# ==============================================================================
cv = TimeSeriesFold(
    initial_train_size = end_validation,
    steps              = 7,
    refit              = False
)

# The validation partition is now used as part of the initial fit
# Epocs is set to the value idenfied in the previous training with early stopping
forecaster.set_fit_kwargs({"epochs": 10, "batch_size": 128})

metrics, predictions = backtesting_forecaster_multiseries(
                          forecaster        = forecaster,
                          series            = series[levels],
                          exog              = exog_features,
                          cv                = cv,
                          levels            = levels,
                          metric            = 'mean_absolute_error',
                          suppress_warnings = False,
                      )

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

Using 'torch' backend with device: cuda

Epoch 1/10
5/5 ━━━━━━━━━━━━━━━━━━━━ 4s 689ms/step - loss: 0.0408
Epoch 2/10
5/5 ━━━━━━━━━━━━━━━━━━━━ 4s 709ms/step - loss: 0.0390
Epoch 3/10
5/5 ━━━━━━━━━━━━━━━━━━━━ 3s 645ms/step - loss: 0.0383
Epoch 4/10
5/5 ━━━━━━━━━━━━━━━━━━━━ 4s 739ms/step - loss: 0.0374
Epoch 5/10
5/5 ━━━━━━━━━━━━━━━━━━━━ 3s 608ms/step - loss: 0.0367
Epoch 6/10
5/5 ━━━━━━━━━━━━━━━━━━━━ 3s 597ms/step - loss: 0.0362
Epoch 7/10
5/5 ━━━━━━━━━━━━━━━━━━━━ 3s 647ms/step - loss: 0.0359
Epoch 8/10
5/5 ━━━━━━━━━━━━━━━━━━━━ 3s 578ms/step - loss: 0.0354
Epoch 9/10
5/5 ━━━━━━━━━━━━━━━━━━━━ 3s 728ms/step - loss: 0.0353
Epoch 10/10
5/5 ━━━━━━━━━━━━━━━━━━━━ 3s 638ms/step - loss: 0.0351

  0%|          | 0/9 [00:00<?, ?it/s]

# 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()
sum_bias_all_buildings = errors_all_buildings['error'].sum().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: 730
Sum of absolute errors for all buildings (x 10,000): 7027
Bias (x 10,000): -589

# Predictions
# ==============================================================================
predictions

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

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

for i, building in enumerate(selected_buildings):
    series.loc[predictions.index, building].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();

Session information

import session_info
session_info.show(html=False)

-----
feature_engine      1.9.3
keras               3.12.0
matplotlib          3.10.8
numpy               2.3.4
pandas              2.3.3
session_info        v1.0.1
skforecast          0.19.0
sklearn             1.7.2
torch               2.8.0
-----
IPython             9.7.0
jupyter_client      8.6.3
jupyter_core        5.9.1
-----
Python 3.13.9 | packaged by conda-forge | (main, Oct 22 2025, 23:12:41) [MSC v.1944 64 bit (AMD64)]
Windows-11-10.0.26100-SP0
-----
Session information updated at 2025-12-01 11:28

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 with deep learning by Joaquín Amat Rodrigo, Javier Escobar Ortiz and Fernando Carazo, available under a CC BY-NC-SA 4.0 at https://www.cienciadedatos.net/documentos/py68-forecasting-at-scale-deep-learning.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. (2025). skforecast (v0.19.0). Zenodo. https://doi.org/10.5281/zenodo.8382788

APA:

Amat Rodrigo, J., & Escobar Ortiz, J. (2025). skforecast (Version 0.19.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.19.0}, month = {11}, year = {2025}, 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.