Predicción (forecasting) de la demanda eléctrica con Python

Predicción (forecasting) de la demanda eléctrica con Python

Joaquín Amat Rodrigo
Abril, 2021 (última actualización Marzo 2022)

Introducción


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 intervalo futuro ($t_{+1}$), ..., ($t_{+n}$)) o un punto alejado en el tiempo ($t_{+n}$). Existen varias estrategias que permiten generar este tipo de predicciones múltiples, la librería skforecast recoge las siguientes:

  • Recursive multi-step forecasting: dado que, para predecir el momento $t_{n}$ se necesita el valor de $t_{n-1}$, y $t_{n-1}$ se desconoce, se sigue un proceso recursivo en el que, cada nueva predicción, hace uso de la predicción anterior.

  • Direct multi-step forecasting: consiste en entrenar un modelo distinto para cada valor futuro (step) que se quiere predecir.

En este documento se muestra un ejemplo de cómo utilizar métodos de forecasting para predecir la demanda eléctrica a nivel horario. En concreto, se hace uso de skforecast, una librería que contiene las clases y funciones necesarias para adaptar cualquier modelo de regresión de scikit-learn a problemas de forecasting.

Para casos de uso más detallados visitar skforecast-examples.

Caso de uso


Se dispone de una serie temporal con la demanda eléctrica (MW) del estado de Victoria (Australia) desde 2011-12-31 al 2014-12-31. Se pretende generar un modelo de forecasting capaz de predecir la demanda energética del día siguiente a nivel horario.

Librerías


Las librerías utilizadas en este documento son:

In [1]:
# Tratamiento de datos
# ==============================================================================
import numpy as np
import pandas as pd

# Gráficos
# ==============================================================================
import matplotlib.pyplot as plt
import seaborn as sns
%matplotlib inline
from statsmodels.graphics.tsaplots import plot_acf
from statsmodels.graphics.tsaplots import plot_pacf
plt.style.use('fivethirtyeight')

# Modelado y Forecasting
# ==============================================================================
from sklearn.linear_model import Ridge
from lightgbm import LGBMRegressor
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import mean_absolute_error
from skforecast.ForecasterAutoreg import ForecasterAutoreg
from skforecast.ForecasterAutoregMultiOutput import ForecasterAutoregMultiOutput
from skforecast.model_selection import grid_search_forecaster
from skforecast.model_selection import backtesting_forecaster

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

Datos


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 del estado de Victoria.

  • Holiday: indicador si el día es festivo (vacaciones).

In [2]:
# Descarga de datos
# ==============================================================================
url = ('https://raw.githubusercontent.com/JoaquinAmatRodrigo/skforecast/master/' +
       'data/vic_elec.csv')
datos = pd.read_csv(url, sep=',')
datos.info()
<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 se ha almacenado 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()

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()
Out[4]:
True
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 de tipo temporal 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 [8]:
# 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[8]:
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 [9]:
# Separación datos train-val-test
# ==============================================================================
datos = datos.loc['2012-01-01 00:00:00': '2014-12-30 23:00:00']
fin_train = '2013-12-31 23:59:00'
fin_validacion = '2014-11-30 23:59:00'
datos_train = datos.loc[: fin_train, :]
datos_val   = datos.loc[fin_train:fin_validacion, :]
datos_test  = datos.loc[fin_validacion:, :]

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


Cuando se quiere generar un modelo de forecasting, es importante representar los valores de la serie temporal. Esto permite identificar patrones tales como tendencias y estacionalidad.

Serie temporal completa

In [10]:
# Gráfico serie temporal
# ==============================================================================
fig, ax = plt.subplots(figsize=(12, 4))
datos_train.Demand.plot(ax=ax, label='entrenamiento', linewidth=1)
datos_val.Demand.plot(ax=ax, label='validación', linewidth=1)
datos_test.Demand.plot(ax=ax, label='test', linewidth=1)
ax.set_title('Demanda eléctrica')
ax.legend();

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 [11]:
# Gráfico serie temporal con zoom
# ==============================================================================
zoom = ('2013-05-01 14:00:00','2013-06-01 14:00:00')

fig = plt.figure(figsize=(12, 6))
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=2)

main_ax.set_title(f'Demanda electricidad: {datos.index.min()}, {datos.index.max()}', fontsize=14)
zoom_ax.set_title(f'Demanda electricidad: {zoom}', fontsize=14)
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.

Estacionalidad anual, semanal y diaria

