Modelos globales: Multi-series forecasting con Python y Skforecast

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

Modelos de forecasting globales: modelado de múltiples series temporales con machine learning

Joaquín Amat Rodrigo, Javier Escobar Ortiz
Octubre 2022 (última actualización Abril 2024)

Modelos de forecasting globales

En los modelos de series individuales (Local Forecasting Model), una única serie temporal se modela como una combinación lineal o no lineal de sus valores pasados (lags) y, opcionalmente, variables exógenas. Aunque este método proporciona una comprensión exhaustiva de cada serie, su escalabilidad puede verse comprometida cuando se dispone de un gran número de series.

Los modelos globales multiserie (Global Forecasting Model) consisten en un único modelo que tiene en cuenta varias series temporales de forma simultánea. Intentan captar los patrones básicos que rigen las series, mitigando así el ruido potencial que pueda introducir cada serie. Este enfoque es eficiente desde el punto de vista computacional, es fácil de mantener y puede conseguir mayor generalización, aunque potencialmente a costa de sacrificar cierta resolución individual. Se pueden distinguir dos estrategias de modelos de previsión global:

Múltiples series temporales independientes

En esta situación, cada serie temporal es independiente de las demás o, dicho de otro modo, los valores pasados de una serie no se utilizan como predictores de las otras series. ¿Por qué es útil entonces modelar todo junto? Aunque las series no dependen unas de otras, pueden seguir el mismo patrón intrínseco en cuanto a sus valores pasados y futuros. Por ejemplo, en una misma tienda, las ventas de los productos A y B pueden no estar relacionadas, pero siguen la misma dinámica, la de la tienda.

Para predecir los siguientes n steps, se sigue una estrategia recurisva, recursive multi-step forecasting.

Múltiples series temporales dependientes

Todas las series se modelan teniendo en cuenta que cada serie temporal depende no sólo de sus valores pasados, sino también de los valores pasados de las demás series. Se espera que el modelo no sólo aprenda la información de cada serie por separado, sino que también las relacione. Por ejemplo, las mediciones realizadas por todos los sensores (caudal, temperatura, presión...) instalados en una máquina industrial como un compresor. Series temporales multivariantes user guide.

🖉 Nota

La clase ForecasterAutoregMultiSeries y ForecasterAutoregMultiSeriesCustom permite modelar series temporales independientes. API Reference La clase ForecasterAutoregMultivariate permite modelar series temporales dependientes. API Reference

Ventajas y limitaciones

Los modelos de forecasting globales no siempre superan al modelado individual de cada serie. Cuál de funciona mejor depende en gran medida de las características del caso de uso al que se aplican. Sin embargo, conviene tener en cuenta la siguiente heurística:

Ventajas de los modelos multiseries

  • Es más fácil mantener y controlar un solo modelo que varios.

  • Dado que todas las series temporales se combinan durante el entrenamiento, cuando las series sean cortas (pocos datos) el modelo tendrá una mayor capacidad de aprendizaje al disponer de más observaciones.

  • Al combinar múltiples series temporales, el modelo puede aprender patrones más generalizables.

Desventajas de los modelos multiseries

  • Si las series no siguen la misma dinámica interna, el modelo puede aprender un patrón que no represente a ninguna de ellas.

  • Las series pueden enmascararse unas a otras, por lo que el modelo puede no predecirlas todas con el mismo rendimiento.

  • Es más exigente desde el punto de vista computacional (tiempo y recursos) entrenar y realizar backtesting de un modelo grande que de varios pequeños.

Caso 1: Series con la misma longitud

El objetivo de este estudio es comparar los resultados obtenidos por un modelo multiserie frente a la utilización de un modelo diferente para cada serie.

Los datos se han obtenido de Store Item Demand Forecasting Challenge. Este dataset contiene 913.000 transacciones de ventas desde el 2013-01-01 hasta el 2017-12-31 para 50 productos (SKU) en 10 tiendas. El objetivo es predecir las ventas de los próximos 7 días de 50 artículos diferentes en una tienda utilizando el historial disponible de 5 años.

Librerías

In [1]:
# Librerías
# ======================================================================================
import sys
import os
import warnings
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from tqdm.notebook import tqdm
from statsmodels.graphics.tsaplots import plot_acf
import sklearn
from sklearn.preprocessing import StandardScaler
from sklearn.ensemble import  HistGradientBoostingRegressor
import skforecast
from skforecast.ForecasterAutoregMultiSeries import ForecasterAutoregMultiSeries
from skforecast.ForecasterAutoreg import ForecasterAutoreg
from skforecast.model_selection import backtesting_forecaster
from skforecast.model_selection import grid_search_forecaster
from skforecast.model_selection_multiseries import backtesting_forecaster_multiseries
from skforecast.model_selection_multiseries import grid_search_forecaster_multiseries
from skforecast.plot import set_dark_theme
from skforecast.preprocessing import series_long_to_dict
from skforecast.preprocessing import exog_long_to_dict
warnings.filterwarnings('once')

print(f"Versión skforecasy: {skforecast.__version__}")
print(f"Versión scikit-learn: {sklearn.__version__}")
Versión skforecasy: 0.12.0
Versión scikit-learn: 1.4.2

Datos

In [2]:
# Lectura de datos
# ======================================================================================
data = pd.read_csv('./train_stores_kaggle.csv')
display(data)
print(f"Shape: {data.shape}")
date store item sales
0 2013-01-01 1 1 13
1 2013-01-02 1 1 11
2 2013-01-03 1 1 14
3 2013-01-04 1 1 13
4 2013-01-05 1 1 10
... ... ... ... ...
912995 2017-12-27 10 50 63
912996 2017-12-28 10 50 59
912997 2017-12-29 10 50 74
912998 2017-12-30 10 50 62
912999 2017-12-31 10 50 82

