Reducir el impacto del Covid en modelos de forecasting

Joaquin Amat Rodrigo, Javier Escobar OrtizDiciembre, 2022 (última actualización Noviembre 2024)

Más sobre forecasting en: cienciadedatos.net

Introducción¶

Si bien en muchos casos reales de forecasting se dispone de datos históricos, no todos son fiables. Algunos ejemplos de estos escenarios son:

Sensores IoT: dentro del Internet de las cosas, los sensores se encargan de capturan los datos del mundo físico. A menudo los sensores se despliegan o instalan en entornos hostiles, lo que inevitablemente implica que los sensores sean propensos a fallos, mal funcionamiento y rápido deterioro, provocando que el sensor produzca lecturas inusuales y erróneas.
Paradas en instalaciones: cada cierto periodo de funcionamiento, las fábricas necesitan ser paradas para actividades de reparación, revisión o mantenimiento. Estos eventos provocan que la producción se detenga, generando un hueco en los datos.
Pandemia (Covid-19): la pandemia de Covid 19 cambió significativamente el comportamiento de la población, impactando directamente en muchas series temporales como la producción, las ventas y el transporte.

La presencia de valores no fiables o no representativos en el historial de datos es un problema importante, ya que dificulta el aprendizaje de los modelos. Para la mayoría de los algoritmos de forecasting, eliminar esa parte de los datos no es una opción, ya que requieren que la serie temporal sea completa. Una solución alternativa es reducir el peso de las observaciones afectadas durante el entrenamiento del modelo. Este documento cómo skforecast facilita la aplicación de esta estrategia.

✎ Nota

En los ejemplos siguientes, una parte de la serie temporal se excluye del entrenamiento del modelo dándole un peso de cero. Sin embargo, el uso de pesos no se limita a incluir o excluir observaciones, sino a equilibrar el grado de influencia de cada observación en el modelo de forecasting. Por ejemplo, una observación con un peso de 10 tiene 10 veces más impacto en el entrenamiento del modelo que una observación con un peso de 1.

⚠ Warning

En la mayoría de las implementaciones de gradient boosting (LightGBM, XGBoost, CatBoost), las observaciones con peso cero se ignoran al calcular los gradientes y hessianos. Sin embargo, los valores de esas observaciones todavía se consideran al construir los histogramas de las variables. Por lo tanto, el modelo resultante puede diferir del modelo entrenado sin las muestras con peso cero. Vea más detalles en este issue.

Librerías¶

# Procesamiento de datos
# ==============================================================================
import pandas as pd
import numpy as np
from skforecast.datasets import fetch_dataset

# Gráficos
# ==============================================================================
import matplotlib.pyplot as plt
from skforecast.plot import set_dark_theme

# Modelado y Forecasting
# ==============================================================================
import skforecast
import sklearn
from sklearn.linear_model import Ridge
from skforecast.recursive import ForecasterRecursive
from skforecast.model_selection import TimeSeriesFold
from skforecast.model_selection import backtesting_forecaster

# Configuracion
# ==============================================================================
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

Covid-19¶

Durante el periodo de confinamiento impuesto como consecuencia de la pandemia de covid-19, el comportamiento de la población se vio alterado. Un ejemplo de ello se puede ver en el uso del servicio de alquiler de bicicletas en la ciudad de Madrid (España).

Datos¶

# Descarga de datos
# ==============================================================================
data = fetch_dataset('bicimad')
data.head()

bicimad
-------
This dataset contains the daily users of the bicycle rental service (BiciMad) in
the city of Madrid (Spain) from 2014-06-23 to 2022-09-30.
The original data was obtained from: Portal de datos abiertos del Ayuntamiento
de Madrid https://datos.madrid.es/portal/site/egob
Shape of the dataset: (3022, 1)

# División datos entrenamiento, validation y test
# ==============================================================================
data = data.loc['2020-01-01': '2021-12-31']
end_train = '2021-06-01'
data_train = data.loc[: end_train, :]
data_test  = data.loc[end_train:, :]

print(f"Fechas train : {data_train.index.min()} --- {data_train.index.max()}  (n={len(data_train)})")
print(f"Fechas test  : {data_test.index.min()} --- {data_test.index.max()}  (n={len(data_test)})")

