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 Mayo 2021)

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.

En este documento se muestra un ejemplo de cómo utilizar métodos de forecasting para predecir la demanda eléctrica a nivel horario. Para ello, se hace uso de skforecast, una sencilla librería de Python que permite adaptar cualquier regresor de scikit-learn a problemas de forecasting.

Multi-Step Time Series Forecasting


Cuando se trabaja con series temporales, raramente se quiere predecir solo el siguiente elemento de la serie ($t_{+1}$), sino todo un intervalo futuro o un punto alejado en el tiempo ($t_{+n}$). A cada paso de predicción se le conoce como step. Existen varias estrategias que permiten generar este tipo de predicciones múltiples:

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, es necesario hacer predicciones recursivas en las que, cada nueva predicción, se basa en la predicción anterior. A este proceso se le conoce como recursive forecasting o recursive multi-step forecasting.

forecasting-python

La principal adaptación que se necesita hacer para aplicar modelos de scikit learn a problemas de recursive multi-step forecasting es transformar la serie temporal en un matriz en la que, cada valor, está asociado a la ventana temporal (lags) que le precede.

forecasting-python

Transformación de una serie temporal en una matriz de 5 lags y un vector con el valor de la serie que sigue a cada fila de la matriz.

Direct multi-step forecasting

Esta estrategia consiste en entrenar un modelo distinto para cada step. Por ejemplo, si se quieren predecir los siguientes 5 valores de una serie temporal, se entrenan 5 modelos distintos, uno para cada step. Como resultado, las predicciones son independientes unas de otras. Esta estrategia tiene un coste computacional más elevado ya que requiere entrenar múltiples modelos.

Multiple output forecasting

Determinados modelos son capaces de predecir de forma simultánea varios valores de una secuencia (one-shot). Un ejemplo de modelo con esta capacidad son las redes neuronales LSTM.

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 [2]:
# Tratamiento de datos
# ==============================================================================
import numpy as np
import pandas as pd

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

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

Además de las anteriores, se utiliza 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 siguiendo la estrategia de recursive multi-step forecasting o direct multi-step forecasting. Puede instalarse directamente desde Github con el comando:

!pip install git+https://github.com/JoaquinAmatRodrigo/skforecast@v0.1.8.1

Última versión (inestable):

!pip install git+https://github.com/JoaquinAmatRodrigo/skforecast#master

In [3]:
# Modelado y Forecasting
# ==============================================================================
#!pip install git+https://github.com/JoaquinAmatRodrigo/skforecast@v0.1.8.1

from sklearn.linear_model import LinearRegression
from sklearn.linear_model import Lasso
from sklearn.linear_model import Ridge
from sklearn.ensemble import RandomForestRegressor
from xgboost import XGBRegressor

from sklearn.metrics import mean_squared_error
from sklearn.metrics import mean_absolute_error
from sklearn.metrics import mean_absolute_percentage_error

from skforecast.ForecasterAutoreg import ForecasterAutoreg
from skforecast.ForecasterAutoregMultiOutput import ForecasterAutoregMultiOutput
from skforecast.model_selection import grid_search_forecaster
from skforecast.model_selection import time_series_spliter
from skforecast.model_selection import cv_forecaster
from skforecast.model_selection import backtesting_forecaster
from skforecast.model_selection import backtesting_forecaster_intervals

Datos


Los datos empleados en este documento se han obtenido del paquete de R tsibble. 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 [4]:
# 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 [5]:
# 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 [6]:
# 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[6]:
True
In [7]:
# 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 [10]:
# 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[10]:
Demand Temperature Holiday
Time
2011-12-31 14:00:00 4323.095350 21.225 True
2011-12-31 15:00:00 3963.264688 20.625 True
2011-12-31 16:00:00 3950.913495 20.325 True
2011-12-31 17:00:00 3627.860675 19.850 True
2011-12-31 18:00:00 3396.251676 19.025 True
... ... ... ...
2014-12-31 09:00:00 4069.625550 21.600 False
2014-12-31 10:00:00 3909.230704 20.300 False
2014-12-31 11:00:00 3900.600901 19.650 False
2014-12-31 12:00:00 3758.236494 18.100 False
2014-12-31 13:00:00 3785.650720 17.200 False

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 evaluar la capacidad predictiva del modelo, se separan los últimos 20 días como conjunto de test.