913000 rows × 4 columns

Shape: (913000, 4)
In [3]:
# Preparación datos
# ======================================================================================
selected_store = 2 # Seleccionar una tienda
selected_items = data.item.unique()
data = data[(data['store'] == selected_store) & (data['item'].isin(selected_items))].copy()
data['date'] = pd.to_datetime(data['date'], format='%Y-%m-%d')
data = pd.pivot_table(
           data    = data,
           values  = 'sales',
           index   = 'date',
           columns = 'item'
       )
data.columns.name = None
data.columns = [f"item_{col}" for col in data.columns]
data = data.asfreq('1D')
data = data.sort_index()
data.head(4)
Out[3]:
item_1 item_2 item_3 item_4 item_5 item_6 item_7 item_8 item_9 item_10 ... item_41 item_42 item_43 item_44 item_45 item_46 item_47 item_48 item_49 item_50
date
2013-01-01 12.0 41.0 19.0 21.0 4.0 34.0 39.0 49.0 28.0 51.0 ... 11.0 25.0 36.0 12.0 45.0 43.0 12.0 45.0 29.0 43.0
2013-01-02 16.0 33.0 32.0 14.0 6.0 40.0 47.0 42.0 21.0 56.0 ... 19.0 21.0 35.0 25.0 50.0 52.0 13.0 37.0 25.0 57.0
2013-01-03 16.0 46.0 26.0 12.0 12.0 41.0 43.0 46.0 29.0 46.0 ... 23.0 20.0 52.0 18.0 56.0 30.0 5.0 45.0 30.0 45.0
2013-01-04 20.0 50.0 34.0 17.0 16.0 41.0 44.0 55.0 32.0 56.0 ... 15.0 28.0 50.0 24.0 57.0 46.0 19.0 32.0 20.0 45.0

4 rows × 50 columns

El dataset se divide en 3 particiones: una para el entrenamiento, otra para la validación y otra para test.

In [4]:
# Separación datos train-validation-test
# ======================================================================================
end_train = '2016-05-31 23:59:00'
end_val = '2017-05-31 23:59:00'

data_train = data.loc[:end_train, :].copy()
data_val   = data.loc[end_train:end_val, :].copy()
data_test  = data.loc[end_val:, :].copy()
print(f"Fechas train      : {data_train.index.min()} --- {data_train.index.max()}  (n={len(data_train)})")
print(f"Fechas validación : {data_val.index.min()} --- {data_val.index.max()}  (n={len(data_val)})")
print(f"Fechas test       : {data_test.index.min()} --- {data_test.index.max()}  (n={len(data_test)})")
Fechas train      : 2013-01-01 00:00:00 --- 2016-05-31 00:00:00  (n=1247)
Fechas validación : 2016-06-01 00:00:00 --- 2017-05-31 00:00:00  (n=365)
Fechas test       : 2017-06-01 00:00:00 --- 2017-12-31 00:00:00  (n=214)

Se representan cuatro de las series para comprender sus tendencias y patrones. Se recomienda al lector que grafique más de ellas para comprenderlas en profundidad.

In [5]:
# Gráfico series temporales
# ======================================================================================
set_dark_theme()
fig, axs = plt.subplots(4, 1, figsize=(7, 5), sharex=True)
data.iloc[:, :4].plot(
    legend   = True,
    subplots = True, 
    title    = 'Ventas de la tienda 2',
    ax       = axs, 
)
for ax in axs:
    ax.axvline(pd.to_datetime(end_train) , color='white', linestyle='--', linewidth=1.5)
    ax.axvline(pd.to_datetime(end_val) , color='white', linestyle='--', linewidth=1.5)
fig.tight_layout();
In [6]:
# Gráficos de autocorrelación
# ======================================================================================
fig, axes = plt.subplots(nrows=4, ncols=1, figsize=(7, 5), sharex=True)
axes = axes.flat
for i, col in enumerate(data.columns[:4]):
    plot_acf(data[col], ax=axes[i], lags=7*5)
    axes[i].set_ylim(-1, 1.1)
    axes[i].set_title(f'{col}', fontsize=10)
fig.tight_layout()
plt.show()

Los gráficos de autocorrelación muestran una clara asociación entre las ventas de un día y las del mismo día una semana antes. Este tipo de correlación es un indicio de que los modelos autorregresivos pueden funcionar bien. También hay una estacionalidad semanal común entre las series. Cuanto más similar sea la dinámica entre las series, más probable será que el modelo multiseries aprenda patrones útiles.

Forecasting individual para cada item

Se entrena un modelo diferente (Gradient Boosting Machine) para cada artículo de la tienda y se estima su error medio absoluto mediante backtesting.

In [7]:
# Entrenar y realizar backtesting de un modelo para cada item
# ======================================================================================
items = []
mae_values = []
predictions = {}

for i, item in enumerate(tqdm(data.columns)):
    # Definir el forecaster
    forecaster = ForecasterAutoreg(
                     regressor     = HistGradientBoostingRegressor(random_state=123),
                     lags          = 14,
                     transformer_y = StandardScaler()
                 )
    # Backtesting forecaster
    metric, preds = backtesting_forecaster(
                        forecaster         = forecaster,
                        y                  = data[item],
                        initial_train_size = len(data_train) + len(data_val),
                        steps              = 7,
                        metric             = 'mean_absolute_error',
                        refit              = False,
                        fixed_train_size   = False,
                        verbose            = False,
                        show_progress      = False
                    )
    items.append(item)
    mae_values.append(metric)
    predictions[item] = preds

# Resultados
uni_series_mae = pd.Series(
                     data  = mae_values,
                     index = items,
                     name  = 'uni_series_mae'
                 )

Modelo global: Forecasting multiseries

