Predicción (forecasting) de la demanda energética con machine learning

Si te gusta  Skforecast , ayúdanos dándonos una estrella en  GitHub! ⭐️

Predicción (forecasting) de la demanda energética con machine learning

Joaquín Amat Rodrigo, Javier Escobar Ortiz
Febrero, 2021 (última actualización diciembre 2023)

Introducción

La predicción de la demanda energética desempeña un papel fundamental en la gestión y planificación de los recursos necesarios para la generación, distribución y utilización de la energía. Predecir la demanda de energía es una tarea compleja en la que influyen factores como los patrones meteorológicos, las condiciones económicas y el comportamiento de la sociedad. Este documento muestra cómo utilizar modelos de machine learning para predecir la demanda de energía.

Series temporales y forecasting

Una serie temporal (time series) es una sucesión de datos ordenados cronológicamente, espaciados a intervalos iguales o desiguales. El proceso de forecasting consiste en predecir el valor futuro de una serie temporal, bien modelando la serie únicamente en función de su comportamiento pasado (autorregresivo) o empleando otras variables externas.


Cuando se trabaja con series temporales, raramente se quiere predecir solo el siguiente elemento de la serie ($t_{+1}$), sino todo un horizonte futuro ($t_{+1}$), ..., ($t_{+n}$)) o un punto alejado en el tiempo ($t_{+n}$). Existen varias estrategias que permiten generar este tipo de predicciones, la librería skforecast recoge las siguientes para series temporales univariantes:

  • Forecasting multi-step recursivo: este método consiste en utilizar propias predicciones del modelo como valores de entrada para predecir el siguiente valor. Por ejemplo, para predecir los 5 valores siguientes de una serie temporal, se entrena un modelo para predecir el siguiente valor ($t_{+1}$), y se utiliza este valor para predecir el siguiente ($t_{+2}$), y así sucesivamente. Todo este proceso se automatiza con las clases ForecasterAutoreg y ForecasterAutoregCustom.
  • Forecasting multi-step directo: este método consiste en entrenar un modelo diferente para cada valor futuro (step) del horizonte de predicción. Por ejemplo, para predecir los 5 siguientes valores de una serie temporal, se entrenan 5 modelos diferentes, uno para cada step. De este modo, las predicciones son independientes entre sí. Todo este proceso se automatiza en la clase ForecasterAutoregDirect.

✎ Nota

Two other great examples of how to use gradient boosting for time series forecasting are: Otros dos ejemplos de cómo utilizar machine learning (*gradient boosting*) para forecasting de series temporales son:

Librerías

Las librerías utilizadas en este documento son:

In [1]:
# Tratamiento de datos
# ==============================================================================
import numpy as np
import pandas as pd
from astral.sun import sun
from astral import LocationInfo
from skforecast.datasets import fetch_dataset

# Gráficos
# ==============================================================================
import matplotlib.pyplot as plt
from statsmodels.graphics.tsaplots import plot_acf
from statsmodels.graphics.tsaplots import plot_pacf
import plotly.graph_objects as go
import plotly.io as pio
import plotly.offline as poff
pio.templates.default = "seaborn"
poff.init_notebook_mode(connected=True)
plt.style.use('seaborn-v0_8-darkgrid')

# Modelado y Forecasting
# ==============================================================================
from lightgbm import LGBMRegressor
from sklearn.preprocessing import PolynomialFeatures
from sklearn.feature_selection import RFECV
from skforecast.ForecasterBaseline import ForecasterEquivalentDate
from skforecast.ForecasterAutoreg import ForecasterAutoreg
from skforecast.ForecasterAutoregDirect import ForecasterAutoregDirect
from skforecast.model_selection import bayesian_search_forecaster
from skforecast.model_selection import backtesting_forecaster
import shap
shap.initjs()

# Configuración warnings
# ==============================================================================
import warnings
warnings.filterwarnings('once')

Datos

Se dispone de una serie temporal de la demanda de electricidad (MW) para el estado de Victoria (Australia) desde 2011-12-31 hasta 2014-12-31. Los datos empleados en este documento se han obtenido del paquete de R tsibbledata. El set de datos contiene 5 columnas y 52608 registros completos. La información de cada columna es:

  • Time: fecha y hora del registro.
  • Date: fecha del registro
  • Demand: demanda de electricidad (MW).
  • Temperature: temperatura en Melbourne, capital de Victoria.
  • Holiday: indicador si el día es festivo (vacaciones).
In [2]:
# Descarga de datos
# ==============================================================================
datos = fetch_dataset(name='vic_electricity', raw=True)
datos.info()
vic_electricity
---------------
Half-hourly electricity demand for Victoria, Australia
O'Hara-Wild M, Hyndman R, Wang E, Godahewa R (2022).tsibbledata: Diverse
Datasets for 'tsibble'. https://tsibbledata.tidyverts.org/,
https://github.com/tidyverts/tsibbledata/.
https://tsibbledata.tidyverts.org/reference/vic_elec.html
Shape of the dataset: (52608, 5)
<class 'pandas.core.frame.DataFrame'>
RangeIndex: 52608 entries, 0 to 52607
Data columns (total 5 columns):
 #   Column       Non-Null Count  Dtype  
---  ------       --------------  -----  
 0   Time         52608 non-null  object 
 1   Demand       52608 non-null  float64
 2   Temperature  52608 non-null  float64
 3   Date         52608 non-null  object 
 4   Holiday      52608 non-null  bool   
dtypes: bool(1), float64(2), object(2)
memory usage: 1.7+ MB

La columna Time está almacenada como string. Para convertirla en datetime, se emplea la función pd.to_datetime(). Una vez en formato datetime, y para hacer uso de las funcionalidades de pandas, se establece como índice. Además, dado que los datos se han registrado cada 30 minutos, se indica la frecuencia '30min'.

In [3]:
# Conversión del formato fecha
# ==============================================================================
datos['Time'] = pd.to_datetime(datos['Time'], format='%Y-%m-%dT%H:%M:%SZ')
datos = datos.set_index('Time')
datos = datos.asfreq('30min')
datos = datos.sort_index()
datos.head(2)
Out[3]:
Demand Temperature Date Holiday
Time
2011-12-31 13:00:00 4382.825174 21.40 2012-01-01 True
2011-12-31 13:30:00 4263.365526 21.05 2012-01-01 True

Uno de los primeros análisis que hay que realizar al trabajar con series temporales es verificar si la serie está completa.

In [4]:
# Verificar que un índice temporal está completo
# ==============================================================================
(datos.index == pd.date_range(start=datos.index.min(),
                              end=datos.index.max(),
                              freq=datos.index.freq)).all()

print(f"Filas con valores ausentes: {datos.isnull().any(axis=1).mean()}")
Filas con valores ausentes: 0.0
In [5]:
# Completar huecos en un índice temporal
# ==============================================================================
# datos.asfreq(freq='30min', fill_value=np.nan)

Aunque los datos se encuentran en intervalos de 30 minutos, el objetivo es crear un modelo capaz de predecir la demanda eléctrica a nivel horario, por lo que se tienen que agregar los datos. Este tipo de transformación es muy sencillas si se combina el índice DatetimeIndex de pandas y su método resample().

Es muy importante utilizar correctamente los argumentos closed='left' y label='right' para no introducir en el entrenamiento información a futuro (leakage)). Supóngase que se dispone de valores para las 10:10, 10:30, 10:45, 11:00, 11:12 y 11:30. Si se quiere obtener el promedio horario, el valor asignado a las 11:00 debe calcularse utilizando los valores de las 10:10, 10:30 y 10:45; y el de las 12:00, con el valor de las 11:00, 11:12 y 11:30.

Diagrama de agregado de datos temporales sin incluir información de futuro.

Para el valor promedio de las 11:00 no se incluye el valor puntual de las 11:00 por que, en la realidad, en ese momento exacto no se dispone todavía del valor.

In [6]:
# Agregado en intervalos de 1H
# ==============================================================================
# Se elimina la columna Date para que no genere error al agregar. La columna Holiday
# no genera error ya que es booleana y se trata como 0-1.
datos = datos.drop(columns='Date')
datos = datos.resample(rule='H', closed='left', label ='right').mean()
datos
Out[6]:
Demand Temperature Holiday
Time
2011-12-31 14:00:00 4323.095350 21.225 1.0
2011-12-31 15:00:00 3963.264688 20.625 1.0
2011-12-31 16:00:00 3950.913495 20.325 1.0
2011-12-31 17:00:00 3627.860675 19.850 1.0
2011-12-31 18:00:00 3396.251676 19.025 1.0
... ... ... ...
2014-12-31 09:00:00 4069.625550 21.600 0.0
2014-12-31 10:00:00 3909.230704 20.300 0.0
2014-12-31 11:00:00 3900.600901 19.650 0.0
2014-12-31 12:00:00 3758.236494 18.100 0.0
2014-12-31 13:00:00 3785.650720 17.200 0.0

26304 rows × 3 columns

El set de datos empieza el 2011-12-31 14:00:00 y termina el 2014-12-31 13:00:00. Se descartan los primeros 10 y los últimos 13 registros para que empiece el 2012-01-01 00:00:00 y termine el 2014-12-30 23:00:00. Además, para poder optimizar los hiperparámetros del modelo y evaluar su capacidad predictiva, se dividen los datos en 3 conjuntos, uno de entrenamiento, uno de validación y otro de test.

In [7]:
# Separación datos train-val-test
# ==============================================================================
datos = datos.loc['2012-01-01 00:00:00': '2014-12-30 23:00:00'].copy()
fin_train = '2013-12-31 23:59:00'
fin_validacion = '2014-11-30 23:59:00'
datos_train = datos.loc[: fin_train, :].copy()
datos_val   = datos.loc[fin_train:fin_validacion, :].copy()
datos_test  = datos.loc[fin_validacion:, :].copy()

print(f"Fechas train      : {datos_train.index.min()} --- {datos_train.index.max()}  (n={len(datos_train)})")
print(f"Fechas validacion : {datos_val.index.min()} --- {datos_val.index.max()}  (n={len(datos_val)})")
print(f"Fechas test       : {datos_test.index.min()} --- {datos_test.index.max()}  (n={len(datos_test)})")
Fechas train      : 2012-01-01 00:00:00 --- 2013-12-31 23:00:00  (n=17544)
Fechas validacion : 2014-01-01 00:00:00 --- 2014-11-30 23:00:00  (n=8016)
Fechas test       : 2014-12-01 00:00:00 --- 2014-12-30 23:00:00  (n=720)

Exploración gráfica

La exploración gráfica de series temporales es una forma eficaz de identificar tendencias, patrones y estacionalidad. Esto, a su vez, ayuda a orientar la selección del modelo de forecasting más adecuado.