Fechas train : 2020-01-01 00:00:00 --- 2021-06-01 00:00:00  (n=518)
Fechas test  : 2021-06-01 00:00:00 --- 2021-12-31 00:00:00  (n=214)

# Time series plot
# ==============================================================================
set_dark_theme()
fig, ax = plt.subplots(figsize=(8, 3))
data_train.users.plot(ax=ax, label='train', linewidth=1)
data_test.users.plot(ax=ax, label='test', linewidth=1)
ax.axvspan(
    pd.to_datetime('2020-03-16'),
    pd.to_datetime('2020-04-21'), 
    label="Covid-19 confinamiento",
    color="red",
    alpha=0.3
)

ax.axvspan(
    pd.to_datetime('2020-04-21'),
    pd.to_datetime('2020-05-31'), 
    label="Tiempo de recuperación",
    color="white",
    alpha=0.3
)

ax.set_title('Número de usuarios diarios')
ax.legend();

Modelado con toda la serie¶

Se crea un forecaster sin tener en cuenta el periodo de confinamiento.

# Crear un objeto ForecasterRecursive
# ==============================================================================
forecaster = ForecasterRecursive(
                 regressor = Ridge(),
                 lags      = 21,
             )   
forecaster

Una vez creado el modelo, se realiza un proceso de backtesting para simular el comportamiento del forecaster si hubiera predicho el conjunto de test en bloques de 7 días.

# Backtesting
# ==============================================================================
cv = TimeSeriesFold(
        steps              = 7,
        initial_train_size = len(data.loc[:end_train]),
        fixed_train_size   = False,
        refit              = False,    
)
metric, predictions_backtest = backtesting_forecaster(
                                   forecaster = forecaster,
                                   y          = data.users,
                                   cv         = cv,
                                   metric     = 'mean_absolute_error',
                                   verbose    = False
                               )
metric

Excluir parte de la serie¶

Para minimizar la influencia en el modelo de estas fechas, se crea una función que asigna pesos siguiendo las reglas:

Peso de 0 si la fecha del índice es:
- Dentro del periodo de confinamiento (2020-03-16 a 2020-04-21).
- Dentro del periodo de recuperación (2020-04-21 a 2020-05-31).
- 21 días después del periodo de recuperación para evitar incluir valores impactados como lags (2020-05-31 a 2020-06-21).
Peso de 1 en otro caso.

Si una observación tiene un peso de 0, no tiene influencia durante el entrenamiento del modelo.

# Función de para crear pesos
# ==============================================================================
def custom_weights(index):
    """
    Devuelve 0 si el índice está entre 2020-03-16 y 2020-06-21. Devuelve 1 en caso contrario.
    """
    weights = np.where((index >= '2020-03-16') & (index <= '2020-06-21'), 0, 1)
    
    return weights

Se crea de nuevo un ForecasterAutoreg pero esta vez incluyendo la función custom_weights.

# Crear un objeto ForecasterAutoreg
# ==============================================================================
forecaster = ForecasterRecursive(
                 regressor   = Ridge(random_state=123),
                 lags        = 21,
                 weight_func = custom_weights
             )

# Backtesting
# ==============================================================================
metric, predictions_backtest = backtesting_forecaster(
                                   forecaster = forecaster,
                                   y          = data.users,
                                   cv         = cv,
                                   metric     = 'mean_absolute_error',
                                   verbose    = False
                               )

metric

Dando un peso de 0 al periodo de confinamiento (excluyéndolo del entrenamiento del modelo) mejora ligeramente el rendimiento del forecasting.

Parada de central eléctrica¶

Las central eléctricas son instalaciones muy complejas que requieren un alto nivel de mantenimiento. Es común que, cada cierto periodo de funcionamiento, la planta tenga que ser detenida para actividades de reparación, revisión o mantenimiento. Estos eventos provocan que la producción se detenga, generando un hueco en los datos.

Datos¶

# Descarga de datos
# ==============================================================================
url = ('https://raw.githubusercontent.com/JoaquinAmatRodrigo/skforecast/master/'
       'data/energy_production_shutdown.csv')
data = pd.read_csv(url, sep=',')