Se entrena un único modelo con todas las series simultáneamente. En el método predict, y en el proceso de backtesting se indica el level para se desea realizar las predicciones. Si level es None, entonces se predicen todas las series disponibles.

In [8]:
# Entrenamiento y backtesting con un único modelo para todos los items
# ======================================================================================
items = list(data.columns)

# Definir el forecaster
forecaster_ms = ForecasterAutoregMultiSeries(
                    regressor          = HistGradientBoostingRegressor(random_state=8523),
                    lags               = 14,
                    encoding           = 'ordinal',
                    transformer_series = StandardScaler(),
                )
# Backtesting forecaster para todos los items
multi_series_mae, predictions_ms = backtesting_forecaster_multiseries(
                                       forecaster         = forecaster_ms,
                                       series             = data,
                                       levels             = items,
                                       steps              = 7,
                                       metric             = 'mean_absolute_error',
                                       initial_train_size = len(data_train) + len(data_val),
                                       refit              = False,
                                       fixed_train_size   = False,
                                       verbose            = False,
                                       show_progress      = True  
                                   )
# Resultados
display(multi_series_mae.head(3))
print('')
display(predictions_ms.head(3))
levels mean_absolute_error
0 item_1 5.564468
1 item_2 9.348777
2 item_3 7.405529

item_1 item_2 item_3 item_4 item_5 item_6 item_7 item_8 item_9 item_10 ... item_41 item_42 item_43 item_44 item_45 item_46 item_47 item_48 item_49 item_50
2017-06-01 34.974911 91.828608 61.151453 35.780062 30.521633 94.265236 92.233415 125.498303 86.283537 120.629623 ... 36.486671 58.575248 78.914890 44.869470 125.695034 97.118028 37.863895 78.212716 48.474518 109.882483
2017-06-02 38.281849 98.901712 61.956171 37.042634 30.804250 100.879910 97.938244 136.358715 91.309997 123.439280 ... 40.501463 58.849345 84.768281 51.738468 139.721728 103.052595 41.185425 83.371411 52.687136 114.615478
2017-06-03 39.945658 101.468483 66.491690 40.452524 33.131752 108.150494 102.250414 145.422706 94.955882 135.510597 ... 41.601258 65.569756 90.529982 52.617657 141.301919 107.370885 41.126271 89.913940 47.741813 118.617881

3 rows × 50 columns

Comparación

In [9]:
# Diferencia de la métrica de backtesting para cada item
# ======================================================================================
multi_series_mae = multi_series_mae.set_index('levels')
multi_series_mae.columns = ['multi_series_mae']
results = pd.concat((uni_series_mae, multi_series_mae), axis = 1)
results['improvement'] = results.eval('uni_series_mae - multi_series_mae')
results['improvement_(%)'] = 100 * results.eval('(uni_series_mae - multi_series_mae) / uni_series_mae')
results = results.round(2)
results.style.bar(subset=['improvement_(%)'], align='mid', color=['#d65f5f', '#5fba7d'])
Out[9]:
  uni_series_mae multi_series_mae improvement improvement_(%)
item_1 6.190000 5.560000 0.630000 10.120000
item_2 9.850000 9.350000 0.500000 5.050000
item_3 8.660000 7.410000 1.260000 14.520000
item_4 5.430000 4.960000 0.470000 8.640000
item_5 5.000000 4.630000 0.370000 7.320000
item_6 10.360000 9.960000 0.400000 3.860000
item_7 10.040000 9.900000 0.150000 1.450000
item_8 11.300000 10.360000 0.940000 8.330000
item_9 9.450000 8.720000 0.730000 7.760000
item_10 11.570000 10.480000 1.090000 9.430000
item_11 11.220000 10.240000 0.980000 8.780000
item_12 12.070000 10.870000 1.190000 9.880000
item_13 11.950000 11.510000 0.440000 3.680000
item_14 10.170000 9.490000 0.680000 6.660000
item_15 12.870000 11.410000 1.470000 11.410000
item_16 6.090000 5.910000 0.180000 2.950000
item_17 7.690000 7.230000 0.460000 5.940000
item_18 12.550000 11.970000 0.580000 4.610000
item_19 7.880000 7.400000 0.480000 6.130000
item_20 8.370000 7.920000 0.450000 5.380000
item_21 8.520000 8.060000 0.460000 5.430000
item_22 11.850000 10.650000 1.190000 10.070000
item_23 7.410000 6.750000 0.650000 8.800000
item_24 10.710000 9.980000 0.730000 6.840000
item_25 12.520000 11.790000 0.730000 5.860000
item_26 9.110000 8.560000 0.550000 6.010000
item_27 5.520000 5.160000 0.350000 6.430000
item_28 13.170000 12.030000 1.130000 8.620000
item_29 10.710000 10.190000 0.530000 4.910000
item_30 8.410000 7.900000 0.510000 6.110000
item_31 10.720000 10.050000 0.670000 6.220000
item_32 9.550000 9.210000 0.340000 3.540000
item_33 9.420000 9.240000 0.180000 1.880000
item_34 6.540000 5.830000 0.710000 10.910000
item_35 10.770000 10.170000 0.600000 5.530000
item_36 11.550000 10.790000 0.760000 6.570000
item_37 6.530000 6.150000 0.370000 5.680000
item_38 12.060000 11.220000 0.840000 6.980000
item_39 8.050000 7.370000 0.680000 8.480000
item_40 7.220000 6.530000 0.690000 9.510000
item_41 5.690000 5.250000 0.440000 7.770000
item_42 7.290000 6.920000 0.370000 5.060000
item_43 8.650000 8.610000 0.040000 0.500000
item_44 6.530000 6.330000 0.200000 3.000000
item_45 12.760000 11.920000 0.850000 6.640000
item_46 10.050000 9.800000 0.250000 2.470000
item_47 5.220000 4.960000 0.270000 5.090000
item_48 9.100000 8.160000 0.940000 10.320000
item_49 6.170000 6.040000 0.130000 2.120000
item_50 12.000000 10.740000 1.250000 10.440000
In [10]:
# Mejora media de todos los items
# ======================================================================================
results[['improvement', 'improvement_(%)']].agg(['mean', 'min', 'max'])
Out[10]:
improvement improvement_(%)
mean 0.6172 6.5938
min 0.0400 0.5000
max 1.4700 14.5200
In [11]:
# Número de series con mejora positiva y negativa
# ======================================================================================
pd.Series(np.where(results['improvement_(%)'] < 0, 'negative', 'positive')).value_counts()
Out[11]:
positive    50
Name: count, dtype: int64

