Más sobre forecasting en: cienciadedatos.net

• Modelos ARIMA y SARIMAX con Python
• Forecasting de series temporales con machine learning
• Forecasting de series temporales con gradient boosting: XGBoost, LightGBM y CatBoost
• Forecasting de series temporales con XGBoost
• Modelos de forecasting globales: forecasting multi-serie
• Modelos de forecasting globales: análisis comparativo de modelos de una y múltiples series
• Forecasting probabilístico
• Forecasting con deep learning
• Forecasting de la demanda energética con machine learning
• Forecasting de visitas web con machine learning
• Forecasting de demanda intermitente
• Modelar series temporales con tendencia utilizando modelos de árboles
• Predicción del precio de Bitcoin con Python
• Stacking ensemble de modelos de machine learning para mejorar el forecasting
• Modelos de forecasting interpretables
• Reducir el impacto del Covid en modelos de forecasting
• Forecasting de series temporales con valores ausentes

Introducción¶

El modelado de forecasting global consiste en crear un único modelo de forecasting que considera todas las series temporales simultáneamente. Este enfoque intenta captar los patrones subyacentes comunes a todas las series, lo que reduce el impacto del ruido presente en las series individuales. Ofrece eficiencia computacional, facilidad de mantenimiento y una generalización robusta entre series temporales. El forecasting global presupone que las series temporales con un comportamiento similar pueden beneficiarse de ser modeladas conjuntamente. Cuando se trabaja con cientos o miles de series, un análisis de clustering inicial puede ayudar a maximizar el rendimiento del modelo, ya que permite identificar grupos de series con un comportamiento similar. Los métodos de clustering utilizados en este documento, basados en características de las series temporales y en medidas de distancia elástica, se describen en detalle en el apartado Clustering de series temporales.

Si después de realizar el análisis de clustering se identifican grupos distintos de series, se pueden utilizar dos estrategias:

• Modelar todas las series juntas, añadiendo una variable que indique a qué cluster pertenece cada serie.

• Construir varios modelos globales, cada uno adaptado a un grupo específico de series.

Para comprender mejor los beneficios del clustering, este documento se centra en examinar y comparar los resultados obtenidos al predecir el consumo energético de más de mil edificios (en el apartado de modelado se utiliza un subconjunto aleatorio de 600 edificios para mantener tiempos de ejecución razonables). Aunque el uso principal de cada edificio está disponible en el conjunto de datos, este puede no reflejar grupos con patrones similares de consumo energético, por lo que se crean grupos adicionales mediante métodos de clustering. Se realizan un total de 4 experimentos:

• Modelar todos los edificios juntos con una estrategia de modelo global único.

• Modelar grupos de edificios según su uso principal (un modelo global por uso principal).

• Modelar grupos de edificios según el clustering basado en características de las series temporales (un modelo global por cluster).

• Modelar grupos de edificios según el clustering basado en Dynamic Time Warping (DTW) (un modelo global por cluster).

Librerías¶

Librerías utilizadas en este documento:

# Manipulación de datos
# ==============================================================================
import numpy as np
import pandas as pd
from datetime import datetime
from skforecast.datasets import fetch_dataset

# Gráficos
# ==============================================================================
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
)

# Extracción de características
# ==============================================================================
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

# Configuración de warnings
# ==============================================================================
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

Datos¶

Los datos utilizados en este documento han sido obtenidos de la competición de Kaggle
Addison Howard, Chris Balbach, Clayton Miller, Jeff Haberl, Krishnan Gowri, Sohier Dane. (2019). ASHRAE - Great Energy Predictor III. Kaggle.

Se utilizan tres archivos para crear el conjunto de datos de modelado:

• weather_train.csv y weather_test.csv: Estos archivos contienen datos meteorológicos para cada edificio, incluyendo temperatura del aire exterior, temperatura del punto de rocío, humedad relativa y otros parámetros climáticos. Los datos meteorológicos son cruciales para comprender el impacto de las condiciones externas en el consumo energético de los edificios.

• building_metadata.csv: Este archivo proporciona metadatos para cada edificio en el conjunto de datos, como el tipo de edificio, uso principal, superficie en pies cuadrados, número de pisos y año de construcción. Esta información ayuda a entender las características de los edificios y su posible influencia en los patrones de consumo energético.