In [11]:
# Separación datos train-test
# ==============================================================================
dias_test = 20
datos = datos.loc['2012-01-01 00:00:00': '2014-12-30 23:00:00']
fecha_fin_train   = (datos.index[-1] - pd.Timedelta(dias_test, unit='D'))
fecha_inicio_test = fecha_fin_train + pd.Timedelta(1, unit='H')

datos_train = datos.loc[: fecha_fin_train, :]
datos_test  = datos.loc[fecha_inicio_test:, :]

print(f"Fechas train: {datos_train.index.min()} --- {datos_train.index.max()}")
print(f"Fechas test : {datos_test.index.min()} --- {datos_test.index.max()}")
Fechas train: 2012-01-01 00:00:00 --- 2014-12-10 23:00:00
Fechas test : 2014-12-11 00:00:00 --- 2014-12-30 23:00:00

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 [12]:
# Gráfico serie temporal
# ==============================================================================
fig, ax = plt.subplots(figsize=(12, 4))
datos_train.Demand.plot(ax=ax, label='train', linewidth=0.5)
datos_test.Demand.plot(ax=ax, label='test', linewidth=0.5)
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 [13]:
# 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 mensual y semanal

In [14]:
# Gráfico boxplot para estacionalidad mensual
# ==============================================================================
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('');
In [15]:
# 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 confirma que hay una estacionalidad mensual, 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.

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 = Ridge(normalize=True),
                lags      = 24
             )

forecaster.fit(y=datos_train.Demand)
forecaster
Out[18]:
=======================ForecasterAutoreg=======================
Regressor: Ridge(normalize=True)
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]
Exogenous variable: False
Parameters: {'alpha': 1.0, 'copy_X': True, 'fit_intercept': True, 'max_iter': None, 'normalize': True, 'random_state': None, 'solver': 'auto', 'tol': 0.001}

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-12-10 23:00 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_train.Demand),
                            steps      = 24,
                            metric     = 'neg_mean_absolute_error',
                            verbose    = True
                        )

# Se añade el índice temporal a las predicciones
predicciones = pd.Series(data=predicciones, index=datos_test.index)
Number of observations used for training: 25800
Number of observations used for testing: 480
    Number of folds: 20
    Number of steps per fold: 24
In [20]:
# Gráfico
# ==============================================================================
fig, ax = plt.subplots(figsize=(9, 4))
datos_test.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: [374.39861662]

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 con validación por backtesting. Este proceso consiste en entrenar un modelo con cada combinación de hiperparámetros y lags, y evaluar su capacidad predictiva mediante backtesting. Dentro del proceso de Grid Search, se entrenan los modelos hasta final del 2013 ('2013-12-31 23:00:00').

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

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

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