Gráfico de la serie temporal

Serie temporal completa

In [8]:
# Gráfico interactivo de la serie temporal
# ==============================================================================
fig = go.Figure()
fig.add_trace(go.Scatter(x=datos_train.index, y=datos_train['Demand'], mode='lines', name='Train'))
fig.add_trace(go.Scatter(x=datos_val.index, y=datos_val['Demand'], mode='lines', name='Validation'))
fig.add_trace(go.Scatter(x=datos_test.index, y=datos_test['Demand'], mode='lines', name='Test'))
fig.update_layout(
    title  = 'Demanda eléctrica horaria',
    xaxis_title="Fecha",
    yaxis_title="Demanda (MWh)",
    legend_title="Partición:",
    width=850,
    height=400,
    margin=dict(l=20, r=20, t=35, b=20),
    legend=dict(
        orientation="h",
        yanchor="top",
        y=1,
        xanchor="left",
        x=0.001
    )
)
#fig.update_xaxes(rangeslider_visible=True)
fig.show()

El gráfico anterior muestra que la demanda eléctrica tiene estacionalidad anual. Se observa un incremento centrado en el mes de Julio y picos de demanda muy acentuados entre enero y marzo.

Sección de la serie temporal

Debido a la varianza de la serie temporal, no es posible apreciar con un solo gráfico el posible patrón intradiario.

In [9]:
# Gráfico serie temporal con zoom
# ==============================================================================
zoom = ('2013-05-01 14:00:00','2013-06-01 14:00:00')
fig = plt.figure(figsize=(8, 4))
grid = plt.GridSpec(nrows=8, ncols=1, hspace=0.6, wspace=0)
main_ax = fig.add_subplot(grid[1:3, :])
zoom_ax = fig.add_subplot(grid[5:, :])
datos.Demand.plot(ax=main_ax, c='black', alpha=0.5, linewidth=0.5)
min_y = min(datos.Demand)
max_y = max(datos.Demand)
main_ax.fill_between(zoom, min_y, max_y, facecolor='blue', alpha=0.5, zorder=0)
main_ax.set_xlabel('')
datos.loc[zoom[0]: zoom[1]].Demand.plot(ax=zoom_ax, color='blue', linewidth=1)
main_ax.set_title(f'Demanda electricidad: {datos.index.min()}, {datos.index.max()}', fontsize=10)
zoom_ax.set_title(f'Demanda electricidad: {zoom}', fontsize=10)
zoom_ax.set_xlabel('')
plt.subplots_adjust(hspace=1)

Al aplicar zoom sobre la serie temporal, se hace patente una clara estacionalidad semanal, con consumos más elevados durante la semana laboral (lunes a viernes) y menor en los fines de semana. Se observa también que existe una clara correlación entre el consumo de un día con el del día anterior.

Gráfico de estacionalidad

Estacionalidad anual, semanal y diaria

In [10]:
# Distribución de la demanda por mes
# ==============================================================================
fig, ax = plt.subplots(figsize=(5.5, 2.5))
datos['mes'] = datos.index.month
datos.boxplot(column='Demand', by='mes', ax=ax,)
datos.groupby('mes')['Demand'].median().plot(style='o-', linewidth=0.8, ax=ax)
ax.set_ylabel('Demand')
ax.set_title('Distribución demanda por mes')
fig.suptitle('');

Se observa que hay una estacionalidad anual, con valores de demanda (mediana) superiores en los meses de Junio, Julio y Agosto, y con elevados picos de demanda en los meses de Noviembre, Diciembre, Enero, Febrero y Marzo.

In [11]:
# Distribución de la demanda por día de la semana
# ==============================================================================
fig, ax = plt.subplots(figsize=(5.5, 2.5))
datos['dia_semana'] = datos.index.day_of_week + 1
datos.boxplot(column='Demand', by='dia_semana', ax=ax)
datos.groupby('dia_semana')['Demand'].median().plot(style='o-', linewidth=0.8, ax=ax)
ax.set_ylabel('Demand')
ax.set_title('Distribución demanda por día de la semana')
fig.suptitle('');

Se aprecia una estacionalidad semanal, con valores de demanda inferiores durante el fin de semana.

In [12]:
# Distribución de la demanda por hora del día
# ==============================================================================
fig, ax = plt.subplots(figsize=(5.5, 2.5))
datos['hora_dia'] = datos.index.hour + 1
datos.boxplot(column='Demand', by='hora_dia', ax=ax)
datos.groupby('hora_dia')['Demand'].median().plot(style='o-', linewidth=0.8, ax=ax)
ax.set_ylabel('Demand')
ax.set_title('Distribución demanda por hora del día')
fig.suptitle('');

También existe una estacionalidad diaria, la demanda se reduce entre las 16 y las 21 horas.

Gráficos de autocorrelación

Los gráficos de autocorrelación muestran la correlación entre una serie temporal y sus valores pasados. Son una herramienta útil para identificar el orden de un modelo autorregresivo, es decir, los valores pasados (lags) que se deben incluir en el modelo.

La función de autocorrelación (ACF) mide la correlación entre una serie temporal y sus valores pasados. La función de autocorrelación parcial (PACF) mide la correlación entre una serie temporal y sus valores pasados, pero solo después de eliminar las variaciones explicadas por los valores pasados intermedios.

In [13]:
# Gráfico autocorrelación
# ==============================================================================
fig, ax = plt.subplots(figsize=(5, 2))
plot_acf(datos.Demand, ax=ax, lags=60)
plt.show()
In [14]:
# Gráfico autocorrelación parcial
# ==============================================================================
fig, ax = plt.subplots(figsize=(5, 2))
plot_pacf(datos.Demand, ax=ax, lags=60)
plt.show()

Los gráficos de autocorrelación y autocorrelación parcial muestran una clara asociación entre la demanda de una hora y las horas anteriores, así como entre la demanda de una hora y la demanda de esa misma hora los días anteriores. Este tipo de correlación, es un indicativo de que los modelos autorregresivos pueden funcionar bien.

Baseline

Al enfrentarse a un problema de forecasting, es recomendable disponer de un modelo de referencia (baseline). Suele tratarse de un modelo muy sencillo que puede utilizarse como referencia para evaluar si merece la pena aplicar modelos más complejos.

Skforecast permite crear fácilmente un modelo de referencia con su clase ForecasterEquivalentDate. Este modelo, también conocido como Seasonal Naive Forecasting, simplemente devuelve el valor observado en el mismo periodo de la temporada anterior (por ejemplo, el mismo día laboral de la semana anterior, la misma hora del día anterior, etc.).

A partir del análisis exploratorio realizado, el modelo de referencia será el que prediga cada hora utilizando el valor de la misma hora del día anterior.

✎ Nota

En las siguientes celdas de código, se entrena un modelo *baseline* y se evalúa su capacidad predictiva mediante un proceso de *backtesting*. Si este concepto es nuevo para ti, no te preocupes, se explicará en detalle a lo largo del documento. Por ahora, basta con saber que el proceso de *backtesting* consiste en entrenar el modelo con una cierta cantidad de datos y evaluar su capacidad predictiva con los datos que el modelo no ha visto. La métrica del error se utilizará como referencia para comparar la capacidad predictiva de los modelos más complejos que se implementarán a lo largo del documento.
In [15]:
#Crear un baseline: valor de la misma hora del día anterior
# ==============================================================================
forecaster = ForecasterEquivalentDate(
                 offset    = pd.DateOffset(days=1),
                 n_offsets = 1
             )

# Entremaiento del forecaster
# ==============================================================================
forecaster.fit(y=datos.loc[:fin_validacion, 'Demand'])
forecaster
Out[15]:
======================== 
ForecasterEquivalentDate 
======================== 
Offset: <DateOffset: days=1> 
Number of offsets: 1 
Aggregation function: mean 
Window size: 24 
Training range: [Timestamp('2012-01-01 00:00:00'), Timestamp('2014-11-30 23:00:00')] 
Training index type: DatetimeIndex 
Training index frequency: H 
Creation date: 2023-12-10 19:09:54 
Last fit date: 2023-12-10 19:09:54 
Skforecast version: 0.11.0 
Python version: 3.11.4 
Forecaster id: None 
In [16]:
# Backtesting
# ==============================================================================
metric, predictions = backtesting_forecaster(
                          forecaster         = forecaster,
                          y                  = datos['Demand'],
                          steps              = 24,
                          metric             = 'mean_absolute_error',
                          initial_train_size = len(datos.loc[:fin_validacion]),
                          refit              = False,
                          n_jobs             = 'auto',
                          verbose            = False,
                          show_progress      = True
                      )

print(f"Backtest error (MAE): {metric}")
Backtest error (MAE): 308.3752715958334

Modelo autoregresivo recursivo

Se entrena un modelo autorregresivo recursivo (ForecasterAutoreg) con gradient boosting LGBMRegressor como regresor. Se utiliza una ventana temporal de 24 horas (24 lags) para predecir la demanda de la hora siguiente. Esto significa que los valores de demanda de las 24 horas anteriores se utilizan como predictores. Los hiperparámetros del regresor se dejan en sus valores por defecto.

In [17]:
# Crear el forecaster
# ==============================================================================
forecaster = ForecasterAutoreg(
                 regressor = LGBMRegressor(random_state=15926, verbose=-1),
                 lags      = 24
             )
             
# Entrena el forecaster
# ==============================================================================
forecaster.fit(y=datos.loc[:fin_validacion, 'Demand'])
forecaster
Out[17]:
================= 
ForecasterAutoreg 
================= 
Regressor: LGBMRegressor(random_state=15926, verbose=-1) 
Lags: [ 1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24] 
Transformer for y: None 
Transformer for exog: None 
Window size: 24 
Weight function included: False 
Differentiation order: None 
Exogenous included: False 
Type of exogenous variable: None 
Exogenous variables names: None 
Training range: [Timestamp('2012-01-01 00:00:00'), Timestamp('2014-11-30 23:00:00')] 
Training index type: DatetimeIndex 
Training index frequency: H 
Regressor parameters: {'boosting_type': 'gbdt', 'class_weight': None, 'colsample_bytree': 1.0, 'importance_type': 'split', 'learning_rate': 0.1, 'max_depth': -1, 'min_child_samples': 20, 'min_child_weight': 0.001, 'min_split_gain': 0.0, 'n_estimators': 100, 'n_jobs': None, 'num_leaves': 31, 'objective': None, 'random_state': 15926, 'reg_alpha': 0.0, 'reg_lambda': 0.0, 'subsample': 1.0, 'subsample_for_bin': 200000, 'subsample_freq': 0, 'verbose': -1} 
fit_kwargs: {} 
Creation date: 2023-12-10 19:09:55 
Last fit date: 2023-12-10 19:09:56 
Skforecast version: 0.11.0 
Python version: 3.11.4 
Forecaster id: None 