• train.csv: El conjunto de datos de entrenamiento contiene la variable objetivo, es decir, los datos de consumo energético de cada edificio, junto con las marcas de tiempo de las mediciones. También incluye los identificadores de edificios y condiciones climáticas para vincular la información entre los distintos conjuntos de datos.

Los tres archivos han sido preprocesados para eliminar edificios con menos del 85% de valores distintos de NaN o cero, utilizar únicamente el medidor de electricidad y agregar los datos a una frecuencia diaria.

# Carga de datos preprocesados
# ==============================================================================
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                                                  │
╰──────────────────────────────────────────────────────────────────────────────────╯

# Asegurar que el índice de cada serie temporal está completo, sin huecos
# ==============================================================================
data = (
    data
    .groupby('building_id')
    .apply(lambda group: group.asfreq('D'), include_groups=False)
    .reset_index()
    .set_index('timestamp')
)

# Imputar valores ausentes de air_temperature y wind_speed con forward y backward fill
# ==============================================================================
# La imputación debe hacerse de forma independiente para cada edificio
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)

Análisis exploratorio¶

Uso designado de los edificios¶

Una de las variables clave asociadas a cada edificio es su uso designado. Esta característica puede desempeñar un papel crucial en la influencia del patrón de consumo energético, ya que distintos usos pueden impactar significativamente tanto en la cantidad como en el momento del consumo de energía.

# Número de edificios y tipos de edificio según su uso principal
# ==============================================================================
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

Para ciertas categorías de uso principal, el número de edificios en el conjunto de datos es limitado. Para simplificar el análisis, las categorías con menos de 50 edificios se agrupan en la categoría "Other".

# Los tipos de edificio (primary use) con menos de 50 edificios se agrupan como "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

A continuación, se crean dos gráficos: el primero muestra el consumo energético de un edificio seleccionado aleatoriamente dentro de cada categoría, y el segundo muestra todas las series temporales disponibles de cada categoría (una línea gris por edificio), con el consumo medio de la categoría en azul.

# Serie temporal de un edificio seleccionado al azar por grupo
# ==============================================================================
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("")
    # Notación científica para el eje y
    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()

# Consumo energético por tipo de edificio (una línea gris por edificio)
# ==============================================================================
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("")
    # Notación científica para el eje y
    axs[i].ticklabel_format(axis='y', style='sci', scilimits=(0, 0))
    axs[i].title.set_size(9)

    # Limitar el eje a 5 veces el valor medio máximo para mejorar la visualización
    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()

El gráfico revela una variabilidad significativa en los patrones de consumo entre edificios con el mismo propósito. Esto sugiere que puede haber margen para mejorar los criterios por los cuales se agrupan los edificios.

# Muestreo de 600 edificios
# ==============================================================================
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 de series temporales¶

La idea detrás del modelado de múltiples series al mismo tiempo es capturar los principales patrones que rigen las series, reduciendo así el impacto del posible ruido presente en cada una. Esto significa que las series con comportamientos similares pueden beneficiarse al ser modeladas juntas. Una forma de identificar posibles grupos de series es realizar un estudio de agrupación (clustering) antes del modelado. Si el análisis de clustering identifica grupos claramente diferenciados, es apropiado modelar cada grupo por separado.

El clustering es una técnica de análisis no supervisado que agrupa un conjunto de observaciones en clusters que contienen observaciones consideradas homogéneas, mientras que las observaciones en distintos clusters se consideran heterogéneas. Los algoritmos de agrupación de series temporales pueden dividirse en dos grupos: aquellos que utilizan una transformación para crear características antes de la agrupación (agrupación basada en características) y aquellos que trabajan directamente sobre las series temporales (medidas de distancia elástica).

• Agrupación basada en características: Se extraen características que describen aspectos estructurales de cada serie temporal, y luego estas características se introducen en algoritmos de clustering convencionales. Estas características se obtienen mediante operaciones estadísticas que capturan mejor los rasgos subyacentes, como tendencia, estacionalidad, periodicidad, correlación serial, asimetría, curtosis, caos, no linealidad y autosimilitud.

• Medidas de distancia elástica: Este enfoque trabaja directamente sobre las series temporales, ajustándolas o "realineándolas" en comparación unas con otras. La medida más conocida de esta familia es Dynamic Time Warping (DTW).