El modelo global logra una mejora media del 6,6% en comparación con el uso de un modelo individual para cada serie. Para todas las series, el error de predicción evaluado mediante backtesting es menor cuando se utiliza el modelo global. Este caso de uso demuestra que un modelo multiserie puede tener ventajas sobre varios modelos individuales a la hora de predecir series temporales que siguen una dinámica similar. Además de las mejoras potenciales en la predicción, también es importante tener en cuenta el beneficio de tener un único modelo que mantener y la velocidad de entrenamiento y predicción.

⚠ Warning

Esta comparación se ha realizado sin optimizar los hiperparámetros del modelo. Consulte la sección Ajuste de hiperparámetros y selección de lags para comprobar que las conclusiones se mantienen cuando los modelos se ajustan con la mejor combinación de hiperparámetros y lags.

Caso 2: Series de diferente longitud y distintas variables exógenas

En escenarios en los que se tienen que modelar múltiples series, es habitual que las series tengan longitudes distintas debido a diferencias en los momentos de inicio del registro de los datos. Para hacer frente a este escenario, la clase ForecasterAutoregMultiSeries permite modelizar simultáneamente series temporales de distintas longitudes y con distintas variables exógenas.

  • Cuando las series modeladas tienen longitudes diferentes, deben almacenarse en un diccionario Python. Las claves del diccionario son los nombres de las series y los valores son las propias series. Todas las series deben ser del tipo pandas.Series, tener un índice datetime y tener la misma frecuencia.

Series values Permitido
[NaN, NaN, NaN, NaN, 4, 5, 6, 7, 8, 9] ✔️
[0, 1, 2, 3, 4, 5, 6, 7, 8, NaN] ✔️
[0, 1, 2, 3, 4, NaN, 6, 7, 8, 9] ✔️
[NaN, NaN, 2, 3, 4, NaN, 6, 7, 8, 9] ✔️


  • Cuando se utilizan variables exógenas diferentes para cada serie o cuando las variables exógenas son las mismas pero tienen valores diferentes para cada serie, deben almacenarse en un diccionario. Las claves del diccionario son los nombres de las series y los valores son las propias variables exógenas. Todas las variables exógenas deben ser de tipo pandas.DataFrame o pandas.Series.

Datos

Los datos de este ejemplo están almacenados en "formato largo" en un único DataFrame. La columna series_id identifica la serie a la que pertenece cada observación. La columna timestamp contiene la fecha de la observación, y la columna value contiene el valor de la serie en esa fecha. Cada serie temporal tiene una longitud diferente.

Las variables exógenas se almacenan en un DataFrame separado, también en "formato largo". La columna series_id identifica la serie a la que pertenece cada observación. La columna timestamp contiene la fecha de la observación, y las columnas restantes contienen los valores de las variables exógenas en esa fecha.

In [12]:
# Lectura series con diferentes longitudes y variables exógenas
# ==============================================================================
series = pd.read_csv(
    'https://raw.githubusercontent.com/JoaquinAmatRodrigo/skforecast-datasets/main/data/demo_multi_series.csv'
)
exog = pd.read_csv(
    'https://raw.githubusercontent.com/JoaquinAmatRodrigo/skforecast-datasets/main/data/demo_multi_series_exog.csv'
)
series['timestamp'] = pd.to_datetime(series['timestamp'])
exog['timestamp'] = pd.to_datetime(exog['timestamp'])
display(series.head())
display(exog.head())
series_id timestamp value
0 id_1000 2016-01-01 1012.500694
1 id_1000 2016-01-02 1158.500099
2 id_1000 2016-01-03 983.000099
3 id_1000 2016-01-04 1675.750496
4 id_1000 2016-01-05 1586.250694
series_id timestamp sin_day_of_week cos_day_of_week air_temperature wind_speed
0 id_1000 2016-01-01 -0.433884 -0.900969 6.416639 4.040115
1 id_1000 2016-01-02 -0.974928 -0.222521 6.366474 4.530395
2 id_1000 2016-01-03 -0.781831 0.623490 6.555272 3.273064
3 id_1000 2016-01-04 0.000000 1.000000 6.704778 4.865404
4 id_1000 2016-01-05 0.781831 0.623490 2.392998 5.228913

Cuando las series tienen longitudes diferentes, los datos deben transformarse en un diccionario. Las claves del diccionario son los nombres de las series y los valores son las propias series. Para ello, se utiliza la función series_long_to_dict, que toma el DataFrame en "formato largo" y devuelve un diccionario de series.

Del mismo modo, cuando las variables exógenas son diferentes (valores o variables) para cada serie, los datos deben transformarse en un diccionario. Las claves del diccionario son los nombres de las series y los valores son las propias variables exógenas. Se utiliza la función exog_long_to_dict, que toma el DataFrame en "formato largo" y devuelve un diccionario de variables exógenas.

In [13]:
# Transform series and exog to dictionaries
# ==============================================================================
series_dict = series_long_to_dict(
    data = series,
    series_id = 'series_id',
    index = 'timestamp',
    values = 'value',
    freq = 'D'
)

