Modelos globales: Multi-series forecasting con Python y Skforecast

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

Multi-series forecasting con Python y Skforecast

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

Forecasting multiseries

En el forecasting de series temporales univariantes, una única serie temporal se modela como una combinación lineal o no lineal de sus lags. Es decir, los valores pasados de la serie se utilizan para predecir su comportamiento futuro. En el forecasting multiserie, dos o más series temporales se modelan conjuntamente mediante un único modelo. Se pueden distinguir dos estrategias:

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.

  Note

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

El forecasting multiserie no siempre superan al modelado de la serie individual. Cuál de ellas 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 de uso

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 numpy as np
import pandas as pd
from tqdm import tqdm
import matplotlib.pyplot as plt
plt.style.use('seaborn-v0_8-darkgrid')
from statsmodels.graphics.tsaplots import plot_acf

from sklearn.preprocessing import StandardScaler
from sklearn.ensemble import  HistGradientBoostingRegressor

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

Datos

In [2]:
# Descarga 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 del dato
# ======================================================================================
selected_store = 2 # Seleccionar una tienda
selected_items = data.item.unique() # Todos los items
#selected_items = [1, 2, 3, 4, 5] # Selección de items para reducir el tiempo de ejecución

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 41 19 21 4 34 39 49 28 51 ... 11 25 36 12 45 43 12 45 29 43
2013-01-02 16 33 32 14 6 40 47 42 21 56 ... 19 21 35 25 50 52 13 37 25 57
2013-01-03 16 46 26 12 12 41 43 46 29 46 ... 23 20 52 18 56 30 5 45 30 45
2013-01-04 20 50 34 17 16 41 44 55 32 56 ... 15 28 50 24 57 46 19 32 20 45

4 rows × 50 columns

⚠ Warning

El modelado de 50 *items* puede requerir un tiempo computacional significativo. Siéntase libre de seleccionar sólo un subconjunto de *items* para acelerar la ejecución.

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 dibujan cuatro de las series para comprender sus tendencias y patrones. Se recomienda encarecidamente al lector que grafique más de ellas para comprenderlas en profundidad.