resultados_grid = grid_search_forecaster(
                        forecaster  = forecaster,
                        y           = datos_train.Demand,
                        param_grid  = param_grid,
                        lags_grid   = lags_grid,
                        steps       = 24,
                        metric      = 'neg_mean_absolute_error',
                        method      = 'backtesting',
                        initial_train_size = len(datos_train.Demand[:'2013-12-31 23:00:00']),
                        return_best = True,
                        verbose     = False
                  )
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:00<00:05,  1.65it/s]
loop param_grid:  20%|██        | 2/10 [00:01<00:04,  1.66it/s]
loop param_grid:  30%|███       | 3/10 [00:01<00:04,  1.65it/s]
loop param_grid:  40%|████      | 4/10 [00:02<00:03,  1.66it/s]
loop param_grid:  50%|█████     | 5/10 [00:03<00:03,  1.63it/s]
loop param_grid:  60%|██████    | 6/10 [00:03<00:02,  1.67it/s]
loop param_grid:  70%|███████   | 7/10 [00:04<00:01,  1.63it/s]
loop param_grid:  80%|████████  | 8/10 [00:04<00:01,  1.64it/s]
loop param_grid:  90%|█████████ | 9/10 [00:05<00:00,  1.63it/s]
loop param_grid: 100%|██████████| 10/10 [00:06<00:00,  1.64it/s]
loop lags_grid:  33%|███▎      | 1/3 [00:06<00:12,  6.09s/it]   
loop param_grid:   0%|          | 0/10 [00:00<?, ?it/s]
loop param_grid:  10%|█         | 1/10 [00:00<00:05,  1.71it/s]
loop param_grid:  20%|██        | 2/10 [00:01<00:04,  1.67it/s]
loop param_grid:  30%|███       | 3/10 [00:01<00:04,  1.72it/s]
loop param_grid:  40%|████      | 4/10 [00:02<00:03,  1.68it/s]
loop param_grid:  50%|█████     | 5/10 [00:02<00:03,  1.66it/s]
loop param_grid:  60%|██████    | 6/10 [00:03<00:02,  1.64it/s]
loop param_grid:  70%|███████   | 7/10 [00:04<00:01,  1.68it/s]
loop param_grid:  80%|████████  | 8/10 [00:04<00:01,  1.66it/s]
loop param_grid:  90%|█████████ | 9/10 [00:05<00:00,  1.64it/s]
loop param_grid: 100%|██████████| 10/10 [00:06<00:00,  1.66it/s]
loop lags_grid:  67%|██████▋   | 2/3 [00:12<00:06,  6.04s/it]   
loop param_grid:   0%|          | 0/10 [00:00<?, ?it/s]
loop param_grid:  10%|█         | 1/10 [00:00<00:05,  1.76it/s]
loop param_grid:  20%|██        | 2/10 [00:01<00:04,  1.78it/s]
loop param_grid:  30%|███       | 3/10 [00:01<00:04,  1.72it/s]
loop param_grid:  40%|████      | 4/10 [00:02<00:03,  1.73it/s]
loop param_grid:  50%|█████     | 5/10 [00:02<00:02,  1.68it/s]
loop param_grid:  60%|██████    | 6/10 [00:03<00:02,  1.67it/s]
loop param_grid:  70%|███████   | 7/10 [00:04<00:01,  1.65it/s]
loop param_grid:  80%|████████  | 8/10 [00:04<00:01,  1.65it/s]
loop param_grid:  90%|█████████ | 9/10 [00:05<00:00,  1.66it/s]
loop param_grid: 100%|██████████| 10/10 [00:05<00:00,  1.63it/s]
loop lags_grid: 100%|██████████| 3/3 [00:18<00:00,  6.03s/it]   
2021-05-16 19:30:30,337 root       INFO  Refitting `forecaster` using the best found parameters and the whole data set: 
lags: [ 1  2  3 23 24 25 47 48 49] 
params: {'alpha': 0.004641588833612777}