exog_dict = exog_long_to_dict(
    data = exog,
    series_id = 'series_id',
    index = 'timestamp',
    freq = 'D'
)

Algunas variables exógenas se omiten en las series 1 y 3 para ilustrar que se pueden utilizar diferentes variables exógenas para cada serie.

In [14]:
exog_dict['id_1000'] = exog_dict['id_1000'].drop(columns=['air_temperature', 'wind_speed'])
exog_dict['id_1003'] = exog_dict['id_1003'].drop(columns=['cos_day_of_week'])
In [15]:
# Particiones de entrenamiento y test
# ==============================================================================
end_train = '2016-07-31 23:59:00'
series_dict_train = {k: v.loc[: end_train,] for k, v in series_dict.items()}
exog_dict_train   = {k: v.loc[: end_train,] for k, v in exog_dict.items()}
series_dict_test  = {k: v.loc[end_train:,] for k, v in series_dict.items()}
exog_dict_test    = {k: v.loc[end_train:,] for k, v in exog_dict.items()}
In [16]:
# Gráfico series
# ==============================================================================
set_dark_theme()
colors = plt.rcParams['axes.prop_cycle'].by_key()['color']
fig, axs = plt.subplots(5, 1, figsize=(8, 4), sharex=True)
for i, s in enumerate(series_dict.values()):
    axs[i].plot(s, label=s.name, color=colors[i])
    axs[i].legend(loc='upper right', fontsize=8)
    axs[i].tick_params(axis='both', labelsize=8)
    axs[i].axvline(pd.to_datetime(end_train) , color='white', linestyle='--', linewidth=1)
In [17]:
# Descripción de cada serie
# ==============================================================================
for k in series_dict.keys():
    print(f"{k}:")
    try:
        print(
            f"\tTrain: len={len(series_dict_train[k])}, {series_dict_train[k].index[0]}"
            f" --- {series_dict_train[k].index[-1]}"
        )
    except:
        print(f"\tTrain: len=0")
    try:
        print(
            f"\tTest : len={len(series_dict_test[k])}, {series_dict_test[k].index[0]}"
            f" --- {series_dict_test[k].index[-1]}"
        )
    except:
        print(f"\tTest : len=0")
id_1000:
	Train: len=213, 2016-01-01 00:00:00 --- 2016-07-31 00:00:00
	Test : len=153, 2016-08-01 00:00:00 --- 2016-12-31 00:00:00
id_1001:
	Train: len=30, 2016-07-02 00:00:00 --- 2016-07-31 00:00:00
	Test : len=153, 2016-08-01 00:00:00 --- 2016-12-31 00:00:00
id_1002:
	Train: len=183, 2016-01-01 00:00:00 --- 2016-07-01 00:00:00
	Test : len=0
id_1003:
	Train: len=213, 2016-01-01 00:00:00 --- 2016-07-31 00:00:00
	Test : len=153, 2016-08-01 00:00:00 --- 2016-12-31 00:00:00
id_1004:
	Train: len=91, 2016-05-02 00:00:00 --- 2016-07-31 00:00:00
	Test : len=31, 2016-08-01 00:00:00 --- 2016-08-31 00:00:00
In [18]:
# Variables exógenas de cada serie
# ==============================================================================
for k in series_dict.keys():
    print(f"{k}:")
    try:
        print(f"\t{exog_dict[k].columns.to_list()}")
    except:
        print(f"\tNo variables exógenas")
id_1000:
	['sin_day_of_week', 'cos_day_of_week']
id_1001:
	['sin_day_of_week', 'cos_day_of_week', 'air_temperature', 'wind_speed']
id_1002:
	['sin_day_of_week', 'cos_day_of_week', 'air_temperature', 'wind_speed']
id_1003:
	['sin_day_of_week', 'air_temperature', 'wind_speed']
id_1004:
	['sin_day_of_week', 'cos_day_of_week', 'air_temperature', 'wind_speed']

Modelo global: Forecaster multiseries

In [19]:
# Fit forecaster
# ==============================================================================
regressor = HistGradientBoostingRegressor(random_state=123, max_depth=5)
forecaster = ForecasterAutoregMultiSeries(
                regressor          = regressor,
                lags               = 14,
                encoding           = "ordinal",
                dropna_from_series = False
            )
forecaster.fit(
    series = series_dict_train,
    exog   = exog_dict_train,
    suppress_warnings = True
)
forecaster
Out[19]:
============================ 
ForecasterAutoregMultiSeries 
============================ 
Regressor: HistGradientBoostingRegressor(max_depth=5, random_state=123) 
Lags: [ 1  2  3  4  5  6  7  8  9 10 11 12 13 14] 
Transformer for series: StandardScaler() 
Transformer for exog: None 
Series encoding: ordinal 
Window size: 14 
Series levels (names): ['id_1000', 'id_1001', 'id_1002', 'id_1003', 'id_1004'] 
Series weights: None 
Weight function included: False 
Differentiation order: None 
Exogenous included: True 
Type of exogenous variable: <class 'dict'> 
Exogenous variables names: ['sin_day_of_week', 'cos_day_of_week', 'air_temperature', 'wind_speed'] 
Training range: 'id_1000': ['2016-01-01', '2016-07-31'], 'id_1001': ['2016-07-02', '2016-07-31'], 'id_1002': ['2016-01-01', '2016-07-01'], ... 
Training index type: DatetimeIndex 
Training index frequency: D 
Regressor parameters: categorical_features: warn, early_stopping: auto, interaction_cst: None, l2_regularization: 0.0, learning_rate: 0.1, ... 
fit_kwargs: {} 
Creation date: 2024-05-06 11:53:32 
Last fit date: 2024-05-06 11:53:32 
Skforecast version: 0.12.0 
Python version: 3.11.8 
Forecaster id: None 
In [20]:
predicciones = forecaster.predict(steps=5, exog=exog_dict_test, suppress_warnings=True)
predicciones
Out[20]:
id_1000 id_1001 id_1003 id_1004
2016-08-01 1335.591458 2616.151907 2795.226332 8002.839878
2016-08-02 1472.518234 2416.633586 2426.971310 8557.739814
2016-08-03 1427.982938 2381.580647 2307.496916 8659.333218
2016-08-04 1382.228393 2348.445322 2158.193360 8701.647616
2016-08-05 1251.755669 2069.264652 2096.527610 9083.129771

