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
• Global Forecasting Models: Multi-series forecasting
• Global Forecasting Models: Comparative Analysis of Single and Multi-Series Forecasting Modeling
• 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¶

Global forecasting modeling involves the creation of a single forecasting model that considers all time series simultaneously. This approach attempts to capture the underlying patterns common to all series, thereby reducing the impact of noise present in individual series. It offers computational efficiency, ease of maintenance and robust generalization across time series. Global forecasting assumes that time series with similar behavior can benefit from being modeled together. When working with hundreds or thousands of series, an initial clustering analysis can help maximize model performance by identifying groups of series that behave similarly. The clustering methods used in this document, based on time series features and on elastic distance measures, are described in detail in the section Clustering time series.

If after performing the clustering analysis, distinct groups of series are identified, two strategies can be used:

• Model all series together, while adding a feature that indicates the cluster to which each series belongs.

• Build several global models, with each model tailored to a specific group of series.

To help understand the benefits of clustering, this document focuses on examining and comparing the results obtained when predicting the energy consumption of over a thousand buildings (a random subset of 600 buildings is used in the modeling section to keep run times reasonable). Although the primary use of each building is available in the dataset, it may not reflect groups with similar patterns of energy use, so additional groups are created using clustering methods. A total of 4 experiments are performed:

• Modeling all buildings together with a unique global model strategy.

• Modeling groups of buildings based on their primary use (a global model per primary use).

• Modeling groups of buildings based on time series feature clustering (a global model per cluster).

• Modeling groups of buildings based on Dynamic Time Warping (DTW) clustering (a global model per cluster).

Libraries¶

Libraries used in this document.

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

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

# Forecasting
# ==============================================================================
import lightgbm
from lightgbm import LGBMRegressor
import skforecast
from skforecast.recursive import ForecasterRecursiveMultiSeries
from skforecast.model_selection import (
    TimeSeriesFold,
    backtesting_forecaster_multiseries
)

from skforecast.preprocessing import (
    CalendarFeatures,
    reshape_series_long_to_dict,
    reshape_exog_long_to_dict
)

# Feature extraction
# ==============================================================================
import tsfresh
from tsfresh import extract_features
from tsfresh.feature_extraction.settings import ComprehensiveFCParameters, from_columns
from tsfresh.utilities.dataframe_functions import impute

# Clustering
# ==============================================================================
import sklearn
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
from sklearn.cluster import KMeans
import sktime
from sktime.clustering.k_means import TimeSeriesKMeans

# 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 tsfresh: {tsfresh.__version__}")
print(f"{color}Version sktime: {sktime.__version__}")

Version skforecast: 0.25.0
Version scikit-learn: 1.7.2
Version lightgbm: 4.7.0
Version tsfresh: 0.21.2
Version sktime: 1.1.0

Data¶

Data used in this document has been obtained from the Kaggle competition Addison Howard, Chris Balbach, Clayton Miller, Jeff Haberl, Krishnan Gowri, Sohier Dane. (2019). ASHRAE - Great Energy Predictor III. Kaggle.

Three files are used to create the modeling data set:

• weather_train.csv and weather_test.csv: These files contain weather-related data for each building, including outdoor air temperature, dew point temperature, relative humidity, and other weather parameters. The weather data is crucial for understanding the impact of external conditions on building energy usage.

• building_metadata.csv: This file provides metadata for each building in the dataset, such as building type, primary use, square footage, floor count, and year built. This information helps in understanding the characteristics of the buildings and their potential influence on energy consumption patterns.

• train.csv: the train dataset contains the target variable, i.e., the energy consumption data for each building, along with the timestamps for the energy consumption readings. It also includes the corresponding building and weather IDs to link the information across different datasets.

The three files have been preprocessed to: remove buildings in which fewer than 85% of the values are valid (neither NaN nor zero), keep only the electricity meter, and aggregate the data to a daily frequency.

# Load preprocessed data
# ==============================================================================
data = fetch_dataset('ashrae_daily')
data.head()

╭────────────────────────────────── ashrae_daily ──────────────────────────────────╮
│ Description:                                                                     │
│ Daily energy consumption data from the ASHRAE competition with building metadata │
│ and weather data.                                                                │
│                                                                                  │
│ Source:                                                                          │
│ Kaggle competition Addison Howard, Chris Balbach, Clayton Miller, Jeff Haberl,   │
│ Krishnan Gowri, Sohier Dane. (2019). ASHRAE - Great Energy Predictor III.        │
│ Kaggle. https://www.kaggle.com/c/ashrae-energy-prediction/overview               │
│                                                                                  │
│ URL:                                                                             │
│ https://huggingface.co/datasets/skforecast/ashrae_daily/resolve/main/ashrae_dail │
│ y.parquet                                                                        │
│                                                                                  │
│ Shape: 444324 rows x 10 columns                                                  │
╰──────────────────────────────────────────────────────────────────────────────────╯

# Ensure all time series indexes are complete without gaps
# ==============================================================================
data = (
    data
    .groupby('building_id')
    .apply(lambda group: group.asfreq('D'), include_groups=False)
    .reset_index()
    .set_index('timestamp')
)

# Fill missing values of air_temperature and wind_speed using forward and backward fill
# ==============================================================================
# Imputation must be done separately for each building
data = data.sort_values(by=['building_id', 'timestamp'])
cols_to_impute = ['air_temperature', 'wind_speed']
data[cols_to_impute] = (
    data
    .groupby('building_id')[cols_to_impute]
    .transform(lambda x: x.ffill().bfill())
)
data = data.sort_index()

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 --- 2016-12-31 00:00:00  (n_days=365)

Exploratory data analysis¶

Building primary use¶

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 and type of buildings based on primary use
# ==============================================================================
n_building = data['building_id'].nunique()
n_type_building = data['primary_use'].nunique()

print(f"Number of buildings: {n_building}")
print(f"Number of building types: {n_type_building}")
display(data.drop_duplicates(subset=['building_id'])['primary_use'].value_counts())

Number of buildings: 1214
Number of building types: 16

primary_use
Education                        463
Office                           228
Entertainment/public assembly    157
Public services                  150
Lodging/residential              112
Other                             19
Healthcare                        19
Parking                           13
Warehouse/storage                 12
Manufacturing/industrial          10
Services                           9
Food sales and service             5
Technology/science                 5
Retail                             5
Utility                            4
Religious worship                  3
Name: count, dtype: int64