Para una revisión más detallada sobre clustering de series temporales, consulta A review and evaluation of elastic distance functions for time series clustering.

Clustering basado en características¶

Extracción de características¶

tsfresh es una potente librería de Python para la ingeniería de características (feature engineering) a partir de series temporales y datos secuenciales, que incluye medidas estadísticas, coeficientes de Fourier y otras características de los dominios del tiempo y de la frecuencia. Proporciona un enfoque sistemático para automatizar el cálculo de las características y seleccionar las más informativas.

Para comenzar, se utiliza la configuración predeterminada de tsfresh, que calcula todas las características disponibles. La forma más sencilla de acceder a esta configuración es instanciar la clase ComprehensiveFCParameters, que devuelve un diccionario que asigna el nombre de cada característica (una cadena de texto) a una lista de diccionarios con los parámetros utilizados al calcular dicha característica (o None si la característica no tiene parámetros).

# Características y configuración por defecto de 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

Muchas de estas características se calculan utilizando diferentes valores para sus argumentos.

# Configuración por defecto para "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}]

Para acceder a la vista detallada de cada característica y los valores de los parámetros incluidos en la configuración predeterminada, utiliza el siguiente código:

# Configuración por defecto para todas las características extraídas
# ==============================================================================
# from pprint import pprint
# pprint(default_features)

Una vez que se ha definido la configuración de las características, el siguiente paso es extraerlas de la serie temporal. Para ello, se utiliza la función extract_features() de tsfresh. Esta función recibe como entrada la serie temporal y la configuración de las características a extraer. La salida es un dataframe con las características extraídas.

# Extracción de características
# ==============================================================================
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)

Como resultado del proceso de extracción, se han creado 783 características para cada serie temporal (building_id en este caso de uso). El dataframe devuelto tiene como índice la columna especificada en el argumento column_id de extract_features.

La extracción predeterminada de tsfresh genera una gran cantidad de características. Sin embargo, solo algunas de ellas pueden ser de interés en cada caso de uso. Para seleccionar las más relevantes, tsfresh incluye un proceso de selección automatizado basado en pruebas de hipótesis, conocido como FRESH (FeatuRe Extraction based on Scalable Hypothesis tests). En este proceso, las características se evalúan de manera individual e independiente para determinar su relevancia en la predicción del objetivo en estudio.

En este caso, dado que el objetivo es identificar grupos de edificios con patrones de consumo de energía similares, independientemente del tipo de edificio, no se utiliza el proceso de selección automatizado. En su lugar, se aplica la tercera estrategia descrita en el aviso anterior: la dimensionalidad de la matriz de características se reduce con PCA antes del clustering.

Algunas características no pueden calcularse para ciertas series temporales y devuelven valores ausentes (NaN) o infinitos. Dado que la mayoría de los algoritmos de clustering no admiten valores ausentes, es necesario tratarlos. En este caso, ya han sido reemplazados durante la extracción por la función impute pasada al argumento impute_function de extract_features(): los valores -inf y +inf se reemplazan por el mínimo y el máximo de cada característica, y los valores NaN por su mediana. Por lo tanto, el siguiente paso, que elimina cualquier característica con valores ausentes, actúa solo como medida de seguridad y no se elimina ninguna característica.

# Eliminar características con valores ausentes
# ==============================================================================
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

Una vez que la matriz final de características está lista, puede ser útil crear un nuevo diccionario para almacenar las características finales y los parámetros utilizados para calcularlas. Esto se puede hacer fácilmente utilizando la función from_columns.

# Diccionario con las características seleccionadas y su configuración
# ==============================================================================
ts_features_info = from_columns(ts_features)
# pprint(ts_features_info['meter_reading'])

K-means clustering¶

Se utiliza el método de clustering K-means para agrupar los edificios. Dado que el clustering se ve negativamente afectado por la alta dimensionalidad, y como se han creado varios cientos de características para cada edificio, se emplea PCA para reducir la dimensionalidad de los datos antes de aplicar K-means. Como PCA es sensible a la escala de las variables, las características se estandarizan primero para que todas tengan media 0 y desviación estándar 1.