Backtesting

Para obtener una estimación robusta de la capacidad predictiva del modelo, se realiza un proceso de backtesting. El proceso de backtesting consiste en generar una predicción para cada observación del conjunto de test, siguiendo el mismo procedimiento que se seguiría si el modelo estuviese en producción, y finalmente comparar el valor predicho con el valor real.

El proceso de backtesting se aplica mediante la función backtesting_forecaster(). Para este caso de uso, la simulación se lleva a cabo de la siguiente manera: el modelo se entrena con datos de 2012-01-01 00:00 a 2014-11-30 23:59, y luego predice las siguientes 24 horas cada día a las 23:59. La métrica de error utilizada es el error absoluto medio (MAE).

  Note

A partir de la versión skforecast 0.9.0, todas las funciones de backtesting y búsqueda de hiperparámetros del módulo model_selection incluyen el parámetro n_jobs que permite su paralelización y, por tanto, mejora su rendimiento. Los beneficios de la paralelización dependen de varios factores, incluyendo el regresor seleccionado, el número de entrenamientos a realizar y el volumen de datos involucrados. Cuando el parámetro n_jobs se establece en 'auto', el nivel de paralelización se selecciona automáticamente basándose en reglas heurísticas que pretenden elegir la mejor configuración posible para cada escenario.
In [18]:
# Backtesting
# ==============================================================================
metrica, predicciones = backtesting_forecaster(
                            forecaster         = forecaster,
                            y                  = datos['Demand'],
                            steps              = 24,
                            metric             = 'mean_absolute_error',
                            initial_train_size = len(datos.loc[:fin_validacion]),
                            refit              = False,
                            n_jobs             = 'auto',
                            verbose            = True,
                            show_progress      = True
                        )
Information of backtesting process
----------------------------------
Number of observations used for initial training: 25560
Number of observations used for backtesting: 720
    Number of folds: 30
    Number of steps per fold: 24
    Number of steps to exclude from the end of each train set before test (gap): 0

Fold: 0
    Training:   2012-01-01 00:00:00 -- 2014-11-30 23:00:00  (n=25560)
    Validation: 2014-12-01 00:00:00 -- 2014-12-01 23:00:00  (n=24)
Fold: 1
    Training:   2012-01-01 00:00:00 -- 2014-11-30 23:00:00  (n=25560)
    Validation: 2014-12-02 00:00:00 -- 2014-12-02 23:00:00  (n=24)
Fold: 2
    Training:   2012-01-01 00:00:00 -- 2014-11-30 23:00:00  (n=25560)
    Validation: 2014-12-03 00:00:00 -- 2014-12-03 23:00:00  (n=24)
Fold: 3
    Training:   2012-01-01 00:00:00 -- 2014-11-30 23:00:00  (n=25560)
    Validation: 2014-12-04 00:00:00 -- 2014-12-04 23:00:00  (n=24)
Fold: 4
    Training:   2012-01-01 00:00:00 -- 2014-11-30 23:00:00  (n=25560)
    Validation: 2014-12-05 00:00:00 -- 2014-12-05 23:00:00  (n=24)
Fold: 5
    Training:   2012-01-01 00:00:00 -- 2014-11-30 23:00:00  (n=25560)
    Validation: 2014-12-06 00:00:00 -- 2014-12-06 23:00:00  (n=24)
Fold: 6
    Training:   2012-01-01 00:00:00 -- 2014-11-30 23:00:00  (n=25560)
    Validation: 2014-12-07 00:00:00 -- 2014-12-07 23:00:00  (n=24)
Fold: 7
    Training:   2012-01-01 00:00:00 -- 2014-11-30 23:00:00  (n=25560)
    Validation: 2014-12-08 00:00:00 -- 2014-12-08 23:00:00  (n=24)
Fold: 8
    Training:   2012-01-01 00:00:00 -- 2014-11-30 23:00:00  (n=25560)
    Validation: 2014-12-09 00:00:00 -- 2014-12-09 23:00:00  (n=24)
Fold: 9
    Training:   2012-01-01 00:00:00 -- 2014-11-30 23:00:00  (n=25560)
    Validation: 2014-12-10 00:00:00 -- 2014-12-10 23:00:00  (n=24)
Fold: 10
    Training:   2012-01-01 00:00:00 -- 2014-11-30 23:00:00  (n=25560)
    Validation: 2014-12-11 00:00:00 -- 2014-12-11 23:00:00  (n=24)
Fold: 11
    Training:   2012-01-01 00:00:00 -- 2014-11-30 23:00:00  (n=25560)
    Validation: 2014-12-12 00:00:00 -- 2014-12-12 23:00:00  (n=24)
Fold: 12
    Training:   2012-01-01 00:00:00 -- 2014-11-30 23:00:00  (n=25560)
    Validation: 2014-12-13 00:00:00 -- 2014-12-13 23:00:00  (n=24)
Fold: 13
    Training:   2012-01-01 00:00:00 -- 2014-11-30 23:00:00  (n=25560)
    Validation: 2014-12-14 00:00:00 -- 2014-12-14 23:00:00  (n=24)
Fold: 14
    Training:   2012-01-01 00:00:00 -- 2014-11-30 23:00:00  (n=25560)
    Validation: 2014-12-15 00:00:00 -- 2014-12-15 23:00:00  (n=24)
Fold: 15
    Training:   2012-01-01 00:00:00 -- 2014-11-30 23:00:00  (n=25560)
    Validation: 2014-12-16 00:00:00 -- 2014-12-16 23:00:00  (n=24)
Fold: 16
    Training:   2012-01-01 00:00:00 -- 2014-11-30 23:00:00  (n=25560)
    Validation: 2014-12-17 00:00:00 -- 2014-12-17 23:00:00  (n=24)
Fold: 17
    Training:   2012-01-01 00:00:00 -- 2014-11-30 23:00:00  (n=25560)
    Validation: 2014-12-18 00:00:00 -- 2014-12-18 23:00:00  (n=24)
Fold: 18
    Training:   2012-01-01 00:00:00 -- 2014-11-30 23:00:00  (n=25560)
    Validation: 2014-12-19 00:00:00 -- 2014-12-19 23:00:00  (n=24)
Fold: 19
    Training:   2012-01-01 00:00:00 -- 2014-11-30 23:00:00  (n=25560)
    Validation: 2014-12-20 00:00:00 -- 2014-12-20 23:00:00  (n=24)
Fold: 20
    Training:   2012-01-01 00:00:00 -- 2014-11-30 23:00:00  (n=25560)
    Validation: 2014-12-21 00:00:00 -- 2014-12-21 23:00:00  (n=24)
Fold: 21
    Training:   2012-01-01 00:00:00 -- 2014-11-30 23:00:00  (n=25560)
    Validation: 2014-12-22 00:00:00 -- 2014-12-22 23:00:00  (n=24)
Fold: 22
    Training:   2012-01-01 00:00:00 -- 2014-11-30 23:00:00  (n=25560)
    Validation: 2014-12-23 00:00:00 -- 2014-12-23 23:00:00  (n=24)
Fold: 23
    Training:   2012-01-01 00:00:00 -- 2014-11-30 23:00:00  (n=25560)
    Validation: 2014-12-24 00:00:00 -- 2014-12-24 23:00:00  (n=24)
Fold: 24
    Training:   2012-01-01 00:00:00 -- 2014-11-30 23:00:00  (n=25560)
    Validation: 2014-12-25 00:00:00 -- 2014-12-25 23:00:00  (n=24)
Fold: 25
    Training:   2012-01-01 00:00:00 -- 2014-11-30 23:00:00  (n=25560)
    Validation: 2014-12-26 00:00:00 -- 2014-12-26 23:00:00  (n=24)
Fold: 26
    Training:   2012-01-01 00:00:00 -- 2014-11-30 23:00:00  (n=25560)
    Validation: 2014-12-27 00:00:00 -- 2014-12-27 23:00:00  (n=24)
Fold: 27
    Training:   2012-01-01 00:00:00 -- 2014-11-30 23:00:00  (n=25560)
    Validation: 2014-12-28 00:00:00 -- 2014-12-28 23:00:00  (n=24)
Fold: 28
    Training:   2012-01-01 00:00:00 -- 2014-11-30 23:00:00  (n=25560)
    Validation: 2014-12-29 00:00:00 -- 2014-12-29 23:00:00  (n=24)
Fold: 29
    Training:   2012-01-01 00:00:00 -- 2014-11-30 23:00:00  (n=25560)
    Validation: 2014-12-30 00:00:00 -- 2014-12-30 23:00:00  (n=24)

In [19]:
# Gráfico prediccion vs valores reales
# ==============================================================================
fig = go.Figure()
trace1 = go.Scatter(x=datos_test.index, y=datos_test['Demand'], name="test", mode="lines")
trace2 = go.Scatter(x=predicciones.index, y=predicciones['pred'], name="prediction", mode="lines")
fig.add_trace(trace1)
fig.add_trace(trace2)
fig.update_layout(
    title="Real value vs predicted in test data",
    xaxis_title="Date time",
    yaxis_title="Demand",
    width=800,
    height=400,
    margin=dict(l=20, r=20, t=35, b=20),
    legend=dict(
        orientation="h",
        yanchor="top",
        y=1.1,
        xanchor="left",
        x=0.001
    )
)
fig.show()
In [20]:
# Error backtest
# ==============================================================================
print(f'Error de backtest (MAE): {metrica}')
Error de backtest (MAE): 236.38849515592196

El modelo autorregresivo alcanza un MAE inferior al del modelo baseline. Esto significa que el modelo autorregresivo es capaz de predecir la demanda de electricidad del día siguiente con mayor precisión que el modelo utilizado como referencia.

Optimización de hiperparámetros (tuning)

El ForecasterAutoreg entrenado utiliza los primeros 24 lags y un modelo LGMBRegressor con los hiperparámetros por defecto. Sin embargo, no hay ninguna razón por la que estos valores sean los más adecuados. Para encontrar los mejores hiperparámetros, se realiza una Búsqueda Bayesiana utilizando la función bayesian_search_forecaster(). La búsqueda se lleva a cabo utilizando el mismo proceso de backtesting que antes, pero cada vez, el modelo se entrena con diferentes combinaciones de hiperparámetros y lags. Es importante señalar que la búsqueda de hiperparámetros debe realizarse utilizando el conjunto de validación, nunca con los datos de test.