For certain primary use categories, there is a limited number of buildings within the dataset. To streamline the analysis, categories with fewer than 50 buildings are grouped into the "Other" category.

# Types of buildings (primary use) with fewer than 50 buildings are grouped as "Other"
# ==============================================================================
infrequent_categories = (
    data
    .drop_duplicates(subset=['building_id'])['primary_use']
    .value_counts()
    .loc[lambda x: x < 50]
    .index
    .tolist()
)
print("Infrequent categories:")
print("======================")
print('\n'.join(infrequent_categories))

data['primary_use'] = np.where(
    data['primary_use'].isin(infrequent_categories),
    'Other',
    data['primary_use']
)

Infrequent categories:
======================
Other
Healthcare
Parking
Warehouse/storage
Manufacturing/industrial
Services
Food sales and service
Technology/science
Retail
Utility
Religious worship

Next, two plots are created: the first one shows the energy consumption of a randomly selected building within each category, and the second one shows all the available time series of each category (one gray line per building), with the average consumption of the category shown in blue.

# Time series for one randomly selected building per group
# ==============================================================================
set_dark_theme()
fig, axs = plt.subplots(nrows=3, ncols=2, figsize=(8, 5.5), sharex=True, sharey=False)
sample_ids = (
    data
    .groupby('primary_use')['building_id']
    .apply(lambda x: x.sample(1, random_state=333))
    .tolist()
)
axs = axs.flatten()

for i, building_id in enumerate(sample_ids):
    data_sample = data[data['building_id'] == building_id]
    building_type = data_sample['primary_use'].unique()[0]
    data_sample.plot(
        y        = 'meter_reading',
        ax       = axs[i],
        legend   = False,
        title    = f"Building: {building_id}, type: {building_type}",
        fontsize = 8
    )
    axs[i].set_xlabel("")
    axs[i].set_ylabel("")
    # Scientific notation for y axis
    axs[i].ticklabel_format(axis='y', style='sci', scilimits=(0, 0))
    axs[i].title.set_size(9)

fig.suptitle('Energy consumption for 6 random buildings', fontsize=12)
fig.tight_layout()
plt.show()

# Energy consumption by type of building (one gray line per building)
# ==============================================================================
fig, axs = plt.subplots(nrows=3, ncols=2, figsize=(8, 5.5), sharex=True, sharey=False)
axs = axs.flatten()

for i, building_type in enumerate(data['primary_use'].unique()):
    data_sample = data[data['primary_use'] == building_type]
    data_sample = data_sample.pivot_table(
                      index   = 'timestamp',
                      columns = 'building_id',
                      values  = 'meter_reading',
                      aggfunc = 'mean'
                  )
    data_sample.plot(
        legend   = False,
        title    = f"Type: {building_type}",
        color    = 'gray',
        alpha    = 0.2,
        fontsize = 8,
        ax       = axs[i]
    )
    mean_consumption = data_sample.mean(axis=1)
    mean_consumption.plot(
        ax       = axs[i],
        color    = 'blue',
        fontsize = 8
    )
    axs[i].set_xlabel("")
    axs[i].set_ylabel("")
    # Scientific notation for y axis
    axs[i].ticklabel_format(axis='y', style='sci', scilimits=(0, 0))
    axs[i].title.set_size(9)

    # Limit the axis to 5 times the maximum mean value to improve visualization
    axs[i].set_ylim(
        bottom = 0,
        top    = 5 * mean_consumption.max()
    )

fig.tight_layout()
fig.suptitle('Energy consumption by type of building', fontsize=12)
plt.subplots_adjust(top=0.9)
plt.show()

The graph reveals that there is significant variability in consumption patterns among buildings of the same purpose. This suggests that there may be scope for improving the criteria by which buildings are grouped.

# Sample 600 buildings
# ==============================================================================
rng = np.random.default_rng(12345)
buildings = data['building_id'].unique()
buildings_selected = rng.choice(
                         buildings,
                         size    = 600,
                         replace = False
                     )

data = data.query("building_id in @buildings_selected")

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 more detailed review of time series clustering, see A review and evaluation of elastic distance functions for time series clustering.

Clustering based on time series features¶

Feature creation¶

tsfresh is a powerful Python library for feature engineering from time-series and sequential data, including statistical measures, Fourier coefficients, and various other time-domain and frequency-domain features. It provides a systematic approach to automate the calculation of features and select the most informative ones.

To start, the default configuration of tsfresh is used, which calculates all available features. The easiest way to access this configuration is to instantiate the ComprehensiveFCParameters class, which returns a dictionary that maps the name of each feature (a string) to a list of dictionaries with the parameters used when that feature is calculated (or None if the feature has no parameters).

# Default features and settings created by tsfresh
# ==============================================================================
default_features = ComprehensiveFCParameters()
print("Name of the features extracted by tsfresh:")
print("=========================================")
print('\n'.join(default_features))

Name of the features extracted by tsfresh:
=========================================
variance_larger_than_standard_deviation
has_duplicate_max
has_duplicate_min
has_duplicate
sum_values
abs_energy
mean_abs_change
mean_change
mean_second_derivative_central
median
mean
length
standard_deviation
variation_coefficient
variance
skewness
kurtosis
root_mean_square
absolute_sum_of_changes
longest_strike_below_mean
longest_strike_above_mean
count_above_mean
count_below_mean
last_location_of_maximum
first_location_of_maximum
last_location_of_minimum
first_location_of_minimum
percentage_of_reoccurring_values_to_all_values
percentage_of_reoccurring_datapoints_to_all_datapoints
sum_of_reoccurring_values
sum_of_reoccurring_data_points
ratio_value_number_to_time_series_length
sample_entropy
maximum
absolute_maximum
minimum
benford_correlation
time_reversal_asymmetry_statistic
c3
cid_ce
symmetry_looking
large_standard_deviation
quantile
autocorrelation
agg_autocorrelation
partial_autocorrelation
number_cwt_peaks
number_peaks
binned_entropy
index_mass_quantile
cwt_coefficients
spkt_welch_density
ar_coefficient
change_quantiles
fft_coefficient
fft_aggregated
value_count
range_count
approximate_entropy
friedrich_coefficients
max_langevin_fixed_point
linear_trend
agg_linear_trend
augmented_dickey_fuller
number_crossing_m
energy_ratio_by_chunks
ratio_beyond_r_sigma
linear_trend_timewise
count_above
count_below
lempel_ziv_complexity
fourier_entropy
permutation_entropy
query_similarity_count
mean_n_absolute_max