# Escalar características para que todas tengan media 0 y desviación estándar 1
# ==============================================================================
scaler = StandardScaler().set_output(transform="pandas")
ts_features_scaled = scaler.fit_transform(ts_features)
ts_features_scaled.head(2)

# Gráfico de reducción de varianza en función del número de componentes PCA
# ==============================================================================
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();

La dimensionalidad de las 783 características originales se reduce utilizando los primeros 57 componentes principales (que explican el 85% de la varianza).

# PCA con tantos componentes como sea necesario para explicar el 85% de la varianza
# ==============================================================================
pca = PCA(n_components=0.85)
pca_projections = pca.fit_transform(ts_features_scaled)

# Crear un data frame con las proyecciones. Cada columna es un componente
# principal y cada fila es el id del edificio.
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

Uno de los desafíos inherentes al clustering es determinar el número óptimo de clusters, ya que en el aprendizaje no supervisado no existe una verdad de referencia ni un número predefinido de clusters. Se han propuesto varias heurísticas para guiar esta selección y, en este caso, se utiliza el método del codo.

El método del codo consiste en representar la suma total de cuadrados intra-cluster (WSS, within-cluster sum of squares) en función del número de clusters (k). El número óptimo de clusters se identifica típicamente en el "codo" de la curva, donde esta se aplana: más allá de este punto, añadir más clusters ya no reduce sustancialmente la WSS. En scikit-learn, la WSS de un modelo KMeans ajustado se almacena en su atributo inertia_.

# Número óptimo de clusters
# ==============================================================================
# Identificar el número óptimo de clusters mediante el método del codo
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();

El gráfico muestra que la disminución de la inercia se desacelera notablemente a partir de 5 clusters, por lo que los valores de esa región son candidatos razonables.

Sin embargo, el criterio del codo no debería ser la única guía. Además de analizar la evolución de la inercia (intra-varianza), es importante verificar el tamaño de los clusters que se crean. La presencia de clusters pequeños puede indicar sobreajuste o muestras anómalas que no encajan bien en ninguno de los grupos.

# Distribución de tamaños de los clusters
# ==============================================================================
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]

Los tamaños confirman que el valor del codo por sí solo no es suficiente. Con 5 clusters, un grupo contiene 481 de los 600 edificios y tres de los grupos restantes tienen 18 edificios o menos, por lo que, en la práctica, casi todas las series se seguirían modelando juntas. Con 7 clusters la partición es más informativa: tres grupos de un tamaño utilizable (277, 221 y 78 edificios) y varios muy pequeños.

Por este motivo, se utilizan 7 clusters, y los más pequeños (menos de 20 edificios) se combinan en un único grupo denominado "Other". Esto evita entrenar un modelo global con solo unas pocas series, lo que iría en contra del propósito del enfoque multi-serie.

# Entrenar el modelo de clustering con 7 clusters y asignar cada edificio a un 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)
})

# Combinar los clusters pequeños en un único 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')

# Añadir el cluster predicho al data frame con la información de los edificios
# ==============================================================================
data = pd.merge(
           data.reset_index(),  # Para no perder el índice
           clusters,
           on       = 'building_id',
           how      = 'left',
           validate = 'm:1'
       ).set_index('timestamp')

data.head(3)

Una vez que se ha asignado un cluster a cada edificio, es útil examinar las características de los edificios agrupados. Por ejemplo, la distribución del uso principal de los edificios.

# Porcentaje de tipos de edificio por 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

Los resultados sugieren (la tabla debe leerse fila a fila, ya que los porcentajes de cada fila suman el 100%) que el proceso de clustering basado en las características extraídas de las series temporales genera grupos que difieren de los formados por el uso principal del edificio. Ningún cluster está dominado por un único tipo de edificio: por ejemplo, los edificios educativos representan entre el 25% y el 55% de cada cluster, y los de entretenimiento/reunión pública, entre el 12% y el 15%.

Clustering basado en distancia elástica (DTW)¶

Dynamic Time Warping (DTW) es una técnica que mide la similitud entre dos secuencias temporales, que pueden variar en velocidad. En esencia, es una medida de distancia elástica que permite que las series temporales se estiren o compriman para alinearse entre sí de manera óptima. El clustering con DTW implica agrupar las series temporales en función de sus distancias DTW, asegurando que las series temporales dentro del mismo cluster tengan formas y patrones similares, incluso si están desfasadas en el tiempo o tienen longitudes diferentes.

