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.
C:\Users\WINDOWS 11\Documents\.virtualenvs\r-tensorflow\lib\site-packages\tqdm\auto.py:21: TqdmWarning:
IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
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 Rtsibbledata. 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).
# 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'.
Uno de los primeros análisis que hay que realizar al trabajar con series temporales es verificar si la serie está completa.
# Verificar que un índice temporal está completo# ==============================================================================fecha_inicio = datos.index.min()fecha_fin = datos.index.max()date_range_completo = pd.date_range(start=fecha_inicio, end=fecha_fin, freq=datos.index.freq)print(f"Índice completo: {(datos.index == date_range_completo).all()}")print(f"Filas con valores ausentes: {datos.isnull().any(axis=1).mean()}")# Completar huecos en un índice temporal# ==============================================================================# datos.asfreq(freq='30min', fill_value=np.nan)
Índice completo: True
Filas con valores ausentes: 0.0
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.
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.
# Agregado en intervalos de 1H# ==============================================================================# Se elimina la columna Date para que no genere error al agregar.datos = datos.drop(columns='Date')datos = ( datos .resample(rule="h", closed="left", label="right") .agg({"Demand": "mean","Temperature": "mean","Holiday": "mean", }))datos
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.
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
# 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=800, 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.
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# ==============================================================================fig, axs = plt.subplots(2, 2, figsize=(8, 5), sharex=False, sharey=True)axs = axs.ravel()# Distribusión de demanda por mesdatos['month'] = datos.index.monthdatos.boxplot(column='Demand', by='month', ax=axs[0], flierprops={'markersize': 3, 'alpha': 0.3})datos.groupby('month')['Demand'].median().plot(style='o-', linewidth=0.8, ax=axs[0])axs[0].set_ylabel('Demand')axs[0].set_title('Distribución de demanda por mes', fontsize=9)# Distribusión de demanda por día de la semanadatos['week_day'] = datos.index.day_of_week +1datos.boxplot(column='Demand', by='week_day', ax=axs[1], flierprops={'markersize': 3, 'alpha': 0.3})datos.groupby('week_day')['Demand'].median().plot(style='o-', linewidth=0.8, ax=axs[1])axs[1].set_ylabel('Demand')axs[1].set_title('Distribución de demanda por día de la semana', fontsize=9)# Distribusión de demanda por hora del díadatos['hour_day'] = datos.index.hour +1datos.boxplot(column='Demand', by='hour_day', ax=axs[2], flierprops={'markersize': 3, 'alpha': 0.3})datos.groupby('hour_day')['Demand'].median().plot(style='o-', linewidth=0.8, ax=axs[2])axs[2].set_ylabel('Demand')axs[2].set_title('Distribución de demanda por hora del día', fontsize=9)# Distribusión de demanda por día de la semana y hora del díamean_day_hour = datos.groupby(["week_day", "hour_day"])["Demand"].mean()mean_day_hour.plot(ax=axs[3])axs[3].set( title ="Promedio de demanda", xticks = [i *24for i inrange(7)], xticklabels = ["Mon", "Tue", "Wed", "Thu", "Fri", "Sat", "Sun"], xlabel ="Día y hora", ylabel ="Demand")axs[3].title.set_size(10)fig.suptitle("Gráficos de estacionalidad", fontsize=12)fig.tight_layout()plt.show()
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. Se aprecia una estacionalidad semanal, con valores de demanda inferiores durante el fin de semana. También existe una estacionalidad diaria, ya que la demanda se reduce entre las 16:00 y las 21:00 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.
# Top 10 lags con mayor autocorrelación parcial absoluta# ==============================================================================calculate_lag_autocorrelation( data = datos['Demand'], n_lags =60, sort_by ="partial_autocorrelation_abs").head(10)
lag
partial_autocorrelation_abs
partial_autocorrelation
autocorrelation_abs
autocorrelation
0
1
0.949526
0.949526
0.949490
0.949490
1
25
0.761400
-0.761400
0.731620
0.731620
2
2
0.657961
-0.657961
0.836792
0.836792
3
26
0.634868
0.634868
0.622458
0.622458
4
24
0.308085
-0.308085
0.785622
0.785622
5
19
0.291113
0.291113
0.302308
0.302308
6
27
0.281261
-0.281261
0.488351
0.488351
7
21
0.269274
0.269274
0.537173
0.537173
8
20
0.201192
0.201192
0.414709
0.414709
9
9
0.184975
0.184975
0.037677
0.037677
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 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.
# 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
========================
ForecasterEquivalentDate
========================
Offset: <DateOffset: days=1>
Number of offsets: 1
Aggregation function: mean
Window size: 24
Series name: Demand
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: 2025-05-29 11:52:54
Last fit date: 2025-05-29 11:52:54
Skforecast version: 0.16.0
Python version: 3.10.11
Forecaster id: None
El error del modelo base se utiliza como referencia para evaluar si merece la pena aplicar los modelos más complejos.
Modelo autoregresivo recursivo
Se entrena un modelo autorregresivo recursivo (ForecasterRecursive) con un gradient boostingLGBMRegressor como regresor para predecir la demanda de energía de las próximas 24 horas.
Se utilizan como predictores los valores de consumo de las últimas 24 horas (24 lags) así como la media movil de los últimos 3 días. Los hiperparámetros del regresor se dejan en sus valores por defecto.
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).
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.
# Backtesting# ==============================================================================metrica, predicciones = backtesting_forecaster( forecaster = forecaster, y = datos['Demand'], cv = cv, metric ='mean_absolute_error', verbose =False, # False para no mostrar info )
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.
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.
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.
# Variables basadas en el calendario# ==============================================================================features_to_extract = ['month','week','day_of_week','hour']calendar_transformer = DatetimeFeatures( variables ='index', features_to_extract = features_to_extract, drop_original =True,)variables_calendario = calendar_transformer.fit_transform(datos)[features_to_extract]# Variables basadas en la luz solar# ==============================================================================location = LocationInfo( latitude =-37.8, longitude =144.95, timezone ='Australia/Melbourne')sunrise_hour = [ sun(location.observer, date=date, tzinfo=location.timezone)['sunrise']for date in datos.index]sunset_hour = [ sun(location.observer, date=date, tzinfo=location.timezone)['sunset']for date in datos.index]sunrise_hour = pd.Series(sunrise_hour, index=datos.index).dt.round("h").dt.hoursunset_hour = pd.Series(sunset_hour, index=datos.index).dt.round("h").dt.hourvariables_solares = pd.DataFrame({'sunrise_hour': sunrise_hour,'sunset_hour': sunset_hour })variables_solares['daylight_hours'] = ( variables_solares['sunset_hour'] - variables_solares['sunrise_hour'])variables_solares["is_daylight"] = np.where( (datos.index.hour >= variables_solares["sunrise_hour"])& (datos.index.hour < variables_solares["sunset_hour"]),1,0,)# Variables basadas en festivos# ==============================================================================variables_festivos = datos[['Holiday']].astype(int)variables_festivos['holiday_previous_day'] = variables_festivos['Holiday'].shift(24)variables_festivos['holiday_next_day'] = variables_festivos['Holiday'].shift(-24)# Variables basadas en temperatura# ==============================================================================wf_transformer = WindowFeatures( variables = ["Temperature"], window = ["1D", "7D"], functions = ["mean", "max", "min"], freq ="h",)variables_temp = wf_transformer.fit_transform(datos[['Temperature']])# Unión de variables exógenas# ==============================================================================assertall(variables_calendario.index == variables_solares.index)assertall(variables_calendario.index == variables_festivos.index)assertall(variables_calendario.index == variables_temp.index)variables_exogenas = pd.concat([ variables_calendario, variables_solares, variables_temp, variables_festivos ], axis=1)# Debido a la creación de medias móviles, hay valores faltantes al principio# de la serie. Y debido a holiday_next_day hay valores faltantes al final.variables_exogenas = variables_exogenas.iloc[7*24:, :]variables_exogenas = variables_exogenas.iloc[:-24, :]variables_exogenas.head(3)
month
week
day_of_week
hour
sunrise_hour
sunset_hour
daylight_hours
is_daylight
Temperature
Temperature_window_1D_mean
Temperature_window_1D_max
Temperature_window_1D_min
Temperature_window_7D_mean
Temperature_window_7D_max
Temperature_window_7D_min
Holiday
holiday_previous_day
holiday_next_day
Time
2012-01-08 00:00:00
1
1
6
0
6
21
15
0
20.575
24.296875
29.0
19.875
23.254018
39.525
14.35
0
0.0
0.0
2012-01-08 01:00:00
1
1
6
1
6
21
15
0
22.500
24.098958
29.0
19.875
23.215774
39.525
14.35
0
0.0
0.0
2012-01-08 02:00:00
1
1
6
2
6
21
15
0
25.250
23.923958
29.0
19.875
23.173214
39.525
14.35
0
0.0
0.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.
# Codificación cíclica de las variables de calendario y luz solar# ==============================================================================features_to_encode = ["month","week","day_of_week","hour","sunrise_hour","sunset_hour",]max_values = {"month": 12,"week": 52,"day_of_week": 6,"hour": 24,"sunrise_hour": 24,"sunset_hour": 24,}cyclical_encoder = CyclicalFeatures( variables = features_to_encode, max_values = max_values, drop_original =False)variables_exogenas = cyclical_encoder.fit_transform(variables_exogenas)variables_exogenas.head(3)
month
week
day_of_week
hour
sunrise_hour
sunset_hour
daylight_hours
is_daylight
Temperature
Temperature_window_1D_mean
...
week_sin
week_cos
day_of_week_sin
day_of_week_cos
hour_sin
hour_cos
sunrise_hour_sin
sunrise_hour_cos
sunset_hour_sin
sunset_hour_cos
Time
2012-01-08 00:00:00
1
1
6
0
6
21
15
0
20.575
24.296875
...
0.120537
0.992709
-2.449294e-16
1.0
0.000000
1.000000
1.0
6.123234e-17
-0.707107
0.707107
2012-01-08 01:00:00
1
1
6
1
6
21
15
0
22.500
24.098958
...
0.120537
0.992709
-2.449294e-16
1.0
0.258819
0.965926
1.0
6.123234e-17
-0.707107
0.707107
2012-01-08 02:00:00
1
1
6
2
6
21
15
0
25.250
23.923958
...
0.120537
0.992709
-2.449294e-16
1.0
0.500000
0.866025
1.0
6.123234e-17
-0.707107
0.707107
3 rows × 30 columns
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.
# Selección de variables exógenas incluidas en el modelo# ==============================================================================exog_features = []# Columnas que terminan con _seno o _coseno son seleccionadasexog_features.extend(variables_exogenas.filter(regex='_sin$|_cos$').columns.tolist())# Columnas que empiezan con temp_ son seleccionadasexog_features.extend(variables_exogenas.filter(regex='^Temperature_.*').columns.tolist())# Columnas que empiezan con holiday_ son seleccionadasexog_features.extend(variables_exogenas.filter(regex='^Holiday_.*').columns.tolist())# Incluir temperatura y festivosexog_features.extend(['Temperature', 'Holiday'])
# 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 ='inner'# Usar solo las filas que coinciden en ambos DataFrames )datos = datos.astype('float32')# Split data into train-val-testdatos_train = datos.loc[: fin_train, :].copy()datos_val = datos.loc[fin_train:fin_validacion, :].copy()datos_test = datos.loc[fin_validacion:, :].copy()
Se vuelve a evaluar el modelo pero, esta vez, las variables exógenas también se incluyen como predictores.
La inclusión de variables exógenas como predictores mejora la capacidad predictiva del modelo.
Optimización de hiperparámetros (tuning)
El ForecasterRecursive 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.
# Búsqueda bayesiana de hiperparámetros# ==============================================================================forecaster = ForecasterRecursive( regressor = LGBMRegressor(random_state=15926, verbose=-1), lags =24, # Este valor se modifica durante la búsqueda window_features = window_features )# Lags utilizados como predictoreslags_grid = [24, [1, 2, 3, 23, 24, 25, 47, 48, 49]]# Espacio de búsqueda de hiperparámetrosdef search_space(trial): search_space = {'n_estimators' : trial.suggest_int('n_estimators', 300, 1000, step=100),'max_depth' : trial.suggest_int('max_depth', 3, 10),'learning_rate': trial.suggest_float('learning_rate', 0.01, 0.5),'reg_alpha' : trial.suggest_float('reg_alpha', 0, 1),'reg_lambda' : trial.suggest_float('reg_lambda', 0, 1),'lags' : trial.suggest_categorical('lags', lags_grid) } return search_space# Partición de entrenamiento y validacióncv_search = TimeSeriesFold( steps =24, initial_train_size =len(datos[:fin_train]), refit =False, )resultados_busqueda, frozen_trial = bayesian_search_forecaster( forecaster = forecaster, y = datos.loc[:fin_validacion, 'Demand'], exog = datos.loc[:fin_validacion, exog_features], cv = cv_search, metric ='mean_absolute_error', search_space = search_space, n_trials =10, # Aumentar para una búsqueda más exhaustiva return_best =True )
0%| | 0/10 [00:00<?, ?it/s]Best trial: 0. Best value: 129.403: 0%| | 0/10 [00:05<?, ?it/s]Best trial: 0. Best value: 129.403: 10%|█ | 1/10 [00:05<00:53, 5.97s/it]Best trial: 0. Best value: 129.403: 10%|█ | 1/10 [00:12<00:53, 5.97s/it]Best trial: 0. Best value: 129.403: 20%|██ | 2/10 [00:12<00:48, 6.07s/it]Best trial: 0. Best value: 129.403: 20%|██ | 2/10 [00:17<00:48, 6.07s/it]Best trial: 0. Best value: 129.403: 30%|███ | 3/10 [00:17<00:41, 5.93s/it]Best trial: 0. Best value: 129.403: 30%|███ | 3/10 [00:24<00:41, 5.93s/it]Best trial: 0. Best value: 129.403: 40%|████ | 4/10 [00:24<00:36, 6.12s/it]Best trial: 0. Best value: 129.403: 40%|████ | 4/10 [00:30<00:36, 6.12s/it]Best trial: 0. Best value: 129.403: 50%|█████ | 5/10 [00:30<00:30, 6.01s/it]Best trial: 0. Best value: 129.403: 50%|█████ | 5/10 [00:35<00:30, 6.01s/it]Best trial: 0. Best value: 129.403: 60%|██████ | 6/10 [00:35<00:23, 5.93s/it]Best trial: 0. Best value: 129.403: 60%|██████ | 6/10 [00:41<00:23, 5.93s/it]Best trial: 0. Best value: 129.403: 70%|███████ | 7/10 [00:41<00:17, 5.90s/it]Best trial: 0. Best value: 129.403: 70%|███████ | 7/10 [00:47<00:17, 5.90s/it]Best trial: 0. Best value: 129.403: 80%|████████ | 8/10 [00:47<00:11, 5.83s/it]Best trial: 0. Best value: 129.403: 80%|████████ | 8/10 [00:53<00:11, 5.83s/it]Best trial: 0. Best value: 129.403: 90%|█████████ | 9/10 [00:53<00:05, 5.77s/it]Best trial: 9. Best value: 128.29: 90%|█████████ | 9/10 [00:59<00:05, 5.77s/it] Best trial: 9. Best value: 128.29: 100%|██████████| 10/10 [00:59<00:00, 5.86s/it]Best trial: 9. Best value: 128.29: 100%|██████████| 10/10 [00:59<00:00, 5.91s/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: {'n_estimators': 800, 'max_depth': 9, 'learning_rate': 0.050765544282894995, 'reg_alpha': 0.7636828414433382, 'reg_lambda': 0.243666374536874}
Backtesting metric: 128.29025410611806
# Resultados de la búsqueda# ==============================================================================best_params = resultados_busqueda.at[0, 'params']best_params = best_params | {'random_state': 15926, 'verbose': -1}best_lags = resultados_busqueda.at[0, 'lags']resultados_busqueda.head(3)
lags
params
mean_absolute_error
n_estimators
max_depth
learning_rate
reg_alpha
reg_lambda
0
[1, 2, 3, 23, 24, 25, 47, 48, 49]
{'n_estimators': 800, 'max_depth': 9, 'learnin...
128.290254
800.0
9.0
0.050766
0.763683
0.243666
1
[1, 2, 3, 23, 24, 25, 47, 48, 49]
{'n_estimators': 800, 'max_depth': 5, 'learnin...
129.402763
800.0
5.0
0.121157
0.551315
0.719469
2
[1, 2, 3, 23, 24, 25, 47, 48, 49]
{'n_estimators': 600, 'max_depth': 8, 'learnin...
130.719469
600.0
8.0
0.099421
0.175452
0.531551
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.
# Mejor modelo# ==============================================================================forecaster
Una vez identificada la mejor combinación de hiperparámetros utilizando los datos de validación, se evalúa de nuevo la capacidad predictiva del modelo cuando se aplica al conjunto de test.
Tras la optimización de lags e hiperparámetros, se ha conseguido reducir aún más el error de predicción.
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 de forma sencilla con la función select_features(). Dos de los métodos más populares son Recursive Feature Elimination y Sequential Feature Selection.
Recursive feature elimination (RFECV)
-------------------------------------
Total number of records available: 17304
Total number of records used for feature selection: 8652
Number of features available: 96
Lags (n=9)
Window features (n=1)
Exog (n=86)
Number of features selected: 55
Lags (n=9) : [1, 2, 3, 23, 24, 25, 47, 48, 49]
Window features (n=1) : ['roll_mean_72']
Exog (n=45) : ['week_cos', 'day_of_week_sin', 'hour_sin', 'hour_cos', 'sunset_hour_cos', 'poly_month_sin__week_sin', 'poly_month_sin__day_of_week_cos', 'poly_month_sin__hour_sin', 'poly_month_sin__hour_cos', 'poly_month_cos__day_of_week_cos', 'poly_month_cos__hour_cos', 'poly_week_sin__day_of_week_sin', 'poly_week_sin__day_of_week_cos', 'poly_week_sin__hour_sin', 'poly_week_sin__hour_cos', 'poly_week_sin__sunrise_hour_sin', 'poly_week_cos__day_of_week_sin', 'poly_week_cos__day_of_week_cos', 'poly_week_cos__hour_sin', 'poly_week_cos__hour_cos', 'poly_day_of_week_sin__day_of_week_cos', 'poly_day_of_week_sin__hour_sin', 'poly_day_of_week_sin__hour_cos', 'poly_day_of_week_sin__sunrise_hour_sin', 'poly_day_of_week_cos__hour_sin', 'poly_day_of_week_cos__hour_cos', 'poly_day_of_week_cos__sunrise_hour_cos', 'poly_day_of_week_cos__sunset_hour_sin', 'poly_day_of_week_cos__sunset_hour_cos', 'poly_hour_sin__hour_cos', 'poly_hour_sin__sunrise_hour_cos', 'poly_hour_sin__sunset_hour_sin', 'poly_hour_sin__sunset_hour_cos', 'poly_hour_cos__sunrise_hour_sin', 'poly_hour_cos__sunrise_hour_cos', 'poly_hour_cos__sunset_hour_sin', 'poly_hour_cos__sunset_hour_cos', 'Temperature_window_1D_mean', 'Temperature_window_1D_max', 'Temperature_window_1D_min', 'Temperature_window_7D_mean', 'Temperature_window_7D_max', 'Temperature_window_7D_min', 'Temperature', 'Holiday']
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.
El forecater se entrena y reevalúa utilizando solo el conjunto de predictores seleccionados.
# Crear forecaster# ==============================================================================forecaster = ForecasterRecursive( regressor = LGBMRegressor(**best_params), lags = lags_select, window_features = window_features )# Backtesting con los predictores seleccionados y los datos de test# ==============================================================================metrica, predicciones = backtesting_forecaster( forecaster = forecaster, y = datos['Demand'], exog = datos[exog_select], cv = cv, metric ='mean_absolute_error' )metrica
Se ha conseguido reducir el número de predictores sin que el rendimiento del modelo se ve comprometido. Esto permite simplificar el modelo y acelera el entrenamiento. Además, se reduce el riesgo de sobreajuste porque es menos probable que el modelo aprenda ruido de características irrelevantes.
Forecasting probabilístico
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. Skforecast implementa varios métodos para el forecasting probabilístico:
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 [5, 95], lo que significa que los intervalos se calculan con los percentiles 5 y 95, lo que da como resultado una cobertura teórica del 90%.
# 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_select], steps =24, interval = [5, 95], method ='conformal' )predicciones.head()
pred
lower_bound
upper_bound
2014-01-01 00:00:00
3657.497280
3617.677953
3697.316607
2014-01-01 01:00:00
3820.682453
3780.863126
3860.501779
2014-01-01 02:00:00
3958.577124
3911.746620
4005.407628
2014-01-01 03:00:00
3997.820878
3950.990375
4044.651382
2014-01-01 04:00:00
4082.630648
4035.800144
4129.461152
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.
Si además de los residuos, se le pasan las correspondientes predicciones al método set_out_sample_residuals(), entonces los residuos utilizados en el proceso de bootstrapping se seleccionan condicionados al rango de valores de las predicciones. Esto puede ayudar a conseguir intervalos con mayor covertura a la vez que se mantienen lo más estrechos posibles.
# Backtesting sobre los datos de validación para obtener los residuos out-sample# ==============================================================================cv = TimeSeriesFold( steps =24, initial_train_size =len(datos.loc[:fin_train]),)_, predicciones_val = backtesting_forecaster( forecaster = forecaster, y = datos.loc[:fin_validacion, 'Demand'], # Train + Validation exog = datos.loc[:fin_validacion, exog_select], cv = cv, metric ='mean_absolute_error' )
# Almacenar residuos out-sample en el forecaster# ==============================================================================forecaster.set_out_sample_residuals( y_true = datos.loc[predicciones_val.index, 'Demand'], y_pred = predicciones_val['pred'])
A continuación, se ejecuta el proceso de backtesting para estimar los intervalos de predicción en el conjunto de test. Se indica el argumento use_in_sample_residuals en False para que se utilicen los residuos out-sample almacenados previamente y use_binned_residuals para que el los residuos utilizados en el bootstrapping se seleccionen condicionados al rango de valores de las predicciones.
# Backtest con intervalos de predicción para el conjunto de test utilizando los# residuos out-sample# ==============================================================================cv = TimeSeriesFold( steps =24, initial_train_size =len(datos.loc[:fin_validacion]), refit =False,)metrica, predicciones = backtesting_forecaster( forecaster = forecaster, y = datos['Demand'], exog = datos[exog_select], cv = cv, metric ='mean_absolute_error', interval = [5, 95], interval_method ='conformal', use_in_sample_residuals =False, # Se utilizan los residuos out-sample use_binned_residuals =True, # Residuos condicionados a los valores predichos )predicciones.head(5)
# Cobertura del intervalo predicho# ==============================================================================cobertura = calculate_coverage( y_true = datos.loc[fin_validacion:, "Demand"], lower_bound = predicciones["lower_bound"], upper_bound = predicciones["upper_bound"] )area = (predicciones['upper_bound'] - predicciones['lower_bound']).sum()print(f"Área total del intervalo: {round(area, 2)}")print(f"Cobertura del intervalo predicho: {round(100* cobertura, 2)} %")
Área total del intervalo: 336324.75
Cobertura del intervalo predicho: 87.64 %
La cobertura observada de los intervalos es próxima a la cobertura teórica esperada (90%).
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
# Crear y entrenar el forecaster# ==============================================================================forecaster = ForecasterRecursive( regressor = LGBMRegressor(**best_params), lags = lags_select, window_features = window_features )forecaster.fit( y = datos.loc[:fin_validacion, 'Demand'], exog = datos.loc[:fin_validacion, exog_select])
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.
# Importancia de los predictores en el modelo# ==============================================================================feature_importances = forecaster.get_feature_importances()feature_importances.head(10)
feature
importance
0
lag_1
2312
53
Temperature
1257
1
lag_2
1132
5
lag_25
1070
2
lag_3
969
4
lag_24
951
3
lag_23
755
39
poly_hour_sin__hour_cos
690
47
Temperature_window_1D_mean
674
49
Temperature_window_1D_min
663
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.
# 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_select] )display(X_train.head(3))display(y_train.head(3))
lag_1
lag_2
lag_3
lag_23
lag_24
lag_25
lag_47
lag_48
lag_49
roll_mean_72
...
poly_hour_cos__sunset_hour_sin
poly_hour_cos__sunset_hour_cos
Temperature_window_1D_mean
Temperature_window_1D_max
Temperature_window_1D_min
Temperature_window_7D_mean
Temperature_window_7D_max
Temperature_window_7D_min
Temperature
Holiday
Time
2012-01-11 00:00:00
5087.274902
4925.770996
4661.530762
4988.975098
4963.185059
4896.571289
5039.643066
4958.630371
4886.074219
4386.356154
...
-0.707107
0.707107
16.140625
19.575001
13.175
19.321875
29.0
13.175
15.150
0.0
2012-01-11 01:00:00
5107.675781
5087.274902
4925.770996
5006.663086
4988.975098
4963.185059
5090.203125
5039.643066
4958.630371
4398.891398
...
-0.683013
0.683013
16.042707
19.575001
13.175
19.291666
29.0
13.175
15.425
0.0
2012-01-11 02:00:00
5126.337891
5107.675781
5087.274902
4977.723633
5006.663086
4988.975098
5108.837402
5090.203125
5039.643066
4410.526191
...
-0.612372
0.612372
15.939584
19.575001
13.175
19.255507
29.0
13.175
14.800
0.0
3 rows × 55 columns
Time
2012-01-11 00:00:00 5107.675781
2012-01-11 01:00:00 5126.337891
2012-01-11 02:00:00 5103.370605
Freq: h, Name: y, dtype: float32
# Crear SHAP explainer (para modelos basados en árboles)# ==============================================================================shap.initjs()explainer = shap.TreeExplainer(forecaster.regressor)# Se selecciona una muestra del 50% de los datos para acelerar el cálculorng = 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)
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.
Los valores SHAP no solo permiten interpretar el comportamiento general del modelo, sino que también son una herramienta poderosa para analizar predicciones individuales. Esto resulta especialmente útil cuando se quiere entender cómo se ha generado una predicción específica y qué variables han contribuido a ella.
Para llevar a cabo este análisis, es necesario acceder a los valores de los predictores —en este caso, los lags y las variables exógenas— en el momento de la predicción. Esto puede lograrse utilizando el método create_predict_X o bien activando el argumento return_predictors=True en la función backtesting_forecaster.
Supóngase que se quiere entender la predicción obtenida durante el backtesting para la fecha 2014-12-16 05:00:00.
# Backtesting indicando que se devuelvan los predictores# ==============================================================================cv = TimeSeriesFold( steps =24, initial_train_size =len(datos.loc[:fin_validacion]),)metrica, predicciones = backtesting_forecaster( forecaster = forecaster, y = datos['Demand'], exog = datos[exog_select], cv = cv, metric ='mean_absolute_error', return_predictors =True, )
Al indicar return_predictors=True, se obtiene un DataFrame con el valor predicho (‘pred’), la partición en la que se encuentra (‘fold’) y el valor de los lags y exógenas utilizadas para realizar cada predicción.
predicciones.head(4)
pred
fold
lag_1
lag_2
lag_3
lag_23
lag_24
lag_25
lag_47
lag_48
...
poly_hour_cos__sunset_hour_sin
poly_hour_cos__sunset_hour_cos
Temperature_window_1D_mean
Temperature_window_1D_max
Temperature_window_1D_min
Temperature_window_7D_mean
Temperature_window_7D_max
Temperature_window_7D_min
Temperature
Holiday
2014-12-01 00:00:00
5619.156432
0
5634.690430
5519.464844
5250.484863
4326.276855
4102.013672
3851.330322
4056.273926
4058.317627
...
-0.866025
0.500000
28.779167
33.849998
23.6
18.208036
33.849998
11.15
20.000000
0.0
2014-12-01 01:00:00
5597.705961
0
5619.156432
5634.690430
5519.464844
4560.871094
4326.276855
4102.013672
4060.130127
4056.273926
...
-0.836516
0.482963
28.556250
33.849998
20.0
18.224703
33.849998
11.15
20.350000
0.0
2014-12-01 02:00:00
5622.181666
0
5597.705961
5619.156432
5634.690430
4782.341309
4560.871094
4326.276855
4071.491943
4060.130127
...
-0.750000
0.433013
28.193750
33.849998
20.0
18.242857
33.849998
11.15
22.299999
0.0
2014-12-01 03:00:00
5637.505193
0
5622.181666
5597.705961
5619.156432
5034.609863
4782.341309
4560.871094
4146.851074
4071.491943
...
-0.612372
0.353553
27.802084
33.849998
20.0
18.273512
33.849998
11.15
23.850000
0.0
4 rows × 57 columns
# Waterfall para una predicción concreta# ==============================================================================predicciones = predicciones.astype(datos[exog_select].dtypes) # Asegurar que los tipos son igualesiloc_predicted_date = predicciones.index.get_loc('2014-12-16 05:00:00')shap_values_single = explainer(predicciones.iloc[:, 2:])shap.plots.waterfall(shap_values_single[iloc_predicted_date], show=False)fig, ax = plt.gcf(), plt.gca()fig.set_size_inches(8, 3.5)ax_list = fig.axesax = ax_list[0]ax.tick_params(labelsize=8)ax.setplt.show()
# Force plot para una predicción concreta# ==============================================================================shap.force_plot( base_value = shap_values_single.base_values[iloc_predicted_date], shap_values = shap_values_single.values[iloc_predicted_date], features = predicciones.iloc[iloc_predicted_date, 2:],)
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.
Forecaster direct multi-step
Los modelos ForecasterRecursive y ForecasterRecursiveCustom 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.
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 DD 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+1D+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:
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).
Se almacenan solo las predicciones del día siguiente, es decir, las 24 últimas.
Se añaden los datos de test hasta las 11:00 del día siguiente.
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.
# Final de initial_train_size + 12 posiciones# ==============================================================================datos.iloc[:len(datos.loc[:fin_validacion]) +12].tail(2)
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 librería skforecast.