Many of these features are calculated using different values for their arguments.

# Default configuration for feature "partial_autocorrelation"
# ==============================================================================
default_features['partial_autocorrelation']

[{'lag': 0},
 {'lag': 1},
 {'lag': 2},
 {'lag': 3},
 {'lag': 4},
 {'lag': 5},
 {'lag': 6},
 {'lag': 7},
 {'lag': 8},
 {'lag': 9}]

To access the detailed view of each feature and the parameter values included in the default configuration, use the following code:

# Default configuration for all extracted features
# ==============================================================================
# from pprint import pprint
# pprint(default_features)

Once the configuration of features has been defined, the next step is to extract them from the time series. The function extract_features() of tsfresh is used for this purpose. This function receives as input the time series and the configuration of the features to be extracted. The output is a dataframe with the extracted features.

# Feature extraction
# ==============================================================================
ts_features = extract_features(
    timeseries_container  = data[['building_id', 'meter_reading']].reset_index(),
    column_id             = "building_id",
    column_sort           = "timestamp",
    column_value          = "meter_reading",
    default_fc_parameters = default_features,
    impute_function       = impute,
    n_jobs                = 4
)

print("Shape of ts_features:", ts_features.shape)
ts_features.head(3)

Shape of ts_features: (600, 783)

As a result of the extraction process, 783 features have been created for each time series (building_id in this use case). The returned dataframe has as its index the column specified in the column_id argument of extract_features.

The default extraction of tsfresh produces a huge number of features. However, only a few of these may be of interest in each use case. To select the most relevant ones, tsfresh includes an automated selection process based on hypothesis tests, known as FRESH (FeatuRe Extraction based on Scalable Hypothesis tests). In this process, the features are individually and independently evaluated for their significance in predicting the target under investigation.

In this case, since the goal is to identify groups of buildings with similar energy consumption patterns, regardless of the building type, the automated selection process is not used. Instead, the third strategy described in the warning above is applied: the dimensionality of the feature matrix is reduced with PCA before clustering.

Some features cannot be calculated for certain time series and return missing (NaN) or infinite values. Since most clustering algorithms do not allow missing values, these must be handled. In this case, they have already been replaced during the extraction by the impute function passed to the impute_function argument of extract_features(): -inf and +inf values are replaced by the minimum and maximum of each feature, and NaN values by its median. Therefore, the next step, which removes any feature with missing values, acts only as a safeguard and no feature is removed.

# Remove features with missing values
# ==============================================================================
ts_features = ts_features.dropna(axis=1, how='any')
print(f"Number of features after removing those with missing values: {ts_features.shape[1]}")

Number of features after removing those with missing values: 783

Once the final feature matrix is ready, it may be useful to create a new dictionary to store the final features and the parameters used to calculate them. This can be easily done using the from_columns function.

# Dictionary with selected features and their configuration
# ==============================================================================
ts_features_info = from_columns(ts_features)
# pprint(ts_features_info['meter_reading'])

K-means clustering¶

The K-means clustering method is used to group the buildings. Since clustering is known to be negatively affected by high dimensionality, and since several hundred features have been created for each building, PCA is used to reduce the dimensionality of the data before K-means is applied. Because PCA is sensitive to the scale of the variables, the features are first standardized so that all of them have mean 0 and standard deviation 1.

# Scale features so all of them have mean 0 and std 1
# ==============================================================================
scaler = StandardScaler().set_output(transform="pandas")
ts_features_scaled = scaler.fit_transform(ts_features)
ts_features_scaled.head(2)

# Plot variance reduction as a function of the number of PCA components
# ==============================================================================
pca = PCA()
pca.fit(ts_features_scaled)
cumulative_variance = np.cumsum(pca.explained_variance_ratio_)
n_components = np.argmax(cumulative_variance > 0.85) + 1

fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(6, 2.5))
ax.plot(np.arange(1, len(cumulative_variance) + 1), cumulative_variance)
ax.set_title('Variance explained by PCA components')
ax.set_xlabel('Number of components')
ax.set_ylabel('Cumulative explained variance')
ax.axhline(y=0.85, color='red', linestyle='--')
ax.axvline(x=n_components, color='red', linestyle='--')
ax.text(
    x=n_components + 10,
    y=0.4,
    s=f"n_components = {n_components}",
    color='red',
    verticalalignment='center'
)

plt.show();

The dimensionality of the original 783 features is reduced by using the first 57 principal components (explaining 85% of the variance).

# PCA with as many components as necessary to explain 85% of the variance
# ==============================================================================
pca = PCA(n_components=0.85)
pca_projections = pca.fit_transform(ts_features_scaled)

# Create a data frame with the projections. Each column is a principal component
# and each row is the id of the building.
pca_projections = pd.DataFrame(
    pca_projections,
    index   = ts_features_scaled.index,
    columns = [f"PC{i}" for i in range(1, pca.n_components_ + 1)]
)

print(f"Number of components selected: {pca.n_components_}")
print(f"Explained variance: {pca.explained_variance_ratio_.sum():.3f}")
pca_projections.head(3)

Number of components selected: 57
Explained variance: 0.852

One of the inherent challenges of clustering is determining the optimal number of clusters, since there is no ground truth or predefined number of clusters in unsupervised learning. Several heuristics have been proposed to guide this selection, and in this case, the elbow method is used.

The elbow method involves plotting the total within-cluster sum of squares (WSS) against the number of clusters (k). The optimal number of clusters is typically identified at the "elbow" of the curve, where the curve flattens: beyond this point, adding more clusters no longer reduces the WSS substantially. In scikit-learn, the WSS of a fitted KMeans model is stored in its inertia_ attribute.

# Optimal number of clusters
# ==============================================================================
# Identify the optimal number of clusters using the elbow method
range_n_clusters = range(1, 30)
inertias = []
size_of_clusters = []