La clase TimeSeriesKMeans de la librería sktime permite la aplicación de clustering K-means con una variedad de métricas de distancia, incluyendo DTW, Euclidiana, ERP, EDR, LCSS, cuadrada, DDTW, WDTW y WDDTW. Muchas de estas métricas son distancias elásticas, lo que hace que el método sea muy adecuado para series temporales.

A diferencia del apartado anterior, aquí el número de clusters no se selecciona con el método del codo: calcular las distancias DTW entre cientos de series lleva varios minutos por ajuste, lo que hace poco práctico explorar un rango amplio de valores. En su lugar, se fijan a priori 4 clusters, el mismo número de grupos utilizado finalmente en el enfoque basado en características (3 clusters más el grupo "Other").

sktime requiere que las series temporales estén estructuradas en un formato largo con un multiíndice. El nivel más externo del índice representa el identificador de la serie, mientras que el nivel más interno corresponde a la fecha.

# Convertir los datos a formato long con un multiíndice: (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-normalizar cada serie (media 0, desviación estándar 1) para que DTW compare
# la forma de los patrones de consumo y no su nivel
# ==============================================================================
grouped = data_long.groupby(level='building_id')['meter_reading']
series_mean = grouped.transform('mean')
series_std = grouped.transform('std').replace(0, 1)  # Evitar la división por cero en series constantes
data_long_scaled = ((data_long['meter_reading'] - series_mean) / series_std).to_frame()
data_long_scaled

# Ajuste del modelo de clustering
# ==============================================================================
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)

# Predicción de clusters
# ==============================================================================
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),
           })

# Tamaño de cada cluster
clusters['cluster_based_on_dtw'].value_counts().sort_index().to_frame(name='Number of series')

En este caso, todos los clusters contienen más de 60 edificios, por lo que, a diferencia del enfoque basado en características, no es necesario combinar clusters pequeños.

# Añadir el cluster predicho al data frame con la información de los edificios
# ==============================================================================
data = pd.merge(
           data.reset_index(),  # Para no perder el índice
           clusters,
           on       = 'building_id',
           how      = 'left',
           validate = 'm:1'
       ).set_index('timestamp')

data.head(3)

# Guardar en disco los datos para modelado para evitar repetir los pasos anteriores
# ==============================================================================
data.to_parquet('data_modelling.parquet')

Modelado y forecasting¶

Después de establecer los tres criterios de agrupación (uso principal del edificio, clustering basado en características de las series temporales y clustering basado en dynamic time warping), se entrena un modelo global multi-serie para cada grupo y se comparan los resultados con los de un único modelo global que incluye todos los edificios. La evaluación que sigue se centra en determinar la eficacia con la que estos modelos predicen el consumo diario durante los últimos cinco meses del año (de agosto a diciembre), que se predicen en horizontes sucesivos de 7 días. Durante esta evaluación, se utilizan tres métricas distintas:

• El valor promedio del Error Absoluto Medio (MAE) para todos los edificios.

• La suma de los errores absolutos, es decir, la suma de las diferencias absolutas entre los valores predichos y el consumo real para todos los edificios y fechas.

• El bias (sesgo), calculado como la suma de los errores (valor predicho menos consumo real) para todos los edificios y fechas. Los valores positivos indican que el modelo sobrestima el consumo, y los negativos, que lo subestima.

Además de los valores rezagados de cada serie temporal, el modelo incluye el día de la semana (codificado con seno-coseno), creado por el forecaster mediante el argumento calendar_features, y las siguientes variables exógenas: temperatura exterior (air_temperature), velocidad del viento (wind_speed) y uso principal del edificio (primary_use). Dado que primary_use es una columna de texto, el forecaster la trata automáticamente como una variable categórica (categorical_features='auto').

Para entrenar los modelos y evaluar su capacidad predictiva, los datos se dividen en dos conjuntos: entrenamiento (de enero a julio) y test (de agosto a diciembre). En este documento no se realiza ninguna búsqueda de hiperparámetros; se utiliza la misma configuración del modelo en todos los experimentos para que las diferencias en los resultados dependan únicamente de cómo se agrupan los edificios.

# Lectura de datos para modelado
# ==============================================================================
data = pd.read_parquet('data_modelling.parquet')