In [21]:
# Búsqueda bayesiana de hiperparámetros
# ==============================================================================
forecaster = ForecasterAutoreg(
                 regressor = LGBMRegressor(random_state=15926, verbose=-1),
                 lags      = 24, # Este valor se modifica durante la búsqueda
             )

# Lags utilizados como predictores
lags_grid = [24, [1, 2, 3, 23, 24, 25, 47, 48, 49]]

# Espacio de búsqueda de hiperparámetros
def search_space(trial):
    search_space  = {
        'n_estimators'  : trial.suggest_int('n_estimators', 600, 1200, step=100),
        'max_depth'     : trial.suggest_int('max_depth', 3, 10, step=1),
        'learning_rate' : trial.suggest_float('learning_rate', 0.01, 0.5),
        'reg_alpha'     : trial.suggest_float('reg_alpha', 0, 1, step=0.1),
        'reg_lambda'    : trial.suggest_float('reg_lambda', 0, 1, step=0.1),
    } 
    return search_space

resultados_busqueda, frozen_trial = bayesian_search_forecaster(
                                        forecaster         = forecaster,
                                        y                  = datos.loc[:fin_validacion, 'Demand'],
                                        steps              = 24,
                                        metric             = 'mean_absolute_error',
                                        search_space       = search_space,
                                        lags_grid          = lags_grid,
                                        initial_train_size = len(datos[:fin_train]),
                                        refit              = False,
                                        n_trials           = 20, # Aumentar para una búsqueda más exhaustiva
                                        random_state       = 123,
                                        return_best        = True,
                                        n_jobs             = 'auto',
                                        verbose            = False,
                                        show_progress      = True
                                    )
Number of models compared: 40,
         20 bayesian search in each lag configuration.
`Forecaster` refitted using the best-found lags and parameters, and the whole data set: 
  Lags: [ 1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24] 
  Parameters: {'n_estimators': 1100, 'max_depth': 10, 'learning_rate': 0.07087975104890648, 'reg_alpha': 0.8, 'reg_lambda': 0.2}
  Backtesting metric: 231.1119966842797