In [12]:
# Gráfico boxplot para estacionalidad anual
# ==============================================================================
fig, ax = plt.subplots(figsize=(7, 3.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 [13]:
# Gráfico boxplot para estacionalidad semanal
# ==============================================================================
fig, ax = plt.subplots(figsize=(7, 3.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 [14]:
# Gráfico boxplot para estacionalidad diaria
# ==============================================================================
fig, ax = plt.subplots(figsize=(9, 3.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.

Días festivos y no festivos

In [15]:
# Grafico violinplot
# ==============================================================================
fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(7, 3.5))
sns.violinplot(
    x       = 'Demand',
    y       = 'Holiday',
    data    = datos.assign(Holiday = datos.Holiday.astype(str)),
    palette = 'tab10',
    ax      = ax
)
ax.set_title('Distribución de la demanda entre festivos y no festivos')
ax.set_xlabel('demanda')
ax.set_ylabel('festivo');

Los días festivos tienden a tener menor consumo.

Gráficos de autocorrelación

In [16]:
# Gráfico autocorrelación
# ==============================================================================
fig, ax = plt.subplots(figsize=(7, 3))
plot_acf(datos.Demand, ax=ax, lags=60)
plt.show()
In [17]:
# Gráfico autocorrelación parcial
# ==============================================================================
fig, ax = plt.subplots(figsize=(7, 3))
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.

Modelo autoregresivo recursivo


Se crea y entrena un modelo autorregresivo recursivo (ForecasterAutoreg) a partir de un modelo de regresión lineal con penalización Ridge y una ventana temporal de 24 lags. Esto último significa que, para cada predicción, se utilizan como predictores la demanda de las 24 horas anteriores.

Entrenamiento del Forecaster

In [18]:
# Crear y entrenar forecaster
# ==============================================================================
forecaster = ForecasterAutoreg(
                regressor = make_pipeline(StandardScaler(), Ridge()),
                lags      = 24
             )

forecaster.fit(y=datos.loc[:fin_validacion, 'Demand'])
forecaster
Out[18]:
================= 
ForecasterAutoreg 
================= 
Regressor: Pipeline(steps=[('standardscaler', StandardScaler()), ('ridge', Ridge())]) 
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] 
Window size: 24 
Included exogenous: 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: {'standardscaler__copy': True, 'standardscaler__with_mean': True, 'standardscaler__with_std': True, 'ridge__alpha': 1.0, 'ridge__copy_X': True, 'ridge__fit_intercept': True, 'ridge__max_iter': None, 'ridge__normalize': 'deprecated', 'ridge__positive': False, 'ridge__random_state': None, 'ridge__solver': 'auto', 'ridge__tol': 0.001} 
Creation date: 2022-04-03 22:09:57 
Last fit date: 2022-04-03 22:09:57 
Skforecast version: 0.4.3 

Predicción (backtest)


Se evalúa el comportamiento que habría tenido el modelo si se hubiese entrenado con los datos desde 2012-01-01 00:00 al 2014-11-30 23:59 y, después, a las 23:59 de cada día, se predijesen las 24 horas siguientes. A este tipo de evaluación se le conoce como backtesting, y puede aplicarse fácilmente con la función backtesting_forecaster(). Esta función devuelve, además de las predicciones, una métrica de error.

In [19]:
# Backtest
# ==============================================================================
metrica, predicciones = backtesting_forecaster(
                            forecaster = forecaster,
                            y          = datos.Demand,
                            initial_train_size = len(datos.loc[:fin_validacion]),
                            fixed_train_size   = False,
                            steps      = 24,
                            metric     = 'mean_absolute_error',
                            refit      = False,
                            verbose    = True
                        )
Information of backtesting process
----------------------------------
Number of observations used for initial training or as initial window: 25560
Number of observations used for backtesting: 720
    Number of folds: 30
    Number of steps per fold: 24

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

In [20]:
# Gráfico
# ==============================================================================
fig, ax = plt.subplots(figsize=(12, 3.5))
datos.loc[predicciones.index, 'Demand'].plot(ax=ax, linewidth=2, label='test')
predicciones.plot(linewidth=2, label='predicción', ax=ax)
ax.set_title('Predicción vs demanda real')
ax.legend();
In [21]:
# Error backtest
# ==============================================================================
print(f'Error backtest: {metrica}')
Error backtest: 289.51915030381366

Optimización de hiperparámetros (tuning)


En el objeto ForecasterAutoreg entrenado, se han utilizado los primeros 24 lags y un modelo Ridge con los hiperparámetros por defecto. Sin embargo, no hay ninguna razón por la que estos valores sean los más adecuados.

Con el objetivo de identificar la mejor combinación de lags e hiperparámetros, se recurre a un Grid Search. Este proceso consiste en entrenar un modelo con cada combinación de hiperparámetros y lags, y evaluar su capacidad predictiva mediante backtesting. En el proceso de búsqueda, es importante evaluar los modelos utilizando únicamente los datos de validación y no incluir los de test, estos se utilizan solo en último lugar para evaluar al modelo final.

In [22]:
# Grid search de hiperparámetros
# ==============================================================================
forecaster = ForecasterAutoreg(
                    regressor = make_pipeline(StandardScaler(), Ridge()),
                    lags      = 24 # Este valor será remplazado en el grid search
             )

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

# Hiperparámetros del regresor
param_grid = {'ridge__alpha': np.logspace(-3, 5, 10)}

resultados_grid = grid_search_forecaster(
                        forecaster  = forecaster,
                        y           = datos.loc[:fin_validacion, 'Demand'],
                        param_grid  = param_grid,
                        lags_grid   = lags_grid,
                        steps       = 24,
                        metric      = 'mean_absolute_error',
                        refit       = False,
                        initial_train_size = len(datos[:fin_train]),
                        fixed_train_size   = False,
                        return_best = True,
                        verbose     = False
                  )
Number of models compared: 30
loop lags_grid:   0%|                                               | 0/3 [00:00<?, ?it/s]
loop param_grid:   0%|                                             | 0/10 [00:00<?, ?it/s]
loop param_grid:  10%|███▋                                 | 1/10 [00:01<00:10,  1.12s/it]
loop param_grid:  20%|███████▍                             | 2/10 [00:02<00:10,  1.31s/it]
loop param_grid:  30%|███████████                          | 3/10 [00:03<00:08,  1.27s/it]
loop param_grid:  40%|██████████████▊                      | 4/10 [00:05<00:07,  1.30s/it]
loop param_grid:  50%|██████████████████▌                  | 5/10 [00:06<00:06,  1.24s/it]
loop param_grid:  60%|██████████████████████▏              | 6/10 [00:07<00:04,  1.23s/it]
loop param_grid:  70%|█████████████████████████▉           | 7/10 [00:08<00:03,  1.28s/it]
loop param_grid:  80%|█████████████████████████████▌       | 8/10 [00:10<00:02,  1.35s/it]
loop param_grid:  90%|█████████████████████████████████▎   | 9/10 [00:11<00:01,  1.28s/it]
loop param_grid: 100%|████████████████████████████████████| 10/10 [00:12<00:00,  1.26s/it]
loop lags_grid:  33%|█████████████                          | 1/3 [00:12<00:25, 12.70s/it]
loop param_grid:   0%|                                             | 0/10 [00:00<?, ?it/s]
loop param_grid:  10%|███▋                                 | 1/10 [00:01<00:11,  1.25s/it]
loop param_grid:  20%|███████▍                             | 2/10 [00:02<00:09,  1.25s/it]
loop param_grid:  30%|███████████                          | 3/10 [00:03<00:08,  1.24s/it]
loop param_grid:  40%|██████████████▊                      | 4/10 [00:04<00:07,  1.25s/it]
loop param_grid:  50%|██████████████████▌                  | 5/10 [00:06<00:06,  1.27s/it]
loop param_grid:  60%|██████████████████████▏              | 6/10 [00:07<00:04,  1.25s/it]
loop param_grid:  70%|█████████████████████████▉           | 7/10 [00:09<00:04,  1.35s/it]
loop param_grid:  80%|█████████████████████████████▌       | 8/10 [00:10<00:02,  1.41s/it]
loop param_grid:  90%|█████████████████████████████████▎   | 9/10 [00:12<00:01,  1.58s/it]
loop param_grid: 100%|████████████████████████████████████| 10/10 [00:14<00:00,  1.59s/it]
loop lags_grid:  67%|██████████████████████████             | 2/3 [00:26<00:13, 13.56s/it]
loop param_grid:   0%|                                             | 0/10 [00:00<?, ?it/s]
loop param_grid:  10%|███▋                                 | 1/10 [00:01<00:11,  1.29s/it]
loop param_grid:  20%|███████▍                             | 2/10 [00:02<00:10,  1.34s/it]
loop param_grid:  30%|███████████                          | 3/10 [00:04<00:09,  1.37s/it]
loop param_grid:  40%|██████████████▊                      | 4/10 [00:05<00:08,  1.36s/it]
loop param_grid:  50%|██████████████████▌                  | 5/10 [00:06<00:06,  1.36s/it]
loop param_grid:  60%|██████████████████████▏              | 6/10 [00:07<00:05,  1.30s/it]
loop param_grid:  70%|█████████████████████████▉           | 7/10 [00:09<00:03,  1.30s/it]
loop param_grid:  80%|█████████████████████████████▌       | 8/10 [00:10<00:02,  1.31s/it]
loop param_grid:  90%|█████████████████████████████████▎   | 9/10 [00:11<00:01,  1.30s/it]
loop param_grid: 100%|████████████████████████████████████| 10/10 [00:13<00:00,  1.29s/it]
loop lags_grid: 100%|███████████████████████████████████████| 3/3 [00:40<00:00, 13.34s/it]
`Forecaster` refitted using the best-found lags and parameters, and the whole data set: 
  Lags: [ 1  2  3 23 24 25 47 48 49] 
  Parameters: {'ridge__alpha': 215.44346900318823}
  Backtesting metric: 257.84317250567267


In [23]:
# Resultados Grid Search
# ==============================================================================
resultados_grid
Out[23]:
lags params metric ridge__alpha
26 [1, 2, 3, 23, 24, 25, 47, 48, 49] {'ridge__alpha': 215.44346900318823} 257.843173 215.443469
25 [1, 2, 3, 23, 24, 25, 47, 48, 49] {'ridge__alpha': 27.825594022071257} 290.555205 27.825594
24 [1, 2, 3, 23, 24, 25, 47, 48, 49] {'ridge__alpha': 3.593813663804626} 306.631981 3.593814
23 [1, 2, 3, 23, 24, 25, 47, 48, 49] {'ridge__alpha': 0.46415888336127775} 309.393349 0.464159
22 [1, 2, 3, 23, 24, 25, 47, 48, 49] {'ridge__alpha': 0.05994842503189409} 309.776084 0.059948
21 [1, 2, 3, 23, 24, 25, 47, 48, 49] {'ridge__alpha': 0.007742636826811269} 309.825962 0.007743
20 [1, 2, 3, 23, 24, 25, 47, 48, 49] {'ridge__alpha': 0.001} 309.832410 0.001000
10 [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14... {'ridge__alpha': 0.001} 325.041129 0.001000
11 [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14... {'ridge__alpha': 0.007742636826811269} 325.043579 0.007743
12 [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14... {'ridge__alpha': 0.05994842503189409} 325.062536 0.059948
13 [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14... {'ridge__alpha': 0.46415888336127775} 325.208755 0.464159
14 [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14... {'ridge__alpha': 3.593813663804626} 326.307375 3.593814
15 [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14... {'ridge__alpha': 27.825594022071257} 333.395125 27.825594
27 [1, 2, 3, 23, 24, 25, 47, 48, 49] {'ridge__alpha': 1668.1005372000557} 356.547658 1668.100537
16 [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14... {'ridge__alpha': 215.44346900318823} 360.841496 215.443469
17 [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14... {'ridge__alpha': 1668.1005372000557} 396.342247 1668.100537
18 [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14... {'ridge__alpha': 12915.496650148827} 421.002019 12915.496650
28 [1, 2, 3, 23, 24, 25, 47, 48, 49] {'ridge__alpha': 12915.496650148827} 443.551888 12915.496650
19 [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14... {'ridge__alpha': 100000.0} 540.659659 100000.000000
29 [1, 2, 3, 23, 24, 25, 47, 48, 49] {'ridge__alpha': 100000.0} 545.502052 100000.000000
7 [1, 2, 3, 4, 5] {'ridge__alpha': 1668.1005372000557} 611.236033 1668.100537
0 [1, 2, 3, 4, 5] {'ridge__alpha': 0.001} 612.352191 0.001000
1 [1, 2, 3, 4, 5] {'ridge__alpha': 0.007742636826811269} 612.352531 0.007743
2 [1, 2, 3, 4, 5] {'ridge__alpha': 0.05994842503189409} 612.355162 0.059948
3 [1, 2, 3, 4, 5] {'ridge__alpha': 0.46415888336127775} 612.375445 0.464159
4 [1, 2, 3, 4, 5] {'ridge__alpha': 3.593813663804626} 612.528081 3.593814
5 [1, 2, 3, 4, 5] {'ridge__alpha': 27.825594022071257} 613.477722 27.825594
6 [1, 2, 3, 4, 5] {'ridge__alpha': 215.44346900318823} 615.109317 215.443469
8 [1, 2, 3, 4, 5] {'ridge__alpha': 12915.496650148827} 625.105850 12915.496650
9 [1, 2, 3, 4, 5] {'ridge__alpha': 100000.0} 681.830571 100000.000000

Los mejores resultados se obtienen si se utilizan los lags [1, 2, 3, 23, 24, 25, 47, 48, 49] y una configuración de Ridge {'alpha': 215.44}. Al indicar return_best = True en la función grid_search_forecaster(), al final del proceso, se reentrena automáticamente el objeto forecaster con la mejor configuración encontrada y el set de datos completo.

In [24]:
forecaster
Out[24]:
================= 
ForecasterAutoreg 
================= 
Regressor: Pipeline(steps=[('standardscaler', StandardScaler()),
                ('ridge', Ridge(alpha=215.44346900318823))]) 
Lags: [ 1  2  3 23 24 25 47 48 49] 
Window size: 49 
Included exogenous: 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: {'standardscaler__copy': True, 'standardscaler__with_mean': True, 'standardscaler__with_std': True, 'ridge__alpha': 215.44346900318823, 'ridge__copy_X': True, 'ridge__fit_intercept': True, 'ridge__max_iter': None, 'ridge__normalize': 'deprecated', 'ridge__positive': False, 'ridge__random_state': None, 'ridge__solver': 'auto', 'ridge__tol': 0.001} 
Creation date: 2022-04-03 22:09:58 
Last fit date: 2022-04-03 22:10:38 
Skforecast version: 0.4.3 

Backtest con el conjunto de test


Una vez identificado y entrenado el mejor modelo, se calcula su error al predecir los datos de test.

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

fig, ax = plt.subplots(figsize=(12, 3.5))
datos.loc[predicciones.index, 'Demand'].plot(linewidth=2, label='test', ax=ax)
predicciones.plot(linewidth=2, label='predicción', ax=ax)
ax.set_title('Predicción vs demanda real')
ax.legend();
In [26]:
# Error backtest
# ==============================================================================
print(f'Error backtest: {metrica}')
Error backtest: 251.93996461683977

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

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 $y$ con una determinada probabilidad. Por ejemplo, es de esperar que el intervalo de predicción (1, 99) contenga el verdadero valor de la predicción con un 98% de 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 el método utilizado en la librería skforecast para los modelos de tipo ForecasterAutoreg y ForecasterAutoregCustom.

In [27]:
# Backtest del conjunto de test con intervalos de predicción
# ==============================================================================
metric, predicciones = backtesting_forecaster(
                            forecaster          = forecaster,
                            y                   = datos.Demand,
                            initial_train_size  = len(datos.Demand[:fin_validacion]),
                            fixed_train_size    = False,
                            steps               = 24,
                            metric              = 'mean_absolute_error',
                            interval            = [10, 90],
                            n_boot              = 500,
                            in_sample_residuals = True,
                            verbose             = False
                       )

print('Métrica backtesting:', metric)
predicciones.head(5)
Métrica backtesting: 251.93996461683977
Out[27]:
pred lower_bound upper_bound
2014-12-01 00:00:00 5727.876561 5598.911947 5849.678211
2014-12-01 01:00:00 5802.876139 5599.270390 5974.888764
2014-12-01 02:00:00 5880.047528 5619.841591 6113.916977
2014-12-01 03:00:00 5953.541637 5657.437015 6240.035968
2014-12-01 04:00:00 6048.740321 5697.793533 6343.068708
In [28]:
# Gráfico
# ==============================================================================
fig, ax = plt.subplots(figsize=(12, 3.5))
datos.loc[predicciones.index, 'Demand'].plot(linewidth=2, label='test', ax=ax)
predicciones['pred'].plot(linewidth=2, label='predicción', ax=ax)
ax.set_title('Predicción vs demanda real')
ax.fill_between(
    predicciones.index,
    predicciones['lower_bound'],
    predicciones['upper_bound'],
    alpha = 0.3,
    color = 'red',
    label = 'Intervalo predicción' 
)
ax.legend();