Los datos se transforman de un data frame en formato largo a un diccionario de series utilizando reshape_series_long_to_dict y reshape_exog_long_to_dict. Aunque no es estrictamente necesario, ya que skforecast admite múltiples formatos de entrada, es el formato recomendado para los modelos multi-serie, ya que permite combinar fácilmente series de diferentes longitudes y con distintos subconjuntos de variables exógenas.

# Transformar series y variables exógenas a diccionarios
# ==============================================================================
series_dict = reshape_series_long_to_dict(
    data      = data.reset_index(),
    series_id = 'building_id',
    index     = 'timestamp',
    values    = 'meter_reading',
    freq      = 'D'
)

# Las variables exógenas incluidas en todos los modelos son: primary_use,
# air_temperature y 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'
)

# Fin del periodo de entrenamiento y series de test
# ==============================================================================
# La división entrenamiento/test utilizada en el modelado la gestiona internamente
# TimeSeriesFold (ver `cv` más abajo), que recibe series_dict/exog_dict completos,
# sin dividir. Las series de test se extraen solo para representarlas frente a
# las predicciones.
end_train = '2016-07-31 23:59:00'
series_dict_test = {k: v.loc[end_train:] for k, v in series_dict.items()}

Se utiliza la misma configuración del forecaster en todos los experimentos, de modo que la única diferencia entre ellos es la forma en que se agrupan los edificios:

• Regresor: un LGBMRegressor con 500 árboles, una profundidad máxima de 10 y una tasa de aprendizaje pequeña (0.01).

• Lags: los últimos 31 valores de cada serie, suficientes para que el modelo capture patrones tanto semanales como mensuales.

• Variables de calendario: el día de la semana, codificado de forma cíclica con CalendarFeatures. La codificación seno-coseno (cíclica) garantiza que el domingo y el lunes estén tan próximos entre sí como cualquier otro par de días consecutivos, algo que no se lograría con una codificación numérica simple (de 0 a 6).

• Codificación de las series: encoding="ordinal_category" añade el identificador de la serie (el id del edificio) como variable categórica, lo que permite que un único modelo global aprenda el comportamiento específico de cada edificio. Las estrategias disponibles se describen en la guía de usuario Global Forecasting Models.

# Definición del forecaster (común a todos los experimentos)
# ==============================================================================
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"
             )

La capacidad predictiva se estima mediante backtesting, utilizando la clase TimeSeriesFold para definir el esquema de validación y la función backtesting_forecaster_multiseries para ejecutarlo:

• initial_train_size = end_train: el modelo se entrena con los datos hasta el 2016-07-31, y los cinco meses restantes (de agosto a diciembre) se utilizan como conjunto de test.

• steps = 7: las predicciones se generan en horizontes sucesivos de 7 días (una semana), lo que da lugar a 22 folds.

• refit = False: el modelo se entrena una sola vez, al inicio del proceso, y no se reentrena a medida que avanza la ventana de backtesting. Esto mantiene asumible, en términos computacionales, la comparación de las cuatro estrategias.

Dado que los datos de entrenamiento siempre preceden a los valores predichos, este esquema respeta el orden temporal de los datos y no se filtra al modelo información del futuro.

# Definición del backtesting
# ==============================================================================
cv = TimeSeriesFold(
        steps              = 7,
        initial_train_size = end_train,
        refit              = False
     )