# Preprocesado de los datos
# ==============================================================================
data['date'] = pd.to_datetime(data['date'], format='%Y-%m-%d')
data = data.set_index('date')
data = data.asfreq('D')
data = data.sort_index()
data.head()

# División datos entrenamiento, validation y test
# ==============================================================================
data = data.loc['2012-01-01 00:00:00': '2014-12-30 23:00:00']
end_train = '2013-12-31 23:59:00'
data_train = data.loc[: end_train, :]
data_test  = data.loc[end_train:, :]

print(f"Dates train : {data_train.index.min()} --- {data_train.index.max()}  (n={len(data_train)})")
print(f"Dates test  : {data_test.index.min()} --- {data_test.index.max()}  (n={len(data_test)})")

Dates train : 2012-01-01 00:00:00 --- 2013-12-31 00:00:00  (n=731)
Dates test  : 2014-01-01 00:00:00 --- 2014-12-30 00:00:00  (n=364)

# Time series plot
# ==============================================================================
fig, ax = plt.subplots(figsize=(8, 3))
data_train.production.plot(ax=ax, label='train', linewidth=1)
data_test.production.plot(ax=ax, label='test', linewidth=1)
ax.axvspan(
    pd.to_datetime('2012-06-01'),
    pd.to_datetime('2012-09-30'), 
    label="Parada",
    color="red",
    alpha=0.1
)
ax.set_title('Producción de energía')
ax.legend();

Modelado con toda la serie¶

Se entrena un modelo incluyendo el periodo de parada de la central eléctrica.

# Crear un ForecasterAutoreg
# ==============================================================================
forecaster = ForecasterRecursive(
                 regressor = Ridge(random_state=123),
                 lags      = 21,
             )

# Backtesting
# ==============================================================================
cv = TimeSeriesFold(
        steps              = 10,
        initial_train_size = len(data.loc[:end_train]),
        refit              = False,    
)
metric, predictions_backtest = backtesting_forecaster(
                                   forecaster = forecaster,
                                   y          = data.production,
                                   cv         = cv,
                                   metric     = 'mean_absolute_error',
                                   verbose    = False
                              )
metric

Excluir parte de la serie¶

La parada de la fábrica tuvo lugar desde 2012-06-01 hasta 2012-09-30. Para minimizar la influencia en el modelo de estas fechas, se crea una función personalizada que asigna un valor de 0 si la fecha del índice está dentro del periodo de parada o 21 días después (lags utilizados por el modelo) y 1 en otro caso. Si una observación tiene un peso de 0, no tiene influencia durante el entrenamiento del modelo.

# Función de pesos
# ==============================================================================
def custom_weights(index):
    """
    Devuelve 0 si el índice está entre 2012-06-01 y 2012-10-21. Devuelve 1 en caso contrario.
    """
    weights = np.where((index >= '2012-06-01') & (index <= '2012-10-21'), 0, 1)
    return weights

# Crear un ForecasterAutoreg
# ==============================================================================
forecaster = ForecasterRecursive(
                 regressor   = Ridge(random_state=123),
                 lags        = 21,
                 weight_func = custom_weights
             )

# Backtesting
# ==============================================================================
metric, predictions_backtest = backtesting_forecaster(
                                   forecaster = forecaster,
                                   y          = data.production,
                                   cv         = cv,
                                   metric     = 'mean_absolute_error',
                                   verbose    = False
                               )
metric

Como en el ejemplo anterior, excluir las observaciones durante el periodo de parada mejora ligeramente el rendimiento del forecasting.

Información de sesión¶

import session_info
session_info.show(html=False)

-----
matplotlib          3.9.2
numpy               2.0.2
pandas              2.2.3
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 15:58

Instrucciones para citar¶

¿Cómo citar este documento?

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

Reducir el impacto del Covid en modelos de forecasting 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/py45-weighted-time-series-forecasting-es.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} }

¿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! 😊

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.

	users
date
2014-06-23	99
2014-06-24	72
2014-06-25	119
2014-06-26	135
2014-06-27	149

	production
date
2012-01-01	375.1
2012-01-02	474.5
2012-01-03	573.9
2012-01-04	539.5
2012-01-05	445.4