In [23]:
# Resultados Grid Search
# ==============================================================================
resultados_grid
Out[23]:
lags params metric
21 [1, 2, 3, 23, 24, 25, 47, 48, 49] {'alpha': 0.004641588833612777} 268.346474
22 [1, 2, 3, 23, 24, 25, 47, 48, 49] {'alpha': 0.021544346900318832} 269.490407
20 [1, 2, 3, 23, 24, 25, 47, 48, 49] {'alpha': 0.001} 295.394793
10 [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14... {'alpha': 0.001} 328.736662
11 [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14... {'alpha': 0.004641588833612777} 342.381520
23 [1, 2, 3, 23, 24, 25, 47, 48, 49] {'alpha': 0.1} 357.465267
12 [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14... {'alpha': 0.021544346900318832} 369.988994
13 [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14... {'alpha': 0.1} 394.608650
14 [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14... {'alpha': 0.46415888336127775} 408.881232
24 [1, 2, 3, 23, 24, 25, 47, 48, 49] {'alpha': 0.46415888336127775} 426.377362
15 [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14... {'alpha': 2.154434690031882} 464.191125
25 [1, 2, 3, 23, 24, 25, 47, 48, 49] {'alpha': 2.154434690031882} 482.242892
26 [1, 2, 3, 23, 24, 25, 47, 48, 49] {'alpha': 10.0} 585.253577
16 [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14... {'alpha': 10.0} 585.399425
3 [1, 2, 3, 4, 5] {'alpha': 0.1} 608.644246
0 [1, 2, 3, 4, 5] {'alpha': 0.001} 611.255382
4 [1, 2, 3, 4, 5] {'alpha': 0.46415888336127775} 612.254957
2 [1, 2, 3, 4, 5] {'alpha': 0.021544346900318832} 612.665889
1 [1, 2, 3, 4, 5] {'alpha': 0.004641588833612777} 612.749912
5 [1, 2, 3, 4, 5] {'alpha': 2.154434690031882} 655.641466
27 [1, 2, 3, 23, 24, 25, 47, 48, 49] {'alpha': 46.41588833612773} 667.312672
17 [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14... {'alpha': 46.41588833612773} 670.024090
6 [1, 2, 3, 4, 5] {'alpha': 10.0} 688.811445
28 [1, 2, 3, 23, 24, 25, 47, 48, 49] {'alpha': 215.44346900318823} 695.267462
18 [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14... {'alpha': 215.44346900318823} 696.141975
7 [1, 2, 3, 4, 5] {'alpha': 46.41588833612773} 700.347947
29 [1, 2, 3, 23, 24, 25, 47, 48, 49] {'alpha': 1000.0} 702.036319
19 [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14... {'alpha': 1000.0} 702.277561
8 [1, 2, 3, 4, 5] {'alpha': 215.44346900318823} 703.186734
9 [1, 2, 3, 4, 5] {'alpha': 1000.0} 703.817946

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': 0.00464}. 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: Ridge(alpha=0.004641588833612777, normalize=True)
Lags: [ 1  2  3 23 24 25 47 48 49]
Exogenous variable: False
Parameters: {'alpha': 0.004641588833612777, 'copy_X': True, 'fit_intercept': True, 'max_iter': None, 'normalize': True, 'random_state': None, 'solver': 'auto', 'tol': 0.001}
In [25]:
# Backtest modelo final
# ==============================================================================
metrica, predicciones = backtesting_forecaster(
                            forecaster = forecaster,
                            y          = datos.Demand,
                            initial_train_size = len(datos_train.Demand),
                            steps      = 24,
                            metric     = 'neg_mean_absolute_error',
                            verbose    = False
                        )

# Se añade el índice a las predicciones
predicciones = pd.Series(data=predicciones, index=datos_test.index)

fig, ax = plt.subplots(figsize=(9, 4))
datos_test.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: [279.52569303]

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

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.

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 ForecasterCustom.

In [27]:
# Backtest con intervalos de predicción
# ==============================================================================
metric, predictions = backtesting_forecaster_intervals(
                            forecaster = forecaster,
                            y          = datos.Demand,
                            initial_train_size = len(datos_train.Demand),
                            steps      = 24,
                            metric     = 'neg_mean_squared_error',
                            interval            = [10, 90],
                            n_boot              = 500,
                            in_sample_residuals = True,
                            verbose             = True
                       )

print('Métrica backtesting:', metric)

# Se añade índice datetime
predictions = pd.DataFrame(data=predictions, index=datos_test.index)
Number of observations used for training: 25800
Number of observations used for testing: 480
    Number of folds: 20
    Number of steps per fold: 24
Métrica backtesting: [132021.90448464]
In [28]:
# Gráfico
# ==============================================================================
fig, ax = plt.subplots(figsize=(9, 4))
datos_test.loc[predicciones.index, 'Demand'].plot(linewidth=2, label='test', ax=ax)
predictions.iloc[:, 0].plot(linewidth=2, label='predicción', ax=ax)
ax.set_title('Predicción vs demanda real')
ax.fill_between(
    predictions.index,
    predictions.iloc[:, 1],
    predictions.iloc[:, 2],
    alpha = 0.2,
    color = 'red',
    label = 'Intervalo predicción' 
)
ax.legend();
In [29]:
# Cobertura del intervalo predicho
# ==============================================================================
dentro_intervalo = np.where(
                     (datos_test.Demand >= predictions.iloc[:, 1]) & \
                     (datos_test.Demand <= predictions.iloc[:, 2]),
                     True,
                     False
                   )

cobertura = dentro_intervalo.mean()
print(f"Cobertura del intervalo predicho: {100 * cobertura}")
Cobertura del intervalo predicho: 79.16666666666666

El intervalo predicho tiene una cobertura inferior a la que cabría esperar (80%). Esto puede deberse al error marcadamente elevado que comete el modelo para los días 21, 24 y 25. Estos días están dentro del periodo vacacional de navidad, que suele caracterizarse por un comportamiento de consumo distinto al resto del mes.

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.

El proceso de backtesting adaptado a este escenario es:

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

  2. Se almacenan solo las predicciones del día siguiente, es decir, de la posición 12 en adelante.

  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:00h de cada día, el modelo tiene acceso a los valores reales de demanda registrados hasta ese momento.

Este proceso puede realizarse fácilmente con el método predict() de un objeto ForecasterAutoreg. Si no se le indica nada, la predicción se inicia después del último valor de entrenamiento, pero, si se le especifica el argumento last_window, utiliza estos valores como punto de partida.

In [30]:
def backtest_predict_next_24h(forecaster, y, hour_init_prediction, exog=None,
                              verbose=False):
    
    '''
    Backtest ForecasterAutoreg object when predicting 24 hours of day D+1
    statring at specific hour of day D.
    
    Parameters
    ----------
    forecaster : ForecasterAutoreg 
        ForecasterAutoreg object already trained.
        
    y : pd.Series with datetime index sorted
        Test time series values. 
        
    exog : pd.Series or pd.Dataframe with datetime index sorted
        Test values of exogen variable. 
    
    hour_init_prediction: int 
        Hour of day D to start predictions of day D+1.


    Returns 
    -------
    predictions: pd.Series
        Value of predictions.

    '''
    
    y = y.sort_index()
    if exog is not None:
        exog = exog.sort_index()
        
    dummy_steps = 24 - (hour_init_prediction + 1)
    steps = dummy_steps + 24
    
    
    
    # First position of `hour_init_prediction` in the series where there is enough
    # previous window to calculate lags.
    for datetime in y.index[y.index.hour == hour_init_prediction]:
        if len(y[:datetime]) >= len(forecaster.last_window):
            datetime_init_backtest = datetime
            print(f"Backtesting starts at day: {datetime_init_backtest}")
            break
    
    days_backtest = np.unique(y[datetime_init_backtest:].index.date)
    days_backtest = pd.to_datetime(days_backtest)
    days_backtest = days_backtest[1:]
    print(f"Days predicted in the backtesting: {days_backtest.strftime('%Y-%m-%d').values}")
    print('')
    backtest_predictions = []
    
    for i, day in enumerate(days_backtest):        
        # Start and end of the last window used to create the lags
        end_window = (day - pd.Timedelta(1, unit='day')).replace(hour=hour_init_prediction)
        start_window = end_window - pd.Timedelta(forecaster.max_lag, unit='hour')
        last_window = y.loc[start_window:end_window]
               
        if exog is None:
            if verbose:
                print(f"Forecasting day {day.strftime('%Y-%m-%d')}")
                print(f"Using window from {start_window} to {end_window}")
                
            pred = forecaster.predict(steps=steps, last_window=last_window)
            
        else:
            start_exog_window = end_window + pd.Timedelta(1, unit='hour')
            end_exog_window   = end_window + pd.Timedelta(steps, unit='hour')
            exog_window = exog.loc[start_exog_window:end_exog_window]
            exog_window = exog_window.to_numpy()
            
            if verbose:
                print(f"Forecasting day {day.strftime('%Y-%m-%d')}")
                print(f"    Using window from {start_window} to {end_window}")
                print(f"    Using exogen variable from {start_exog_window} to {end_exog_window}")
            
            pred = forecaster.predict(steps=steps, last_window=last_window, exog=exog_window)

        # Only store predictions of day D+1
        pred = pred[dummy_steps:]
        backtest_predictions.append(pred)
    
    backtest_predictions = np.concatenate(backtest_predictions)
    # Add datetime index
    backtest_predictions = pd.Series(
                             data  = backtest_predictions,
                             index = pd.date_range(
                                        start = days_backtest[0],
                                        end   = days_backtest[-1].replace(hour=23),
                                        freq  = 'h'
                                    )
                           )
    
    return backtest_predictions
In [31]:
# Backtest
# ==============================================================================
predicciones = backtest_predict_next_24h(
                    forecaster = forecaster,
                    y          = datos_test.Demand,
                    hour_init_prediction = 11,
                    verbose    = False
                )
Backtesting starts at day: 2014-12-13 11:00:00
Days predicted in the backtesting: ['2014-12-14' '2014-12-15' '2014-12-16' '2014-12-17' '2014-12-18'
 '2014-12-19' '2014-12-20' '2014-12-21' '2014-12-22' '2014-12-23'
 '2014-12-24' '2014-12-25' '2014-12-26' '2014-12-27' '2014-12-28'
 '2014-12-29' '2014-12-30']

In [32]:
# Gráfico
# ==============================================================================
fig, ax = plt.subplots(figsize=(9, 4))
# datos_train.plot(ax=ax, label='train')
datos_test.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();