Más sobre forecasting en: cienciadedatos.net
- Forecasting series temporales con machine learning
- Modelos ARIMA y SARIMAX
- Forecasting series temporales con gradient boosting: XGBoost, LightGBM y CatBoost
- Global Forecasting: Multi-series forecasting
- Forecasting de la demanda eléctrica con machine learning
- Forecasting con deep learning
- Forecasting de visitas a página web con machine learning
- Forecasting del precio de Bitcoin
- Forecasting probabilístico
- Forecasting de demanda intermitente
- Reducir el impacto del Covid en modelos de forecasting
- Modelar series temporales con tendencia utilizando modelos de árboles

Introducción
Los modelos de forecasting globales son modelos de forecasting que consideran múltiples 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 buena generalización. 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, la agrupación inicial puede ayudar a maximizar el rendimiento del modelo.
El clustering es una técnica de análisis no supervisada que agrupa un conjunto de observaciones en clústeres que contienen observaciones que se consideran homogéneas, mientras que las observaciones de los clústeres diferentes se consideran heterogéneas. Los algoritmos que agrupan series temporales pueden dividirse en dos grupos: los que utilizan una transformación para extraer características de las series antes de la agrupación (agrupación de series temporales basada en características) y los que trabajan directamente con las series temporales (medidas de distancia elásticas, Dynamic Time Warping).
Agrupación basada en características: Se extraen características que describen aspectos estructurales de cada serie temporal y luego se aplican algoritmos de agrupación convencionales. Estas características se obtienen mediante operaciones estadísticas que capturan mejor los rasgos subyacentes: tendencia, estacionalidad, periodicidad, correlación serial, asimetría, curtosis, no linealidad ...
Medidas de distancia elástica: Este enfoque trabaja directamente sobre las series temporales, ajustándolas o "realineándolas" en comparación con otras. La medida más conocida de esta familia es Dynamic Time Warping (DTW).
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, agregando una variable que indique a qué grupo pertenece cada serie.
Construir varios modelos globales, cada uno adaptado a un grupo específico de series.
Para comprender mejor los beneficios de la agrupación, este documento se centra en examinar y comparar los resultados obtenidos al predecir el consumo energético de más de mil edificios. Aunque el uso principal de cada edificio está disponible en el conjunto de datos, este puede no reflejar patrones similares de consumo energético, por lo que se crean grupos adicionales mediante métodos de clustering. Se realizan un total de cuatro 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 la agrupación basada en características de series temporales (un modelo global por grupo).
Modelar grupos de edificios según la agrupación basada en Dynamic Time Warping (DTW) (un modelo global por grupo).
💡 Tip
Para aprender más sobre los Modelos de Pronóstico Global, consulta los siguientes artículos:Librerías
Librerías utilizadas en este documento:
# 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
plt.style.use('seaborn-v0_8-darkgrid')
# Forecasting
# ==============================================================================
import lightgbm
from lightgbm import LGBMRegressor
from sklearn.preprocessing import StandardScaler
import skforecast
from skforecast.recursive import ForecasterRecursiveMultiSeries
from skforecast.model_selection import TimeSeriesFold
from skforecast.model_selection import backtesting_forecaster_multiseries
# Feature extraction
# ==============================================================================
import tsfresh
from tsfresh import extract_features
from tsfresh.feature_extraction.settings import ComprehensiveFCParameters
from tsfresh.feature_extraction.settings import from_columns
# 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}Versión skforecast: {skforecast.__version__}")
print(f"{color}Versión scikit-learn: {sklearn.__version__}")
print(f"{color}Versión lightgbm: {lightgbm.__version__}")
print(f"{color}Versión tsfresh: {tsfresh.__version__}")
print(f"{color}Versión sktime: {sktime.__version__}")
Versión skforecast: 0.15.1 Versión scikit-learn: 1.5.2 Versión lightgbm: 4.6.0 Versión tsfresh: 0.21.0 Versión sktime: 0.36.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
yweather_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.
# Lectura de datos
# ==============================================================================
data = fetch_dataset('ashrae_daily')
data.head()
ashrae_daily ------------ Daily energy consumption data from the ASHRAE competition with building metadata and weather data. 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 Shape of the dataset: (444324, 10)
building_id | meter_reading | site_id | primary_use | square_feet | air_temperature | dew_temperature | sea_level_pressure | wind_direction | wind_speed | |
---|---|---|---|---|---|---|---|---|---|---|
timestamp | ||||||||||
2016-01-01 | id_105 | 1102.766933 | 1 | Education | 50623 | 3.8 | 2.4 | 1020.9 | 240.0 | 3.1 |
2016-01-01 | id_106 | 17.606200 | 1 | Education | 5374 | 3.8 | 2.4 | 1020.9 | 240.0 | 3.1 |
2016-01-01 | id_107 | 8233.625107 | 1 | Education | 97532 | 3.8 | 2.4 | 1020.9 | 240.0 | 3.1 |
2016-01-01 | id_108 | 4289.478432 | 1 | Education | 81580 | 3.8 | 2.4 | 1020.9 | 240.0 | 3.1 |
2016-01-01 | id_109 | 3911.921989 | 1 | Education | 56995 | 3.8 | 2.4 | 1020.9 | 240.0 | 3.1 |
# Imputación de valores faltantes de air_temperature y wind_speed usando fowrard y backward fill
# ==============================================================================
# La imputación debe hacerse de forma independiente para cada serie temporal
data = data.sort_values(by=['building_id', 'timestamp'])
data['air_temperature'] = data.groupby('building_id')['air_temperature'].ffill().bfill()
data['wind_speed'] = data.groupby('building_id')['wind_speed'].ffill().bfill()
data = data.sort_index()
print(
f"Rage of dates available : {data.index.min()} --- {data.index.max()} "
f"(n_days={(data.index.max() - data.index.min()).days})"
)
Rage of dates available : 2016-01-01 00:00:00 --- 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 tipo de edificios basado en el uso principal
# ==============================================================================
n_building = data['building_id'].nunique()
n_type_building = data['primary_use'].nunique()
range_datetime = [data.index.min(), data.index.max()]
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 "Otros"
.
# Tipos de edificios (primary use) con menos de 50 muestras 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 genera un gráfico que muestra el consumo energético de un edificio seleccionado aleatoriamente dentro de cada categoría, junto con un gráfico de todas las series temporales disponibles para cada categoría.
# Serie temporal para 1 edificio seleccionado al azar por grupo
# ==============================================================================
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()
# 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]
)
data_sample.mean(axis=1).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 * data_sample.mean(axis=1).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.
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 librería de Python para calcular características a partir de datos secuenciales, incluyendo medidas estadísticas, coeficientes de Fourier e información en los dominios del tiempo y la frecuencia.
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 mediante las funciones proporcionadas por la clase ComprehensiveFCParameters
. Esta crea un diccionario que asigna una cadena de texto (el nombre de cada característica) a una lista de diccionarios de parámetros que se utilizarán cuando se llame a la función con ese nombre.
# Configuración por defecto de tsfresh
# ==============================================================================
default_features = ComprehensiveFCParameters()
print("Name of the features extracted by tsfresh:")
print("=========================================")
print('\n'.join(list(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.
# 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 = tsfresh.utilities.dataframe_functions.impute,
n_jobs = 4
)
print("Shape of ts_features:", ts_features.shape)
ts_features.head(3)
Shape of ts_features: (1214, 783)
meter_reading__variance_larger_than_standard_deviation | meter_reading__has_duplicate_max | meter_reading__has_duplicate_min | meter_reading__has_duplicate | meter_reading__sum_values | meter_reading__abs_energy | meter_reading__mean_abs_change | meter_reading__mean_change | meter_reading__mean_second_derivative_central | meter_reading__median | ... | meter_reading__fourier_entropy__bins_5 | meter_reading__fourier_entropy__bins_10 | meter_reading__fourier_entropy__bins_100 | meter_reading__permutation_entropy__dimension_3__tau_1 | meter_reading__permutation_entropy__dimension_4__tau_1 | meter_reading__permutation_entropy__dimension_5__tau_1 | meter_reading__permutation_entropy__dimension_6__tau_1 | meter_reading__permutation_entropy__dimension_7__tau_1 | meter_reading__query_similarity_count__query_None__threshold_0.0 | meter_reading__mean_n_absolute_max__number_of_maxima_7 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
id_1000 | 1.0 | 0.0 | 0.0 | 1.0 | 5.424863e+05 | 8.444178e+08 | 235.994991 | 0.332875 | -0.255837 | 1515.625595 | ... | 0.170467 | 0.260704 | 0.780228 | 1.515622 | 2.520615 | 3.496026 | 4.341908 | 5.045132 | 0.0 | 2144.607384 |
id_1001 | 1.0 | 0.0 | 0.0 | 1.0 | 2.543133e+05 | 4.114736e+08 | 212.298183 | -0.049314 | 0.026097 | 304.500150 | ... | 0.090729 | 0.090729 | 0.700821 | 1.752526 | 3.038890 | 4.327469 | 5.278212 | 5.735413 | 0.0 | 2871.142857 |
id_1002 | 1.0 | 0.0 | 0.0 | 1.0 | 1.409981e+06 | 6.255919e+09 | 960.832035 | 3.226027 | -0.308379 | 3649.003765 | ... | 0.488746 | 0.765782 | 2.650153 | 1.747629 | 3.035159 | 4.431478 | 5.394148 | 5.780678 | 0.0 | 7882.572432 |
3 rows × 783 columns
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 relevantes en cada caso de uso. Para seleccionar las más importantes, tsfresh incluye un proceso de selección automatizado basado en pruebas de hipótesis (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.
⚠ Warning
El proceso de selección utilizado por tsfresh se basa en la significancia de cada característica para predecir con precisión la variable objetivo. Para llevar a cabo este proceso, se requiere una variable objetivo, por ejemplo, el tipo de edificio asociado a una serie temporal determinada. Sin embargo, hay casos en los que no se dispone fácilmente de una variable objetivo. En tales situaciones, se pueden utilizar estrategias alternativas:- En lugar de calcular todas las características estándar, centrarse solo en aquellas que probablemente sean relevantes para la aplicación específica, basándose en el conocimiento experto.
- Excluir características según criterios como baja varianza y alta correlación. Esto ayuda a refinar el conjunto de características a considerar, enfocándose en aquellas que proporcionan información más significativa para el análisis.
- Utilizar técnicas como PCA, t-SNE o autoencoders para reducir la dimensionalidad.
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.
Algunas características creadas pueden tener valores faltantes en algunas de las series temporales. Dado que la mayoría de los algoritmos de clustering no permiten valores faltantes, estos son excluidos.
# Eliminar características con valores faltantes
# ==============================================================================
ts_features = ts_features.dropna(axis=1, how='any')
print(f"Number of features after selection and removing missing values: {ts_features.shape[1]}")
Number of features after selection and removing 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 seleccionadas 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 con PCA
Se utiliza un 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 un PCA para reducir la dimensionalidad de los datos antes de aplicar el k-means.
# 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)
meter_reading__variance_larger_than_standard_deviation | meter_reading__has_duplicate_max | meter_reading__has_duplicate_min | meter_reading__has_duplicate | meter_reading__sum_values | meter_reading__abs_energy | meter_reading__mean_abs_change | meter_reading__mean_change | meter_reading__mean_second_derivative_central | meter_reading__median | ... | meter_reading__fourier_entropy__bins_5 | meter_reading__fourier_entropy__bins_10 | meter_reading__fourier_entropy__bins_100 | meter_reading__permutation_entropy__dimension_3__tau_1 | meter_reading__permutation_entropy__dimension_4__tau_1 | meter_reading__permutation_entropy__dimension_5__tau_1 | meter_reading__permutation_entropy__dimension_6__tau_1 | meter_reading__permutation_entropy__dimension_7__tau_1 | meter_reading__query_similarity_count__query_None__threshold_0.0 | meter_reading__mean_n_absolute_max__number_of_maxima_7 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
id_1000 | 0.028712 | -0.119173 | -0.649265 | 0.933043 | -0.306332 | -0.136331 | -0.144129 | 0.102783 | 0.114077 | -0.301155 | ... | -0.254284 | -0.331247 | -0.741751 | -1.961882 | -1.922583 | -1.994894 | -1.987187 | -1.477824 | 0.0 | -0.303740 |
id_1001 | 0.028712 | -0.119173 | -0.649265 | 0.933043 | -0.406251 | -0.138365 | -0.173691 | 0.010218 | 0.257708 | -0.453833 | ... | -0.744451 | -0.917204 | -0.843649 | 0.680139 | 0.864391 | 0.743247 | 0.601810 | 0.510372 | 0.0 | -0.241759 |
2 rows × 783 columns
# Gráfico de reducción de varianza en función del número de componentes PCA
# ==============================================================================
pca = PCA()
pca.fit(ts_features_scaled)
n_components = np.argmax(np.cumsum(pca.explained_variance_ratio_) > 0.85) + 1
fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(6, 2.5))
ax.plot(np.cumsum(pca.explained_variance_ratio_))
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 68 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)
# 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: 68 Explained variance: 0.851
PC1 | PC2 | PC3 | PC4 | PC5 | PC6 | PC7 | PC8 | PC9 | PC10 | ... | PC59 | PC60 | PC61 | PC62 | PC63 | PC64 | PC65 | PC66 | PC67 | PC68 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
id_1000 | -3.645060 | -1.691436 | 6.209303 | -1.049321 | -0.267084 | 3.060263 | -3.501740 | -3.391411 | 1.606899 | -1.649764 | ... | 1.250706 | 0.045484 | -1.065746 | 0.761409 | 0.175742 | -0.626030 | 1.222963 | 0.521544 | 0.899080 | -0.198204 |
id_1001 | -3.323758 | 4.391786 | -5.026659 | 3.759591 | 4.651089 | -1.628526 | -0.657600 | 3.759905 | -2.050114 | 2.949916 | ... | -0.329982 | 1.179811 | 2.015101 | 1.224391 | -3.262146 | -0.705693 | 1.116380 | 0.101948 | 0.511592 | 0.251012 |
id_1002 | 3.051851 | 1.715388 | 1.847019 | -4.892750 | 2.391237 | 1.259424 | -0.030738 | 0.445865 | 2.549206 | 3.459638 | ... | 1.456168 | -2.346037 | 1.190361 | -0.627113 | -1.784680 | 1.856020 | -0.481244 | 0.719710 | -0.532219 | 3.647449 |
3 rows × 68 columns
Uno de los desafíos inherentes al clustering es determinar el número óptimo de clusters, ya que no existe una verdad fundamental o un número predefinido de clusters en el aprendizaje no supervisado. 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 graficar la suma total de intra-varianza 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 la tasa de disminución de intra-varianza total se estabiliza. Esto sugiere que agregar más clusters más allá de este punto no aporta una mejora sustancial en el rendimiento del clustering.
# Número óptimo de clusters usando 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)
dicrease_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('Intra-cluster variance (inertia)')
dicrease_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 después de 7 clusters, la disminución de la inercia se desacelera, por lo tanto, 7 clusters podrían ser una buena opción.
Además de analizar la evolución de la inercia (intra-varianza), es importante verificar el tamaño de los clusters que se están creando. La presencia de clusters pequeños puede indicar sobreajuste o muestras anómalas que no se ajustan bien a 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): [1214] Size of clusters (n = 2): [1201, 13] Size of clusters (n = 3): [1186, 27, 1] Size of clusters (n = 4): [1140, 4, 4, 66] Size of clusters (n = 5): [1105, 4, 1, 97, 7] Size of clusters (n = 6): [466, 1, 654, 6, 82, 5] Size of clusters (n = 7): [81, 465, 8, 1, 5, 653, 1] Size of clusters (n = 8): [652, 466, 8, 4, 1, 1, 81, 1] Size of clusters (n = 9): [163, 4, 530, 1, 405, 26, 1, 81, 3] Size of clusters (n = 10): [139, 552, 1, 1, 4, 81, 2, 22, 411, 1] Size of clusters (n = 11): [545, 1, 5, 5, 1, 146, 410, 1, 18, 1, 81] Size of clusters (n = 12): [382, 1, 1, 4, 11, 80, 1, 400, 107, 6, 220, 1] Size of clusters (n = 13): [210, 124, 1, 84, 2, 368, 1, 26, 2, 389, 1, 4, 2] Size of clusters (n = 14): [370, 1, 394, 1, 81, 1, 1, 4, 26, 1, 1, 123, 1, 209] Size of clusters (n = 15): [367, 7, 396, 9, 1, 1, 4, 1, 24, 82, 1, 1, 218, 1, 101] Size of clusters (n = 16): [247, 103, 4, 1, 1, 23, 80, 1, 1, 1, 1, 9, 213, 338, 12, 179] Size of clusters (n = 17): [1, 110, 1, 367, 382, 2, 19, 77, 1, 1, 20, 2, 1, 213, 1, 1, 15] Size of clusters (n = 18): [196, 247, 88, 1, 1, 39, 1, 199, 1, 1, 1, 1, 1, 80, 339, 5, 11, 2] Size of clusters (n = 19): [23, 16, 1, 298, 1, 2, 190, 96, 81, 1, 39, 1, 2, 1, 253, 1, 1, 202, 5] Size of clusters (n = 20): [93, 254, 1, 57, 1, 82, 344, 1, 1, 1, 141, 16, 2, 1, 16, 5, 1, 1, 195, 1] Size of clusters (n = 21): [330, 12, 1, 1, 384, 1, 1, 1, 5, 1, 8, 1, 2, 17, 107, 5, 1, 80, 1, 152, 103] Size of clusters (n = 22): [92, 1, 1, 1, 245, 9, 1, 14, 24, 79, 2, 1, 1, 199, 330, 7, 1, 199, 1, 1, 4, 1] Size of clusters (n = 23): [242, 1, 1, 194, 1, 14, 1, 1, 300, 1, 8, 1, 21, 105, 5, 1, 79, 1, 39, 1, 1, 1, 195] Size of clusters (n = 24): [8, 71, 98, 4, 1, 305, 1, 1, 16, 107, 12, 242, 1, 1, 1, 1, 1, 1, 194, 1, 1, 25, 114, 7] Size of clusters (n = 25): [197, 2, 1, 1, 5, 1, 102, 1, 1, 134, 1, 1, 75, 12, 242, 1, 58, 1, 1, 2, 1, 1, 342, 15, 16] Size of clusters (n = 26): [77, 1, 141, 1, 16, 234, 101, 1, 16, 1, 4, 1, 243, 57, 1, 1, 1, 1, 7, 1, 116, 188, 1, 1, 1, 1] Size of clusters (n = 27): [16, 192, 1, 1, 1, 139, 1, 248, 1, 1, 13, 77, 1, 1, 1, 4, 1, 95, 7, 1, 57, 1, 311, 1, 40, 1, 1] Size of clusters (n = 28): [76, 223, 1, 6, 4, 141, 1, 200, 240, 16, 1, 1, 1, 1, 8, 1, 40, 1, 93, 1, 7, 1, 1, 1, 1, 89, 57, 1] Size of clusters (n = 29): [1, 134, 1, 1, 1, 8, 135, 75, 216, 18, 97, 1, 1, 57, 1, 1, 78, 1, 1, 1, 1, 1, 282, 1, 2, 1, 4, 7, 86]
Cuando se crean 7 clusters, el tamaño de los clusters no está bien equilibrado. Antes de modelar las series temporales, los clusters más pequeños (menos de 20 observaciones) se combinan en un solo cluster denominado "otros".
# 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)
kmeans.fit(pca_projections)
clusters = kmeans.predict(pca_projections)
clusters = pd.DataFrame({
'building_id': pca_projections.index,
'cluster_base_on_features': clusters.astype(str)
})
# Merge cluster minor clusters into a single cluster
threshold = 20
cluster_size = clusters['cluster_base_on_features'].value_counts()
cluster_size = cluster_size[cluster_size < threshold].index.tolist()
clusters['cluster_base_on_features'] = np.where(
clusters['cluster_base_on_features'].isin(cluster_size),
'Other',
clusters['cluster_base_on_features']
)
display(clusters.head(3))
clusters['cluster_base_on_features'].value_counts().to_frame(name='Number of series')
building_id | cluster_base_on_features | |
---|---|---|
0 | id_1000 | 5 |
1 | id_1001 | 1 |
2 | id_1002 | 5 |
Number of series | |
---|---|
cluster_base_on_features | |
5 | 653 |
1 | 465 |
0 | 81 |
Other | 15 |
# Añadir el cluster predicho al dataframe con la información de los edificios
# ==============================================================================
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)
building_id | meter_reading | site_id | primary_use | square_feet | air_temperature | dew_temperature | sea_level_pressure | wind_direction | wind_speed | cluster_base_on_features | |
---|---|---|---|---|---|---|---|---|---|---|---|
timestamp | |||||||||||
2016-01-01 | id_1000 | 1012.500694 | 10 | Education | 121146 | -8.9 | NaN | NaN | NaN | 3.1 | 5 |
2016-01-01 | id_970 | 0.000000 | 9 | Other | 346056 | 7.8 | NaN | NaN | NaN | 3.1 | 5 |
2016-01-01 | id_1066 | 639.270004 | 12 | Education | 55800 | 1.9 | -1.2 | 1016.2 | 200.0 | 5.0 | 5 |
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 edificios por cluster
# ==============================================================================
primary_usage_per_cluster = (
data
.groupby('cluster_base_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
primary_use | Education | Entertainment/public assembly | Lodging/residential | Office | Other | Public services |
---|---|---|---|---|---|---|
cluster_base_on_features | ||||||
0 | 59.3 | 11.1 | 1.2 | 16.0 | 2.5 | 9.9 |
1 | 26.0 | 15.3 | 18.7 | 9.2 | 11.6 | 19.1 |
5 | 43.8 | 11.6 | 3.7 | 25.4 | 7.4 | 8.1 |
Other | 53.3 | 6.7 | 0.0 | 40.0 | 0.0 | 0.0 |
Los resultados sugieren (la tabla debe leerse horizontalmente) que el proceso de clustering basado en las características extraídas de las series temporales genera grupos que difieren de aquellos formados por el propósito principal del edificio.
Clustering basado en Dynamic Time Warping (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 de manera óptima. El clustering con DTW implica agrupar los datos de series temporales en función de sus distancias dinámicamente ajustadas, 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.
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
meter_reading | ||
---|---|---|
building_id | timestamp | |
id_1000 | 2016-01-01 | 1012.500694 |
2016-01-02 | 1158.500099 | |
2016-01-03 | 983.000099 | |
2016-01-04 | 1675.750496 | |
2016-01-05 | 1586.250694 | |
... | ... | ... |
id_999 | 2016-12-27 | 2627.250000 |
2016-12-28 | 2667.250000 | |
2016-12-29 | 2495.000000 | |
2016-12-30 | 2059.500000 | |
2016-12-31 | 1899.500000 |
444324 rows × 1 columns
⚠ Warning
La siguiente celda tiene un tiempo de ejecución de aproximadamente 10 minutos.# Fit clustering
# ==============================================================================
model = TimeSeriesKMeans(
n_clusters = 4,
metric = "dtw",
max_iter = 10,
random_state = 123456
)
model.fit(data_long)
TimeSeriesKMeans(max_iter=10, n_clusters=4, random_state=123456)Please rerun this cell to show the HTML repr or trust the notebook.
TimeSeriesKMeans(max_iter=10, n_clusters=4, random_state=123456)
# Predocción de clusters
# ==============================================================================
clusters = model.predict(data_long)
clusters = pd.DataFrame({
'building_id': data_long.index.get_level_values('building_id').drop_duplicates(),
'cluster_base_on_dtw': clusters.astype(str),
})
# Size of each cluster
clusters['cluster_base_on_dtw'].value_counts().sort_index().to_frame(name='Number of series')
Number of series | |
---|---|
cluster_base_on_dtw | |
0 | 22 |
1 | 717 |
2 | 347 |
3 | 128 |
# Añadir el cluster predicho al dataframe con la información de los edificios
# ==============================================================================
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)
building_id | meter_reading | site_id | primary_use | square_feet | air_temperature | dew_temperature | sea_level_pressure | wind_direction | wind_speed | cluster_base_on_features | cluster_base_on_dtw | |
---|---|---|---|---|---|---|---|---|---|---|---|---|
timestamp | ||||||||||||
2016-01-01 | id_1000 | 1012.500694 | 10 | Education | 121146 | -8.9 | NaN | NaN | NaN | 3.1 | 5 | 1 |
2016-01-01 | id_970 | 0.000000 | 9 | Other | 346056 | 7.8 | NaN | NaN | NaN | 3.1 | 5 | 1 |
2016-01-01 | id_1066 | 639.270004 | 12 | Education | 55800 | 1.9 | -1.2 | 1016.2 | 200.0 | 5.0 | 5 | 1 |
✎ Note
Sktime ofrece algoritmos de clustering adicionales, como TimeSeriesKMeansTslearn, TimeSeriesKMedoids y TimeSeriesKShapes, que valen la pena explorar. Otras excelentes bibliotecas para clustering de series temporales son: + DTAIDistance + Tslearn + AEON# 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 agrupamiento (uso del edificio, clustering basado en características de series temporales y clustering basado en dynamic time warping), se entrenan modelos globales multiserie. La evaluación posterior se centra en determinar qué tan efectivos son estos modelos para predecir los datos diarios de las dos últimas semanas. Durante esta evaluación, se utilizan tres métricas distintas:
El valor promedio del Error Absoluto Medio (MAE) para todos los edificios.
La suma absoluta de los errores, es decir, la desviación absoluta entre el valor predicho y el consumo real para todos los edificios.
El bias para las predicciones sumado sobre todos los edificios.
Además de los valores rezagados de las propias series temporales, se incluyen como variables exógenas el día de la semana (codificado en seno-coseno), la temperatura exterior y la velocidad del viento.
✎ Note
Para una explicación más detallada sobre la validación de modelos de series temporales, se recomienda a los lectores consultar la guía del usuario de Backtesting. Para más información sobre variables del calendario y codificación cíclica, visita Variables del calendario y Variables cíclicas en series temporales.Para entrenar los modelos, buscar los hiperparámetros óptimos y evaluar su rendimiento predictivo, los datos se dividen en tres conjuntos separados: entrenamiento, validación y prueba.
# Lectura de datos para modelado
# ==============================================================================
data = pd.read_parquet('data_modelling.parquet')
# División de los datos en entrenamiento, validación y test
# ==============================================================================
end_train = '2016-07-31 23:59:00'
end_validation = '2016-09-30 23:59:00'
data_train = data.loc[: end_train, :].copy()
data_val = data.loc[end_train:end_validation, :].copy()
data_test = data.loc[end_validation:, :].copy()
print(
f"Rage of dates available : {data.index.min()} --- {data.index.max()} "
f"(n_days={(data.index.max() - data.index.min()).days})"
)
print(
f" Dates for training : {data_train.index.min()} --- {data_train.index.max()} "
f"(n_days={(data_train.index.max() - data_train.index.min()).days})"
)
print(
f" Dates for validation : {data_val.index.min()} --- {data_val.index.max()} "
f"(n_days={(data_val.index.max() - data_val.index.min()).days})"
)
print(
f" Dates for test : {data_test.index.min()} --- {data_test.index.max()} "
f"(n_days={(data_test.index.max() - data_test.index.min()).days})"
)
Rage of dates available : 2016-01-01 00:00:00 --- 2016-12-31 00:00:00 (n_days=365) Dates for training : 2016-01-01 00:00:00 --- 2016-07-31 00:00:00 (n_days=212) Dates for validation : 2016-08-01 00:00:00 --- 2016-09-30 00:00:00 (n_days=60) Dates for test : 2016-10-01 00:00:00 --- 2016-12-31 00:00:00 (n_days=91)
# Tabla de resultados
# ==============================================================================
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})
# Función para añadir información del día de la semana
# ==============================================================================
def day_of_week_cyclical_encoding(data):
"""
Calculate cyclical encoding for the day of week using a
DataFrame with a datetime index.
"""
data['day_of_week'] = data.index.dayofweek
data['sin_day_of_week'] = np.sin(2*np.pi*data['day_of_week']/7)
data['cos_day_of_week'] = np.cos(2*np.pi*data['day_of_week']/7)
return data
# Variables exógenas incluidas en el modelo
# ==============================================================================
exog_features = ['sin_day_of_week', 'cos_day_of_week', 'air_temperature', 'wind_speed']
Modelo global para todos los edificios
Se entrena y prueba un modelo global para todos los edificios utilizando la clase ForecasterRecursiveMultiSeries
de skforecast. Este forecaster permite al usuario pasar las series temporales en varios formatos, incluyendo un DataFrame con las series temporales organizadas como columnas.
Para más información sobre cómo utilizar series de diferentes longitudes o variables exógenas diferentes para cada serie, consulta Global Forecasting Models.
# Transformar los datos para tener una columna por edificio
# ==============================================================================
data_pivot = data.pivot_table(
index = 'timestamp',
columns = 'building_id',
values = 'meter_reading',
aggfunc = 'mean'
)
data_pivot.columns = data_pivot.columns.astype(str)
data_pivot = data_pivot.asfreq('D').sort_index()
# Add calendar features
data_pivot = day_of_week_cyclical_encoding(data_pivot)
# Add temperature and wind speed as exogenous variables
data_pivot = data_pivot.merge(
data[['air_temperature', 'wind_speed']].resample('D').mean(),
left_index = True,
right_index = True,
how = 'left',
validate = '1:m'
)
data_pivot.head(3)
id_1000 | id_1001 | id_1002 | id_1003 | id_1004 | id_1007 | id_1008 | id_1009 | id_1010 | id_1013 | ... | id_993 | id_995 | id_996 | id_998 | id_999 | day_of_week | sin_day_of_week | cos_day_of_week | air_temperature | wind_speed | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
timestamp | |||||||||||||||||||||
2016-01-01 | 1012.500694 | 142.999700 | 3114.500198 | 2294.750893 | 7188.003021 | 682.750099 | 2007.500893 | 507.630098 | 109.000300 | 996.750298 | ... | 27809.0 | 585.0 | 2643.0 | 1757.750298 | 2401.250595 | 4 | -0.433884 | -0.900969 | 6.416639 | 4.040115 |
2016-01-02 | 1158.500099 | 141.000801 | 4110.000000 | 1750.000198 | 8494.001007 | 661.750000 | 2025.250793 | 593.370001 | 111.000000 | 947.750504 | ... | 29067.0 | 639.0 | 2771.0 | 1756.000496 | 2244.370499 | 5 | -0.974928 | -0.222521 | 6.366474 | 4.530395 |
2016-01-03 | 983.000099 | 137.000300 | 2965.000000 | 1455.750492 | 7130.001007 | 713.250099 | 1922.000496 | 561.630098 | 111.000301 | 940.180004 | ... | 28937.0 | 605.0 | 2705.0 | 1716.000893 | 1878.750893 | 6 | -0.781831 | 0.623490 | 6.555272 | 3.273064 |
3 rows × 1219 columns
# Forecaster multiseries para modelar todos los edificios a la vez
# ==============================================================================
exog_features = ['sin_day_of_week', 'cos_day_of_week', 'air_temperature', 'wind_speed']
params_lgbm = {
'n_estimators': 500,
'learning_rate': 0.01,
'max_depth': 10,
'random_state': 8520,
'verbose': -1
}
forecaster = ForecasterRecursiveMultiSeries(
regressor = LGBMRegressor(**params_lgbm),
lags = 31,
encoding = 'ordinal_category'
)
start = datetime.now()
cv = TimeSeriesFold(
steps = 7,
initial_train_size = len(data_pivot.loc[:end_validation, :]),
refit = False
)
metric, predictions = backtesting_forecaster_multiseries(
forecaster = forecaster,
series = data_pivot.filter(regex='^id_'), # Select only buildings
exog = data_pivot[exog_features],
cv = cv,
metric = 'mean_absolute_error',
add_aggregated_metric = False,
)
end = datetime.now()
predictions = predictions.pivot(columns="level", values="pred")
mean_metric_all_buildings = metric['mean_absolute_error'].mean()
errors_all_buildings = predictions - data_pivot.loc[predictions.index, predictions.columns]
sum_abs_errors_all_buildings = errors_all_buildings.abs().sum().sum()
sum_bias_all_buildings = errors_all_buildings.sum().sum()
table_results.loc['Global model', 'mae'] = mean_metric_all_buildings
table_results.loc['Global model', 'abs_error'] = sum_abs_errors_all_buildings
table_results.loc['Global model', 'bias'] = sum_bias_all_buildings
table_results.loc['Global model', 'elapsed_time'] = end - start
print("")
print(f"Average mean absolute error for all buildings: {mean_metric_all_buildings:.0f}")
print(f"Sum of absolute errors for all buildings (x 10,000): {sum_abs_errors_all_buildings/10000:.0f}")
print(f"Bias (x 10,000): {sum_bias_all_buildings/10000:.0f}")
0%| | 0/14 [00:00<?, ?it/s]
Average mean absolute error for all buildings: 419 Sum of absolute errors for all buildings (x 10,000): 4684 Bias (x 10,000): 1205
# Gráfico de predicciones vs valor real para 2 edificios seleccionados al azar
# ==============================================================================
rng = np.random.default_rng(147)
selected_buildings = rng.choice(data['building_id'].unique(), size=2, replace=False)
fig, axs = plt.subplots(2, 1, figsize=(6, 4), sharex=True)
axs = axs.flatten()
for i, building in enumerate(selected_buildings):
data_test.query("building_id == @building")['meter_reading'].plot(ax=axs[i], label='test')
predictions[building].plot(ax=axs[i], label='One model per building')
axs[i].set_title(f"Building {building}", fontsize=10)
axs[i].set_xlabel("")
axs[i].legend()
fig.tight_layout()
plt.show();
Modelo global por uso del edificio
# Forecaster multi-series para los edificios agrupados por uso principal
# ==============================================================================
predictions_all_buildings = []
metrics_all_buildings = []
params_lgbm = {
'n_estimators': 500,
'learning_rate': 0.01,
'max_depth': 10,
'random_state': 8520,
'verbose': -1
}
forecaster = ForecasterRecursiveMultiSeries(
regressor = LGBMRegressor(**params_lgbm),
lags = 31,
encoding = "ordinal_category"
)
start = datetime.now()
for primary_usage in data['primary_use'].unique():
print(
f"Training and testing model for primary usage: {primary_usage} "
f"(n = {data[data['primary_use'] == primary_usage]['building_id'].nunique()})"
)
# Create a subset based on primary use clusters
data_subset = data[data['primary_use'] == primary_usage]
data_subset = data_subset.pivot_table(
index = 'timestamp',
columns = 'building_id',
values = 'meter_reading',
aggfunc = 'mean'
)
data_subset.columns = data_subset.columns.astype(str)
data_subset = data_subset.asfreq('D').sort_index()
# Add calendar features
data_subset = day_of_week_cyclical_encoding(data_subset)
data_subset = data_subset.merge(
data[['air_temperature', 'wind_speed']].resample('D').mean(),
left_index = True,
right_index = True,
how = 'left',
validate = '1:m'
)
cv = TimeSeriesFold(
steps = 7,
initial_train_size = len(data_subset.loc[:end_validation, :]),
refit = False
)
metric, predictions = backtesting_forecaster_multiseries(
forecaster = forecaster,
series = data_subset.filter(regex='^id_'),
exog = data_subset[['sin_day_of_week', 'cos_day_of_week']],
cv = cv,
metric = 'mean_absolute_error',
add_aggregated_metric = False,
show_progress = False
)
predictions_all_buildings.append(predictions)
metrics_all_buildings.append(metric)
end = datetime.now()
predictions_all_buildings = pd.concat(predictions_all_buildings, axis=0)
predictions_all_buildings = predictions_all_buildings.pivot(columns="level", values="pred")
metrics_all_buildings = pd.concat(metrics_all_buildings, axis=0)
errors_all_buildings = (
predictions_all_buildings - data_pivot.loc[predictions_all_buildings.index, predictions_all_buildings.columns]
)
mean_metric_all_buildings = metrics_all_buildings['mean_absolute_error'].mean()
sum_abs_errors_all_buildings = errors_all_buildings.abs().sum().sum()
sum_bias_all_buildings = errors_all_buildings.sum().sum()
table_results.loc['Global model per primary usage', 'mae'] = mean_metric_all_buildings
table_results.loc['Global model per primary usage', 'abs_error'] = sum_abs_errors_all_buildings
table_results.loc['Global model per primary usage', 'bias'] = sum_bias_all_buildings
table_results.loc['Global model per primary usage', 'elapsed_time'] = end - start
print("")
print(f"Average mean absolute error for all buildings: {mean_metric_all_buildings:.0f}")
print(f"Sum of absolute errors for all buildings (x 10,000): {sum_abs_errors_all_buildings/10000:.0f}")
print(f"Bias (x 10,000): {sum_bias_all_buildings/10000:.0f}")
Training and testing model for primary usage: Education (n = 463) Training and testing model for primary usage: Other (n = 104) Training and testing model for primary usage: Entertainment/public assembly (n = 157) Training and testing model for primary usage: Office (n = 228) Training and testing model for primary usage: Lodging/residential (n = 112) Training and testing model for primary usage: Public services (n = 150) Average mean absolute error for all buildings: 405 Sum of absolute errors for all buildings (x 10,000): 4523 Bias (x 10,000): 1149
# Añadir predicciones al gráfico ya existente (no mostrar el gráfico todavía)
# ==============================================================================
for i, building in enumerate(selected_buildings):
predictions_all_buildings[building].plot(ax=axs[i], label='Global model per primary usage')
axs[i].legend()
Modelo global por cluster (características)
Se entrena y evalúa un modelo global para cada cluster basado en las características de las series temporales. Para ello, se utilizan las mismas funciones y clases que en el caso anterior, pero con las series temporales agrupadas por cluster.
# Forecaster multiserie para los edificios agrupados por características de las series
# ==============================================================================
predictions_all_buildings = []
metrics_all_buildings = []
params_lgbm = {
'n_estimators': 500,
'learning_rate': 0.01,
'max_depth': 10,
'random_state': 8520,
'verbose': -1
}
forecaster = ForecasterRecursiveMultiSeries(
regressor = LGBMRegressor(**params_lgbm),
lags = 31,
encoding = "ordinal_category"
)
start = datetime.now()
for cluster in data['cluster_base_on_features'].unique():
print(
f"Training and testing model for cluster: {cluster} "
f"(n = {data[data['cluster_base_on_features'] == cluster]['building_id'].nunique()})"
)
# Create subset based on features clusters
data_subset = data[data['cluster_base_on_features'] == cluster]
data_subset = data_subset.pivot_table(
index = 'timestamp',
columns = 'building_id',
values = 'meter_reading',
aggfunc = 'mean'
)
data_subset.columns = data_subset.columns.astype(str)
data_subset = data_subset.asfreq('D').sort_index()
# Add calendar features
data_subset = day_of_week_cyclical_encoding(data_subset)
data_subset = data_subset.merge(
data[['air_temperature', 'wind_speed']].resample('D').mean(),
left_index = True,
right_index = True,
how = 'left',
validate = '1:m'
)
cv = TimeSeriesFold(
steps = 7,
initial_train_size = len(data_subset.loc[:end_validation, :]),
refit = False
)
metric, predictions = backtesting_forecaster_multiseries(
forecaster = forecaster,
series = data_subset.filter(regex='^id_'),
exog = data_subset[['sin_day_of_week', 'cos_day_of_week']],
cv = cv,
metric = 'mean_absolute_error',
add_aggregated_metric = False,
show_progress = False
)
predictions_all_buildings.append(predictions)
metrics_all_buildings.append(metric)
end = datetime.now()
predictions_all_buildings = pd.concat(predictions_all_buildings, axis=0)
predictions_all_buildings = predictions_all_buildings.pivot(columns="level", values="pred")
metrics_all_buildings = pd.concat(metrics_all_buildings, axis=0)
errors_all_buildings = (
predictions_all_buildings - data_pivot.loc[predictions_all_buildings.index, predictions_all_buildings.columns]
)
mean_metric_all_buildings = metrics_all_buildings['mean_absolute_error'].mean()
sum_abs_errors_all_buildings = errors_all_buildings.abs().sum().sum()
sum_bias_all_buildings = errors_all_buildings.sum().sum()
table_results.loc['Global model per cluster (features)', 'mae'] = mean_metric_all_buildings
table_results.loc['Global model per cluster (features)', 'abs_error'] = sum_abs_errors_all_buildings
table_results.loc['Global model per cluster (features)', 'bias'] = sum_bias_all_buildings
table_results.loc['Global model per cluster (features)', 'elapsed_time'] = end - start
print("")
print(f"Average mean absolute error for all buildings: {mean_metric_all_buildings:.0f}")
print(f"Sum of absolute errors for all buildings (x 10,000): {sum_abs_errors_all_buildings/10000:.0f}")
print(f"Bias (x 10,000): {sum_bias_all_buildings/10000:.0f}")
Training and testing model for cluster: 5 (n = 653) Training and testing model for cluster: 0 (n = 81) Training and testing model for cluster: 1 (n = 465) Training and testing model for cluster: Other (n = 15) Average mean absolute error for all buildings: 383 Sum of absolute errors for all buildings (x 10,000): 4274 Bias (x 10,000): 1072
# Añaadir predicciones al gráfico ya existente (no mostrar el gráfico todavía)
# ==============================================================================
for i, building in enumerate(selected_buildings):
predictions_all_buildings[building].plot(ax=axs[i], label='Global model per cluster (features)')
axs[i].legend()
Modelo global por cluster (DTW)
Se entrena y evalúa un modelo global para cada cluster basado en la distancia elástica (DTW).
# Forecaster multiserie para los edificios agrupados por DTW
# ==============================================================================
predictions_all_buildings = []
metrics_all_buildings = []
params_lgbm = {
'n_estimators': 500,
'learning_rate': 0.01,
'max_depth': 10,
'random_state': 8520,
'verbose': -1
}
forecaster = ForecasterRecursiveMultiSeries(
regressor = LGBMRegressor(**params_lgbm),
lags = 31,
encoding = 'ordinal_category'
)
start = datetime.now()
for cluster in data['cluster_base_on_dtw'].unique():
print(
f"Training and testing model for cluster: {cluster} "
f"(n = {data[data['cluster_base_on_dtw'] == cluster]['building_id'].nunique()})"
)
# Create subset based on DTW clusters
data_subset = data[data['cluster_base_on_dtw'] == cluster]
data_subset = data_subset.pivot_table(
index = 'timestamp',
columns = 'building_id',
values = 'meter_reading',
aggfunc = 'mean'
)
data_subset.columns = data_subset.columns.astype(str)
data_subset = data_subset.asfreq('D').sort_index()
# Add calendar features
data_subset = day_of_week_cyclical_encoding(data_subset)
data_subset = data_subset.merge(
data[['air_temperature', 'wind_speed']].resample('D').mean(),
left_index = True,
right_index = True,
how = 'left',
validate = '1:m'
)
cv = TimeSeriesFold(
steps = 7,
initial_train_size = len(data_subset.loc[:end_validation, :]),
refit = False
)
metric, predictions = backtesting_forecaster_multiseries(
forecaster = forecaster,
series = data_subset.filter(regex='^id_'),
exog = data_subset[['sin_day_of_week', 'cos_day_of_week']],
cv = cv,
metric = 'mean_absolute_error',
add_aggregated_metric = False,
show_progress = False
)
predictions_all_buildings.append(predictions)
metrics_all_buildings.append(metric)
end = datetime.now()
predictions_all_buildings = pd.concat(predictions_all_buildings, axis=0)
predictions_all_buildings = predictions_all_buildings.pivot(columns="level", values="pred")
metrics_all_buildings = pd.concat(metrics_all_buildings, axis=0)
errors_all_buildings = (
predictions_all_buildings - data_pivot.loc[predictions_all_buildings.index, predictions_all_buildings.columns]
)
mean_metric_all_buildings = metrics_all_buildings['mean_absolute_error'].mean()
sum_abs_errors_all_buildings = errors_all_buildings.abs().sum().sum()
sum_bias_all_buildings = errors_all_buildings.sum().sum()
table_results.loc['Global model per cluster (DTW)', 'mae'] = mean_metric_all_buildings
table_results.loc['Global model per cluster (DTW)', 'abs_error'] = sum_abs_errors_all_buildings
table_results.loc['Global model per cluster (DTW)', 'bias'] = sum_bias_all_buildings
table_results.loc['Global model per cluster (DTW)', 'elapsed_time'] = end - start
print("")
print(f"Average mean absolute error for all buildings: {mean_metric_all_buildings:.0f}")
print(f"Sum of absolute errors for all buildings (x 10,000): {sum_abs_errors_all_buildings/10000:.0f}")
print(f"Bias (x 10,000): {sum_bias_all_buildings/10000:.0f}")
Training and testing model for cluster: 1 (n = 717) Training and testing model for cluster: 2 (n = 347) Training and testing model for cluster: 3 (n = 128) Training and testing model for cluster: 0 (n = 22) Average mean absolute error for all buildings: 377 Sum of absolute errors for all buildings (x 10,000): 4207 Bias (x 10,000): 1044
# Añaadir predicciones al gráfico ya existente (no mostrar el gráfico todavía)
# ==============================================================================
for i, building in enumerate(selected_buildings):
predictions_all_buildings[building].plot(ax=axs[i], label='Global model per cluster (DTW)')
axs[i].legend()
Resultados
# Tabla de resultados
# ==============================================================================
table_results['elapsed_time'] = table_results['elapsed_time'].astype(str).str[:7]
table_results.style.highlight_min(axis=0, color='green').format(precision=0)
mae | abs_error | bias | elapsed_time | |
---|---|---|---|---|
model | ||||
Global model | 419 | 46842049 | 12052318 | 0:00:33 |
Global model per primary usage | 405 | 45233247 | 11487379 | 0:00:31 |
Global model per cluster (features) | 383 | 42744087 | 10720536 | 0:00:32 |
Global model per cluster (DTW) | 377 | 42073831 | 10441483 | 0:00:32 |
# 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.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. Al agrupar series temporales con patrones similares, se puede reducir el impacto del ruido, lo que lleva a predicciones más precisas. Los resultados de este estudio demuestran que la elección del algoritmo de clustering puede afectar significativamente el rendimiento del modelo de forecasting. En este caso, la distancia dynamic time warping (DTW) superó al método de clustering basado en características y al uso principal del edificio como criterio de agrupación.
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 del usuario de skforecast para variables del calendario y luz solar que suelen afectar el consumo de energía.
Optimizar lags e hiperparámetros: utiliza búsqueda diferentes estrategias de búsqueda para encontrar la mejor configuración del modelo.
Probar otros algoritmos de machine learning.
Información de sesión
import session_info
session_info.show(html=False)
----- lightgbm 4.6.0 matplotlib 3.10.1 numpy 2.0.2 pandas 2.2.3 session_info 1.0.0 skforecast 0.15.1 sklearn 1.5.2 sktime 0.36.0 tsfresh 0.21.0 ----- IPython 9.0.2 jupyter_client 8.6.3 jupyter_core 5.7.2 notebook 6.4.12 ----- Python 3.12.9 | packaged by Anaconda, Inc. | (main, Feb 6 2025, 18:56:27) [GCC 11.2.0] Linux-5.15.0-1077-aws-x86_64-with-glibc2.31 ----- Session information updated at 2025-03-28 15:29
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://www.cienciadedatos.net/documentos/py64-clustering-series-temporales-forecasting.ipynb
¿Cómo citar skforecast?
Si utilizas skforecast, te agradeceríamos mucho que lo cites. ¡Muchas gracias!
Zenodo:
Amat Rodrigo, Joaquin, & Escobar Ortiz, Javier. (2025). skforecast (v0.15.1). Zenodo. https://doi.org/10.5281/zenodo.8382788
APA:
Amat Rodrigo, J., & Escobar Ortiz, J. (2025). skforecast (Version 0.15.1) [Computer software]. https://doi.org/10.5281/zenodo.8382788
BibTeX:
@software{skforecast, author = {Amat Rodrigo, Joaquin and Escobar Ortiz, Javier}, title = {skforecast}, version = {0.15.1}, month = {03}, year = {2025}, 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.