Backtesting

In [21]:
# Backtesting
# ==============================================================================
forecaster = ForecasterAutoregMultiSeries(
    regressor=regressor, lags=14, encoding="ordinal", dropna_from_series=False
)

metrics_levels, backtest_predictions = backtesting_forecaster_multiseries(
    forecaster=forecaster,
    series=series_dict,
    exog=exog_dict,
    steps=24,
    metric="mean_absolute_error",
    initial_train_size=len(series_dict_train["id_1000"]),
    fixed_train_size=True,
    gap=0,
    allow_incomplete_fold=True,
    refit=False,
    n_jobs="auto",
    verbose=True,
    show_progress=True,
    suppress_warnings=True,
)

display(metrics_levels)
display(backtest_predictions)
Information of backtesting process
----------------------------------
Number of observations used for initial training: 213
Number of observations used for backtesting: 153
    Number of folds: 7
    Number of steps per fold: 24
    Number of steps to exclude from the end of each train set before test (gap): 0
    Last fold only includes 9 observations.

Fold: 0
    Training:   2016-01-01 00:00:00 -- 2016-07-31 00:00:00  (n=213)
    Validation: 2016-08-01 00:00:00 -- 2016-08-24 00:00:00  (n=24)
Fold: 1
    Training:   2016-01-01 00:00:00 -- 2016-07-31 00:00:00  (n=213)
    Validation: 2016-08-25 00:00:00 -- 2016-09-17 00:00:00  (n=24)
Fold: 2
    Training:   2016-01-01 00:00:00 -- 2016-07-31 00:00:00  (n=213)
    Validation: 2016-09-18 00:00:00 -- 2016-10-11 00:00:00  (n=24)
Fold: 3
    Training:   2016-01-01 00:00:00 -- 2016-07-31 00:00:00  (n=213)
    Validation: 2016-10-12 00:00:00 -- 2016-11-04 00:00:00  (n=24)
Fold: 4
    Training:   2016-01-01 00:00:00 -- 2016-07-31 00:00:00  (n=213)
    Validation: 2016-11-05 00:00:00 -- 2016-11-28 00:00:00  (n=24)
Fold: 5
    Training:   2016-01-01 00:00:00 -- 2016-07-31 00:00:00  (n=213)
    Validation: 2016-11-29 00:00:00 -- 2016-12-22 00:00:00  (n=24)
Fold: 6
    Training:   2016-01-01 00:00:00 -- 2016-07-31 00:00:00  (n=213)
    Validation: 2016-12-23 00:00:00 -- 2016-12-31 00:00:00  (n=9)

levels mean_absolute_error
0 id_1000 174.456354
1 id_1001 997.630352
2 id_1002 NaN
3 id_1003 315.758865
4 id_1004 910.152611
id_1000 id_1001 id_1003 id_1004
2016-08-01 1335.591458 2616.151907 2795.226332 8002.839878
2016-08-02 1472.518234 2416.633586 2426.971310 8557.739814
2016-08-03 1427.982938 2381.580647 2307.496916 8659.333218
2016-08-04 1382.228393 2348.445322 2158.193360 8701.647616
2016-08-05 1251.755669 2069.264652 2096.527610 9083.129771
... ... ... ... ...
2016-12-27 1690.936215 1681.918483 1904.545843 NaN
2016-12-28 1606.220869 1588.118897 1835.349567 NaN
2016-12-29 1469.959950 1675.787962 1743.484785 NaN
2016-12-30 1388.803603 1527.865827 1918.761755 NaN
2016-12-31 1272.047660 1294.140529 1937.245496 NaN

153 rows × 4 columns

In [22]:
# Gráfico predicciones backtesting
# ==============================================================================
colors = plt.rcParams['axes.prop_cycle'].by_key()['color']
fig, axs = plt.subplots(5, 1, figsize=(8, 4), sharex=True)
for i, s in enumerate(series_dict.keys()):
    axs[i].plot(series_dict[s], label=series_dict[s].name, color=colors[i])
    axs[i].axvline(pd.to_datetime(end_train) , color='white', linestyle='--', linewidth=1)
    try:
        axs[i].plot(backtest_predictions[s], label='prediction', color="white")
    except:
        pass
    axs[i].legend(loc='upper left', fontsize=8)
    axs[i].tick_params(axis='both', labelsize=8)

Al permitir la modelización de series temporales de diferentes longitudes y con diferentes variables exógenas, la clase ForecasterAutoregMultiSeries proporciona una herramienta flexible y potente para utilizar toda la información disponible para entrenar los modelos.

Caso 3: Optimización de hiperparámetros (tuning) y selección de lags

En caso 1, la comparación entre forecaster se ha realizado sin optimizar los hiperparámetros de los regresores. Para una comparación justa, se utiliza una estrategia de grid search con el fin de seleccionar la mejor configuración para cada forecaster. Véase más información en hyperparameter tuning and lags selection.

⚠ Warning