for n_clusters in range_n_clusters:
    kmeans = KMeans(
                 n_clusters   = n_clusters, 
                 n_init       = 20, 
                 random_state = 963852
             )
    kmeans.fit(pca_projections)
    inertias.append(kmeans.inertia_)
    size_of_clusters.append(np.unique(kmeans.labels_, return_counts=True)[1].tolist())

inertias = pd.Series(inertias, index=range_n_clusters)
decrease_in_inertia = inertias.pct_change() * -100

fig, axs = plt.subplots(1, 2, figsize=(9, 2.5))
inertias.plot(marker='o', ax=axs[0])
axs[0].xaxis.set_major_locator(ticker.MaxNLocator(integer=True))
axs[0].set_title("Cluster inertia vs number of clusters")
axs[0].set_xlabel('Number of clusters')
axs[0].set_ylabel('Within-cluster sum of squares (inertia)')

decrease_in_inertia.plot(kind='bar', ax=axs[1])
axs[1].set_title("Decrease in inertia")
axs[1].set_xlabel('Number of clusters')
axs[1].set_ylabel('Decrease in inertia (%)')
fig.tight_layout()
plt.show();

The chart shows that the decrease in inertia slows down noticeably after 5 clusters, so values in that region are reasonable candidates.

However, the elbow criterion should not be the only guide. In addition to analyzing the evolution of the inertia (intra-variance), it is important to check the size of the clusters being created. The presence of small clusters may indicate overfitting or anomalous samples that do not fit well into any of the groups.

# Distribution of cluster sizes
# ==============================================================================
for n_clusters, sizes in zip(range_n_clusters, size_of_clusters):
    print(f"Size of clusters (n = {n_clusters}): {sizes}")

Size of clusters (n = 1): [600]
Size of clusters (n = 2): [578, 22]
Size of clusters (n = 3): [22, 575, 3]
Size of clusters (n = 4): [31, 562, 3, 4]
Size of clusters (n = 5): [481, 3, 95, 3, 18]
Size of clusters (n = 6): [293, 15, 1, 64, 3, 224]
Size of clusters (n = 7): [277, 78, 221, 3, 3, 17, 1]
Size of clusters (n = 8): [60, 207, 3, 1, 30, 4, 12, 283]
Size of clusters (n = 9): [194, 3, 47, 298, 3, 1, 1, 7, 46]
Size of clusters (n = 10): [42, 188, 1, 3, 4, 200, 10, 1, 45, 106]
Size of clusters (n = 11): [98, 193, 185, 1, 2, 16, 45, 1, 1, 3, 55]
Size of clusters (n = 12): [45, 264, 1, 3, 1, 73, 1, 194, 1, 1, 15, 1]
Size of clusters (n = 13): [46, 280, 1, 1, 67, 1, 3, 184, 1, 5, 1, 9, 1]
Size of clusters (n = 14): [56, 5, 1, 1, 100, 1, 44, 189, 1, 4, 1, 2, 4, 191]
Size of clusters (n = 15): [187, 1, 44, 12, 190, 2, 1, 1, 2, 98, 3, 1, 3, 54, 1]
Size of clusters (n = 16): [93, 4, 1, 25, 1, 3, 1, 235, 185, 1, 45, 1, 1, 2, 1, 1]
Size of clusters (n = 17): [200, 2, 1, 3, 8, 45, 1, 188, 76, 1, 1, 1, 1, 1, 3, 1, 67]
Size of clusters (n = 18): [46, 1, 187, 1, 1, 4, 168, 105, 30, 1, 44, 1, 1, 3, 1, 2, 3, 1]
Size of clusters (n = 19): [90, 6, 1, 1, 45, 134, 2, 49, 1, 1, 150, 1, 1, 3, 1, 1, 1, 81, 31]
Size of clusters (n = 20): [104, 1, 30, 1, 137, 1, 173, 2, 2, 1, 1, 10, 1, 1, 96, 1, 1, 33, 1, 3]
Size of clusters (n = 21): [118, 2, 1, 32, 153, 1, 1, 1, 46, 3, 1, 10, 5, 1, 96, 1, 124, 1, 1, 1, 1]
Size of clusters (n = 22): [114, 1, 1, 38, 1, 1, 113, 24, 103, 46, 2, 1, 1, 3, 1, 1, 1, 1, 1, 113, 2, 31]
Size of clusters (n = 23): [10, 115, 109, 1, 1, 1, 1, 4, 1, 41, 1, 33, 1, 1, 76, 3, 1, 1, 40, 3, 1, 154, 1]
Size of clusters (n = 24): [95, 112, 1, 1, 45, 1, 1, 4, 1, 31, 1, 2, 10, 1, 124, 1, 37, 31, 1, 1, 2, 1, 1, 95]
Size of clusters (n = 25): [9, 128, 1, 93, 1, 1, 45, 176, 10, 2, 1, 1, 1, 3, 40, 1, 1, 32, 1, 1, 1, 2, 1, 1, 47]
Size of clusters (n = 26): [176, 1, 1, 10, 112, 2, 1, 66, 1, 1, 1, 1, 9, 1, 1, 1, 1, 1, 41, 1, 103, 33, 32, 1, 1, 1]
Size of clusters (n = 27): [38, 1, 9, 70, 145, 1, 1, 10, 1, 1, 1, 2, 1, 42, 1, 2, 10, 84, 1, 55, 86, 1, 1, 1, 1, 2, 32]
Size of clusters (n = 28): [40, 96, 31, 1, 1, 4, 2, 1, 1, 41, 105, 1, 1, 93, 114, 1, 1, 46, 1, 1, 1, 1, 3, 9, 1, 1, 1, 1]
Size of clusters (n = 29): [64, 1, 24, 110, 1, 1, 1, 2, 30, 3, 31, 1, 1, 1, 27, 62, 1, 1, 9, 81, 1, 1, 1, 1, 1, 1, 40, 1, 101]

The sizes confirm that the elbow value alone is not enough. With 5 clusters, one group contains 481 of the 600 buildings and three of the remaining groups have 18 buildings or fewer, so in practice almost all the series would still be modeled together. With 7 clusters the partition is more informative: three groups of a usable size (277, 221 and 78 buildings) plus several very small ones.

For this reason, 7 clusters are used, and the smallest ones (fewer than 20 buildings) are combined into a single group named "Other". This avoids training a global model with only a handful of series, which would defeat the purpose of the multi-series approach.

