Forecasting de demanda intermitente con skforecast

Forecasting demanda intermitente con skforecast

Joaquín Amat Rodrigo, Javier Escobar Ortiz
April, 2023

Introducción


El forecasting de la demanda intermitente es un método estadístico utilizado para predecir la demanda de productos que tienen patrones de venta esporádicos o irregulares. Este tipo de productos se caracteriza por periodos de gran demanda seguidos de periodos de escasa o nula demanda. La previsión de la demanda intermitente se utiliza en las industrias manufacturera, minorista y sanitaria para gestionar los niveles de inventario, optimizar los programas de producción y reducir las roturas de stock y los costes por exceso de inventario.

La demanda intermitente regular se refiere a patrones de demanda predecibles con intervalos conocidos entre periodos de demanda. Este tipo de demanda intermitente se produce con cierta regularidad, como un producto que tiene una demanda estacional o un producto que se pide todos los meses en una fecha determinada. La previsión de la demanda intermitente regular puede realizarse mediante métodos estadísticos, como el método de Croston, que es una técnica habitual para prever la demanda en estas situaciones. Sin embargo, la mayoría de implementaciones no suelen permitir incorporar variables exógenas, a pesar de que pueden contener información muy relevante sobre los intervalos de demanda. Por esta razón, el forecasting basado en modelos de machine learning puede ser una alternativa interesante.

Por otro lado, la demanda intermitente irregular se refiere a patrones de demanda impredecibles, sin intervalos conocidos entre periodos de demanda. Este tipo de demanda intermitente es aleatoria e imprevisible, como un producto que se pide sólo unas pocas veces al año y en cantidades variables. Prever la demanda intermitente irregular es un reto porque no hay un patrón predecible para la demanda. Los métodos de previsión tradicionales pueden no funcionar eficazmente en estas situaciones y puede ser necesario recurrir a métodos de previsión alternativos, como el bootstrapping o la simulación.

En resumen, la demanda intermitente regular tiene un patrón predecible, mientras que la demanda intermitente irregular no lo tiene. La previsión de la demanda intermitente regular es más fácil que la previsión de la demanda intermitente irregular debido a la previsibilidad del patrón de demanda.

Este documento muestra cómo se puede utilizar la librería Python skforecast para predecir escenarios de demanda intermitente regular. Utilizando esta librería, el modelo de machine learning puede centrarse en aprender a predecir la demanda durante los periodos de interés, evitando la influencia de los periodos en los que no hay demanda.

Librerías

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

# Gráficos
# ==============================================================================
import matplotlib.pyplot as plt
import plotly.graph_objects as go
import plotly.express as px
import plotly.io as pio
pio.templates.default = "seaborn"

# Modelado y Forecasting
# ==============================================================================
from lightgbm import LGBMRegressor
from sklearn.preprocessing import FunctionTransformer
from sklearn.metrics import mean_absolute_error
from skforecast.ForecasterAutoreg import ForecasterAutoreg
from skforecast.ForecasterAutoregCustom import ForecasterAutoregCustom
from skforecast.model_selection import grid_search_forecaster
from skforecast.model_selection import backtesting_forecaster

Datos


Los datos utilizados en este ejemplo representan el número de usuarios que visitaron una tienda durante su horario de apertura de lunes a viernes, entre las 7:00 y las 20:00 horas. Por lo tanto, cualquier predicción fuera de este periodo no es útil y puede ignorarse o fijarse en 0.

In [2]:
# DEscarga de datos
# ======================================================================================
url = ('https://raw.githubusercontent.com/JoaquinAmatRodrigo/Estadistica-machine'
       '-learning-python/master/data/intermittent_demand.csv')
data = pd.read_csv(url, sep=',')
data['date_time'] = pd.to_datetime(data['date_time'], format='%Y-%m-%d %H:%M:%S')
data = data.set_index('date_time')
data = data.asfreq('H')
data = data.sort_index()
data.head(3)
Out[2]:
users
date_time
2011-01-01 00:00:00 0.0
2011-01-01 01:00:00 0.0
2011-01-01 02:00:00 0.0
In [3]:
# División train-valalidación-test
# ======================================================================================
end_train = '2012-03-31 23:59:00'
end_validation = '2012-08-31 23:59:00'
data_train = data.loc[: end_train, :]
data_val   = data.loc[end_train:end_validation, :]
data_test  = data.loc[end_validation:, :]

