Forecasting de demanda intermitente con skforecast

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

Forecasting demanda intermitente con skforecast

Joaquín Amat Rodrigo, Javier Escobar Ortiz
April, 2023 (última actualización Noviembre 2024)

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

# Gráficos
# ==============================================================================
import plotly.graph_objects as go
import plotly.express as px
import plotly.io as pio
import plotly.offline as poff
pio.templates.default = "seaborn"
poff.init_notebook_mode(connected=True)

# Modelado y Forecasting
# ==============================================================================
import skforecast
import sklearn
from lightgbm import LGBMRegressor
from sklearn.metrics import mean_absolute_error
from skforecast.recursive import ForecasterRecursive
from skforecast.model_selection import TimeSeriesFold
from skforecast.model_selection import grid_search_forecaster
from skforecast.model_selection import backtesting_forecaster

# Warnings
# ==============================================================================
import warnings
warnings.filterwarnings('once')

color = '\033[1m\033[38;5;208m' 
print(f"{color}Version skforecast: {skforecast.__version__}")
print(f"{color}Version scikit-learn: {sklearn.__version__}")
print(f"{color}Version pandas: {pd.__version__}")
print(f"{color}Version numpy: {np.__version__}")
Version skforecast: 0.14.0
Version scikit-learn: 1.5.1
Version pandas: 2.2.3
Version numpy: 2.0.2

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 [ ]:
# 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[ ]:
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  = 750,
    height = 350,
    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 = ForecasterRecursive(
                 regressor = LGBMRegressor(
                                 learning_rate = 0.1,
                                 max_depth     = 5,
                                 n_estimators  = 500,
                                 random_state  = 123,
                                 verbose       = -1
                             ),
                 lags = 24
             )
forecaster
Out[8]:

ForecasterRecursive

General Information
  • Regressor: LGBMRegressor(max_depth=5, n_estimators=500, random_state=123, verbose=-1)
  • 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]
  • Window features: None
  • Window size: 24
  • Exogenous included: False
  • Weight function included: False
  • Differentiation order: None
  • Creation date: 2024-11-04 16:20:14
  • Last fit date: None
  • Skforecast version: 0.14.0
  • Python version: 3.12.5
  • Forecaster id: None
Exogenous Variables
    None
Data Transformations
  • Transformer for y: None
  • Transformer for exog: None
Training Information
  • Training range: Not fitted
  • Training index type: Not fitted
  • Training index frequency: Not fitted
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': None, 'num_leaves': 31, 'objective': None, 'random_state': 123, 'reg_alpha': 0.0, 'reg_lambda': 0.0, 'subsample': 1.0, 'subsample_for_bin': 200000, 'subsample_freq': 0, 'verbose': -1}
Fit Kwargs
    {}

🛈 API Reference    🗎 User Guide

In [9]:
# Backtesting del periodo de test
# ==============================================================================
cv = TimeSeriesFold(
        steps              = 36,
        initial_train_size = len(data.loc[:end_validation]),
        refit              = False
     )
metric, predictions = backtesting_forecaster(
                          forecaster    = forecaster,
                          y             = data['users'],
                          cv            = cv,
                          metric        = custom_metric,
                          verbose       = False,
                          show_progress = True
                      )
metric
Out[9]:
custom_metric
0 108.045136
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  = 750,
    height = 350,
    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.astype(int)

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'],
                        cv            = cv,
                        metric        = custom_metric,
                        verbose       = False,
                        show_progress = True
                    )

metric
Out[14]:
custom_metric
0 52.93065
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 = ForecasterRecursive(
                 regressor = LGBMRegressor(
                                 learning_rate = 0.1,
                                 max_depth     = 5,
                                 n_estimators  = 500,
                                 random_state  = 123,
                                 verbose       = -1
                             ),
                 lags = 24
             )

# Lags
lags_grid = [24, 48, 72]

# Hiperparameters
param_grid = {
    "n_estimators": [50, 100, 500],
    "max_depth": [5, 10, 15],
    "learning_rate": [0.01, 0.1],
}

# Folds
cv_search = TimeSeriesFold(
                steps              = 36,
                initial_train_size = len(data.loc[:end_train]),
                refit              = False
            )

results_grid = grid_search_forecaster(
                   forecaster  = forecaster,
                   y           = data.loc[:end_validation, 'users'],
                   exog        = data.loc[:end_validation, 'is_closed'],
                   cv          = cv_search,
                   param_grid  = param_grid,
                   lags_grid   = lags_grid,
                   metric      = custom_metric,
                   return_best = True,
                   verbose     = False
               )
`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'],
                          cv            = cv,
                          metric        = custom_metric,
                          verbose       = False,
                          show_progress = True
                      )
metric
Out[18]:
custom_metric
0 39.861555
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  = 750,
    height = 350,
    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            4.5.0
numpy               2.0.2
pandas              2.2.3
plotly              5.24.1
session_info        1.0.0
skforecast          0.14.0
sklearn             1.5.1
-----
IPython             8.27.0
jupyter_client      8.6.3
jupyter_core        5.7.2
notebook            6.4.12
-----
Python 3.12.5 | packaged by Anaconda, Inc. | (main, Sep 12 2024, 18:27:27) [GCC 11.2.0]
Linux-5.15.0-1071-aws-x86_64-with-glibc2.31
-----
Session information updated at 2024-11-04 16:23

Instrucciones para citar

¿Cómo citar este documento?

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

Forecasting demanda intermitente con skforecast por Joaquín Amat Rodrigo y Javier Escobar Ortiz, disponible bajo una licencia Attribution-NonCommercial-ShareAlike 4.0 International (CC BY-NC-SA 4.0 DEED) en https://www.cienciadedatos.net/documentos/py48-forecasting-demanda-intermitente.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.14.0). Zenodo. https://doi.org/10.5281/zenodo.8382788

APA:

Amat Rodrigo, J., & Escobar Ortiz, J. (2024). skforecast (Version 0.14.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.14.0}, month = {11}, year = {2024}, license = {BSD-3-Clause}, url = {https://skforecast.org/}, doi = {10.5281/zenodo.8382788} }

%%html


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

Se permite:

  • Compartir: copiar y redistribuir el material en cualquier medio o formato.

  • Adaptar: remezclar, transformar y crear a partir del material.

Bajo los siguientes términos:

  • Atribución: Debes otorgar el crédito adecuado, proporcionar un enlace a la licencia e indicar si se realizaron cambios. Puedes hacerlo de cualquier manera razonable, pero no de una forma que sugiera que el licenciante te respalda o respalda tu uso.

  • NoComercial: No puedes utilizar el material para fines comerciales.

  • CompartirIgual: Si remezclas, transformas o creas a partir del material, debes distribuir tus contribuciones bajo la misma licencia que el original.