La sección siguiente puede requerir un tiempo de ejecución considerable (alrededor de 45 minutos). Siéntase libre de seleccionar sólo un subconjunto de items para acelerar la ejecución.
In [23]:
# Búsqueda de hiperparámetros y backtesting de un modelo para cada item
# ======================================================================================
with open(os.devnull, 'w') as devnull: # Ocultar prints
    sys.stdout = devnull
    items = []
    mae_values  = []

    lags_grid = [7, 14]
    param_grid = {
        'max_iter': [100, 500],
        'max_depth': [5, 10],
        'learning_rate': [0.01, 0.1]
    }

    for item in tqdm(data.columns):

        forecaster = ForecasterAutoreg(
                        regressor     = HistGradientBoostingRegressor(random_state=123),
                        lags          = 14,
                        transformer_y = StandardScaler()
                    )

        results_grid = grid_search_forecaster(
                        forecaster         = forecaster,
                        y                  = data.loc[:end_val, item],
                        lags_grid          = lags_grid,
                        param_grid         = param_grid,
                        steps              = 7,
                        metric             = 'mean_absolute_error',
                        initial_train_size = len(data_train),
                        refit              = False,
                        fixed_train_size   = False,
                        return_best        = True,
                        verbose            = False,
                        show_progress      = False
                    )

        metric, preds = backtesting_forecaster(
                            forecaster         = forecaster,
                            y                  = data[item],
                            initial_train_size = len(data_train) + len(data_val),
                            steps              = 7,
                            metric             = 'mean_absolute_error',
                            refit              = False,
                            fixed_train_size   = False,
                            verbose            = False,
                            show_progress      = False
                        )

        items.append(item)
        mae_values.append(metric)

    uni_series_mae = pd.Series(
                        data  = mae_values,
                        index = items,
                        name  = 'uni_series_mae'
                    )
    
sys.stdout = sys.__stdout__
In [24]:
# Búsqueda de hiperparámetros y backtesting para un modelo multiserie
# ======================================================================================
lags_grid = [7, 14]
param_grid = {
    'max_iter': [100, 500],
    'max_depth': [5, 10],
    'learning_rate': [0.01, 0.1]
}

forecaster_ms = ForecasterAutoregMultiSeries(
                    regressor          = HistGradientBoostingRegressor(random_state=123),
                    lags               = 14,
                    transformer_series = StandardScaler(),
                    encoding           = 'ordinal'
                )

results_grid_ms = grid_search_forecaster_multiseries(
                      forecaster         = forecaster_ms,
                      series             = data.loc[:end_val, :],
                      levels             = None, # Si es None se seleccionan todos los niveles
                      lags_grid          = lags_grid,
                      param_grid         = param_grid,
                      steps              = 7,
                      metric             = 'mean_absolute_error',
                      initial_train_size = len(data_train),
                      refit              = False,
                      fixed_train_size   = False,
                      return_best        = True,
                      verbose            = False
                  )      

multi_series_mae, predictions_ms = backtesting_forecaster_multiseries(
                                       forecaster         = forecaster_ms,
                                       series             = data,
                                       levels             = None, # Si es None se seleccionan todos los niveles
                                       steps              = 7,
                                       metric             = 'mean_absolute_error',
                                       initial_train_size = len(data_train) + len(data_val),
                                       refit              = False,
                                       fixed_train_size   = False,
                                       verbose            = False
                                   )
16 models compared for 50 level(s). Number of iterations: 16.
`Forecaster` refitted using the best-found lags and parameters, and the whole data set: 
  Lags: [ 1  2  3  4  5  6  7  8  9 10 11 12 13 14] 
  Parameters: {'learning_rate': 0.1, 'max_depth': 5, 'max_iter': 100}
  Backtesting metric: 8.27973886640365
  Levels: ['item_1', 'item_2', 'item_3', 'item_4', 'item_5', 'item_6', 'item_7', 'item_8', 'item_9', 'item_10', 'item_11', 'item_12', 'item_13', 'item_14', 'item_15', 'item_16', 'item_17', 'item_18', 'item_19', 'item_20', 'item_21', 'item_22', 'item_23', 'item_24', 'item_25', 'item_26', 'item_27', 'item_28', 'item_29', 'item_30', 'item_31', 'item_32', 'item_33', 'item_34', 'item_35', 'item_36', 'item_37', 'item_38', 'item_39', 'item_40', 'item_41', 'item_42', 'item_43', 'item_44', 'item_45', 'item_46', 'item_47', 'item_48', 'item_49', 'item_50']

-----
matplotlib          3.8.4
numpy               1.26.4
pandas              2.1.4
session_info        1.0.0
skforecast          0.12.0
sklearn             1.4.2
statsmodels         0.14.1
tqdm                4.66.2
-----
IPython             8.22.2
jupyter_client      8.6.1
jupyter_core        5.7.2
notebook            6.4.12
-----
Python 3.11.8 (main, Feb 26 2024, 21:39:34) [GCC 11.2.0]
Linux-5.15.0-1058-aws-x86_64-with-glibc2.31
-----
Session information updated at 2024-05-06 13:01
In [25]:
# Diferencia de la métrica de backtesting para cada item
# ======================================================================================
multi_series_mae = multi_series_mae.set_index('levels')
multi_series_mae.columns = ['multi_series_mae']
results = pd.concat((uni_series_mae, multi_series_mae), axis = 1)
results['improvement'] = results.eval('uni_series_mae - multi_series_mae')
results['improvement_(%)'] = 100 * results.eval('(uni_series_mae - multi_series_mae) / uni_series_mae')
results = results.round(2)

# Mejora media de todos los items
# ======================================================================================
results[['improvement', 'improvement_(%)']].agg(['mean', 'min', 'max'])
Out[25]:
improvement improvement_(%)
mean 0.4424 4.767
min -0.0200 -0.260
max 1.1800 12.820
In [26]:
# Número de series con mejora positiva y negativa
# ======================================================================================
pd.Series(np.where(results['improvement_(%)'] < 0, 'negative', 'positive')).value_counts()
Out[26]:
positive    48
negative     2
Name: count, dtype: int64