In [22]:
# Resultados de la búsqueda
# ==============================================================================
resultados_busqueda.head(10)
Out[22]:
lags params mean_absolute_error n_estimators max_depth learning_rate reg_alpha reg_lambda
15 [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14... {'n_estimators': 1100, 'max_depth': 10, 'learn... 231.111997 1100.0 10.0 0.070880 0.8 0.2
11 [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14... {'n_estimators': 1200, 'max_depth': 8, 'learni... 233.187901 1200.0 8.0 0.071589 1.0 0.0
14 [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14... {'n_estimators': 1200, 'max_depth': 10, 'learn... 233.293696 1200.0 10.0 0.074394 0.8 0.2
17 [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14... {'n_estimators': 1100, 'max_depth': 9, 'learni... 234.490443 1100.0 9.0 0.071958 0.7 0.2
16 [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14... {'n_estimators': 1100, 'max_depth': 9, 'learni... 239.766898 1100.0 9.0 0.167380 0.8 0.2
10 [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14... {'n_estimators': 1200, 'max_depth': 8, 'learni... 242.800676 1200.0 8.0 0.038751 0.9 0.0
0 [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14... {'n_estimators': 1000, 'max_depth': 5, 'learni... 242.914116 1000.0 5.0 0.121157 0.6 0.7
6 [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14... {'n_estimators': 600, 'max_depth': 6, 'learnin... 243.588252 600.0 6.0 0.221123 0.5 0.4
18 [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14... {'n_estimators': 900, 'max_depth': 7, 'learnin... 243.741637 900.0 7.0 0.158435 1.0 0.1
3 [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14... {'n_estimators': 1100, 'max_depth': 4, 'learni... 244.666519 1100.0 4.0 0.095971 0.5 0.5

Al indicar return_best = True, el objeto forecaster se actualiza automáticamente con la mejor configuración encontrada y se entrena con el conjunto de datos completo. Este modelo final puede utilizarse para obtener predicciones sobre nuevos datos.

In [23]:
# Mejor modelo
# ==============================================================================
forecaster
Out[23]:
================= 
ForecasterAutoreg 
================= 
Regressor: LGBMRegressor(learning_rate=0.07087975104890648, max_depth=10,
              n_estimators=1100, random_state=15926, reg_alpha=0.8,
              reg_lambda=0.2, verbose=-1) 
Lags: [ 1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24] 
Transformer for y: None 
Transformer for exog: None 
Window size: 24 
Weight function included: False 
Differentiation order: None 
Exogenous included: False 
Type of exogenous variable: None 
Exogenous variables names: None 
Training range: [Timestamp('2012-01-01 00:00:00'), Timestamp('2014-11-30 23:00:00')] 
Training index type: DatetimeIndex 
Training index frequency: H 
Regressor parameters: {'boosting_type': 'gbdt', 'class_weight': None, 'colsample_bytree': 1.0, 'importance_type': 'split', 'learning_rate': 0.07087975104890648, 'max_depth': 10, 'min_child_samples': 20, 'min_child_weight': 0.001, 'min_split_gain': 0.0, 'n_estimators': 1100, 'n_jobs': None, 'num_leaves': 31, 'objective': None, 'random_state': 15926, 'reg_alpha': 0.8, 'reg_lambda': 0.2, 'subsample': 1.0, 'subsample_for_bin': 200000, 'subsample_freq': 0, 'verbose': -1} 
fit_kwargs: {} 
Creation date: 2023-12-10 19:10:00 
Last fit date: 2023-12-10 19:19:33 
Skforecast version: 0.11.0 
Python version: 3.11.4 
Forecaster id: None 

Backtesting con datos de test

Una vez identificada la mejor combinación de hiperparámetros utilizando los datos de validación, se evalúa la capacidad predictiva del modelo cuando se aplica al conjunto de test. Se recomienda revisar la documentación de la función backtesting_forecaster() para comprender mejor sus capacidades. Esto ayudará a utilizar todo su potencial para analizar la capacidad predictiva del modelo.

In [25]:
# Backtest modelo final
# ==============================================================================
metrica, predicciones = backtesting_forecaster(
                            forecaster         = forecaster,
                            y                  = datos['Demand'],
                            steps              = 24,
                            metric             = 'mean_absolute_error',
                            initial_train_size = len(datos[:fin_validacion]),
                            refit              = False,
                            n_jobs             = 'auto',
                            verbose            = False,
                            show_progress      = True
                        )

print(f"Error de backtest (MAE): {metrica}")
predicciones.head()
Error de backtest (MAE): 210.06355860742855
Out[25]:
pred
2014-12-01 00:00:00 5586.288619
2014-12-01 01:00:00 5592.875437
2014-12-01 02:00:00 5567.013856
2014-12-01 03:00:00 5553.220511
2014-12-01 04:00:00 5521.734716

Tras la optimización de lags e hiperparámetros, se ha conseguido reducir el error de predicción de forma notable.

Forecasting con variables exógenas

Hasta ahora, sólo se han utilizado como predictores los valores pasados (lags) de la serie temporal. Sin embargo, es posible incluir otras variables como predictores. Estas variables se conocen como variables exógenas (features) y su uso puede mejorar la capacidad predictiva del modelo. Un punto muy importante que hay que tener en cuenta es que los valores de las variables exógenas deben conocerse en el momento de la predicción.

Ejemplos habituales de variables exógenas son aquellas obtenidas del calendario, como el día de la semana, el mes, el año o los días festivos. Las variables meteorológicas como la temperatura, la humedad y el viento también entran en esta categoría, al igual que las variables económicas como la inflación y los tipos de interés.

Warning

Las variables exógenas deben conocerse en el momento de la predicción. Por ejemplo, si se utiliza la temperatura como variable exógena, el valor de la temperatura para la hora siguiente debe conocerse en el momento de la previsión. Si no se conoce el valor de la temperatura, la predicción no será posible.

Las variables meteorológicas deben utilizarse con precaución. Cuando el modelo se pone en producción, las condiciones meteorológicas futuras no se conocen, sino que son predicciones realizadas por los servicios meteorológicos. Al tratarse de predicciones, introducen errores en el modelo de previsión. Como consecuencia, es probable que las predicciones del modelo empeoren. Una forma de anticiparse a este problema, y conocer (no evitar) el rendimiento esperado del modelo, es utilizar las previsiones meteorológicas disponibles en el momento en que se entrena el modelo, en lugar de las condiciones reales registradas.

Variables exógenas

A continuación, se crean variables exógenas basadas en información del calendario, las horas de salida y puesta del sol, la temperatura y los días festivos. Estas nuevas variables se añaden a los conjuntos de entrenamiento, validación y test, y se utilizan como predictores en el modelo autorregresivo.

In [26]:
# Variables basadas en el calendario
# ==============================================================================
variables_calendario = pd.DataFrame(index=datos.index)
variables_calendario['mes'] = variables_calendario.index.month
variables_calendario['semana_anyo'] = variables_calendario.index.isocalendar().week
variables_calendario['dia_semana'] = variables_calendario.index.day_of_week + 1
variables_calendario['hora_dia'] = variables_calendario.index.hour + 1

# Variables basadas en la luz solar
# ==============================================================================
location = LocationInfo(
    "Melbourne",
    "Australia",
    latitude=-37.8,
    longitude=144.95,
    timezone='Australia/Melbourne'
)
hora_amanecer = [
    sun(location.observer, date=date, tzinfo=location.timezone)['sunrise'].hour
    for date in datos.index
]
hora_anochecer = [
    sun(location.observer, date=date, tzinfo=location.timezone)['sunset'].hour
    for date in datos.index
]
variables_solares = pd.DataFrame({
                         'hora_amanecer': hora_amanecer,
                         'hora_anochecer': hora_anochecer}, 
                         index = datos.index
                     )
variables_solares['horas_luz_solar'] = (
    variables_solares['hora_anochecer'] - variables_solares['hora_amanecer']
)
variables_solares['es_de_dia'] = np.where(
                                    (datos.index.hour >= variables_solares['hora_amanecer']) & \
                                    (datos.index.hour < variables_solares['hora_anochecer']),
                                    1,
                                    0
                                )

# Variables basadas en festivos
# ==============================================================================
variables_festivos = datos[['Holiday']].astype(int)
variables_festivos['holiday_dia_anterior'] = variables_festivos['Holiday'].shift(24)
variables_festivos['holiday_dia_siguiente'] = variables_festivos['Holiday'].shift(-24)

# Temperature features
# ==============================================================================
variables_temp = datos[['Temperature']].copy()
variables_temp['temp_roll_mean_1_dia'] = variables_temp['Temperature'].rolling(24, closed='left').mean()
variables_temp['temp_roll_mean_7_dia'] = variables_temp['Temperature'].rolling(24*7, closed='left').mean()
variables_temp['temp_roll_max_1_dia'] = variables_temp['Temperature'].rolling(24, closed='left').max()
variables_temp['temp_roll_min_1_dia'] = variables_temp['Temperature'].rolling(24, closed='left').min()
variables_temp['temp_roll_max_7_dia'] = variables_temp['Temperature'].rolling(24*7, closed='left').max()
variables_temp['temp_roll_min_7_dia'] = variables_temp['Temperature'].rolling(24*7, closed='left').min()


# Merge all exogenous variables
# ==============================================================================
variables_exogenas = pd.concat([
                         variables_calendario,
                         variables_solares,
                         variables_temp,
                         variables_festivos
                     ], axis=1)

variables_exogenas.head(4)
Out[26]:
mes semana_anyo dia_semana hora_dia hora_amanecer hora_anochecer horas_luz_solar es_de_dia Temperature temp_roll_mean_1_dia temp_roll_mean_7_dia temp_roll_max_1_dia temp_roll_min_1_dia temp_roll_max_7_dia temp_roll_min_7_dia Holiday holiday_dia_anterior holiday_dia_siguiente
Time
2012-01-01 00:00:00 1 52 7 1 6 20 14 0 27.000 NaN NaN NaN NaN NaN NaN 1 NaN 1.0
2012-01-01 01:00:00 1 52 7 2 6 20 14 0 29.650 NaN NaN NaN NaN NaN NaN 1 NaN 1.0
2012-01-01 02:00:00 1 52 7 3 6 20 14 0 31.650 NaN NaN NaN NaN NaN NaN 1 NaN 1.0
2012-01-01 03:00:00 1 52 7 4 6 20 14 0 32.675 NaN NaN NaN NaN NaN NaN 1 NaN 1.0

Algunos aspectos del calendario, como las horas o los días, son cíclicos. Por ejemplo, la hora del día va de 0 a 23 horas. Aunque se interpreta como una variable continua, la hora 23:00 sólo dista una hora de las 00:00. Lo mismo ocurre con los meses del año, ya que diciembre sólo dista un mes de enero. El uso de funciones trigonométricas como seno y coseno permite representar patrones cíclicos y evitar incoherencias en la representación de los datos. Este enfoque se conoce como codificación cíclica y puede mejorar significativamente la capacidad predictiva de los modelos.

In [27]:
# Codificación cíclica de las variables de calendario y luz solar
# ==============================================================================
def codificacion_ciclica(datos: pd.Series, longitud_ciclo: int) -> pd.DataFrame:
    """
    Codifica una variable cíclica con dos nuevas variables: seno y coseno.
    Se asume que el valor mínimo de la variable es 0. El valor máximo de la
    variable se pasa como argumento.
      
    Parameters
    ----------
    datos : pd.Series
        Serie con la variable a codificar.
    longitud_ciclo : int
        La longitud del ciclo. Por ejemplo, 12 para meses, 24 para horas, etc.
        Este valor se utiliza para calcular el ángulo del seno y coseno.

    Returns
    -------
    resultado : pd.DataFrame
        Dataframe con las dos nuevas características: seno y coseno.

    """

    seno = np.sin(2 * np.pi * datos/longitud_ciclo)
    coseno = np.cos(2 * np.pi * datos/longitud_ciclo)
    resultado =  pd.DataFrame({
                  f"{datos.name}_seno": seno,
                  f"{datos.name}_coseno": coseno
              })

    return resultado


mes_encoded = codificacion_ciclica(variables_exogenas['mes'], longitud_ciclo=12)
semana_anyo_encoded = codificacion_ciclica(variables_exogenas['semana_anyo'], longitud_ciclo=52)
dia_semana_encoded = codificacion_ciclica(variables_exogenas['dia_semana'], longitud_ciclo=7)
hora_dia_encoded = codificacion_ciclica(variables_exogenas['hora_dia'], longitud_ciclo=24)
hora_amanecer_encoded = codificacion_ciclica(variables_exogenas['hora_amanecer'], longitud_ciclo=24)
hora_anochecer_encoded = codificacion_ciclica(variables_exogenas['hora_anochecer'], longitud_ciclo=24)

variables_ciclicas = pd.concat([
                            mes_encoded,
                            semana_anyo_encoded,
                            dia_semana_encoded,
                            hora_dia_encoded,
                            hora_amanecer_encoded,
                            hora_anochecer_encoded
                        ], axis=1)  
variables_exogenas = pd.concat([variables_exogenas, variables_ciclicas], axis=1)
variables_exogenas.head(3)
Out[27]:
mes semana_anyo dia_semana hora_dia hora_amanecer hora_anochecer horas_luz_solar es_de_dia Temperature temp_roll_mean_1_dia ... semana_anyo_seno semana_anyo_coseno dia_semana_seno dia_semana_coseno hora_dia_seno hora_dia_coseno hora_amanecer_seno hora_amanecer_coseno hora_anochecer_seno hora_anochecer_coseno
Time
2012-01-01 00:00:00 1 52 7 1 6 20 14 0 27.00 NaN ... 0.0 1.0 -2.449294e-16 1.0 0.258819 0.965926 1.0 6.123234e-17 -0.866025 0.5
2012-01-01 01:00:00 1 52 7 2 6 20 14 0 29.65 NaN ... 0.0 1.0 -2.449294e-16 1.0 0.500000 0.866025 1.0 6.123234e-17 -0.866025 0.5
2012-01-01 02:00:00 1 52 7 3 6 20 14 0 31.65 NaN ... 0.0 1.0 -2.449294e-16 1.0 0.707107 0.707107 1.0 6.123234e-17 -0.866025 0.5

3 rows × 30 columns

✎ Nota

Con fines ilustrativos, se han creado manualmente las variables del calendario y su codificación cíclica. Sin embargo, para un enfoque más eficiente y automatizado, se recomienda utiliar las clases DatetimeFeatures y CyclicalFeatures disponibles en la librería Feature-engine y pasarlas directamente al argumento transformer_exog del forecaster.

En muchos casos, las variables exógenas no son independientes. Más bien, su efecto sobre la variable objetivo depende del valor de otras variables. Por ejemplo, el efecto de la temperatura en la demanda de electricidad depende de la hora del día. La interacción entre las variables exógenas puede captarse mediante nuevas variables que se obtienen multiplicando entre sí las variables existentes. Estas interacciones se obtienen fácilmente con la clase PolynomialFeatures de scikit-learn.

In [28]:
# Interacción entre variables exógenas
# ==============================================================================
transformer_poly = PolynomialFeatures(
                       degree           = 2,
                       interaction_only = True,
                       include_bias     = False
                   ).set_output(transform="pandas")
poly_cols = [
    'mes_seno',
    'mes_coseno',
    'semana_anyo_seno',
    'semana_anyo_coseno',
    'dia_semana_seno',
    'dia_semana_coseno',
    'hora_dia_seno',
    'hora_dia_coseno',
    'hora_amanecer_seno',
    'hora_amanecer_coseno',
    'hora_anochecer_seno',
    'hora_anochecer_coseno',
    'horas_luz_solar',
    'es_de_dia',
    'holiday_dia_anterior',
    'holiday_dia_siguiente',
    'temp_roll_mean_1_dia',
    'temp_roll_mean_7_dia',
    'temp_roll_max_1_dia',
    'temp_roll_min_1_dia',
    'temp_roll_max_7_dia',
    'temp_roll_min_7_dia',
    'Temperature',
    'Holiday'
]

variables_poly = transformer_poly.fit_transform(variables_exogenas[poly_cols].dropna())
variables_poly = variables_poly.drop(columns=poly_cols)
variables_poly.columns = [f"poly_{col}" for col in variables_poly.columns]
variables_poly.columns = variables_poly.columns.str.replace(" ", "__")
variables_exogenas = pd.concat([variables_exogenas, variables_poly], axis=1)
variables_exogenas.head(3)
Out[28]:
mes semana_anyo dia_semana hora_dia hora_amanecer hora_anochecer horas_luz_solar es_de_dia Temperature temp_roll_mean_1_dia ... poly_temp_roll_min_1_dia__temp_roll_max_7_dia poly_temp_roll_min_1_dia__temp_roll_min_7_dia poly_temp_roll_min_1_dia__Temperature poly_temp_roll_min_1_dia__Holiday poly_temp_roll_max_7_dia__temp_roll_min_7_dia poly_temp_roll_max_7_dia__Temperature poly_temp_roll_max_7_dia__Holiday poly_temp_roll_min_7_dia__Temperature poly_temp_roll_min_7_dia__Holiday poly_Temperature__Holiday
Time
2012-01-01 00:00:00 1 52 7 1 6 20 14 0 27.00 NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
2012-01-01 01:00:00 1 52 7 2 6 20 14 0 29.65 NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
2012-01-01 02:00:00 1 52 7 3 6 20 14 0 31.65 NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN

3 rows × 306 columns

In [29]:
# Selección de variables exógenas incluidas en el modelo
# ==============================================================================
exog_features = []

# Columnas que terminan con _seno o _coseno son seleccionadas
exog_features.extend(variables_exogenas.filter(regex='_seno$|_coseno$').columns.tolist())

# Columnas que empiezan con temp_ son seleccionadas
exog_features.extend(variables_exogenas.filter(regex='^temp_.*').columns.tolist())

# Columnas que empiezan con holiday_ son seleccionadas
exog_features.extend(variables_exogenas.filter(regex='^holiday_.*').columns.tolist())

# Incluir temperatura y festivos
exog_features.extend(['Temperature', 'Holiday'])
In [30]:
# Combinar variables exógenas seleccionadas con la serie temporal en un único DataFrame
# ==============================================================================
datos = datos[['Demand']].merge(
           variables_exogenas[exog_features],
           left_index=True,
           right_index=True,
           how='left'
       )
# Debido a la creación de medias móviles, existen valores ausentes al principio
# de la serie. Debido a holiday_dia_siguiente existen valores ausentes al final.
datos = datos.dropna()
datos = datos.astype('float32')

# Split data into train-val-test
datos_train = datos.loc[: fin_train, :].copy()
datos_val   = datos.loc[fin_train:fin_validacion, :].copy()
datos_test  = datos.loc[fin_validacion:, :].copy()

Los hiperparámetros y lags identificados como óptimos en la sección anterior se utilizan de nuevo para entrenar el modelo, pero esta vez, las variables exógenas también se incluyen como predictores.

In [31]:
# Crear forecaster
# ==============================================================================
params = {
    'n_estimators': 1100,
    'max_depth': 10,
    'learning_rate': 0.0709,
    'reg_alpha': 0.8,
    'reg_lambda': 0.2,
    'random_state': 15926,
    'verbose': -1
}

forecaster = ForecasterAutoreg(
                 regressor = LGBMRegressor(**params),
                 lags      = 24
             )
In [32]:
# Backtesting
# ==============================================================================
metrica, predicciones = backtesting_forecaster(
                            forecaster         = forecaster,
                            y                  = datos['Demand'],
                            exog               = datos[exog_features],
                            steps              = 24,
                            metric             = 'mean_absolute_error',
                            initial_train_size = len(datos[:fin_validacion]),
                            refit              = False,
                            n_jobs             = 'auto',
                            verbose            = False,
                            show_progress      = True
                       )

print(f"Error de backtest (MAE): {metrica:.2f}")
predicciones.head()
Error de backtest (MAE): 132.77
Out[32]:
pred
2014-12-01 00:00:00 5552.413345
2014-12-01 01:00:00 5489.721561
2014-12-01 02:00:00 5465.588132
2014-12-01 03:00:00 5451.111073
2014-12-01 04:00:00 5525.933964

La inclusión de variables exógenas como predictores mejora aún más la capacidad predictiva del modelo.

Intervalos de predicción

Un intervalo de predicción define el intervalo dentro del cual es de esperar que se encuentre el verdadero valor de la variable respuesta con una determinada probabilidad. Rob J Hyndman y George Athanasopoulos, listan en su libro Forecasting: Principles and Practice mútiples formas de estimar intervalos de predicción, la mayoría los cuales requieren que los resudios (errores) del modelo se distribuyan de forma normal. Cuando no se puede asumir esta propiedad, se puede recurrir a bootstrapping, que solo asume que los residuos no están correlacionados. Este es uno de los métodos disponibles en skforecast. Se pude encontrar una explicación más detallada de los intervalos de predicción en Probabilistic forecasting: prediction intervals and prediction distribution.

El siguiente código muestra cómo generar intervalos de predicción con un modelo autorregresivo. La función prediction_interval() se utiliza para estimar los intervalos para cada step predicho. Después, se utiliza la función backtesting_forecaster() para generar los intervalos de predicción de todo el conjunto de test. El argumento interval se utiliza para especificar la probabilidad de cobertura deseada de los intervalos. En este caso, interval se establece en [10, 90], lo que significa que los intervalos se calculan con los percentiles 10 y 90, lo que da como resultado una cobertura teórica del 80%. El argumento n_boot se utiliza para especificar el número de iteraciones de bootstraping que se utilizan para estimar los intervalos. Cuanto mayor sea el número de muestras, más precisos serán los intervalos de predicción, pero más tiempo tardará el proceso.

In [33]:
# Crear y entrenar el forecaster
# ==============================================================================
params = {
    'n_estimators': 1100,
    'max_depth': 10,
    'learning_rate': 0.0709,
    'reg_alpha': 0.8,
    'reg_lambda': 0.2,
    'random_state': 15926,
    'verbose': -1
}

forecaster = ForecasterAutoreg(
                 regressor = LGBMRegressor(**params),
                 lags      = 24
             )

forecaster.fit(
    y    = datos.loc[:fin_train, 'Demand'],
    exog = datos.loc[:fin_train, exog_features]
)
In [34]:
# Predecir intervalos
# ==============================================================================
# Dado que el modelo ha sido entrenado con variables exógenas, estas deben ser
# proporcionadas también en la predicción.
predicciones = forecaster.predict_interval(
                  exog     = datos.loc[fin_train:, exog_features],
                  steps    = 24,
                  interval = [10, 90], 
              )
predicciones.head()
Out[34]:
pred lower_bound upper_bound
2014-01-01 00:00:00 3600.929322 3572.643772 3624.916926
2014-01-01 01:00:00 3700.916873 3665.079895 3725.453531
2014-01-01 02:00:00 3758.736410 3727.163598 3807.674444
2014-01-01 03:00:00 3751.756941 3718.551059 3812.145448
2014-01-01 04:00:00 3750.495804 3717.478448 3827.383248

Por defecto, los intervalos se calculan utilizando los residuos in-sample (residuos del conjunto de entrenamiento). Sin embargo, esto puede dar lugar a intervalos demasiado estrechos (demasiado optimistas). Para evitarlo, se utiliza el método set_out_sample_residuals() para almacenar residuos out-sample calculados mediante backtesting con un conjunto de validación

In [35]:
# Backtesting sobre los datos de validación para obtener los residuos out-sample
# ==============================================================================
_, predicciones_val = backtesting_forecaster(
                         forecaster         = forecaster,
                         y                  = datos.loc[:fin_validacion, 'Demand'], # Train + Validation
                         exog               = datos.loc[:fin_validacion, exog_features],
                         steps              = 24,
                         metric             = 'mean_absolute_error',
                         initial_train_size = len(datos.loc[:fin_train]),
                         refit              = False,
                         n_jobs             = 'auto',
                         verbose            = False,
                         show_progress      = True
                     )
residuos = predicciones_val['pred'] - datos.loc[predicciones_val.index, 'Demand']
residuos = residuos.dropna().to_numpy()

# Almacenar residuos out-sample en el forecaster
# ==============================================================================
forecaster.set_out_sample_residuals(residuals=residuos)

A continuación, se repite el proceso de backtesting para generar los intervalos de predicción en el conjunto de test. Se indica el argumento in_sample_residuals en False para que se utilicen los residuos out-sample almacenados previamente.

In [36]:
# Backtest con intervalos de predicción para el conjunto de test utilizando los
# residuos out-sample
# ==============================================================================
metrica, predicciones = backtesting_forecaster(
                            forecaster          = forecaster,
                            y                   = datos['Demand'],
                            exog                = datos[exog_features],
                            steps               = 24,
                            metric              = 'mean_absolute_error',
                            initial_train_size  = len(datos.loc[:fin_validacion]),
                            refit               = False,
                            interval            = [10, 90],
                            n_boot              = 1000,
                            in_sample_residuals = False, # Se utilizan los residuos out-sample
                            n_jobs              = 'auto',
                            verbose             = False,
                            show_progress       = True
                        )
predicciones.head(5)
Out[36]:
pred lower_bound upper_bound
2014-12-01 00:00:00 5552.413345 5393.982951 5751.242281
2014-12-01 01:00:00 5489.721561 5254.015126 5840.891039
2014-12-01 02:00:00 5465.588132 5196.604384 5974.429540
2014-12-01 03:00:00 5451.111073 5168.466441 6054.435283
2014-12-01 04:00:00 5525.933964 5220.894376 6250.380375
In [38]:
# Gráfico intervalos de predicción vs valores reales
# ==============================================================================
fig = go.Figure([
    go.Scatter(
        name='Predicción',
        x=predicciones.index,
        y=predicciones['pred'],
        mode='lines',
    ),
    go.Scatter(
        name='Valor real',
        x=datos_test.index,
        y=datos_test['Demand'],
        mode='lines',
    ),
    go.Scatter(
        name='Upper Bound',
        x=predicciones.index,
        y=predicciones['upper_bound'],
        mode='lines',
        marker=dict(color="#444"),
        line=dict(width=0),
        showlegend=False
    ),
    go.Scatter(
        name='Lower Bound',
        x=predicciones.index,
        y=predicciones['lower_bound'],
        marker=dict(color="#444"),
        line=dict(width=0),
        mode='lines',
        fillcolor='rgba(68, 68, 68, 0.3)',
        fill='tonexty',
        showlegend=False
    )
])
fig.update_layout(
    title="Predicción con intervalos de confianza vs valor real",
    xaxis_title="Fecha",
    yaxis_title="Demanda (MWh)",
    width=800,
    height=400,
    margin=dict(l=20, r=20, t=35, b=20),
    hovermode="x",
    legend=dict(
        orientation="h",
        yanchor="top",
        y=1.1,
        xanchor="left",
        x=0.001
    )
)
fig.show()
In [39]:
# Cobertura del intervalo predicho
# ==============================================================================
dentro_intervalo = np.where(
                       (datos.loc[fin_validacion:, 'Demand'] >= predicciones['lower_bound']) & \
                       (datos.loc[fin_validacion:, 'Demand'] <= predicciones['upper_bound']),
                       True,
                       False
                   )

cobertura = dentro_intervalo.mean()
print(f"Cobertura del intervalo predicho: {round(100*cobertura, 2)} %")
Cobertura del intervalo predicho: 93.39 %

Una cobertura empírica superior al 80% esperado sugiere dos posibilidades: o bien los intervalos se estiman de forma conservadora, lo que da lugar a rangos más amplios de lo necesario, o bien el modelo capta la variabilidad de los datos de forma más eficaz de lo previsto. Ampliar el conjunto de entrenamiento para incluir más datos puede ayudar a distinguir entre estas dos posibilidades.

Explicabilidad del modelo

Debido a la naturaleza compleja de muchos de los actuales modelos de machine learning, a menudo funcionan como cajas negras, lo que dificulta entender por qué han hecho una predicción u otra. Las técnicas de explicabilidad pretenden desmitificar estos modelos, proporcionando información sobre su funcionamiento interno y ayudando a generar confianza, mejorar la transparencia y cumplir los requisitos normativos en diversos ámbitos. Mejorar la explicabilidad de los modelos no sólo ayuda a comprender su comportamiento, sino también a identificar sesgos, mejorar su rendimiento y permitir a las partes interesadas tomar decisiones más informadas basadas en los conocimientos del machine learning

Skforecast es compatible con algunos de los métodos de explicabilidad más populares: model-specific feature importances, SHAP values, and partial dependence plots.

In [40]:
# Crear y entrenar el forecaster
# ==============================================================================
params = {
    'n_estimators': 1100,
    'max_depth': 10,
    'learning_rate': 0.0709,
    'reg_alpha': 0.8,
    'reg_lambda': 0.2,
    'random_state': 15926,
    'verbose': -1
}
forecaster = ForecasterAutoreg(
                 regressor = LGBMRegressor(**params),
                 lags      = 24
             )
forecaster.fit(
    y    = datos.loc[:fin_validacion, 'Demand'],
    exog = datos.loc[:fin_validacion, exog_features]
)

Explicabilidad propia del modelo

La explicabilidad propia del modelo se refiere a la capacidad inherente de un modelo para ser interpretado y comprendido sin la necesidad de técnicas adicionales.

In [41]:
# Importancia de los predictores en el modelo
# ==============================================================================
feature_importances = forecaster.get_feature_importances()
feature_importances.sort_values(by='importance', ascending=False).head(10)
Out[41]:
feature importance
0 lag_1 2429
110 Temperature 1300
1 lag_2 1219
23 lag_24 1182
105 temp_roll_min_1_dia 821
22 lag_23 814
2 lag_3 788
102 temp_roll_mean_1_dia 727
104 temp_roll_max_1_dia 718
21 lag_22 659

Warning

El método get_feature_importances() sólo devuelve valores si el regresor del forecaster tiene el atributo coef_ o feature_importances_, que son los estándares en scikit-learn.

Shap values

Los valores SHAP (SHapley Additive exPlanations) son un método muy utilizado para explicar los modelos de machine learning, ya que ayudan a comprender cómo influyen las variables y los valores en las predicciones de forma visual y cuantitativa.

Se puede obtener un análisis SHAP a partir de modelos skforecast con sólo dos elementos:

  • El regresor interno del forecaster.

  • Las matrices de entrenamiento creadas a partir de la serie temporal y variables exógenas, utilizadas para ajustar el pronosticador.

Aprovechando estos dos componentes, los usuarios pueden crear explicaciones interpretables para sus modelos de skforecast. Estas explicaciones pueden utilizarse para verificar la fiabilidad del modelo, identificar los factores más significativos que contribuyen a las predicciones y comprender mejor la relación subyacente entre las variables de entrada y la variable objetivo.

In [42]:
# Matrices de entrenamiento utilizadas por el forecaster para ajustar el regresor interno
# ==============================================================================
X_train, y_train = forecaster.create_train_X_y(
                       y    = datos_train['Demand'],
                       exog = datos_train[exog_features]
                   )

display(X_train.head(3))
display(y_train.head(3))
lag_1 lag_2 lag_3 lag_4 lag_5 lag_6 lag_7 lag_8 lag_9 lag_10 ... temp_roll_mean_1_dia temp_roll_mean_7_dia temp_roll_max_1_dia temp_roll_min_1_dia temp_roll_max_7_dia temp_roll_min_7_dia holiday_dia_anterior holiday_dia_siguiente Temperature Holiday
Time
2012-01-09 00:00:00 4886.074219 4664.689453 4390.567871 3873.921631 3498.672607 3348.528809 3458.686523 3739.490967 3630.717285 3897.311035 ... 20.281250 22.389137 27.075001 14.8 39.525002 14.35 0.0 0.0 19.200001 0.0
2012-01-09 01:00:00 4958.630371 4886.074219 4664.689453 4390.567871 3873.921631 3498.672607 3348.528809 3458.686523 3739.490967 3630.717285 ... 20.223959 22.308035 27.075001 14.8 39.525002 14.35 0.0 0.0 20.525000 0.0
2012-01-09 02:00:00 5039.643066 4958.630371 4886.074219 4664.689453 4390.567871 3873.921631 3498.672607 3348.528809 3458.686523 3739.490967 ... 20.141666 22.224852 27.075001 14.8 39.525002 14.35 0.0 0.0 21.100000 0.0

3 rows × 112 columns

Time
2012-01-09 00:00:00    4958.630371
2012-01-09 01:00:00    5039.643066
2012-01-09 02:00:00    5090.203125
Freq: H, Name: y, dtype: float32
In [43]:
# Crear SHAP explainer (para modelos basados en árboles)
# ==============================================================================
explainer = shap.TreeExplainer(forecaster.regressor)

# Se selecciona una muestra del 50% de los datos para acelerar el cálculo
rng = np.random.default_rng(seed=785412)
sample = rng.choice(X_train.index, size=int(len(X_train)*0.5), replace=False)
X_train_sample = X_train.loc[sample, :]

shap_values = explainer.shap_values(X_train_sample)

✎ Note

La librería Shap cuenta con varios Explainers, cada uno diseñado para un tipo de modelo diferente. El shap.TreeExplainer explainer se utiliza para modelos basados en árboles, como el LGBMRegressor utilizado en este ejemplo. Para más información, consultar la documentación de SHAP.
In [44]:
# Shap summary plot (top 10)
# ==============================================================================
shap.summary_plot(shap_values, X_train_sample, max_display=10, show=False)
fig, ax = plt.gcf(), plt.gca()
ax.set_title("SHAP Summary plot")
ax.tick_params(labelsize=8)
fig.set_size_inches(8, 4.5)
In [45]:
# Force plot para la primera observación de la muestra
# ==============================================================================
shap.force_plot(explainer.expected_value, shap_values[0,:], X_train_sample.iloc[0,:])
Out[45]:
Visualization omitted, Javascript library not loaded!
Have you run `initjs()` in this notebook? If this notebook was from another user you must also trust this notebook (File -> Trust notebook). If you are viewing this notebook on github the Javascript has been stripped for security. If you are using JupyterLab this error is because a JupyterLab extension has not yet been written.

Selección de predictores

La selección de predictores (feature selection) es el proceso de identificar un subconjunto de predictores relevantes para su uso en la creación del modelo. Es un paso importante en el proceso de machine learning, ya que puede ayudar a reducir el sobreajuste, mejorar la precisión del modelo y reducir el tiempo de entrenamiento. Dado que los regresores subyacentes de skforecast siguen la API de scikit-learn, es posible utilizar los métodos de selección de predictores disponibles en scikit-learn. Dos de los métodos más populares son Recursive Feature Elimination y Sequential Feature Selection.

💡 Tip

La selección de predictores es una herramienta potente para mejorar el rendimiento de los modelos de machine learning. Sin embargo, es computacionalmente costosa y puede requerir mucho tiempo. Dado que el objetivo es encontrar el mejor subconjunto de predictores, no el mejor modelo, no es necesario utilizar todos los datos disponibles ni un modelo muy complejo. En su lugar, se recomienda utilizar un pequeño subconjunto de los datos y un modelo sencillo. Una vez identificados los mejores predictores, el modelo puede entrenarse utilizando todo el conjunto de datos y una configuración más compleja. Por ejemplo, en este caso de uso, el mejor modelo es un LGMBRegressor con 1100 árboles y una profundidad máxima de 10. Sin embargo, para la selección predictores, sólo se utilizan 100 árboles y una profundidad máxima de 5.

In [46]:
# Crear forecaster
# ==============================================================================
forecaster = ForecasterAutoreg(
                 regressor = LGBMRegressor(random_state=123, verbose=-1),
                 lags = 24
             )

# Matrices de entrenamiento utilizadas por el forecaster para ajustar el regresor interno
# ==============================================================================
X_train, y_train = forecaster.create_train_X_y(
                       y    = datos_train['Demand'],
                       exog = datos_train[exog_features]
                   )

El RFECV de scikit-learn empieza entrenando un modelo con todos los predictores disponibles y calculando la importancia de cada uno en base a los atributos como coef_ o feature_importances_. A continuación, se elimina el predictor menos importante y se realiza una validación cruzada para calcular el rendimiento del modelo con los predictores restantes. Este proceso se repite hasta que la eliminación de predictores adicionales no mejora la metrica de rendimiento elegida o se alcanza el min_features_to_select.

El resultado final es un subconjunto de predictores que idealmente equilibra la simplicidad del modelo y su capacidad predictiva, determinada por el proceso de validación cruzada.

In [47]:
# Eliminación recursiva de predictores con validación cruzada
# ==============================================================================
# Muestrear el 50% de los datos para acelerar el cálculo
rng = np.random.default_rng(seed=785412)
sample = rng.choice(X_train.index, size=int(len(X_train)*0.5), replace=False)
X_train_sample = X_train.loc[sample, :]
y_train_sample = y_train.loc[sample]

print(f"Número de predictores disponibles: {X_train.shape[1]}") 
print(f"Número de observaciones disponibles: {X_train.shape[0]}")
print(f"Número de observaciones utilizadas para selección de predictores: {X_train_sample.shape[0]}")

regressor = LGBMRegressor(
                n_estimators = 100,
                max_depth = 5,
                random_state = 15926,
                verbose = -1
            )
selector = RFECV(
    estimator              = regressor,
    step                   = 1,
    cv                     = 3,
    min_features_to_select = 25,
    n_jobs                 = -1
)
selector.fit(X_train_sample, y_train_sample)
predictores_seleccionados = selector.get_feature_names_out()

print("\nEliminación recursiva de predictores")
print("-------------------------------------")
print(f"Número de predictores seleccionados: {selector.n_features_}")
print(f"Predictores seleccionados: \n {predictores_seleccionados}")
Número de predictores disponibles: 112
Número de observaciones disponibles: 17352
Número de observaciones utilizadas para selección de predictores: 8676

Eliminación recursiva de predictores
-------------------------------------
Número de predictores seleccionados: 74
Predictores seleccionados: 
 ['lag_1' 'lag_2' 'lag_3' 'lag_4' 'lag_5' 'lag_6' 'lag_7' 'lag_8' 'lag_9'
 'lag_10' 'lag_11' 'lag_12' 'lag_13' 'lag_14' 'lag_15' 'lag_16' 'lag_17'
 'lag_18' 'lag_19' 'lag_20' 'lag_21' 'lag_22' 'lag_23' 'lag_24'
 'mes_coseno' 'semana_anyo_seno' 'semana_anyo_coseno' 'dia_semana_seno'
 'hora_dia_seno' 'hora_dia_coseno' 'hora_anochecer_coseno'
 'poly_mes_seno__semana_anyo_seno' 'poly_mes_seno__semana_anyo_coseno'
 'poly_mes_seno__dia_semana_seno' 'poly_mes_seno__dia_semana_coseno'
 'poly_mes_seno__hora_dia_seno' 'poly_mes_seno__hora_dia_coseno'
 'poly_mes_coseno__semana_anyo_seno' 'poly_mes_coseno__hora_dia_seno'
 'poly_mes_coseno__hora_dia_coseno'
 'poly_semana_anyo_seno__dia_semana_seno'
 'poly_semana_anyo_seno__dia_semana_coseno'
 'poly_semana_anyo_seno__hora_dia_seno'
 'poly_semana_anyo_seno__hora_dia_coseno'
 'poly_semana_anyo_coseno__dia_semana_seno'
 'poly_semana_anyo_coseno__dia_semana_coseno'
 'poly_semana_anyo_coseno__hora_dia_seno'
 'poly_semana_anyo_coseno__hora_dia_coseno'
 'poly_semana_anyo_coseno__hora_amanecer_seno'
 'poly_dia_semana_seno__dia_semana_coseno'
 'poly_dia_semana_seno__hora_dia_seno'
 'poly_dia_semana_seno__hora_dia_coseno'
 'poly_dia_semana_seno__hora_anochecer_seno'
 'poly_dia_semana_coseno__hora_dia_seno'
 'poly_dia_semana_coseno__hora_dia_coseno'
 'poly_hora_dia_seno__hora_dia_coseno'
 'poly_hora_dia_seno__hora_amanecer_seno'
 'poly_hora_dia_seno__hora_amanecer_coseno'
 'poly_hora_dia_seno__hora_anochecer_seno'
 'poly_hora_dia_seno__hora_anochecer_coseno'
 'poly_hora_dia_coseno__hora_amanecer_seno'
 'poly_hora_dia_coseno__hora_amanecer_coseno'
 'poly_hora_dia_coseno__hora_anochecer_seno'
 'poly_hora_dia_coseno__hora_anochecer_coseno'
 'poly_hora_amanecer_seno__hora_anochecer_seno' 'temp_roll_mean_1_dia'
 'temp_roll_mean_7_dia' 'temp_roll_max_1_dia' 'temp_roll_min_1_dia'
 'temp_roll_max_7_dia' 'temp_roll_min_7_dia' 'holiday_dia_anterior'
 'Temperature' 'Holiday']

El forecater se entrena y reevalúa utilizando el conjunto de predictores seleccionados.

In [48]:
# Crear forecaster
# ==============================================================================
params = {
    'n_estimators': 1100,
    'max_depth': 10,
    'learning_rate': 0.0709,
    'reg_alpha': 0.8,
    'reg_lambda': 0.2,
    'random_state': 15926,
    'verbose': -1
}

lags_seleccionados = [
    int(predictor.replace("lag_", ""))
    for predictor in predictores_seleccionados
    if predictor.startswith("lag_")
]
exog_seleccionadas = [
    predictor
    for predictor in predictores_seleccionados
    if not predictor.startswith("lag_")
]

forecaster = ForecasterAutoreg(
                 regressor = LGBMRegressor(**params),
                 lags      = lags_seleccionados,
             )

# Backtesting con los predictores seleccionados y los datos de test
# ==============================================================================
metrica, predicciones = backtesting_forecaster(
                            forecaster         = forecaster,
                            y                  = datos['Demand'],
                            exog               = datos[exog_seleccionadas],
                            steps              = 24,
                            metric             = 'mean_absolute_error',
                            initial_train_size = len(datos[:fin_validacion]),
                            refit              = False,
                            n_jobs             = 'auto',
                            verbose            = False,
                            show_progress      = True
                        )

print(f"Error de backtest (MAE): {metrica:.2f}")
Error de backtest (MAE): 134.38

El rendimiento del modelo no se ve comprometido, sólo se observa un pequeño aumento del MAE. El número de características se reduce de 112 a 74, lo que simplifica y acelera el entrenamiento del modelo. Además, se reduce el riesgo de sobreajuste porque es menos probable que el modelo aprenda ruido de características irrelevantes.

Como era de esperar, al aumentar el horizonte de predicción de 24 a 36 horas, también lo hace el error.

Forecaster direct multi-step

Los modelos ForecasterAutoreg y ForecasterAutoregCustom siguen una estrategia recursiva en la que cada nueva predicción se basa en las anteriores. Otra estrategia para predecir múltiples valores futuros consiste en entrenar un modelo diferente para cada paso (step) a predecir. Esto se conoce como direct multi-step forecasting y, aunque es más costosa computacionalmente que la estrategia recursiva debido a la necesidad de entrenar múltiples modelos, puede dar mejores resultados.

In [49]:
# Forecaster con le método direct
# ==============================================================================
forecaster = ForecasterAutoregDirect(
                 regressor        = LGBMRegressor(**params),
                 steps            = 24, # horizonte de predicción
                 lags             = lags_seleccionados
             )

# Backtesting
# ==============================================================================
metrica, predicciones = backtesting_forecaster(
                          forecaster         = forecaster,
                          y                  = datos['Demand'],
                          exog               = datos[exog_seleccionadas],
                          steps              = 24,
                          metric             = 'mean_absolute_error',
                          initial_train_size = len(datos[:fin_validacion]),
                          refit              = False,
                          n_jobs             = 'auto',
                          verbose            = False,
                          show_progress      = True
                      )

print(f"Error de backtest (MAE): {metrica:.2f}")
Error de backtest (MAE): 125.51

El modelo direct multi-step supera la capacidad de predicción del modelo recursivo, reduciendo aún más el error medio absoluto. Sin embargo, es importante tener en cuenta su mayor coste computacional a la hora de evaluar si merece la pena aplicarlo.

Predicción diaria anticipada

En el apartado anterior, se evaluó el modelo asumiendo que las predicciones del día siguiente se ejecutan justo al final del día anterior. En la práctica, esto no resulta muy útil ya que, para las primeras horas del día, apenas se dispone de anticipación.

Supóngase ahora que, para poder tener suficiente margen de acción, a las 11:00 horas de cada día se tienen que generar las predicciones del día siguiente. Es decir, a las 11:00 del dia $D$ se tienen que predecir las horas [12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23] de ese mismo día, y las horas [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23] del día $D+1$. Esto implica que se tienen que predecir un total de 36 horas a futuro aunque solo se almacenen las 24 últimas.

Este tipo de evaluación puede realizarse mediante la función backtesting_forecaster() y su parámetro gap. Además, el parámetro allow_incomplete_fold permite excluir la última iteración del análisis si no tiene el mismo tamaño que el número de steps requeridos. El proceso adaptado a este escenario se ejecuta diariamente y consta de los siguientes pasos:

  1. A las 11:00h del primer día del conjunto de test, se predicen las 36 horas siguientes (las 12 horas que quedan del día más las 24 horas del día siguiente).

  2. Se almacenan solo las predicciones del día siguiente, es decir, las 24 últimas.

  3. Se añaden los datos de test hasta las 11:00 del día siguiente.

  4. Se repite el proceso.

De esta forma, a las 11:00 de cada día, el modelo tiene acceso a los valores reales de demanda registrados hasta ese momento.


Warning

Se deben tener en cuenta las siguientes consideraciones:
  • Los datos de entrenamiento deben terminar donde empieza el gap. En este caso, el initial_train_size debe ampliarse en 12 posiciones para que el primer punto a predecir sea el 2014-12-01 12:00:00.
  • En este ejemplo, aunque sólo se almacenan las últimas 24 predicciones (steps) para la evaluación del modelo, el número total de pasos predichos en cada iteración es de 36 (steps + gap).
In [50]:
# Final de initial_train_size + 12 posiciones
# ==============================================================================
datos.iloc[:len(datos.loc[:fin_validacion])+12].tail(2)
Out[50]:
Demand mes_seno mes_coseno semana_anyo_seno semana_anyo_coseno dia_semana_seno dia_semana_coseno hora_dia_seno hora_dia_coseno hora_amanecer_seno ... temp_roll_mean_1_dia temp_roll_mean_7_dia temp_roll_max_1_dia temp_roll_min_1_dia temp_roll_max_7_dia temp_roll_min_7_dia holiday_dia_anterior holiday_dia_siguiente Temperature Holiday
Time
2014-12-01 10:00:00 5084.011230 -2.449294e-16 1.0 -0.354605 0.935016 0.781832 0.62349 2.588190e-01 -0.965926 0.965926 ... 25.589582 18.573214 30.950001 20.0 33.849998 11.15 0.0 0.0 19.90 0.0
2014-12-01 11:00:00 4851.066895 -2.449294e-16 1.0 -0.354605 0.935016 0.781832 0.62349 1.224647e-16 -1.000000 0.965926 ... 25.129168 18.601191 29.500000 19.9 33.849998 11.15 0.0 0.0 19.35 0.0

2 rows × 89 columns

In [55]:
# Forecaster
# ==============================================================================
forecaster = ForecasterAutoreg(
                 regressor = LGBMRegressor(**params),
                 lags      = lags_seleccionados,
             )

# Backtesting con gap
# ==============================================================================
metrica, predicciones = backtesting_forecaster(
                            forecaster            = forecaster,
                            y                     = datos['Demand'],
                            steps                 = 24,
                            metric                = 'mean_absolute_error',
                            initial_train_size    = len(datos.loc[:fin_validacion])+12,
                            gap                   = 12,
                            allow_incomplete_fold = True,
                            refit                 = False,
                            n_jobs                = 'auto',
                            verbose               = False,
                            show_progress         = True  
                        )
print(f"Error de backtest (MAE): {metrica:.2f}")
Error de backtest (MAE): 308.74

Como era de esperar, el error aumenta a medida que el horizonte de previsión pasa de 24 a 36 horas.

Conclusiones

El uso de modelos de gradient boosting ha demostrado ser una potente herramienta para predecir la demanda de energía. Una de las principales ventajas de estos modelos es que pueden incorporar fácilmente variables exógenas, lo que puede mejorar significativamente el poder predictivo del modelo. Además, el uso de técnicas de explicabilidad puede ayudar a comprender visual y cuantitativamente cómo afectan las variables y los valores a las predicciones. Todas estas cuestiones se abordan fácilmente con la biblioteca skforecast.

Forecaster Exogenous Variables MAE backtest
ForecasterEquivalentDate (Baseline) False 308.4
ForecasterAutoreg False 210.1
ForecasterAutoreg True 132.8
ForecasterAutoreg True (Feature Selection) 134.4
ForecasterAutoregDirect True (Feature Selection) 125.5

Información de sesión

In [56]:
import session_info
session_info.show(html=False)
-----
astral              3.2
lightgbm            4.1.0
matplotlib          3.7.2
numpy               1.25.0
optuna              3.2.0
pandas              2.1.4
plotly              5.16.1
session_info        1.0.0
shap                0.43.0
skforecast          0.11.0
sklearn             1.3.2
statsmodels         0.14.0
-----
IPython             8.14.0
jupyter_client      8.3.0
jupyter_core        5.3.1
-----
Python 3.11.4 | packaged by Anaconda, Inc. | (main, Jul  5 2023, 13:47:18) [MSC v.1916 64 bit (AMD64)]
Windows-10-10.0.19045-SP0
-----
Session information updated at 2023-12-10 19:59

Bibliografía

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.

Instrucciones para citar

Cómo citar este documento

Si utiliza este documento o parte de él, por favor, cite la fuente, ¡gracias!

Predicción (forecasting) de la demanda energética con machine learning por Joaquín Amat Rodrigo y Javier Escobar Ortiz, disponible con licencia Attribution-NonCommercial-ShareAlike 4.0 International (CC BY-NC-SA 4.0 DEED) at https://www.cienciadedatos.net/documentos/py29-forecasting-demanda-energia-electrica-python.html

Cómo citar skforecast

Si utiliza skforecast, le agradecemos que cite el software publicado.

Zenodo:

Amat Rodrigo, Joaquin, & Escobar Ortiz, Javier. (2023). skforecast (v0.10.0). Zenodo. https://doi.org/10.5281/zenodo.8382788

APA:

Amat Rodrigo, J., & Escobar Ortiz, J. (2023). skforecast (Version 0.10.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.10.0}, month = {9}, year = {2023}, 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! 😊


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