# Tabla de resultados de todos los modelos
# ==============================================================================
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.
    """
    # Combinar las predicciones con los valores reales
    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}"
    )

Modelo global multi-serie para todos los edificios¶

Se entrena y evalúa un modelo global para todos los edificios utilizando la clase ForecasterRecursiveMultiSeries de skforecast. Este forecaster admite las series temporales en varios formatos; en este caso, se utilizan los diccionarios de series y de variables exógenas creados anteriormente.

Para más información sobre cómo utilizar series de diferentes longitudes o variables exógenas diferentes para cada serie, consulta Series con diferentes longitudes y diferentes variables exógenas.

# Forecaster multi-serie para modelar todos los edificios a la vez
# ==============================================================================
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

# Gráfico de predicciones frente a valor real para 2 edificios aleatorios
# ==============================================================================
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();

Modelo global multi-serie por uso principal¶

Se entrena y evalúa un modelo global para cada uso principal de los edificios.

# Forecasters multi-serie para los edificios agrupados por uso principal
# ==============================================================================
predictions_all_buildings = []
metrics_all_buildings = []
start = datetime.now()

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

    # Crear un subconjunto según el uso principal
    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

# Añadir predicciones al gráfico ya existente (sin mostrar el gráfico todavía)
# ==============================================================================
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()

Modelo global multi-serie por cluster (características)¶

Se entrena y evalúa un modelo global para cada cluster basado en características de las series temporales.

# Forecasters multi-serie para los edificios agrupados por características de las series
# ==============================================================================
predictions_all_buildings = []
metrics_all_buildings = []
start = datetime.now()

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

    # Crear un subconjunto según los clusters basados en características
    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

# Añadir predicciones al gráfico ya existente (sin mostrar el gráfico todavía)
# ==============================================================================
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()

Modelo global multi-serie por cluster (distancia elástica DTW)¶

Se entrena y evalúa un modelo global para cada cluster basado en la distancia elástica (DTW).

# Forecasters multi-serie para los edificios agrupados por DTW
# ==============================================================================
predictions_all_buildings = []
metrics_all_buildings = []
start = datetime.now()

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

    # Crear un subconjunto según los clusters DTW
    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

# Añadir predicciones al gráfico ya existente (sin mostrar el gráfico todavía)
# ==============================================================================
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()

Resultados¶

# Tabla de resultados
# ==============================================================================
def highlight_best(column):
    # El mejor bias es el más cercano a cero, no el más negativo;
    # para mae y abs_error, el mejor valor es el menor.
    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)

Solo la agrupación basada en características de las series temporales supera claramente al modelo global único que incluye todos los edificios: el MAE promedio baja de 456 a 418, una mejora de alrededor del 8%, y el error absoluto total disminuye en la misma proporción (de 41,8 a 38,3 millones). El bias también se reduce sustancialmente, de 8,8 a 5,9 millones. Dado que el bias es positivo en todos los casos, todos los modelos tienden a sobrestimar el consumo energético de los edificios durante el periodo de test.

En cambio, ni agrupar los edificios por su uso principal ni agruparlos mediante clustering con DTW mejora la precisión del modelo global único: el MAE promedio es de 457 y 460, respectivamente, frente a 456, y el error absoluto total también es ligeramente mayor (41,9 y 42,2 millones frente a 41,8 millones). Ambas estrategias reducen el bias (6,1 y 6,4 millones), pero no la magnitud de los errores.

• En el caso del uso principal, esto es coherente con el análisis exploratorio: los edificios destinados a un mismo uso pueden tener patrones de consumo muy diferentes, por lo que el uso principal no es un buen criterio para decidir qué series deben modelarse juntas.

• En el caso de DTW, las series se z-normalizaron antes del clustering, por lo que los edificios se agrupan según la forma de su patrón de consumo, independientemente de su nivel de consumo. Una posible explicación de la falta de mejora es que edificios con perfiles normalizados similares pero magnitudes muy diferentes acaban en el mismo grupo. En cambio, varias de las características extraídas con tsfresh dependen del nivel de la serie (por ejemplo, mean, sum_values o standard_deviation), por lo que los clusters basados en características combinan información tanto de la forma como del nivel de consumo.

Por último, cabe señalar que el tiempo total de entrenamiento aumenta cuando se utilizan varios modelos en lugar de uno (aproximadamente 15 segundos para el modelo global único frente a 21-22 segundos para las estrategias con agrupación en esta ejecución), aunque la diferencia es pequeña porque cada modelo se entrena con menos series.

Finalmente, se comparan las predicciones de los cuatro modelos con el consumo real de los dos edificios seleccionados aleatoriamente.

# Gráfico de predicciones frente a valor real de 2 edificios para todos los modelos
# ==============================================================================
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

Conclusión¶

Este documento ha mostrado que el clustering de series temporales puede ser una herramienta valiosa para mejorar los modelos de forecasting. Agrupar los edificios mediante clusters basados en características extraídas de las series (tsfresh) y entrenar un modelo global para cada grupo redujo el error promedio en torno a un 8% en comparación con un único modelo global entrenado con todos los edificios a la vez. Una explicación plausible es que cada modelo puede centrarse en los patrones compartidos por un conjunto de series más homogéneo, aunque este análisis no aísla el mecanismo exacto responsable de la mejora.

Los resultados también muestran que no todos los criterios de agrupación son útiles. Agrupar los edificios por su uso principal, un atributo disponible en los metadatos y que a priori puede parecer informativo, no mejoró la precisión del modelo global único. Tampoco lo hizo el clustering de las series z-normalizadas con DTW, a pesar de ser un método diseñado específicamente para identificar series con formas similares. Esto pone de manifiesto que el criterio de similitud utilizado para formar los grupos, incluidas decisiones de preprocesamiento como la normalización, tiene un impacto directo en que el clustering resulte útil o no. Por lo tanto, la utilidad de un clustering debe validarse con las métricas de forecasting, y no darse por supuesta únicamente a partir de la calidad de los clusters.

Estos resultados deben interpretarse con cierta cautela: los clusters se obtuvieron utilizando la serie observada completa, incluido el periodo de test (ver el aviso en el apartado Clustering de series temporales), la evaluación se basa en una única partición entrenamiento/test y solo se utilizó un subconjunto aleatorio de 600 edificios. Sería necesario repetir el análisis en condiciones diferentes para confirmar las conclusiones.

Por último, conviene recordar que el clustering añade una capa adicional de complejidad: más modelos que entrenar, monitorizar y mantener. Este coste debe sopesarse frente a la mejora en precisión que aporta.

Posibles mejoras

Este análisis ha proporcionado ideas interesantes sobre la efectividad de combinar clustering y modelos de forecasting. Los siguientes pasos posibles podrían incluir:

• Revisar manualmente los edificios con alto error. Identificar si hay algún grupo para el cual el modelo no está funcionando bien.

• Añadir más variables exógenas: consulta la guía de usuario de skforecast Calendar Features para variables del calendario y de luz solar que suelen afectar al consumo de energía.

• Optimizar lags e hiperparámetros: utiliza búsqueda Grid, Random o Bayesiana para encontrar la mejor configuración del modelo.

• Probar otros algoritmos de machine learning.

• Construir los clusters utilizando únicamente los datos disponibles hasta el final del periodo de entrenamiento, para evitar el sesgo de anticipación (lookahead bias) descrito en el apartado Clustering de series temporales.

• Comparar los clusters DTW obtenidos con y sin z-normalización, para evaluar en qué medida el nivel de consumo influye en la agrupación de los edificios y en la precisión de los modelos resultantes.

• Utilizar otros criterios para seleccionar el número de clusters, como el coeficiente de silueta.

Información de sesión¶

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:54

Instrucciones para citar¶

¿Cómo citar este documento?

Si utilizas este documento o alguna parte de él, te agradecemos que lo cites. ¡Muchas gracias!

Clustering de series temporales para mejorar modelos de forecasting por Joaquín Amat Rodrigo y Javier Escobar Ortiz, disponible bajo una licencia Attribution-NonCommercial-ShareAlike 4.0 International (CC BY-NC-SA 4.0 DEED) en https://cienciadedatos.net/documentos/py64-clustering-series-temporales-forecasting.html

¿Cómo citar skforecast?

Si utilizas skforecast, te agradeceríamos mucho que lo cites. ¡Muchas gracias!

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} }

¿Te ha gustado el artículo? Tu ayuda es importante

Tu contribución me ayudará a seguir generando contenido divulgativo gratuito. ¡Muchísimas gracias! 😊

Este documento creado por Joaquín Amat Rodrigo y Javier Escobar Ortiz tiene licencia Attribution-NonCommercial-ShareAlike 4.0 International.

Se permite:

• Compartir: copiar y redistribuir el material en cualquier medio o formato.

• Adaptar: remezclar, transformar y crear a partir del material.

Bajo los siguientes términos:

• Atribución: Debes otorgar el crédito adecuado, proporcionar un enlace a la licencia e indicar si se realizaron cambios. Puedes hacerlo de cualquier manera razonable, pero no de una forma que sugiera que el licenciante te respalda o respalda tu uso.

• No-Comercial: No puedes utilizar el material para fines comerciales.

• Compartir-Igual: Si remezclas, transformas o creas a partir del material, debes distribuir tus contribuciones bajo la misma licencia que el original.