print(f"Fechas train      : {data_train.index.min()} --- {data_train.index.max()}  (n={len(data_train)})")
print(f"Fechas validation : {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      : 2011-01-01 00:00:00 --- 2012-03-31 23:00:00  (n=10944)
Fechas validation : 2012-04-01 00:00:00 --- 2012-08-31 23:00:00  (n=3672)
Fechas test       : 2012-09-01 00:00:00 --- 2012-12-31 23:00:00  (n=2928)
In [4]:
# Gráfico time series
# ======================================================================================
fig = go.Figure()
trace1 = go.Scatter(x=data_train.index, y=data_train['users'], name="train", mode="lines")
trace2 = go.Scatter(x=data_val.index, y=data_val['users'], name="validation", mode="lines")
trace3 = go.Scatter(x=data_test.index, y=data_test['users'], name="test", mode="lines")
fig.add_trace(trace1)
fig.add_trace(trace2)
fig.add_trace(trace3)
fig.update_layout(
    title="Serie temporal de usuarios",
    xaxis_title="Date time",
    yaxis_title="Usuarios",
    width  = 800,
    height = 400,
    margin=dict(l=20, r=20, t=50, b=20),
    legend=dict(
        orientation="h",
        yanchor="top",
        y=1.1,
        xanchor="left",
        x=0.001
    )
)
fig.show()
In [5]:
# Boxplot para la estacionalidad semanal
# ==============================================================================
data['week_day'] = data.index.day_of_week + 1
fig = px.box(
        data,
        x="week_day",
        y="users",
        title = 'Distribusión de usuarios por día de la semana',
        width=600,
        height=300
     )
median_values = data.groupby('week_day')['users'].median()
fig.add_trace(
    go.Scatter(
        x=median_values.index,
        y=median_values.values,
        mode='lines+markers',
        line=dict(color='blue', dash='dash'),
        showlegend=False
    )
)
fig.update_layout(margin=dict(l=20, r=20, t=35, b=20))
fig.show()
In [6]:
# Boxplot para la estacionalidad intra-diaria
# ==============================================================================
data['hour_day'] = data.index.hour + 1
fig = px.box(
        data,
        x="hour_day",
        y="users",
        title = 'Distribución de usuarios por hora del día',
        width=600,
        height=300
     )
median_values = data.groupby('hour_day')['users'].median()
fig.add_trace(
    go.Scatter(
        x=median_values.index,
        y=median_values.values,
        mode='lines+markers',
        line=dict(color='blue', dash='dash'),
        showlegend=False
    )
)
fig.update_layout(margin=dict(l=20, r=20, t=35, b=20))
fig.show()

Definir una métrica personalizada para evaluar el modelo

Para evaluar con adecuadamente el rendimiento del modelo, es crucial definir una métrica que refleje fielmente el escenario en el que se utilizará el modelo. En este caso, el rendimiento del modelo debe optimizarse durante los días laborables de 9:00 a 20:00.

Afortunadamente, skforecast ofrece la flexibilidad de utilizar funciones personalizadas como métricas para backtesting y búsqueda de hiperparámetros.

In [7]:
def custom_metric(y_true, y_pred):
    """
    Calcular el error medio absoluto utilizando sólo los valores predichos para
    los días laborables de 9:00 a 20:00.
    """

    day_of_week = y_true.index.day_of_week
    hour_of_day = y_true.index.hour
    mask = day_of_week.isin([0, 1, 2, 3, 4]) | ((hour_of_day > 7) | (hour_of_day < 20))
    metric = mean_absolute_error(y_true[mask], y_pred[mask])
    
    return metric

Forecasting

In [8]:
# Crear forecaster
# ==============================================================================
forecaster = ForecasterAutoreg(
                 regressor = LGBMRegressor(
                                 learning_rate = 0.1,
                                 max_depth     = 5,
                                 n_estimators  = 500,
                                 random_state  = 123,
                             ),
                 lags = 24
             )
forecaster
Out[8]:
================= 
ForecasterAutoreg 
================= 
Regressor: LGBMRegressor(max_depth=5, n_estimators=500, random_state=123) 
Lags: [ 1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24] 
Transformer for y: None 
Transformer for exog: None 
Window size: 24 
Weight function included: False 
Exogenous included: False 
Type of exogenous variable: None 
Exogenous variables names: None 
Training range: None 
Training index type: None 
Training index frequency: None 
Regressor parameters: {'boosting_type': 'gbdt', 'class_weight': None, 'colsample_bytree': 1.0, 'importance_type': 'split', 'learning_rate': 0.1, 'max_depth': 5, 'min_child_samples': 20, 'min_child_weight': 0.001, 'min_split_gain': 0.0, 'n_estimators': 500, 'n_jobs': -1, 'num_leaves': 31, 'objective': None, 'random_state': 123, 'reg_alpha': 0.0, 'reg_lambda': 0.0, 'silent': 'warn', 'subsample': 1.0, 'subsample_for_bin': 200000, 'subsample_freq': 0} 
Creation date: 2023-04-30 12:40:48 
Last fit date: None 
Skforecast version: 0.7.0 
Python version: 3.10.9 
Forecaster id: None 
In [9]:
# Backtesting del periodo de test
# ==============================================================================
metric, predictions = backtesting_forecaster(
                          forecaster         = forecaster,
                          y                  = data['users'],
                          initial_train_size = len(data.loc[:end_validation]),
                          fixed_train_size   = False,
                          steps              = 36,
                          refit              = False,
                          metric             = custom_metric,
                          verbose            = False # True para ver información de las iteraciones
                      )

print(f"Backtest error: {metric}")
Backtest error: 108.04513580389987
In [10]:
# Remplazar por cero las predicciones para las horas cerradas
# ==============================================================================
hour = data_test.index.hour
day_of_week = data_test.index.day_of_week
closed_hours = (hour < 7) | (hour > 20)
closed_days = day_of_week.isin([5, 6])
is_closed = (closed_hours) | (closed_days)
predictions[is_closed] = 0
In [11]:
# Gráfico real vs predicción en el periodo de test
# ======================================================================================
fig = go.Figure()
trace1 = go.Scatter(x=data_test.index, y=data_test['users'], name="test", mode="lines")
trace2 = go.Scatter(x=predictions.index, y=predictions['pred'], name="prediction", mode="lines")
fig.add_trace(trace1)
fig.add_trace(trace2)
fig.update_layout(
    title="Valores reales vs predicciones en el periodo de test",
    xaxis_title="Date time",
    yaxis_title="Users",
    width  = 800,
    height = 400,
    margin=dict(l=20, r=20, t=50, b=20),
    legend=dict(
        orientation="h",
        yanchor="top",
        y=1.1,
        xanchor="left",
        x=0.001
    )
)
fig.show()

Al analizar las predicciones, queda claro que el modelo tiene dificultades para captar con precisión las pautas de los horarios de apertura y cierre de las tiendas. Además, la influencia de los lags lleva a subestimar las primeras horas del día y los días siguientes a los de cierre.

In [12]:
# Distribución de los residuos por día de la semana y hora del día
# ==============================================================================
residuals = (predictions['pred'] - data_test['users']).to_frame('residuals')
residuals['week_day'] = residuals.index.day_of_week + 1
residuals['hour_day'] = residuals.index.hour + 1

fig = px.box(
        residuals,
        x="week_day",
        y="residuals",
        title = 'Distribución de los residuos por día de la semana',
        width=600,
        height=300
     )
fig.update_layout(margin=dict(l=20, r=20, t=35, b=20))
fig.show()

fig = px.box(
        residuals,
        x="hour_day",
        y="residuals",
        title = 'Distribución de los residuos por hora del día',
        width=600,
        height=300
     )
fig.update_layout(margin=dict(l=20, r=20, t=35, b=20))
fig.show()

Informar al modelo cuando la tienda está cerrada


Las variables exógenas pueden utilizarse en un modelo de forecasting para proporcionar información adicional y mejorar la capacidad del modelo para detectar patrones. Este enfoque ofrece la ventaja de incorporar factores externos que podrían influir en la exactitud de la previsión, lo que conduce a un modelo más fiable y preciso. El argumento exog facilita enormemente la integración de variables exógenas en el forecaster.

In [13]:
# Crear variable exogena para indicar si el local está cerrado
# ==============================================================================
hour = data.index.hour
day_of_week = data.index.day_of_week
closed_hours = (hour < 7) | (hour > 20)
closed_days = day_of_week.isin([5, 6])
is_closed = (closed_hours) | (closed_days)
data['is_closed'] = is_closed

end_train = '2012-03-31 23:59:00'
end_validation = '2012-08-31 23:59:00'
data_train = data.loc[: end_train, :]
data_val   = data.loc[end_train:end_validation, :]
data_test  = data.loc[end_validation:, :]
In [14]:
# Backtesting del periodo de test
# ==============================================================================
metric, predictions = backtesting_forecaster(
    forecaster         = forecaster,
    y                  = data['users'],
    exog               = data['is_closed'],
    initial_train_size = len(data.loc[:end_validation]),
    fixed_train_size   = False,
    steps              = 36,
    refit              = False,
    metric             = custom_metric,
    verbose            = False # True para ver información de las iteraciones
)

print(f"Backtest error: {metric}")
Backtest error: 52.93064967792484
In [15]:
# Remplazar por cero las predicciones para las horas cerradas
# ==============================================================================
hour = data_test.index.hour
day_of_week = data_test.index.day_of_week
closed_hours = (hour < 7) | (hour > 20)
closed_days = day_of_week.isin([5, 6])
is_closed = (closed_hours) | (closed_days)
predictions[is_closed] = 0
In [16]:
# Grafico de real vs predicción en el periodo de test
# ======================================================================================
fig = go.Figure()
trace1 = go.Scatter(x=data_test.index, y=data_test['users'], name="test", mode="lines")
trace2 = go.Scatter(x=predictions.index, y=predictions['pred'], name="prediction", mode="lines")
fig.add_trace(trace1)
fig.add_trace(trace2)
fig.update_layout(
    title="Valores reales vs predicciones en el periodo de test",
    xaxis_title="Date time",
    yaxis_title="Users",
    width  = 800,
    height = 400,
    margin=dict(l=20, r=20, t=50, b=20),
    legend=dict(
        orientation="h",
        yanchor="top",
        y=1.1,
        xanchor="left",
        x=0.001
    )
)
fig.show()

La incorporación de la variable exógena al modelo ha permitido una mejora significativa de las predicciones, reduciéndose el error a la mitad. Esto pone de relevancia la importancia de utilizar factores externos para mejorar el rendimiento del modelo de forecasting.

Model tunning


La optimización de hiperparámetros es un aspecto crítico del desarrollo de machine learning. Los hiperparámetros son valores que no pueden aprenderse a partir de los datos y que el usuario debe establecer antes de entrenar el modelo. Estos hiperparámetros pueden tener un impacto significativo en el rendimiento del modelo, y un ajuste cuidadoso puede mejorar su precisión y generalización a nuevos datos. En el caso de los modelos de forecasting, los lags incluidos en el modelo pueden considerarse un hiperparámetro adicional.

La optimización de hiperparámetros consiste en probar sistemáticamente distintos valores o combinaciones de hiperparámetros (incluidos los lags) para encontrar la configuración óptima que ofrezca los mejores resultados. La librería skforecast ofrece varias estrategias de ajuste de hiperparámetros, como la búsqueda en cuadrícula (grid), la búsqueda aleatoria y la búsqueda bayesiana, para identificar la combinación óptima de lags e hiperparámetros que logre el mejor rendimiento.

In [17]:
# Grid search de hiperparámetros y lags
# ==============================================================================
forecaster = ForecasterAutoreg(
                 regressor = LGBMRegressor(
                                 learning_rate = 0.1,
                                 max_depth     = 5,
                                 n_estimators  = 500,
                                 random_state  = 123,
                             ),
                 lags = 24
             )

# Lags utilizados
lags_grid = [24, 48, 72]

# Hiperparámetros del regresor
param_grid = {'n_estimators': [50, 100, 500],
              'max_depth': [5, 10, 15],
              'learning_rate': [0.01, 0.1]}

results_grid = grid_search_forecaster(
                   forecaster         = forecaster,
                   y                  = data.loc[:end_validation, 'users'],
                   exog               = data.loc[:end_validation, 'is_closed'],
                   param_grid         = param_grid,
                   lags_grid          = lags_grid,
                   steps              = 36,
                   refit              = False,
                   metric             = custom_metric,
                   initial_train_size = len(data.loc[:end_train]),
                   fixed_train_size   = False,
                   return_best        = True,
                   verbose            = False
               )
Number of models compared: 54.
loop lags_grid: 100%|███████████████████████████████████████| 3/3 [01:46<00:00, 35.38s/it]
`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 22 23 24
 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48
 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72] 
  Parameters: {'learning_rate': 0.1, 'max_depth': 15, 'n_estimators': 100}
  Backtesting metric: 39.30052275988041

In [18]:
# Backtesting periodo de test
# ==============================================================================
metric, predictions = backtesting_forecaster(
                          forecaster         = forecaster,
                          y                  = data['users'],
                          exog               = data['is_closed'],
                          initial_train_size = len(data.loc[:end_validation]),
                          fixed_train_size   = False,
                          steps              = 36,
                          refit              = False,
                          metric             = custom_metric,
                          verbose            = False 
                      )

print(f"Backtest error: {metric}")
Backtest error: 39.86155496691617
In [19]:
# Remplazar por cero las predicciones para las horas cerradas
# ==============================================================================
hour = data_test.index.hour
day_of_week = data_test.index.day_of_week
closed_hours = (hour < 7) | (hour > 20)
closed_days = day_of_week.isin([5, 6])
is_closed = (closed_hours) | (closed_days)
predictions[is_closed] = 0
In [20]:
# Gráfico de valores reales vs predicciones en el periodo de test
# ======================================================================================
fig = go.Figure()
trace1 = go.Scatter(x=data_test.index, y=data_test['users'], name="test", mode="lines")
trace2 = go.Scatter(x=predictions.index, y=predictions['pred'], name="prediction", mode="lines")
fig.add_trace(trace1)
fig.add_trace(trace2)
fig.update_layout(
    title="Valores reales vs predicciones en el periodo de test",
    xaxis_title="Date time",
    yaxis_title="Users",
    width  = 800,
    height = 400,
    margin=dict(l=20, r=20, t=50, b=20),
    legend=dict(
        orientation="h",
        yanchor="top",
        y=1.1,
        xanchor="left",
        x=0.001
    )
)
fig.show()

Se ha enconrado un conjunto de hiperparámetros mediante la búsqueda en cuadrícula que ha permitido reducir aun más de los errores de predicción.

Información de sesión

In [21]:
import session_info
session_info.show(html=False)
-----
lightgbm            3.3.5
matplotlib          3.6.3
numpy               1.23.5
pandas              1.5.3
plotly              5.13.1
session_info        1.0.0
skforecast          0.7.0
sklearn             1.2.2
-----
IPython             8.10.0
jupyter_client      7.3.4
jupyter_core        5.2.0
-----
Python 3.10.9 (main, Jan 11 2023, 15:21:40) [GCC 11.2.0]
Linux-5.15.0-1035-aws-x86_64-with-glibc2.31
-----
Session information updated at 2023-04-30 12:42

%%html

¿Cómo citar este documento?

Forecasting de demanda intermitente con skforecast by Joaquín Amat Rodrigo and Javier Escobar Ortiz, available under a CC BY-NC-SA 4.0 at https://www.cienciadedatos.net/documentos/py48-forecasting-demanda-intermitente.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.