# Train clustering model with 7 clusters and assign each building to a cluster
# ==============================================================================
kmeans = KMeans(n_clusters=7, n_init=20, random_state=963852)
cluster_labels = kmeans.fit_predict(pca_projections)
clusters = pd.DataFrame({
    'building_id': pca_projections.index,
    'cluster_based_on_features': cluster_labels.astype(str)
})

# Merge minor clusters into a single cluster
min_cluster_size = 20
cluster_size = clusters['cluster_based_on_features'].value_counts()
small_clusters = cluster_size[cluster_size < min_cluster_size].index.tolist()
clusters['cluster_based_on_features'] = np.where(
    clusters['cluster_based_on_features'].isin(small_clusters),
    'Other',
    clusters['cluster_based_on_features']
)
display(clusters.head(3))
clusters['cluster_based_on_features'].value_counts().to_frame(name='Number of series')

# Add the predicted cluster to the data frame with the building information
# ==============================================================================
data = pd.merge(
           data.reset_index(),  # To avoid losing the index
           clusters,
           on       = 'building_id',
           how      = 'left',
           validate = 'm:1'
       ).set_index('timestamp')

data.head(3)

Once each building has been assigned to a cluster, it is useful to examine the characteristics of the grouped buildings. For example, the distribution of the primary use of the buildings.

# Percentage of building types per cluster
# ==============================================================================
primary_usage_per_cluster = (
    data
    .drop_duplicates(subset='building_id')
    .groupby('cluster_based_on_features')['primary_use']
    .value_counts(normalize=True)
    .unstack()
    .fillna(0)
)

primary_usage_per_cluster = (100 * primary_usage_per_cluster).round(1)
primary_usage_per_cluster

The results suggest (the table must be read row by row, since the percentages of each row add up to 100%) that the clustering process based on features extracted from the time series generates groups that differ from those formed by the main purpose of the building. No cluster is dominated by a single type of building: for example, educational buildings account for between 25% and 55% of every cluster, and entertainment/public assembly buildings for between 12% and 15%.

Clustering based on elastic distance (DTW)¶

Dynamic Time Warping (DTW) is a technique that measures the similarity between two temporal sequences, which may vary in speed. In essence, it's an elastic distance measure that allows the time series to be stretched or compressed to align with each other optimally. Clustering with DTW involves grouping time series data based on their dynamic time-warped distances, ensuring that time series within the same cluster have similar shapes and patterns, even if they are out of phase in time or have different lengths.

The TimeSeriesKMeans class from the sktime library enables the application of K-means clustering with a variety of distance metrics, including DTW, Euclidean, ERP, EDR, LCSS, squared, DDTW, WDTW, and WDDTW. Many of these metrics are elastic distances, making the method well-suited for time series data.

Unlike the previous section, the number of clusters is not selected here with the elbow method: computing DTW distances between hundreds of series takes several minutes per fit, which makes it impractical to explore a wide range of values. Instead, 4 clusters are set a priori, the same number of groups finally used in the feature-based approach (3 clusters plus the "Other" group).

sktime requires time series to be structured in a long format with a multi-index. The outermost index level represents the series ID, while the innermost level corresponds to the datetime.

# Convert data to long format with a multiindex: (building_id, timestamp)
# ==============================================================================
data_long = (
    data
    .reset_index()
    .loc[:, ['timestamp', 'building_id', 'meter_reading']]
    .set_index(['building_id', 'timestamp'])
    .sort_index(ascending=True)
)
data_long

# Z-normalize each series (mean 0, standard deviation 1) so that DTW compares
# the shape of the consumption patterns rather than their level
# ==============================================================================
grouped = data_long.groupby(level='building_id')['meter_reading']
series_mean = grouped.transform('mean')
series_std = grouped.transform('std').replace(0, 1)  # Avoid division by zero in constant series
data_long_scaled = ((data_long['meter_reading'] - series_mean) / series_std).to_frame()
data_long_scaled

# Fit clustering model
# ==============================================================================
model = TimeSeriesKMeans(
            n_clusters   = 4,
            metric       = "dtw",
            max_iter     = 10,
            random_state = 123456
        )

model.fit(data_long_scaled)

TimeSeriesKMeans(max_iter=10, n_clusters=4, random_state=123456)

TimeSeriesKMeans(max_iter=10, n_clusters=4, random_state=123456)

# Cluster prediction
# ==============================================================================
cluster_labels = model.predict(data_long_scaled)
clusters = pd.DataFrame({
               'building_id': data_long_scaled.index.get_level_values('building_id').drop_duplicates(),
               'cluster_based_on_dtw': cluster_labels.astype(str),
           })

# Size of each cluster
clusters['cluster_based_on_dtw'].value_counts().sort_index().to_frame(name='Number of series')

In this case, all the clusters contain more than 60 buildings, so, unlike in the feature-based approach, it is not necessary to merge small clusters.

# Add the predicted cluster to the data frame with the building information
# ==============================================================================
data = pd.merge(
           data.reset_index(),  # To avoid losing the index
           clusters,
           on       = 'building_id',
           how      = 'left',
           validate = 'm:1'
       ).set_index('timestamp')

data.head(3)

# Save data for modeling to avoid repeating the previous steps
# ==============================================================================
data.to_parquet('data_modelling.parquet')

Modeling and forecasting¶

After establishing the three grouping criteria (building primary use, clustering based on time series features, and clustering based on dynamic time warping), a multi-series global model is trained for each group, and the results are compared with those of a single global model that includes all the buildings. The evaluation that follows focuses on determining how effectively these models can forecast daily consumption during the last five months of the year (August to December), which are predicted in successive 7-day horizons. During this assessment, three distinct metrics are used:

• The average value of the Mean Absolute Error (MAE) across all buildings.

• The sum of the absolute errors, i.e., the sum of the absolute differences between the predicted values and the actual consumption over all buildings and dates.

• The bias, calculated as the sum of the errors (predicted value minus actual consumption) over all buildings and dates. Positive values indicate that the model overestimates consumption, and negative values that it underestimates it.

In addition to the lagged values of each time series, the model includes the day of the week (sine-cosine encoded), created by the forecaster through the calendar_features argument, and the following exogenous variables: outdoor temperature (air_temperature), wind speed (wind_speed) and primary use of the building (primary_use). Since primary_use is a text column, the forecaster automatically treats it as a categorical feature (categorical_features='auto').