In [5]:
# Gráfico series temporales
# ======================================================================================
fig, ax = plt.subplots(figsize=(8, 5))
data.iloc[:, :4].plot(
    legend   = True,
    subplots = True, 
    sharex   = True,
    title    = 'Ventas de la tienda 2',
    ax       = ax, 
);
/tmp/ipykernel_36202/2053242801.py:4: UserWarning: To output multiple subplots, the figure containing the passed axes is being cleared.
  data.iloc[:, :4].plot(
In [6]:
# Gráficos de autocorrelación
# ======================================================================================
fig, axes = plt.subplots(nrows=4, ncols=1, figsize=(8, 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_title(f'{col}')
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 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'
                 )
100%|██████████| 50/50 [00:32<00:00,  1.52it/s]

Forecasting multiseries

Se entrena un único modelo multiserie para predecir las ventas de cada producto en los próximos 7 días. En este caso, se debe indicar a qué item, level, se desea realizar el backtesting.

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

# Definir el forecaster
forecaster_ms = ForecasterAutoregMultiSeries(
                    regressor          = HistGradientBoostingRegressor(random_state=123),
                    lags               = 14,
                    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.579400
1 item_2 9.269760
2 item_3 7.410133

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.970517 90.425181 61.410945 35.095362 30.546009 92.451404 91.585215 127.189217 87.655798 120.136830 ... 37.011384 58.206884 78.447951 43.877367 125.384112 97.990755 37.572624 78.323208 49.073243 110.726256
2017-06-02 38.519069 99.603758 62.314616 37.184234 31.385269 100.206850 94.790210 135.543159 91.844634 124.148161 ... 41.811890 59.124785 85.434846 50.532889 137.563561 101.337474 41.599957 84.519961 52.524124 114.854633
2017-06-03 37.490938 101.565408 65.009803 40.614632 33.810371 108.394930 103.483348 143.192527 95.888009 133.022574 ... 42.815077 66.180392 87.478442 52.633649 140.414283 108.508473 41.503211 89.098035 47.387566 117.525369

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.580000 0.610000 9.880000
item_2 9.850000 9.270000 0.580000 5.850000
item_3 8.660000 7.410000 1.250000 14.470000
item_4 5.430000 4.980000 0.440000 8.170000
item_5 5.000000 4.660000 0.340000 6.820000
item_6 10.360000 10.040000 0.310000 3.020000
item_7 10.040000 9.800000 0.240000 2.390000
item_8 11.300000 10.420000 0.880000 7.800000
item_9 9.450000 8.770000 0.680000 7.240000
item_10 11.570000 10.570000 1.010000 8.700000
item_11 11.220000 10.320000 0.910000 8.080000
item_12 12.070000 10.980000 1.090000 9.010000
item_13 11.950000 11.370000 0.580000 4.870000
item_14 10.170000 9.600000 0.570000 5.650000
item_15 12.870000 11.380000 1.490000 11.600000
item_16 6.090000 6.030000 0.070000 1.110000
item_17 7.690000 7.290000 0.400000 5.170000
item_18 12.550000 12.010000 0.540000 4.300000
item_19 7.880000 7.330000 0.550000 6.930000
item_20 8.370000 7.870000 0.500000 5.970000
item_21 8.520000 8.100000 0.420000 4.960000
item_22 11.850000 10.790000 1.060000 8.940000
item_23 7.410000 6.720000 0.680000 9.230000
item_24 10.710000 9.950000 0.760000 7.080000
item_25 12.520000 11.740000 0.780000 6.240000
item_26 9.110000 8.590000 0.520000 5.760000
item_27 5.520000 5.150000 0.360000 6.610000
item_28 13.170000 11.870000 1.300000 9.860000
item_29 10.710000 10.380000 0.340000 3.140000
item_30 8.410000 7.820000 0.600000 7.110000
item_31 10.720000 10.050000 0.670000 6.210000
item_32 9.550000 9.230000 0.320000 3.300000
item_33 9.420000 9.300000 0.120000 1.280000
item_34 6.540000 5.860000 0.680000 10.430000
item_35 10.770000 10.230000 0.540000 4.980000
item_36 11.550000 10.680000 0.870000 7.520000
item_37 6.530000 6.090000 0.430000 6.660000
item_38 12.060000 11.270000 0.790000 6.510000
item_39 8.050000 7.270000 0.780000 9.730000
item_40 7.220000 6.530000 0.690000 9.560000
item_41 5.690000 5.250000 0.440000 7.800000
item_42 7.290000 6.980000 0.310000 4.290000
item_43 8.650000 8.480000 0.170000 1.970000
item_44 6.530000 6.390000 0.140000 2.110000
item_45 12.760000 11.780000 0.980000 7.710000
item_46 10.050000 9.700000 0.350000 3.480000
item_47 5.220000 4.990000 0.240000 4.510000
item_48 9.100000 8.200000 0.900000 9.900000
item_49 6.170000 5.930000 0.240000 3.810000
item_50 12.000000 10.650000 1.340000 11.180000
In [10]:
# Mejora media de todos los items
# ======================================================================================
results[['improvement', 'improvement_(%)']].agg(['mean', 'min', 'max'])
Out[10]:
improvement improvement_(%)
mean 0.6172 6.578
min 0.0700 1.110
max 1.4900 14.470
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
In [12]:
# Gráfico del item con la máxima mejora
# ======================================================================================
fig, ax=plt.subplots(figsize=(8, 4))
data_test['item_3'].tail(60).plot(ax=ax)
predictions['item_3'].tail(60).plot(ax=ax)
predictions_ms['item_3'].tail(60).plot(ax=ax)
ax.legend(['real', 'predicción forecaster individual', 'predicción multiseries forecaster'])
ax.set_xlabel('')
ax.set_title('Forecasting for Item 3');
In [13]:
# Gráfico del item con la menor mejora
# ======================================================================================
fig, ax=plt.subplots(figsize=(8, 4))
data_test['item_16'].tail(60).plot(ax=ax)
predictions['item_16'].tail(60).plot(ax=ax)
predictions_ms['item_16'].tail(60).plot(ax=ax)
ax.legend(['real', 'predicción forecaster individual', 'predicción multiseries forecaster'])
ax.set_xlabel('')
ax.set_title('Forecasting Item 16');

Impacto de la longitud de la serie

Si una serie temporal tiene pocas observaciones, la cantidad de información disponible para que el modelo aprenda es limitada. Este es un problema común en casos reales en los que no hay muchos datos históricos disponibles. Los forecasters multiseries no multivariantes combinan todas las series durante el entrenamiento, por lo que pueden acceder a más datos.

En esta sección, se realiza la misma comparación entre el forecasting de una sola serie y el forecasting multiserie, pero, esta vez, la longitud de las series es mucho menor.

In [14]:
# Separación datos train-validation-test
# ======================================================================================
start_train = '2017-01-01 00:00:00'
end_train = '2017-05-01 00:00:00'
end_val = '2017-07-31 23:59:00'
end_test = '2017-09-30 23:59:00'

data = data.loc[start_train:, :].copy()
data_train = data.loc[:end_train, :].copy()
data_val   = data.loc[end_train:end_val, :].copy()
data_test  = data.loc[end_val:end_test, :].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      : 2017-01-01 00:00:00 --- 2017-05-01 00:00:00  (n=121)
Fechas validación : 2017-05-01 00:00:00 --- 2017-07-31 00:00:00  (n=92)
Fechas test       : 2017-08-01 00:00:00 --- 2017-09-30 00:00:00  (n=61)
In [15]:
# Entrenar y realizar backtesting de un modelo para cada item
# ======================================================================================
items = []
mae_values = []
predictions = {}

for i, item in enumerate(tqdm(data.columns)):

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

    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

uni_series_mae = pd.Series(
                     data  = mae_values,
                     index = items,
                     name  = 'uni_series_mae'
                 )
100%|██████████| 50/50 [00:18<00:00,  2.75it/s]
In [16]:
# Entrenar y realizar backtesting con un único modelo para todos los items
# ======================================================================================

# Definir el forecaster
forecaster_ms = ForecasterAutoregMultiSeries(
                    regressor          = HistGradientBoostingRegressor(random_state=123),
                    lags               = 14,
                    transformer_series = StandardScaler(),
                )

# Backtesting forecaster para todos los items
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
                                   )

# Results
display(multi_series_mae.head(3))
print('')
display(predictions_ms.head(3))
levels mean_absolute_error
0 item_1 5.916400
1 item_2 10.241195
2 item_3 7.937301

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-08-02 38.799183 96.508038 59.663414 38.651256 32.065031 92.758344 97.089159 128.457251 88.622465 119.606434 ... 34.575309 59.960024 85.529356 52.144640 134.024475 96.528340 35.984520 88.372264 50.788241 121.528017
2017-08-03 38.712972 110.883112 65.769765 36.440864 36.865280 97.905376 108.625140 137.329917 93.654646 126.784431 ... 39.612689 63.967339 87.927004 58.093810 143.851475 109.246742 39.443198 96.031760 54.062914 123.135648
2017-08-04 40.099963 116.237826 69.827511 39.391681 38.751137 107.745055 115.040061 146.882289 103.897172 138.875529 ... 40.061436 65.427221 96.417850 62.312794 163.098872 120.418589 42.338213 98.744629 60.984447 130.095691

3 rows × 50 columns

In [17]:
# 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

# Mejora media de todos los items
# ======================================================================================
results[['improvement', 'improvement_(%)']].agg(['mean', 'min', 'max'])
Out[17]:
improvement improvement_(%)
mean 0.88 8.169
min -0.19 -3.390
max 2.25 18.440
In [18]:
# Número de series con mejora positiva y negativa
# ======================================================================================
pd.Series(np.where(results['improvement_(%)'] < 0, 'negative', 'positive')).value_counts()
Out[18]:
positive    49
negative     1
Name: count, dtype: int64

La mejora media ha pasado del 6.6 al 8.2%. La ventaja de utilizar un forecaster multiserie parece crecer a medida que se reduce la longitud de las series disponibles.

Optimización de hiperparámetros (tuning)

En las secciones anteriores, 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 [19]:
# Ocultar progress bar tqdm
# ======================================================================================
from tqdm import tqdm
from functools import partialmethod
tqdm.__init__ = partialmethod(tqdm.__init__, disable=True)
In [ ]:
# Búsqueda de hiperparámetros y backtesting de un modelo para cada item
# ======================================================================================
items = []
mae_values  = []

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

for i, item in enumerate(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'
                 )
In [21]:
# Búsqueda de hiperparámetros y backtesting para un modelo multiserie
# ======================================================================================
lags_grid = [7, 14, 21]
param_grid = {
    'max_iter': [100, 500],
    'max_depth': [3, 5, 10, None],
    'learning_rate': [0.01, 0.1]
}

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

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
                                   )
48 models compared for 50 level(s). Number of iterations: 48.
`Forecaster` refitted using the best-found lags and parameters, and the whole data set: 
  Lags: [ 1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21]
  Parameters: {'learning_rate': 0.01, 'max_depth': 3, 'max_iter': 500}
  Backtesting metric: 10.79552245163781
  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']

In [22]:
# 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[22]:
  uni_series_mae multi_series_mae improvement improvement_(%)
item_1 6.740000 6.070000 0.670000 9.940000
item_2 12.590000 11.290000 1.290000 10.280000
item_3 8.330000 8.030000 0.300000 3.610000
item_4 5.870000 5.810000 0.060000 0.980000
item_5 5.590000 5.230000 0.370000 6.560000
item_6 11.410000 11.260000 0.150000 1.320000
item_7 11.340000 11.270000 0.070000 0.660000
item_8 14.340000 13.910000 0.420000 2.940000
item_9 11.940000 10.740000 1.190000 9.980000
item_10 14.790000 13.630000 1.160000 7.850000
item_11 12.730000 13.090000 -0.360000 -2.850000
item_12 14.430000 14.110000 0.310000 2.180000
item_13 15.600000 13.920000 1.670000 10.740000
item_14 11.770000 10.840000 0.930000 7.890000
item_15 15.750000 15.210000 0.540000 3.410000
item_16 6.920000 6.830000 0.090000 1.320000
item_17 8.310000 7.890000 0.420000 5.030000
item_18 16.290000 15.190000 1.090000 6.720000
item_19 8.620000 8.290000 0.320000 3.750000
item_20 9.630000 9.320000 0.310000 3.210000
item_21 9.100000 8.570000 0.530000 5.870000
item_22 14.300000 13.830000 0.470000 3.300000
item_23 7.650000 7.510000 0.140000 1.830000
item_24 14.080000 12.450000 1.630000 11.580000
item_25 14.260000 15.180000 -0.920000 -6.460000
item_26 11.090000 10.680000 0.410000 3.720000
item_27 5.830000 5.390000 0.440000 7.530000
item_28 14.990000 14.720000 0.270000 1.830000
item_29 12.550000 12.630000 -0.080000 -0.650000
item_30 10.300000 8.770000 1.530000 14.860000
item_31 11.900000 11.610000 0.290000 2.420000
item_32 10.190000 10.620000 -0.430000 -4.250000
item_33 12.710000 12.230000 0.480000 3.780000
item_34 6.420000 6.130000 0.300000 4.660000
item_35 12.580000 12.540000 0.040000 0.290000
item_36 14.640000 13.670000 0.970000 6.660000
item_37 7.190000 6.780000 0.410000 5.670000
item_38 15.660000 14.650000 1.010000 6.480000
item_39 8.620000 8.730000 -0.110000 -1.220000
item_40 8.260000 7.220000 1.030000 12.500000
item_41 6.450000 5.770000 0.680000 10.510000
item_42 8.860000 7.840000 1.020000 11.540000
item_43 10.910000 10.340000 0.560000 5.180000
item_44 7.960000 6.980000 0.990000 12.430000
item_45 14.170000 15.230000 -1.060000 -7.450000
item_46 12.270000 12.040000 0.230000 1.890000
item_47 5.930000 5.800000 0.130000 2.210000
item_48 9.440000 9.680000 -0.240000 -2.550000
item_49 7.270000 7.260000 0.010000 0.160000
item_50 12.800000 12.600000 0.200000 1.550000
In [23]:
# Mejora media de todos los items
# ======================================================================================
results[['improvement', 'improvement_(%)']].agg(['mean', 'min', 'max'])
Out[23]:
improvement improvement_(%)
mean 0.4386 4.2278
min -1.0600 -7.4500
max 1.6700 14.8600
In [24]:
# Número de series con mejora positiva y negativa
# ======================================================================================
pd.Series(np.where(results['improvement_(%)'] < 0, 'negative', 'positive')).value_counts()
Out[24]:
positive    43
negative     7
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 [25]:
# Pesos en forecasting multiseries
# ======================================================================================

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


# Se puede pasar un diccionario a `weight_func` para aplicar diferentes funciones para cada serie
# Las series que no aparezcan en el dict tienen un peso de 1
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_09_py11/lib/python3.11/site-packages/skforecast/ForecasterAutoregMultiSeries/ForecasterAutoregMultiSeries.py:577: IgnoredArgumentWarning: {'item_44', 'item_8', 'item_11', 'item_5', 'item_37', 'item_42', 'item_25', 'item_13', 'item_14', 'item_24', 'item_36', 'item_48', 'item_18', 'item_3', 'item_19', 'item_45', 'item_49', 'item_31', 'item_23', 'item_10', 'item_12', 'item_16', 'item_7', 'item_50', 'item_46', 'item_34', 'item_39', 'item_32', 'item_38', 'item_6', 'item_41', 'item_47', 'item_4', 'item_28', 'item_21', 'item_26', 'item_29', 'item_43', 'item_35', 'item_30', 'item_17', 'item_22', 'item_9', 'item_20', 'item_27', 'item_33', 'item_15', 'item_40', 'item_2'} not present in `series_weights`. A weight of 1 is given to all their samples. 
 You can suppress this warning using: warnings.simplefilter('ignore', category=IgnoredArgumentWarning)
  warnings.warn(
Out[25]:
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
2018-01-01 20.109874 53.478138 36.348853 23.820464 23.794303 51.894103 55.784459 68.754883 40.716822 60.196889 ... 26.809501 34.880140 44.365278 27.547850 81.444684 53.123010 22.365198 47.706126 27.452497 64.263828
2018-01-02 22.510628 63.946284 40.750709 25.113519 19.989210 58.160237 61.580644 82.784148 48.883716 73.419269 ... 23.802414 39.521469 53.034414 30.260770 94.674182 59.403028 26.274588 50.975857 32.787511 61.725721
2018-01-03 20.182876 66.101162 39.411403 26.318990 19.076542 63.713520 63.230728 84.571132 52.234398 70.178156 ... 19.533864 40.694235 52.909487 32.867555 97.204512 62.760504 26.893860 52.239564 35.155388 71.421888

3 rows × 50 columns

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 [26]:
import session_info
session_info.show(html=False)
-----
matplotlib          3.7.2
numpy               1.24.4
pandas              2.0.3
session_info        1.0.0
skforecast          0.9.0
sklearn             1.3.0
statsmodels         0.14.0
tqdm                4.65.0
-----
IPython             8.14.0
jupyter_client      8.3.0
jupyter_core        5.3.1
notebook            6.5.4
-----
Python 3.11.4 (main, Jul  5 2023, 13:45:01) [GCC 11.2.0]
Linux-5.15.0-1039-aws-x86_64-with-glibc2.31
-----
Session information updated at 2023-07-13 14:32

¿Cómo citar este documento?

Modelos globales: Multi-series forecasting con python y skforecast por Joaquín Amat Rodrigo and Javier Escobar Ortiz, disponible con licencia a CC BY-NC-SA 4.0 en https://www.cienciadedatos.net/documentos/py44-multi-series-forecasting-skforecast-español.html DOI


¿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 contenido, creado por Joaquín Amat Rodrigo y Javier Escobar Ortiz, tiene licencia Attribution-NonCommercial-ShareAlike 4.0 International.