Si te gusta Skforecast , ayúdanos dándonos una estrella en GitHub! ⭐️
Modelos de forecasting globales: Análisis comparativo de modelos de una y múltiples series
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
- Forecasting series temporales con XGBoost
- Global Forecasting: Multi-series forecasting
- Forecasting de la demanda eléctrica con machine learning
- Modelos de forecasting globales: Análisis comparativo de modelos de una y múltiples series
- 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¶
En casos de uso en los que se necesita predecir cientos o miles de series temporales ¿se debe desarrollar un modelo individual para cada serie o un único modelo capaz de predecir todas las series a la vez?
En los modelo de forecasting locales, se crea un modelo independiente para cada serie temporal. Aunque este método proporciona una comprensión exhaustiva de cada serie, su escalabilidad puede verse dificultada por la necesidad de crear y mantener cientos o miles de modelos.
El modelo de forecasting global consiste en crear un único modelo que tenga en cuenta todas las series temporales simultáneamente. Intenta captar los patrones comunes que rigen las series, mitigando así el ruido que pueda introducir cada serie. Este enfoque es eficiente desde el punto de vista computacional, fácil de mantener y puede producir generalizaciones más sólidas, aunque potencialmente a costa de sacrificar algunos aspectos individuales.
Para ayudar a comprender las ventajas de cada estrategia de forecasting, este documento examina y compara los resultados obtenidos al predecir el consumo energético de más de mil edificios utilizando el conjunto de datos ASHRAE - Great Energy Predictor III disponible en Kaggle.
El forecasting con modelos globales parte de la base de que las series que se comportan de forma similar pueden beneficiarse de ser modelizadas conjuntamente. Aunque el uso principal de cada edificio está disponible en el conjunto de datos, puede que no refleje grupos con patrones similares de consumo de energía, por lo que se crean grupos adicionales utilizando métodos de clustering. Se realizan un total de 5 experimentos:
Forecasting individual de cada edificio.
Forecasting de todos los edificios juntos con una única estrategia de modelo global.
Forecasting de grupos de edificios en función de su uso principal (un modelo global por uso principal).
Forecasting de grupos de edificios basada en la agrupación de características de series temporales (un modelo global por agrupación).
Forecasting de grupos de edificios basada en la agrupación por Dynamic Time Warping (DTW) (un modelo global por agrupación).
El consumo de energía de cada edificio se predice semanalmente (resolución diaria) durante 13 semanas siguiendo cada estrategia. La eficacia de cada enfoque se evalúa utilizando varias métricas de rendimiento, como el error absoluto medio (MAE), el error absoluto y el sesgo. El objetivo del estudio es identificar el enfoque más eficaz tanto para las predicciones globales como para un grupo específico de edificios.
✎ Note
Si prefieres una visión general rápida antes de sumergirte en los detalles, considera comenzar con la sección de Conclusiones. Este enfoque te permite adaptar tu lectura a tus intereses y limitaciones de tiempo, y resume nuestros hallazgos e ideas. Después de leer las conclusiones, es posible que encuentres ciertas secciones particularmente relevantes o interesantes. Siéntete libre de navegar directamente a esas partes del artículo para una comprensión más profunda.✎ Note
Este documento se centra en un caso de uso. Para una descripción detallada de las muchas características que skforecast proporciona para la construcción de modelos de forecasting globales, consulta Modelos de forecasting globales: Modelado de múltiples series temporales con machine learning.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
plt.style.use('seaborn-v0_8-darkgrid')
from tqdm.auto import tqdm
# Forecasting
# ==============================================================================
import lightgbm
from lightgbm import LGBMRegressor
from sklearn.preprocessing import StandardScaler
import skforecast
from skforecast.recursive import ForecasterRecursive
from skforecast.model_selection import backtesting_forecaster
from skforecast.recursive import ForecasterRecursiveMultiSeries
from skforecast.model_selection import TimeSeriesFold
from skforecast.model_selection import backtesting_forecaster_multiseries
# Feature engineering
# ==============================================================================
import tsfresh
from tsfresh import extract_features
from tsfresh import select_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 tslearn
from tslearn.clustering import TimeSeriesKMeans
# Warnings
# ==============================================================================
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 tslearn: {tslearn.__version__}")
print(f"{color}Versión pandas: {pd.__version__}")
print(f"{color}Versión numpy: {np.__version__}")
⚠ Warning
En el momento de escribir este documento,tslearn solo es compatible con versiones de numpy inferiores a 2.0. Si tienes una versión superior, puedes reducirla ejecutando el siguiente comando: pip install numpy==1.26.4
Datos¶
Los datos utilizados en este documento se han obtenido 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.
Tres archivos se utilizan para crear el conjunto de datos de modelado:
weather_train.csv y weather_test.csv: Estos archivos contienen datos relacionados con el clima para cada edificio, incluida la temperatura del aire exterior, la temperatura de rocío, la humedad relativa y otros parámetros meteorológicos. Los datos meteorológicos son cruciales para comprender el impacto de las condiciones externas en el uso de energía de los edificios.
building_metadata.csv: este archivo proporciona metadatos para cada edificio en el conjunto de datos, como el tipo de edificio, el uso principal, el metraje cuadrado, el número de pisos y el año de construcción. Esta información ayuda a comprender las características de los edificios y su posible influencia en los patrones de consumo de energía.
train.csv: el conjunto de datos de entrenamiento contiene la variable objetivo, es decir, los datos de consumo de energía para cada edificio, junto con las fechas de las lecturas de consumo de energía. También incluye los identificadores de edificio y clima correspondientes para vincular la información en diferentes conjuntos de datos.
Los tres archivos se han preprocesado para eliminar los edificios con menos del 85% de valores diferentes de NaN o cero, para utilizar solo el medidor de electricidad y para agregar los datos a una frecuencia diaria.
# Lectura de datos
# ==============================================================================
data = fetch_dataset('ashrae_daily')
data.head()
# Imputar valores faltantes de air_temperature y wind_speed usando fowrard y backward fill
# ==============================================================================
# La imputación debe hacerse por separado para cada edificio
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"Rango de fechas disponibles : {data.index.min()} --- {data.index.max()} "
f"(n_días={(data.index.max() - data.index.min()).days})"
)
Análisis exploratorio de datos¶
Uso principal de los edificios¶
Uno de los atributos clave asociados con cada edificio es su uso designado. Esta característica puede desempeñar un papel crucial en la influencia del patrón de consumo de energía, ya que los usos distintos 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"Número de edificios: {n_building}")
print(f"Número de tipos de edificios: {n_type_building}")
display(data.drop_duplicates(subset=['building_id'])['primary_use'].value_counts())
Para algunas categorías de uso principal, hay un número limitado de edificios dentro del conjunto de datos. Para simplificar el análisis, las categorías con menos de 100 edificios se agrupan en la categoría "Other".
# Tipos de edificios (primary use) con menos de 100 muestras se agrupan como "Other".
# ==============================================================================
infrequent_categories = (
data
.drop_duplicates(subset=['building_id'])['primary_use']
.value_counts()
.loc[lambda x: x < 100]
.index
.tolist()
)
print("Categorias poco frecuentes:")
print("===========================")
print('\n'.join(infrequent_categories))
data['primary_use'] = np.where(
data['primary_use'].isin(infrequent_categories),
'Other',
data['primary_use']
)
A continuación, se crea un gráfico que muestra el consumo de energía para un edificio seleccionado al azar dentro de cada categoría respectiva, y un gráfico de todas las series temporales disponibles para cada categoría.
# Series temporales 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"Edificio: {building_id}, tipo: {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('Consumo de energía para 6 edificios aleatorios', fontsize=12)
fig.tight_layout()
plt.show()
# Consumo de energía 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"Tipo: {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("")
axs[i].ticklabel_format(axis='y', style='sci', scilimits=(0,0))
axs[i].title.set_size(9)
# Limitar el eje y a 5 veces el valor medio máximo para mejorar la visualización
axs[i].set_ylim(
bottom = 0,
top = 5 * data_sample.mean(axis=1).max()
)
fig.tight_layout()
fig.suptitle('Consumo de energía por tipo de edificio', fontsize=12)
plt.subplots_adjust(top=0.9)
plt.show()
El gráfico revela que hay una variabilidad notable en los patrones de consumo entre edificios del mismo propósito. Esto sugiere que puede haber margen para mejorar los criterios por los que se agrupan los edificios.
Clustering por los patrones de consumo de energía¶
La idea detrás de modelar múltiples series de forma conjunta es poder capturar los patrones principales que rigen las series, reduciendo así el impacto del ruido que cada serie pueda tener. Esto significa que las series que se comportan de manera similar pueden beneficiarse de ser modelizadas juntas. Una forma de identificar posibles grupos de series es realizar un estudio de clustering antes de modelizar. Si como resultado del clustering se identifican grupos claros, es apropiado modelar cada uno de ellos por separado.
El clustering es una técnica de análisis no supervisado que agrupa un conjunto de observaciones en clústeres que contienen observaciones consideradas homogéneas, mientras que las observaciones en diferentes clústeres se consideran heterogéneas. Los algoritmos que agrupan series temporales se pueden dividir en dos grupos: aquellos que utilizan una transformación para crear variables antes de agrupar (clustering de series temporales basado en características) y aquellos que trabajan directamente en las series temporales (medidas de distancia elástica).
Clustering basado en características de series temporales: Se extraen varaibles que describen las características estructurales de cada serie temporal, y luego estas variables se introducen en algoritmos de clustering. Estas variables se obtienen aplicando operaciones estadísticas que capturan mejor las características subyacentes: tendencia, estacionalidad, periodicidad, correlación serial, asimetría, curtosis, caos, no linealidad y auto-similitud.
Medidas de distancia elástica: Este enfoque trabaja directamente en las series temporales, ajustando o "reajustando" las series en comparación con otras. La más conocida de esta familia de medidas es Dynamic Time Warping (DTW).
Para información más detallada sobre el clustering de series temporales, consulta A review and evaluation of elastic distance functions for time series clustering.
Clustering basado en características de series temporales¶
Extracción de variables¶
tsfresh es una librería de Python para extraer variables de series temporales y datos secuenciales, que incluye medidas estadísticas, coeficientes de Fourier y otras variables en el dominio del tiempo y de la frecuencia. Proporciona un enfoque sistemático para automatizar el cálculo de variables y seleccionar las más informativas.
Para empezar, se utiliza la configuración predeterminada de tsfresh, que calcula todas las variables disponibles. La forma más sencilla de acceder a esta configuración es utilizar las funciones proporcionadas por la clase ComprehensiveFCParameters. Esta función crea un diccionario que asocia el nombre de cada característica a una lista de parámetros que se utilizarán cuando se llame a la función con ese nombre.
# Variables por defecto
# ==============================================================================
default_features = ComprehensiveFCParameters()
print("Nombre de las variables extraídas por tsfresh")
print("=============================================")
print('\n'.join(list(default_features)))
Muchas de las variables se calculan utilizando diferentes valores para sus argumentos.
# Configuración por defecto para "partial_autocorrelation"
# ==============================================================================
default_features['partial_autocorrelation']
Para acceder a la vista detallada de cada característica y los valores de los parámetros incluidos en la configuración predeterminada, utilice el siguiente código:
# Configuración por defecto para todas las variables
# ==============================================================================
# from pprint import pprint
# pprint(default_features)
Una vez que se ha definido la configuración de las variables, el siguiente paso es extraerlas de las series temporales. Para ello, se utiliza la función extract_features() de tsfresh. Esta función recibe como entrada la o las series temporales y la configuración de las variables a extraer. La salida es un dataframe con las variables extraídas.
# Extracción de variables
# ==============================================================================
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("Dimensiones de ts_features:", ts_features.shape)
ts_features.head(3)
Como resultado del proceso de extracción, se han creado 783 variables para cada serie temporal (building_id en este caso). El dataframe devuelto tiene como índice la columna especificada en el argumento column_id de extract_features.
La extracción por defecto de tsfresh genera un gran número de variables. Sin embargo, solo unas pocas 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 test de hipótesis (FeatuRE Extraction based on Scalable Hypothesis tests).
⚠ Warning
El proceso de selección utilizado por tsfresh se basa en la importancia de cada característica para predecir con exactitud la variable objetivo. Para realizar este proceso, se necesita 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 casos, pueden utilizarse estrategias alternativas:- En lugar de calcular todas las variables estándar, calcular sólo las que puedan ser relevantes para la aplicación específica, basándose en los conocimientos de los expertos.
- Excluir variables basadas en criterios como baja varianza y alta correlación. Esto ayuda a refinar el conjunto de variables a considerar, centrándose en aquellas que proporcionan la información más relevante para el análisis.
- Utilizar técnicas como PCA, t-SNE o autocodificadores para reducir la dimensionalidad.
Como la selección de variables es un paso crítico que afecta a la información disponible para los siguientes pasos del análisis, es recomendable entender los parámetros que controlan el comportamiento de select_features().
X: DataFrame con las variables creadas conextract_features. Puede contener variables tanto binarias como de valor real al mismo tiempo.y: Vector objetivo que se necesita para probar qué variables son relevantes. Puede ser binario o de valor real.test_for_binary_target_binary_feature: Prueba que se utilizará cuando el objetivo es binario y la característica es binaria. Actualmente no se usa, valor predeterminado'fisher'.test_for_binary_target_real_feature: Prueba que se utilizará cuando el objetivo es binario y la característica es continua (valor real). Valor predeterminado'mann'.test_for_real_target_binary_feature: Prueba que se utilizará cuando el objetivo es continuo (valor real) y la característica es binaria. Actualmente no se usa, valor predeterminado'mann'.test_for_real_target_real_feature: Prueba que se utilizará cuando el objetivo es continuo (valor real) y la característica es continua (valor real). Actualmente no se usa, valor predeterminado'kendall'.fdr_level: El nivel de FDR que debe respetarse, este es el porcentaje teórico esperado de variables irrelevantes entre todas las variables creadas. Valor predeterminado0.05.hypotheses_independent: ¿Se puede suponer que la significancia de las variables es independiente? Normalmente, esto debe establecerse en False, ya que las variables nunca son independientes (por ejemplo, media y mediana). Valor predeterminadoFalse.n_jobs: Número de procesos a utilizar durante el cálculo del valor p. Valor predeterminado4.show_warnings: Mostrar advertencias durante el cálculo del valor p (necesario para la depuración de los calculadores). Valor predeterminadoFalse.chunksize: El tamaño de un chunk que se envía al proceso de trabajo para la paralelización. Donde un chunk se define como los datos para una característica. Si estableces el chunksize en 10, significa que una tarea es filtrar 10 variables. Si se establece enNone, dependiendo del distribuidor, se utilizan heurísticas para encontrar el tamaño óptimo de chunk. Si recibes excepciones de falta de memoria, puedes probar con el distribuidor dask y un chunksize más pequeño. Valor predeterminadoNone.ml_task: La tarea de machine learning prevista. Puede ser ‘classification’, ‘regression’ o ‘auto’. El valor predeterminado es ‘auto’, lo que significa que la tarea prevista se infiere a partir dey. Siytiene un tipo de datos booleano, entero u objeto, se asume que la tarea es clasificación, de lo contrario, regresión. Valor predeterminado'auto'.multiclass: Si el problema es de clasificación multiclase. Esto modifica la forma en que se seleccionan las variables. La clasificación multiclase requiere que las variables sean estadísticamente significativas para predecirn_significantvariables. Valor predeterminadoFalse.n_significant: El número de clases para las cuales las variables deben ser predictores estadísticamente significativos para ser consideradas ‘relevantes’. Solo especificar cuandomulticlass=True. Valor predeterminado1.
⚠ Warning
El orden del índice devuelto en el dataframe de variables no es el mismo que el orden de las columnas en el dataframe original. Por lo tanto, los datos pasados al argumentoy en select_features deben estar ordenados para garantizar la correcta asociación entre las variables y la variable objetivo.
# Seleccionar caracteristicas relecantes
# ==============================================================================
target = (
data[['building_id', 'primary_use']]
.drop_duplicates()
.set_index('building_id')
.loc[ts_features.index, :]
['primary_use']
)
assert ts_features.index.equals(target.index)
ts_features_selected = select_features(
X = ts_features,
y = target,
fdr_level = 0.001 # Un filtrado muy estricto
)
ts_features_selected.index.name = 'building_id'
print(f"Número de variables antes de la selección: {ts_features.shape[1]}")
print(f"Número de variables tras la selección: {ts_features_selected.shape[1]}")
ts_features_selected.head()
Algunas de las variables creadas pueden tener valores ausentes para algunas de las series temporales. Dado que la mayoría de los algoritmos de clustering no permiten valores ausentes, se excluyen.
# Eliminar variables con valores ausentes
# ==============================================================================
ts_features_selected = ts_features_selected.dropna(axis=1, how='any')
print(
f"Número de variables tras la selección y la eliminación de valores perdidos: "
f"{ts_features_selected.shape[1]}"
)
Una vez que la matriz de variables final esté lista, puede ser útil crear un nuevo diccionario para almacenar las variables seleccionadas y los parámetros utilizados para calcularlas. Esto se puede hacer fácilmente utilizando la función from_columns.
# Diccionario con las variables seleccionadas y su configuración
# ==============================================================================
selected_features_info = from_columns(ts_features_selected)
# pprint(selected_features_info['meter_reading'])
K-means clustering¶
El método de clustering K-means se utiliza para agrupar los edificios. Dado que se sabe que el clustering se ve negativamente afectado por la alta dimensionalidad, y dado que se han creado varias centenas de características para cada edificio, se utiliza un PCA para reducir la dimensionalidad de los datos antes de aplicar el k-means.
# Escalado de variables para que tengan media 0 y desviación estándar 1
# ==============================================================================
scaler = StandardScaler().set_output(transform="pandas")
ts_features_selected_scaled = scaler.fit_transform(ts_features_selected)
ts_features_selected_scaled.head(2)
# Gráfico de la reducción de la varianza en función del número de componentes PCA
# ==============================================================================
pca = PCA()
pca.fit(ts_features_selected_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 457 variables originales se reduce utilizando los primeros 30 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_selected_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_selected_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)
Uno de los retos inherentes del clustering es determinar el número óptimo de clústeres, ya que no hay una verdad absoluta o un número predefinido de clústeres 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 implica trazar la suma total de cuadrados dentro del clúster (WSS) frente al número de clústeres (k). El número óptimo de clústeres se identifica típicamente en el "codo" de la curva, donde la tasa de disminución de WSS comienza a disminuir notablemente. Esto sugiere que añadir más clústeres más allá de este punto proporciona poca mejora en el rendimiento del clustering.
# Número óptimo de clusters (método del codo)
# ==============================================================================
range_n_clusters = range(1, 20)
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(list(np.unique(kmeans.labels_, return_counts=True)[1]))
inertias = pd.Series(inertias, index=range_n_clusters)
dicrease_in_inertia = inertias.pct_change() * -100
fig, axs = plt.subplots(1, 2, figsize=(7, 2.5))
inertias.plot(marker='o', ax=axs[0])
axs[0].xaxis.set_major_locator(ticker.MaxNLocator(integer=True))
axs[0].set_title("Varianza Intra-cluster vs número de clusters")
axs[0].set_xlabel('Número de clusters')
axs[0].set_ylabel('Varianza Intra-cluster (Inercia)')
dicrease_in_inertia.plot(kind='bar', ax=axs[1])
axs[1].set_title("Reducción de la inercia en % vs número de clusters")
axs[1].set_xlabel('Número de clusters')
axs[1].set_ylabel('Reducción de la inercia (%)')
fig.tight_layout()
plt.show();
El gráfico muestra que después de 9 clústeres, la disminución de la inercia se ralentiza, por lo que 9 clústeres pueden ser una buena elección.
Además de analizar la evolución de la inercia (intra-varianza), es importante comprobar el tamaño de los clústeres que se están creando. La presencia de clústeres pequeños puede indicar sobreajuste o muestras anómalas que no encajan bien en ninguno de los grupos.
# Distribución del tamaño de los clusters
# ==============================================================================
for n_clusters, sizes in zip(range_n_clusters, size_of_clusters):
print(f"Tamaño del cluster (n = {n_clusters}): {sizes}")
Cuando se crean 9 clústeres, el tamaño de los clústeres no está bien equilibrado. Antes de modelizar las series temporales, los clústeres más pequeños (menos de 20 observaciones) se combinan en un único clúster denominado "other".
# Entrenamiento de un modelo de clustering con 6 clusters y asignación de cada edificio a un cluster
# ==============================================================================
kmeans = KMeans(n_clusters=9, 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)
})
# Combinar los clusters con menos de 20 edificios en un solo 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']
)
clusters['cluster_base_on_features'].value_counts()
# Añadir el cluster predicho al data frame con la información de los edificios
# ==============================================================================
data = pd.merge(
data.reset_index(), # Para evitar perder el índice
clusters,
on ='building_id',
how ='left',
validate = 'm:1'
).set_index('timestamp')
data.head(3)
Una vez que se han asignado los edificios a un clúster, es útil examinar las variables de los edificios agrupados. Por ejemplo, la distribución del uso principal de los edificios.
# Porcentaje de cada tipo de edificio en cada 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
Los resultados sugieren (la tabla debe leerse horizontalmente) que el proceso de clustering basado en características extraídas de las series temporales genera grupos que difieren de los formados por el uso principal del edificio.
Clustering basado en la distancia elástica (DTW)¶
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 forma óptima. El clustering con DTW implica agrupar datos de series temporales en función de sus distancias DTW, asegurando que las series temporales dentro del mismo clúster tengan formas y patrones similares, incluso si están desfasadas en el tiempo o tienen longitudes diferentes.
Los datos tienen que ser transformados de un formato estrecho (largo) a un formato ancho (pivot) para poder utilizar un StandardScaler de scikit-learn, que escala las series temporales de forma que todas tengan media cero y varianza unitaria. Cada columna representa un edificio.
# Tabla pivotada con series temporales para cada edificio
# ==============================================================================
data_pivot = data.pivot_table(
index = 'timestamp',
columns = 'building_id',
values = 'meter_reading',
aggfunc = 'mean'
)
data_pivot = data_pivot.fillna(value=0)
data_pivot.head(3)
# Scalado de las series
# ==============================================================================
data_pivot = pd.DataFrame(
StandardScaler().fit_transform(data_pivot),
index = data_pivot.index,
columns = data_pivot.columns
)
data_pivot.head(3)
TimeSeriesKMeans de la librería tslearn es un algoritmo de clustering especializado diseñado para datos de series temporales. Amplía el enfoque de clustering K-means tradicional para manejar específicamente las formas y patrones inherentes en los conjuntos de datos de series temporales. Agrupa series temporales similares en clústeres basados en sus variables temporales, teniendo en cuenta aspectos como tendencias, estacionalidad y varios patrones dependientes del tiempo.
⚠ Warning
La siguiente celda tiene un tiempo de ejecución de aproximadamente 15 minutos.# Entrenamieto de un modelo de clustering con 4 clusters y distancia DTW
# ==============================================================================
model = TimeSeriesKMeans(
n_clusters = 4,
metric = "dtw",
max_iter = 100,
random_state = 123456
)
# TimeSeriesKMeans asume que hay una serie temporal por fila, por lo que se
# transpone el data frame.
model.fit(data_pivot.transpose())
# Predicción de los clusters
# ==============================================================================
clusters = model.predict(data_pivot.transpose())
clusters = pd.DataFrame({
'building_id': data_pivot.columns,
'cluster_base_on_dtw': clusters.astype(str),
})
# Tamaño de cada cluster
print("Tamaño de cada cluster:")
clusters['cluster_base_on_dtw'].value_counts().sort_index()
Los 4 cuatro cluster generados por Dynamic Time Warping (DTW) están bien equilibrados en tamaño.
# Añadir el cluster predicho al data frame con la información de los edificios
# ==============================================================================
data = pd.merge(
data.reset_index(), # Para evitar perder el índice
clusters,
on = 'building_id',
how = 'left',
validate = 'm:1'
).set_index('timestamp')
data.head(3)
✎ Nota
DTAIDistance también es una excelente librería de Python para calcular distancias entre series temporales. Está basada en C++ y Cython y es muy rápida. También tiene implementadas muchas métricas de distancia y métodos de clustering. Es una buena alternativa a tslearn.# Guardar los datos para modelado para evitar repetir los pasos anteriores
# ==============================================================================
data.to_parquet('data_modelling.parquet')
Modelado y forecasting¶
After establishing the three grouping criteria (building usage, clustering based on time series features, and clustering based on dynamic time warping), both single and multi-series global models are trained. T The evaluation that follows concentrates on determining how effectively these models can forecast daily data for the final two weeks. During this assessment, three distinct metrics are used:
Después de establecer los tres criterios de agrupación (uso del edificio, clustering basado en características de series temporales y clustering basado en Dynamic Time Warping), se entrenan los modelos. La evaluación se centra en determinar la capacidadd de estos modelos para predecir los datos diarios de las dos últimas semanas. En esta evaluación, se utilizan tres métricas distintas:
El promedio del error absoluto medio (MAE) para todos los edificios.
La suma de los errores absolutos, es decir, la desviación absoluta entre el valor predicho y el consumo real para todos los edificios.
The bias for the predictions summed over all buildings.
El sesgo (bias) para las predicciones para todos los edificios.
Además de los valores rezagados (lags) de la propia serie temporal, se incluyen el día de la semana (codificado en seno-coseno), la temperatura exterior y la velocidad del viento como variables exógenas.
✎ Note
Para una explicación más detallada de la validación de modelos de series temporales, se recomienda consultar la Backtesting user guide. Para más información sobre las variables de calendario y la codificación cíclica, visita Calendar features and Cyclical features in time series.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 test.
# Lectura de los 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"Rango de fechas disponible : {data.index.min()} --- {data.index.max()} "
f"(n_días={(data.index.max() - data.index.min()).days})"
)
print(
f" Fechas entrenamiento : {data_train.index.min()} --- {data_train.index.max()} "
f"(n_días={(data_train.index.max() - data_train.index.min()).days})"
)
print(
f" Fechas validación : {data_val.index.min()} --- {data_val.index.max()} "
f"(n_días={(data_val.index.max() - data_val.index.min()).days})"
)
print(
f" Fechas test : {data_test.index.min()} --- {data_test.index.max()} "
f"(n_días={(data_test.index.max() - data_test.index.min()).days})"
)
⚠ Warning
Si bien se dispone 1214 edificios, para mantener el entrenamiento del modelo dentro de un rango de tiempo razonable, se puede utilizar un subconjunto de, por ejemplo, 600 edificios seleccionados al azar. Se anima al lector a adaptar el número de edificios si es necesario y comprobar si las conclusiones se mantienen.# Muestra 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")
# data_train = data_train.query("building_id in @buildings_selected")
# data_val = data_val.query("building_id in @buildings_selected")
# data_test = data_test.query("building_id in @buildings_selected")
# Tabla de resultados para todos los modelos
# ==============================================================================
table_results = pd.DataFrame(columns=['modelo', 'mae', 'abs_error', 'bias', 'elapsed_time'])
table_results = table_results.set_index('modelo')
table_results = table_results.astype({'mae': float, 'abs_error': float, 'bias': float, 'elapsed_time': object})
# Función para añadir variables de calendario a un DataFrame
# ==============================================================================
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 individual para cada edificio¶
Un modelo de forecasting individual se entrena y evalúa para cada edificio.
# Entrenamiento y predicción de un modelo para cada edificio
# ==============================================================================
predictions_all_buildings = {}
metrics_all_buildings = {}
errors_all_buildings = {}
# Forecaster
params_lgbm = {
'n_estimators': 500,
'learning_rate': 0.01,
'max_depth': 4,
'random_state': 8520,
'verbose': -1
}
forecaster = ForecasterRecursive(
regressor = LGBMRegressor(**params_lgbm),
lags = 31
)
# Entrenar y predecir para cada edificio
start = datetime.now()
for building in tqdm(data['building_id'].unique(), desc='Modelling buildings'):
# Seleccionar los datos del edificio
data_building = data[data['building_id'] == building]
data_building = data_building.asfreq('D').sort_index()
# Añadir variables de calendario
data_building = day_of_week_cyclical_encoding(data_building)
# Backtesting
try:
cv = TimeSeriesFold(
steps = 7,
initial_train_size = len(data_building.loc[:end_validation, :]),
refit = False
)
metric, predictions = backtesting_forecaster(
forecaster = forecaster,
y = data_building['meter_reading'],
exog = data_building[exog_features],
cv = cv,
metric = 'mean_absolute_error',
verbose = False,
show_progress = False
)
predictions_all_buildings[building] = predictions['pred']
metrics_all_buildings[building] = metric.at[0, 'mean_absolute_error']
errors_all_buildings[building] = (
predictions['pred'] - data_building.loc[predictions.index, 'meter_reading']
)
except:
print(f"Error modelling building {building}")
end = datetime.now()
predictions_all_buildings = pd.DataFrame(predictions_all_buildings)
errors_all_buildings = pd.DataFrame(errors_all_buildings)
mean_metric_all_buildings = pd.Series(metrics_all_buildings).mean()
sum_abs_errors_all_buildings = errors_all_buildings.abs().sum(axis=1).sum()
sum_bias_all_buildings = errors_all_buildings.sum().sum()
table_results.loc['One model per building', 'mae'] = mean_metric_all_buildings
table_results.loc['One model per building', 'abs_error'] = sum_abs_errors_all_buildings
table_results.loc['One model per building', 'bias'] = sum_bias_all_buildings
table_results.loc['One model per building', 'elapsed_time'] = end - start
print("")
print(f"Media de errores absolutos para todos los edificios: {mean_metric_all_buildings:.0f}")
print(f"Suma de los errores absolutos para todos los edificios (x 10,000): {sum_abs_errors_all_buildings/10000:.0f}")
print(f"Bias (x 10,000): {sum_bias_all_buildings/10000:.0f}")
# Gráfico con las predicciones y los valores reales 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_all_buildings[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 para todos los edificios¶
Se entrena un modelo global para todos los edificios y se evalúa utilizando la clase ForecasterRecursiveMultiSeries de skforecast. Este modelo permite al usuario pasar las series temporales en varios formatos, incluido un dataframe con las series temporales organizadas como columnas. Para más información, sobre como utilizar series de diferente longitud, o diferentes variables exógenas por serie, consulta la documentación de Global Forecasting Models.
# Modificar el formato de 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()
# Añadir variables de calendario
data_pivot = day_of_week_cyclical_encoding(data_pivot)
# Añadir temperatura y velocidad del viento como variables exógenas
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)
⚠ Warning
Desde la versión0.13.0, la clase ForecasterRecursiveMultiSeries utiliza por defecto encoding='ordinal_category' para codificar los identificadores de las series. Esto crea una nueva columna (_level_skforecast) de tipo pandas category. En consecuencia, los regresores deben ser capaces de manejar variables categóricas. Si los regresores no admiten variables categóricas, el usuario debe establecer la codificación en 'ordinal' o 'onehot' para garantizar la compatibilidad.
Algunos ejemplos de regresores que admiten variables categóricas y cómo habilitarlas son:
-
HistGradientBoostingRegressor:
HistGradientBoostingRegressor(categorical_features="from_dtype")
-
LightGBM:
LGBMRegressorno permite la configuración de variables categóricas durante la inicialización, sino en su métodofit. Por lo tanto, se tiene que indicarfit_kwargs = {'categorical_feature':'auto'}. Este es el comportamiento por defecto deLGBMRegressorsi no se da ninguna indicación.
-
XGBoost:
XGBRegressor(enable_categorical=True)
# Forecaster multi-series 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,
verbose = False,
show_progress = True
)
end = datetime.now()
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"Media de errores absolutos para todos los edificios: {mean_metric_all_buildings:.0f}")
print(f"Suma de los errores absolutos para todos los edificios (x 10,000): {sum_abs_errors_all_buildings/10000:.0f}")
print(f"Bias (x 10,000): {sum_bias_all_buildings/10000:.0f}")
# Añadir las 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')
axs[i].legend()
Modelo global por uso principal¶
# Forecaster multi-series models para 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"Entrenamiento y test del modelo para el uso principal: {primary_usage} "
f"(n = {data[data['primary_use'] == primary_usage]['building_id'].nunique()})"
)
# Seleccionar solo los edificios de un tipo de uso principal
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()
# Añadir variables de calendario
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,
verbose = 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=1)
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"Media de errores absolutos para todos los edificios: {mean_metric_all_buildings:.0f}")
print(f"Suma de los errores absolutos para todos los edificios (x 10,000): {sum_abs_errors_all_buildings/10000:.0f}")
print(f"Bias (x 10,000): {sum_bias_all_buildings/10000:.0f}")
# Añadir las predicciones al gráfico 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 basado en características de las series¶
Se entrena y evalua un modelo global para cada cluster basado en las variables de las series temporales.
# Forecaster multi-series para edificios agrupados por las 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"Entrenamiento y test del modelo para el cluster: {cluster} "
f"(n = {data[data['cluster_base_on_features'] == cluster]['building_id'].nunique()})"
)
# Seleccionar solo los edificios del cluster
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()
# Añadir variables de calendario
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,
verbose = 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=1)
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"Media de errores absolutos para todos los edificios: {mean_metric_all_buildings:.0f}")
print(f"Suma de los errores absolutos para todos los edificios (x 10,000): {sum_abs_errors_all_buildings/10000:.0f}")
print(f"Bias (x 10,000): {sum_bias_all_buildings/10000:.0f}")
# Añadir las predicciones al gráfico 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 basado DTW¶
Se entrenan y evalúan un modelo global para cada clúster basado en la distancia elástica (DTW).
# Forecaster multi-series para 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"Entrenamiento y test del modelo para el cluster: {cluster} "
f"(n = {data[data['cluster_base_on_dtw'] == cluster]['building_id'].nunique()})"
)
# Seleccionar solo los edificios del cluster
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()
# Añadir variables de calendario
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,
verbose = 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=1)
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"Media de errores absolutos para todos los edificios: {mean_metric_all_buildings:.0f}")
print(f"Suma de los errores absolutos para todos los edificios (x 10,000): {sum_abs_errors_all_buildings/10000:.0f}")
print(f"Bias (x 10,000): {sum_bias_all_buildings/10000:.0f}")
# Añadir las predicciones al gráfico 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()
Results¶
# 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)
# Gráfico con las predicciones y los valores reales para 2 edificios seleccionados al azar
# ==============================================================================
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¶
La decisión entre utilizar modelos individuales globales puede influir mucho en el resultado. Este análisis compara estas dos estrategias, destacando el balance entre la capacidad predictiva y la eficiencia computacional.
Capacidad de predicción: Los modelos globales superan a los individuales en términos de error medio absoluto, lo que sugiere que captan mejor los patrones comunes entre las series, dando lugar a mejores predicciones.
Eficiencia computacional: Los modelos globales requieren menos tiempo que la ejecución de múltiples modelos individuales. Esto pone de relieve la ventaja del modelo global en aplicaciones sensibles al tiempo, en las que se valora un entrenamiento y predicción rápidos.
Mejor modelo: para este caso de uso, entre todos lo modelos probados, el "Global model per cluster (features)" muestra el menor error medio absoluto y el menor error absoluto.
Próximos pasos
Este análisis ha proporcionado ideas interesantes sobre la eficacia de los modelos globales frente a los individuales en el forecasting de series temporales. Los posibles próximos pasos podrían incluir:
Revisar los edificios con un error alto. Identifica si hay un grupo para el que el modelo no funciona bien.
Add more exogenous features: see skforecast's user guide Calendar Features for calendar and sunlight features that typically affect energy consumption.
Añadir más variables exógenas: consultar skforecast's user guide Calendar Features para obtener variables de calendario y luz solar que suelen afectar al consumo de energía.
Optimización de lags e hiperparámetros: Utilizar Grid, Random o Bayesian search 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)
Bibliography¶
Hyndman, R.J., & Athanasopoulos, G. (2021) Forecasting: principles and practice, 3rd edition, OTexts: Melbourne, Australia.
Time Series Analysis and Forecasting with ADAM Ivan Svetunkov
Joseph, M. (2022). Modern time series forecasting with Python: Explore industry-ready time series forecasting using modern machine learning and Deep Learning. Packt Publishing.
Time Series Analysis and Forecasting with ADAM Ivan Svetunkov
Holder, C., Middlehurst, M. & Bagnall, A. A review and evaluation of elastic distance functions for time series clustering. Knowl Inf Syst (2023)
Forecasting: theory and practice
Instrucciones para citar¶
¿Cómo citar este documento?
Si utilizas este documento o alguna parte de él, te agradecemos que lo cites. ¡Muchas gracias!
Modelos de forecasting globales: Análisis comparativo de modelos de una y múltiples series 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/py53-modelos-forecasting-globales.html
¿Cómo citar skforecast?
Si utilizas skforecast en tu investigación o publicación, te lo agradeceríamos mucho que lo cites. ¡Muchas gracias!
Zenodo:
Amat Rodrigo, Joaquin, & Escobar Ortiz, Javier. (2024). skforecast (v0.14.0). Zenodo. https://doi.org/10.5281/zenodo.8382788
APA:
Amat Rodrigo, J., & Escobar Ortiz, J. (2024). skforecast (Version 0.14.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.14.0}, month = {11}, year = {2024}, license = {BSD-3-Clause}, url = {https://skforecast.org/}, doi = {10.5281/zenodo.8382788} }
¿Te ha gustado el artículo? Tu ayuda es importante
Mantener un sitio web tiene unos costes elevados, 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.
-
NoComercial: No puedes utilizar el material para fines comerciales.
-
CompartirIgual: Si remezclas, transformas o creas a partir del material, debes distribuir tus contribuciones bajo la misma licencia que el original.