To train the models and evaluate their predictive performance, the data is divided into two sets: training (January to July) and test (August to December). No hyperparameter search is carried out in this document; the same model configuration is used in all the experiments so that the differences in the results depend only on how the buildings are grouped.

# Load data for modeling
# ==============================================================================
data = pd.read_parquet('data_modelling.parquet')

The data is reshaped from a long data frame into a dictionary of series using reshape_series_long_to_dict and reshape_exog_long_to_dict. Although this is not strictly required, since skforecast allows multiple input formats, it is the recommended format for multi-series models, since it allows series of different lengths and with different subsets of exogenous features to be combined easily.

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

# The exogenous variables included in all the models are: primary_use,
# air_temperature and wind_speed
exog_dict = reshape_exog_long_to_dict(
    data      = data[['building_id', 'primary_use', 'air_temperature', 'wind_speed']].reset_index(),
    series_id = 'building_id',
    index     = 'timestamp',
    freq      = 'D'
)

# End of the training period and test series
# ==============================================================================
# The train/test split used for modeling is handled internally by TimeSeriesFold
# (see `cv` below), which receives the full, undivided series_dict/exog_dict.
# The test series are extracted only to plot them against the predictions.
end_train = '2016-07-31 23:59:00'
series_dict_test = {k: v.loc[end_train:] for k, v in series_dict.items()}

The same forecaster configuration is used in all the experiments, so that the only difference between them is the way the buildings are grouped:

• Regressor: a LGBMRegressor with 500 trees, a maximum depth of 10 and a small learning rate (0.01).

• Lags: the last 31 values of each series, enough for the model to capture both weekly and monthly patterns.

• Calendar features: the day of the week, cyclically encoded with CalendarFeatures. The sine-cosine (cyclical) encoding ensures that Sunday and Monday are as close to each other as any other pair of consecutive days, something that a plain numeric encoding (0 to 6) would not achieve.

• Series encoding: encoding="ordinal_category" adds the series identifier (the building id) as a categorical feature, which allows a single global model to learn the specific behavior of each building. The available strategies are described in the Global Forecasting Models user guide.

# Forecaster definition (shared by all the experiments)
# ==============================================================================
calendar_features = CalendarFeatures(features=['day_of_week'], encoding='cyclical')

params_lgbm = {
    'n_estimators': 500,
    'learning_rate': 0.01,
    'max_depth': 10,
    'random_state': 8520,
    'verbose': -1
}
forecaster = ForecasterRecursiveMultiSeries(
                 estimator         = LGBMRegressor(**params_lgbm),
                 lags              = 31,
                 calendar_features = calendar_features,
                 encoding          = "ordinal_category"
             )

The predictive performance is estimated with backtesting, using the class TimeSeriesFold to define the validation scheme and the function backtesting_forecaster_multiseries to run it:

• initial_train_size = end_train: the model is trained with the data up to 2016-07-31, and the remaining five months (August to December) are used as the test set.

• steps = 7: predictions are generated in successive horizons of 7 days (one week), which results in 22 folds.

• refit = False: the model is trained only once, at the beginning of the process, and it is not retrained as the backtesting window advances. This keeps the comparison of the four strategies affordable in computational terms.

Since the training data always precede the predicted values, this scheme respects the temporal order of the data and no information from the future leaks into the model.

# Backtesting definition
# ==============================================================================
cv = TimeSeriesFold(
        steps              = 7,
        initial_train_size = end_train,
        refit              = False
     )

# Table of results for all models
# ==============================================================================
table_results = pd.DataFrame(columns=['model', 'mae', 'abs_error', 'bias', 'elapsed_time'])
table_results = table_results.set_index('model')
table_results = table_results.astype({'mae': float, 'abs_error': float, 'bias': float, 'elapsed_time': object})


def add_results_to_table(table_results, model_name, metrics, predictions, data, elapsed_time):
    """
    Calculate the average MAE, the sum of absolute errors and the bias of the
    backtesting predictions, store them in `table_results` and print them.
    """
    # Combine predictions with real values
    results = predictions.reset_index(names=['timestamp']).merge(
        data[['building_id', 'meter_reading']].reset_index(names=['timestamp']),
        left_on  = ['timestamp', 'level'],
        right_on = ['timestamp', 'building_id'],
        how      = 'inner'
    )
    errors = results['pred'] - results['meter_reading']
    mae_mean = metrics['mean_absolute_error'].mean()
    sum_abs_errors = errors.abs().sum()
    bias = errors.sum()

    table_results.loc[model_name, ['mae', 'abs_error', 'bias', 'elapsed_time']] = [
        mae_mean,
        sum_abs_errors,
        bias,
        elapsed_time
    ]

    print(
        f"\nAverage mean absolute error for all buildings: {mae_mean:.0f}\n"
        f"Sum of absolute errors for all buildings (x 10,000): {sum_abs_errors / 10000:.0f}\n"
        f"Bias (x 10,000): {bias / 10000:.0f}"
    )

Global multi-series model for all buildings¶

A global model for all buildings is trained and tested using the skforecast class ForecasterRecursiveMultiSeries. This forecaster accepts the time series in several formats; here, the dictionaries of series and exogenous variables created above are used.

For more information about using series of different lengths, or different exogenous variables for each series, see Series with different lengths and different exogenous variables.

# Forecaster multi-series to model all buildings at the same time
# ==============================================================================
start = datetime.now()

metrics, predictions = backtesting_forecaster_multiseries(
                           forecaster            = forecaster,
                           series                = series_dict,
                           exog                  = exog_dict,
                           cv                    = cv,
                           metric                = 'mean_absolute_error',
                           add_aggregated_metric = False,
                           verbose               = False,
                           show_progress         = True,
                           suppress_warnings     = True
                       )

end = datetime.now()

add_results_to_table(
    table_results = table_results,
    model_name    = 'Global model',
    metrics       = metrics,
    predictions   = predictions,
    data          = data,
    elapsed_time  = end - start
)

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

Average mean absolute error for all buildings: 456
Sum of absolute errors for all buildings (x 10,000): 4182
Bias (x 10,000): 881

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

fig, axs = plt.subplots(2, 1, figsize=(6, 4), sharex=True)