Tras identificar la combinación de lags e hiperparámetros que logran el mejor rendimiento predictivo para cada forecaster, un número superior de modelos univariantes han logrado una mayor capacidad predictiva al generalizar mejor sus propios datos (un item). Aun así, el modelo multiserie proporciona mejores resultados para la mayoría de los items.

Pesos en forecasting multiseries

Los pesos se utilizan para controlar la influencia que tiene cada observación en el entrenamiento del modelo. ForecasterAutoregMultiseries acepta dos tipos de pesos:

  • series_weights controla la importancia relativa de cada serie. Si una serie tiene el doble de peso que las demás, las observaciones de esa serie influyen el doble en el entrenamiento. Cuanto mayor sea el peso de una serie en relación con las demás, más se centrará el modelo en intentar aprender esa serie.

  • weight_func controla la importancia relativa de cada observación en función del índice. Por ejemplo, una función que asigna un peso menor a ciertas fechas.

Si se indican los dos tipos de pesos, estos se multiplican para crear los pesos finales como se muestra en la figura. El sample_weight resultante no puede contener valores negativos.

Más información sobre weights in multi-series forecasting y weighted time series forecasting con skforecast.

En este ejemplo, item_1 tiene una mayor importancia relativa entre series (pesa 3 veces más que el resto de series), y las observaciones entre '2013-12-01' y '2014-01-31' se consideran no representativas y se les aplica un peso de 0.

In [27]:
# Pesos en forecasting multiseries
# ======================================================================================

# Pesos de cada serie
series_weights = {'item_1': 3.} # Las series que no aparezcan en el dict tienen un peso de 1


# Pesos de cada observación (por índice)
def custom_weights(index):
    """
    Devuelve 0 si el índice está entre '2013-12-01' y '2014-01-31', 1 en caso contrario.
    """
    weights = np.where(
                  (index >= '2013-12-01') & (index <= '2014-01-31'),
                   0,
                   1
              )
    
    return weights

forecaster = ForecasterAutoregMultiSeries(
                 regressor          = HistGradientBoostingRegressor(random_state=123),
                 lags               = 14,
                 transformer_series = StandardScaler(),
                 transformer_exog   = None,
                 weight_func        = custom_weights, 
                 series_weights     = series_weights
             )

forecaster.fit(series=data)
forecaster.predict(steps=7).head(3)
/home/ubuntu/anaconda3/envs/skforecast_12_py11/lib/python3.11/site-packages/sklearn/ensemble/_hist_gradient_boosting/gradient_boosting.py:399: FutureWarning: The categorical_features parameter will change to 'from_dtype' in v1.6. The 'from_dtype' option automatically treats categorical dtypes in a DataFrame as categorical features.
  warnings.warn(
Out[27]:
item_1 item_10 item_11 item_12 item_13 item_14 item_15 item_16 item_17 item_18 ... item_46 item_47 item_48 item_49 item_5 item_50 item_6 item_7 item_8 item_9
2018-01-01 20.583876 61.580282 62.323468 62.581761 79.205140 56.450206 85.846570 22.718565 32.502652 79.377183 ... 53.421299 22.260726 47.144927 28.682841 22.372861 63.199237 53.271716 57.232493 68.655821 41.978625
2018-01-02 22.501538 71.421524 75.349079 75.283254 89.691063 59.121550 90.118244 24.848422 36.265055 87.467801 ... 62.983381 26.943935 51.961949 34.340793 19.502197 63.945234 57.812538 61.734525 82.100161 49.326013
2018-01-03 20.907183 73.185299 76.323151 75.547258 83.861397 63.508950 95.531382 25.331899 33.468072 90.110389 ... 62.890533 25.834892 52.942216 32.730445 18.498621 71.040223 61.123112 63.709698 84.395482 50.310361

3 rows × 50 columns

🖉 Nota

Se puede pasar un diccionario a `weight_func` para aplicar diferentes funciones a cada serie. Si una serie no se presenta en el diccionario, tendrá peso

Conclusiones

Este caso de uso muestra como un modelo multiserie puede presentar ventajas sobre varios modelos individuales cuando se predicen series temporales con una dinámica similar. Más allá de las posibles mejoras en la predicción, también es importante tener en cuenta la ventaja de tener un solo modelo que mantener.

Información de sesión

In [28]:
import session_info
session_info.show(html=False)

Instrucciones para citar

¿Cómo citar este documento?

Si utilizas este documento o alguna parte de él, te agradecemos que lo cites. ¡Muchas gracias!

Modelos de forecasting globales: modelado de múltiples series temporales con machine learning por Joaquín Amat Rodrigo and Javier Escobar Ortiz, disponible con licencia Attribution-NonCommercial-ShareAlike 4.0 International (CC BY-NC-SA 4.0 DEED) en https://www.cienciadedatos.net/documentos/py44-multi-series-forecasting-skforecast-español.html

¿Cómo citar skforecast?

Si utilizas skforecast en tu investigación o publicación, te lo agradeceríamos mucho que lo cites. ¡Muchas gracias!

Zenodo:

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

APA:

Amat Rodrigo, J., & Escobar Ortiz, J. (2024). skforecast (Version 0.12.0) [Computer software]. https://doi.org/10.5281/zenodo.8382788

BibTeX:

@software{skforecast, author = {Amat Rodrigo, Joaquin and Escobar Ortiz, Javier}, title = {skforecast}, version = {0.12.0}, month = {4}, year = {2024}, license = {BSD-3-Clause}, url = {https://skforecast.org/}, doi = {10.5281/zenodo.8382788} }


¿Te ha gustado el artículo? Tu ayuda es importante

Mantener un sitio web tiene unos costes elevados, tu contribución me ayudará a seguir generando contenido divulgativo gratuito. ¡Muchísimas gracias! 😊


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