for i, building in enumerate(selected_buildings):
    series_dict_test[building].plot(ax=axs[i], label='Real value')
    predictions.query("level==@building")['pred'].plot(ax=axs[i], label='Global model')
    axs[i].set_title(f"Building {building}", fontsize=10)
    axs[i].set_xlabel("")
    axs[i].legend()

fig.tight_layout()
plt.show();

Global multi-series model by primary use¶

A global model for each primary use of the buildings is trained and tested.

# Forecaster multi-series models for buildings grouped by primary use
# ==============================================================================
predictions_all_buildings = []
metrics_all_buildings = []
start = datetime.now()

for primary_use in data['primary_use'].unique():

    # Create a subset based on primary use
    building_id_subset = set(data.loc[data['primary_use'] == primary_use, 'building_id'].unique())
    series_dict_subset = {k: v for k, v in series_dict.items() if k in building_id_subset}
    exog_dict_subset   = {k: v for k, v in exog_dict.items() if k in building_id_subset}

    print(
        f"Training and testing model for primary use: {primary_use} "
        f"(n = {len(building_id_subset)})"
    )

    metrics, predictions = backtesting_forecaster_multiseries(
                               forecaster            = forecaster,
                               series                = series_dict_subset,
                               exog                  = exog_dict_subset,
                               cv                    = cv,
                               metric                = 'mean_absolute_error',
                               add_aggregated_metric = False,
                               verbose               = False,
                               show_progress         = True,
                               suppress_warnings     = True
                           )
    predictions_all_buildings.append(predictions)
    metrics_all_buildings.append(metrics)

end = datetime.now()

predictions_all_buildings = pd.concat(predictions_all_buildings, axis=0)
metrics_all_buildings = pd.concat(metrics_all_buildings, axis=0)

add_results_to_table(
    table_results = table_results,
    model_name    = 'Global model per primary use',
    metrics       = metrics_all_buildings,
    predictions   = predictions_all_buildings,
    data          = data,
    elapsed_time  = end - start
)

Training and testing model for primary use: Other (n = 62)

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

Training and testing model for primary use: Education (n = 229)

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

Training and testing model for primary use: Entertainment/public assembly (n = 82)

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

Training and testing model for primary use: Office (n = 108)

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

Training and testing model for primary use: Public services (n = 67)

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

Training and testing model for primary use: Lodging/residential (n = 52)

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

Average mean absolute error for all buildings: 457
Sum of absolute errors for all buildings (x 10,000): 4193
Bias (x 10,000): 608

# Add predictions to the already existing plot (not showing the plot yet)
# ==============================================================================
for i, building in enumerate(selected_buildings):
    predictions_all_buildings.query("level==@building")['pred'].plot(ax=axs[i], label='Global model per primary use')
    axs[i].legend()

Global multi-series model by cluster (features)¶

A global model for each cluster based on time series features is trained and tested.

# Forecaster multi-series models for buildings grouped by time series features
# ==============================================================================
predictions_all_buildings = []
metrics_all_buildings = []
start = datetime.now()

for cluster in data['cluster_based_on_features'].unique():

    # Create a subset based on time series feature clusters
    building_id_subset = set(data.loc[data['cluster_based_on_features'] == cluster, 'building_id'].unique())
    series_dict_subset = {k: v for k, v in series_dict.items() if k in building_id_subset}
    exog_dict_subset   = {k: v for k, v in exog_dict.items() if k in building_id_subset}

    print(
        f"Training and testing model for cluster: {cluster} "
        f"(n = {len(building_id_subset)})"
    )

    metrics, predictions = backtesting_forecaster_multiseries(
                               forecaster            = forecaster,
                               series                = series_dict_subset,
                               exog                  = exog_dict_subset,
                               cv                    = cv,
                               metric                = 'mean_absolute_error',
                               add_aggregated_metric = False,
                               verbose               = False,
                               show_progress         = True,
                               suppress_warnings     = True
                           )
    predictions_all_buildings.append(predictions)
    metrics_all_buildings.append(metrics)

end = datetime.now()

predictions_all_buildings = pd.concat(predictions_all_buildings, axis=0)
metrics_all_buildings = pd.concat(metrics_all_buildings, axis=0)

add_results_to_table(
    table_results = table_results,
    model_name    = 'Global model per cluster (features)',
    metrics       = metrics_all_buildings,
    predictions   = predictions_all_buildings,
    data          = data,
    elapsed_time  = end - start
)

Training and testing model for cluster: 0 (n = 277)

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

Training and testing model for cluster: 1 (n = 78)

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

Training and testing model for cluster: 2 (n = 221)

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

Training and testing model for cluster: Other (n = 24)

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

Average mean absolute error for all buildings: 418
Sum of absolute errors for all buildings (x 10,000): 3834
Bias (x 10,000): 591

# Add predictions to the already existing plot (not showing the plot yet)
# ==============================================================================
for i, building in enumerate(selected_buildings):
    predictions_all_buildings.query("level==@building")['pred'].plot(ax=axs[i], label='Global model per cluster (features)')
    axs[i].legend()

Global multi-series model by cluster (elastic distance DTW)¶

A global model for each cluster based on elastic distance (DTW) is trained and tested.

# Forecaster multi-series models for buildings grouped by DTW
# ==============================================================================
predictions_all_buildings = []
metrics_all_buildings = []
start = datetime.now()

for cluster in data['cluster_based_on_dtw'].unique():

    # Create a subset based on DTW clusters
    building_id_subset = set(data.loc[data['cluster_based_on_dtw'] == cluster, 'building_id'].unique())
    series_dict_subset = {k: v for k, v in series_dict.items() if k in building_id_subset}
    exog_dict_subset   = {k: v for k, v in exog_dict.items() if k in building_id_subset}

    print(
        f"Training and testing model for cluster: {cluster} "
        f"(n = {len(building_id_subset)})"
    )

    metrics, predictions = backtesting_forecaster_multiseries(
                               forecaster            = forecaster,
                               series                = series_dict_subset,
                               exog                  = exog_dict_subset,
                               cv                    = cv,
                               metric                = 'mean_absolute_error',
                               add_aggregated_metric = False,
                               verbose               = False,
                               show_progress         = True,
                               suppress_warnings     = True
                           )
    predictions_all_buildings.append(predictions)
    metrics_all_buildings.append(metrics)

end = datetime.now()

predictions_all_buildings = pd.concat(predictions_all_buildings, axis=0)
metrics_all_buildings = pd.concat(metrics_all_buildings, axis=0)

add_results_to_table(
    table_results = table_results,
    model_name    = 'Global model per cluster (DTW)',
    metrics       = metrics_all_buildings,
    predictions   = predictions_all_buildings,
    data          = data,
    elapsed_time  = end - start
)

Training and testing model for cluster: 1 (n = 303)

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

Training and testing model for cluster: 0 (n = 153)

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

Training and testing model for cluster: 2 (n = 81)

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

Training and testing model for cluster: 3 (n = 63)

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

Average mean absolute error for all buildings: 460
Sum of absolute errors for all buildings (x 10,000): 4225
Bias (x 10,000): 639

# Add predictions to the already existing plot (not showing the plot yet)
# ==============================================================================
for i, building in enumerate(selected_buildings):
    predictions_all_buildings.query("level==@building")['pred'].plot(ax=axs[i], label='Global model per cluster (DTW)')
    axs[i].legend()

Results¶

# Table of results
# ==============================================================================
def highlight_best(column):
    # Bias is best when closest to zero, not when most negative;
    # mae and abs_error are best when smallest.
    best_idx = column.abs().idxmin() if column.name == 'bias' else column.idxmin()
    return ['background-color: green' if idx == best_idx else '' for idx in column.index]

table_results['elapsed_time'] = table_results['elapsed_time'].astype(str).str[:7]
table_results.style.apply(highlight_best, subset=['mae', 'abs_error', 'bias'], axis=0).format(precision=0)

Only the grouping based on time series features clearly outperforms the single global model that includes all the buildings: the average MAE falls from 456 to 418, an improvement of around 8%, and the total absolute error decreases in the same proportion (from 41.8 to 38.3 million). The bias is also substantially reduced, from 8.8 to 5.9 million. Since the bias is positive in all cases, all the models tend to overestimate the energy consumption of the buildings during the test period.

In contrast, neither grouping the buildings by their primary use nor clustering them with DTW improves the accuracy of the single global model: the average MAE is 457 and 460, respectively, versus 456, and the total absolute error is also slightly higher (41.9 and 42.2 million versus 41.8 million). Both strategies reduce the bias (6.1 and 6.4 million), but not the magnitude of the errors.

• For the primary use, this is consistent with the exploratory analysis: buildings devoted to the same use can have very different consumption patterns, so the primary use is not a good criterion for deciding which series should be modeled together.

• For DTW, the series were z-normalized before clustering, so buildings are grouped according to the shape of their consumption pattern, regardless of their consumption level. A possible explanation for the lack of improvement is that buildings with similar normalized profiles but very different magnitudes end up in the same group. In contrast, several of the features extracted with tsfresh depend on the level of the series (for example, mean, sum_values or standard_deviation), so the feature-based clusters combine information about both the shape and the level of consumption.

Finally, note that the total training time increases when several models are used instead of one (approximately 15 seconds for the single global model versus 21-22 seconds for the grouped strategies in this execution), although the difference is small because each model is trained with fewer series.

Finally, the predictions of the four models are compared with the real consumption of the two randomly selected buildings.

# Plot predictions vs real value for 2 random buildings for all models
# ==============================================================================
handles, labels = axs[0].get_legend_handles_labels()
fig.legend(handles, labels, loc='upper center', bbox_to_anchor=(0.5, 0.05), ncol=2)
for ax in axs:
    ax.get_legend().remove()
fig

Conclusion¶

This document has shown that clustering time series can be a valuable tool for improving forecasting models. Grouping the buildings with clusters based on features extracted from the series (tsfresh) and training a global model for each group lowered the average error by around 8% compared with a single global model trained on all the buildings at once. A plausible explanation is that each model can focus on the patterns shared by a more homogeneous set of series, although this analysis does not isolate the exact mechanism behind the improvement.

The results also show that not every grouping criterion is useful. Grouping the buildings by their primary use, an attribute that is available in the metadata and that may seem informative a priori, did not improve the accuracy of the single global model. Neither did clustering the z-normalized series with DTW, despite being a method specifically designed to identify series with similar shapes. This highlights that the similarity criterion used to form the groups, including preprocessing decisions such as normalization, has a direct impact on whether clustering helps. The usefulness of a clustering should therefore be validated with the forecasting metrics, not assumed from the quality of the clusters alone.

These results should be interpreted with some caution: the clusters were obtained using the entire observed series, including the test period (see the warning in the section Clustering time series), the evaluation is based on a single train/test split, and only a random subset of 600 buildings was used. Repeating the analysis under different conditions would be necessary to confirm the conclusions.

Finally, it is worth remembering that clustering adds an additional layer of complexity: more models to train, monitor and maintain. This cost must be weighed against the improvement in accuracy that it provides.

Further research

This analysis has provided interesting insights into the effectiveness of combining clustering and forecasting models. Possible next steps could include:

• Manually review buildings with high error. Identify if there is a group for which the model is not performing well.

• Add more exogenous features: see skforecast's user guide Calendar Features for calendar and sunlight features that typically affect energy consumption.

• Optimize lags and hyperparameters: Use Grid, Random or Bayesian search to find the best model configuration.

• Try other machine learning algorithms.

• Build the clusters using only the data available up to the end of the training period, to avoid the lookahead bias described in the section Clustering time series.

• Compare the DTW clusters obtained with and without z-normalization, to assess how much the level of consumption influences the grouping of the buildings and the accuracy of the resulting models.

• Use other criteria to select the number of clusters, such as the silhouette coefficient.

Session information¶

import session_info
session_info.show(html=False)

-----
lightgbm            4.7.0
matplotlib          3.10.9
numpy               2.4.6
pandas              2.3.3
session_info        v1.0.1
skforecast          0.25.0
sklearn             1.7.2
sktime              1.1.0
tsfresh             0.21.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-15 22:08

Citation¶

How to cite this document

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

Clustering Time Series to Improve Forecasting Models by Joaquín Amat Rodrigo and Javier Escobar Ortiz, available under a CC BY-NC-SA 4.0 at https://www.cienciadedatos.net/documentos/py64-clustering-time-series-forecasting.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.