Más sobre forecasting en: cienciadedatos.net

• Modelos ARIMA y SARIMAX con Python
• Forecasting de series temporales con machine learning
• Forecasting de series temporales con gradient boosting: XGBoost, LightGBM y CatBoost
• Forecasting de series temporales con XGBoost
• Forecasting probabilístico
• Forecasting con deep learning
• Forecasting de la demanda energética con machine learning
• Forecasting de visitas web con machine learning
• Forecasting de demanda intermitente
• Modelar series temporales con tendencia utilizando modelos de árboles
• Predicción del precio de Bitcoin con Python
• Stacking ensemble de modelos de machine learning para mejorar el forecasting
• Modelos de forecasting interpretables
• Reducir el impacto del Covid en modelos de forecasting
• Forecasting de series temporales con valores ausentes

Introducción¶

En escenarios que implican la predicción de cientos o miles de series temporales, surge una decisión crucial: ¿se deben desarrollar modelos individuales para cada serie o se debe utilizar un único modelo para manejarlas todas a la vez?

En la modelación de una sola serie (modelo de forecasting local), se crea un modelo de predicción independiente para cada serie temporal. Aunque este método proporciona una comprensión exhaustiva de cada serie, su escalabilidad puede verse dificultada por la necesidad de crear y mantener cientos o miles de modelos.

La modelización multiserie (modelo de forecasting global) consiste en crear un único modelo predictivo que tenga en cuenta todas las series temporales simultáneamente. Intenta captar los patrones básicos que rigen las series, mitigando así el ruido potencial que pueda introducir cada serie. Este enfoque es eficiente desde el punto de vista computacional, fácil de mantener y puede producir generalizaciones más sólidas, aunque potencialmente a costa de sacrificar algunos conocimientos individuales.

Este documento muestra cómo predecir más de 1,000 series temporales con un único modelo que incluye características exógenas, algunas de las cuales tienen valores diferentes en cada serie.

Librerías¶

# Data management
# ==============================================================================
import numpy as np
import pandas as pd

# Plots
# ==============================================================================
import matplotlib.pyplot as plt
plt.style.use('seaborn-v0_8-darkgrid')

# Forecasting
# ==============================================================================
import skforecast
import lightgbm
from lightgbm import LGBMRegressor
from sklearn.feature_selection import RFECV
from skforecast.recursive import ForecasterRecursiveMultiSeries
from skforecast.model_selection import (
     TimeSeriesFold,
     OneStepAheadFold,
     backtesting_forecaster_multiseries,
     bayesian_search_forecaster_multiseries
)
from skforecast.feature_selection import select_features_multiseries
from skforecast.preprocessing import (
     CalendarFeatures,
     RollingFeatures,
     reshape_series_long_to_dict,
     reshape_exog_long_to_dict
)
from feature_engine.timeseries.forecasting import WindowFeatures
from skforecast.datasets import fetch_dataset

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

color = '\033[1m\033[38;5;208m' 
print(f'{color}Versión skforecast: {skforecast.__version__}')
print(f'{color}Versión lightgbm: {lightgbm.__version__}')

Versión skforecast: 0.25.0
Versión lightgbm: 4.7.0

Datos¶

Los datos utilizados en este documento se han obtenido del proyecto The Building Data Genome Project 2. El conjunto de datos contiene información sobre el consumo energético de más de 1500 edificios. El rango temporal de los datos de las series temporales abarca los dos años completos (2016 y 2017) y la frecuencia es de mediciones horarias de electricidad, agua de calefacción y refrigeración, vapor y contadores de riego. Además, el conjunto de datos incluye información sobre las características de los edificios y las condiciones meteorológicas. Los datos se han agregado a una resolución diaria y solo se ha considerado la electricidad entre las distintas fuentes de energía.

# Descarga de datos
# ==============================================================================
data = fetch_dataset(name='bdg2_daily')
print('Data shape:', data.shape)
data.head(3)

╭─────────────────────────────────── bdg2_daily ───────────────────────────────────╮
│ Description:                                                                     │
│ Daily energy consumption data from the The Building Data Genome Project 2 with   │
│ building metadata and weather data. https://github.com/buds-lab/building-data-   │
│ genome-project-2                                                                 │
│                                                                                  │
│ Source:                                                                          │
│ Miller, C., Kathirgamanathan, A., Picchetti, B. et al. The Building Data Genome  │
│ Project 2, energy meter data from the ASHRAE Great Energy Predictor III          │
│ competition. Sci Data 7, 368 (2020). https://doi.org/10.1038/s41597-020-00712-x  │
│                                                                                  │
│ URL:                                                                             │
│ https://huggingface.co/datasets/skforecast/bdg2_daily/resolve/main/bdg2_daily.pa │
│ rquet                                                                            │
│                                                                                  │
│ Shape: 1153518 rows x 17 columns                                                 │
╰──────────────────────────────────────────────────────────────────────────────────╯

Data shape: (1153518, 17)

# Rango de fechas disponibles
# ==============================================================================
print(
    f'Rango de fechas disponibles : {data.index.min()} --- {data.index.max()}  '
    f'(n_días={(data.index.max() - data.index.min()).days})'
)

Rango de fechas disponibles : 2016-01-01 00:00:00 --- 2017-12-31 00:00:00  (n_días=730)

# Asegurar que el índice de las series esté completo sin huecos intermedios
# ==============================================================================
data = (
    data
    .groupby('building_id')
    .apply(lambda group: group.asfreq('D', fill_value=np.nan), include_groups=False)
    .reset_index(level=0)
)
print('Data shape:', data.shape)

Data shape: (1153518, 17)

# Rango de fechas disponibles por serie
# ==============================================================================
available_dates_per_series = (
    data
    .dropna(subset='meter_reading')
    .reset_index()
    .groupby('building_id')
    .agg(
        min_index=('timestamp', 'min'),
        max_index=('timestamp', 'max'),
        n_values=('timestamp', 'nunique')
    )
)
display(available_dates_per_series)
print(f'Longitudes de las series : {available_dates_per_series.n_values.unique()}')

Longitudes de las series : [731]

Todas las series temporales tienen la misma longitud, comenzando el 1 de enero de 2016 y terminando el 31 de diciembre de 2017. Algunas variables exógenas contienen valores faltantes. Skforecast no requiere que las series temporales tengan la misma longitud, y se permiten valores faltantes siempre que el regresor subyacente pueda manejarlos, que es el caso de LightGBM, XGBoost y HistGradientBoostingRegressor.

# Valores nulos por variable
# ==============================================================================
data.isna().mean().mul(100).round(2)

building_id               0.00
meter_reading             0.00
site_id                   0.00
primaryspaceusage         1.20
sub_primaryspaceusage     1.20
sqm                       0.00
lat                      14.83
lng                      14.83
timezone                  0.00
airTemperature            0.02
cloudCoverage             7.02
dewTemperature            0.03
precipDepth1HR            0.02
precipDepth6HR            0.02
seaLvlPressure            9.56
windDirection             0.02
windSpeed                 0.02
dtype: float64

Variables exógenas¶

Las variables exógenas son variables externas a la serie temporal y pueden utilizarse como predictores para mejorar la predicción. En este caso, las variables exógenas utilizadas son: las características de los edificios, variables de calendario y condiciones meteorológicas.

Características de los edificios¶

Uno de los atributos clave asociados a cada edificio es su uso designado. Este atributo puede desempeñar un papel crucial en el patrón de consumo de energía.

# Número de edificios, tipos y subtipos
# ==============================================================================
print(f"Número de edificios: {data['building_id'].nunique()}")
print(f"Número de tipos de edificios: {data['primaryspaceusage'].nunique()}")
print(f"Número de subtipos de edificios: {data['sub_primaryspaceusage'].nunique()}")

Número de edificios: 1578
Número de tipos de edificios: 16
Número de subtipos de edificios: 104

Algunos tipos y subtipos de edificios aparecen con poca frecuencia en el conjunto de datos. Tipos con menos de 100 edificios y subtipos con menos de 50 edificios se agrupan en la categoría "Other".

# Agregación de categorías infrecuentes
# ==============================================================================
infrequent_types = (
    data
    .drop_duplicates(subset=['building_id'])['primaryspaceusage']
    .value_counts()
    .loc[lambda x: x < 100]
    .index
    .tolist()
)
infrequent_subtypes = (
    data
    .drop_duplicates(subset=['building_id'])['sub_primaryspaceusage']
    .value_counts()
    .loc[lambda x: x < 50]
    .index
    .tolist()
)

data['primaryspaceusage'] = np.where(
    data['primaryspaceusage'].isin(infrequent_types),
    'Other',
    data['primaryspaceusage']
)
data['sub_primaryspaceusage'] = np.where(
    data['sub_primaryspaceusage'].isin(infrequent_subtypes),
    'Other',
    data['sub_primaryspaceusage']
)

display(data.drop_duplicates(subset=['building_id'])['primaryspaceusage'].value_counts(dropna=False))
display(data.drop_duplicates(subset=['building_id'])['sub_primaryspaceusage'].value_counts(dropna=False))

primaryspaceusage
Education                        604
Office                           296
Entertainment/public assembly    203
Public services                  166
Lodging/residential              149
Other                            141
NaN                               19
Name: count, dtype: int64

sub_primaryspaceusage
Other                          612
Office                         295
College Classroom              131
College Laboratory             116
K-12 School                    109
Dormitory                       91
Primary/Secondary Classroom     84
Education                       67
Library                         54
NaN                             19
Name: count, dtype: int64

Variables de calendario¶

# Variables de calendario
# ==============================================================================
features_to_extract = [
    'month',
    'week',
    'day_of_week',
]
calendar_transformer = CalendarFeatures(
                            features = features_to_extract,
                            encoding = 'cyclical',
                            keep_original_columns = True,
                       )
data = calendar_transformer.fit_transform(data)
print(f'Data shape: {data.shape}')
data.head(3)

Data shape: (1153518, 23)

Variables meteorológicas¶

Las variables meteorológicas se han registrado a nivel de localidad, lo que significa que los datos meteorológicos varían según la ubicación del edificio, incluso en el mismo instante de tiempo. En otras palabras, aunque las variables exógenas son consistentes en todas las series, sus valores difieren por ubicación.

# Valores meteorológicos para cada ubicación para una fecha dada
# ==============================================================================
data.loc['2016-01-01'].groupby('site_id', observed=True).agg(
    {
        'airTemperature': 'first',
        'cloudCoverage': 'first',
        'dewTemperature': 'first',
        'precipDepth1HR': 'first',
        'precipDepth6HR': 'first',
        'seaLvlPressure': 'first',
        'windDirection': 'first',
        'windSpeed': 'first',
    }
)

Skforecast permite incluir variables exógenas distintas y/o valores distintos para cada serie dentro del conjunto de datos (detalles proporcionados en la siguiente sección).

Rolling features de variables exógenas¶

# Rolling features
# ==============================================================================
wf_transformer = WindowFeatures(
    variables      = ['airTemperature', 'windSpeed'],
    window         = ['7D', '14D'],
    functions      = ['mean'],
    freq           = 'D',
    missing_values = 'ignore',
    drop_na        = False,
)
data = data.groupby('building_id').apply(wf_transformer.fit_transform, include_groups=False).reset_index(level=0)
print(f'Data shape: {data.shape}')
data.head(3)

Data shape: (1153518, 27)

Variables categóricas¶

Desde la versión 0.22.0, skforecast proporciona un parámetro categorical_features que maneja automáticamente la codificación y configura de forma nativa los estimadores de gradient boosting (XGBoost, LightGBM, CatBoost e HistGradientBoosting), sin necesidad de pipelines de codificación manuales o parámetros específicos del estimador.

Modelado y predicción¶

ForecasterRecursiveMultiSeries permite modelar series temporales de diferentes longitudes y con distintas variables exógenas.

Cuando las series tienen longitudes diferentes, los datos deben transformarse en un diccionario o en un DataFrame con MultiIndex. Si es un diccionario, las claves corresponden a los nombres de las series y los valores son las propias series. Si es un DataFrame con MultiIndex, el primer nivel del índice es el nombre de la serie y el segundo nivel es el índice temporal de tipo DateTime.

Lo mismo aplica para las variables exógenas cuando son distintas para cada serie. Estas variables deben transformarse en un diccionario o en un DataFrame con MultiIndex, donde las claves del diccionario son los nombres de las series y los valores son las variables exógenas correspondientes.

Skforecast facilita esta transformación de los datos a la estructura requerida mediante las funciones reshape_series_long_to_dict y reshape_exog_long_to_dict.

Cuando todas las series tienen la misma longitud y las mismas variables exógenas, no es necesario utilizar diccionarios. En ese caso, las series pueden pasarse como un único DataFrame con cada serie en una columna, y las variables exógenas como un DataFrame con la misma longitud que las series.

# Variables exógenas para el modelo
# ==============================================================================
exog_features = [
    'primaryspaceusage',
    'sub_primaryspaceusage',
    'timezone',
    'sqm',
    'airTemperature',
    'cloudCoverage',
    'dewTemperature',
    'precipDepth1HR',
    'precipDepth6HR',
    'seaLvlPressure',
    'windDirection',
    'windSpeed',
    'day_of_week_sin',
    'day_of_week_cos',
    'week_sin',
    'week_cos',
    'month_sin',
    'month_cos',
    'airTemperature_window_7D_mean',
    'windSpeed_window_7D_mean',
    'airTemperature_window_14D_mean',
    'windSpeed_window_14D_mean',
]

# Transformación de las series y variables exógenas a formato dict
# ==============================================================================
series_dict = reshape_series_long_to_dict(
    data      = data.reset_index(),
    series_id = 'building_id',
    index     = 'timestamp',
    values    = 'meter_reading',
    freq      = 'D'
)

exog_dict = reshape_exog_long_to_dict(
    data      = data[exog_features + ['building_id']].reset_index(),
    series_id = 'building_id',
    index     = 'timestamp',
    freq      = 'D'
)

Para entrenar los modelos, buscar los hiperparámetros óptimos y evaluar su rendimiento predictivo, los datos se dividen en tres conjuntos separados: entrenamiento, validación y test.

# Partición de los datos en entrenamiento, validación y test
# ==============================================================================
data = data.sort_index()
end_train = '2017-08-31 23:59:00'
end_validation = '2017-10-31 23:59:00'
# Solo se guardan las particiones que se utilizan más adelante en el documento
series_dict_train = {k: v.loc[:end_train] for k, v in series_dict.items()}
exog_dict_train = {k: v.loc[:end_train] for k, v in exog_dict.items()}
exog_dict_valid = {k: v.loc[end_train:end_validation] for k, v in exog_dict.items()}

print(
    f'Rango de fechas disponibles : {data.index.min()} --- {data.index.max()} '
    f'(n_días={(data.index.max() - data.index.min()).days})'
)
print(
    f'  Fechas de entrenamiento : {data.loc[: end_train, :].index.min()} --- {data.loc[: end_train, :].index.max()} '
    f'(n_días={(data.loc[: end_train, :].index.max() - data.loc[: end_train, :].index.min()).days})'
)
print(
    f'  Fechas de validación    : {data.loc[end_train:end_validation, :].index.min()} --- {data.loc[end_train:end_validation, :].index.max()} '
    f'(n_días={(data.loc[end_train:end_validation, :].index.max() - data.loc[end_train:end_validation, :].index.min()).days})'
)
print(
    f'  Fechas de test          : {data.loc[end_validation:, :].index.min()} --- {data.loc[end_validation:, :].index.max()} '
    f'(n_días={(data.loc[end_validation:, :].index.max() - data.loc[end_validation:, :].index.min()).days})'
)

Rango de fechas disponibles : 2016-01-01 00:00:00 --- 2017-12-31 00:00:00 (n_días=730)
  Fechas de entrenamiento : 2016-01-01 00:00:00 --- 2017-08-31 00:00:00 (n_días=608)
  Fechas de validación    : 2017-09-01 00:00:00 --- 2017-10-31 00:00:00 (n_días=60)
  Fechas de test          : 2017-11-01 00:00:00 --- 2017-12-31 00:00:00 (n_días=60)

Búsqueda de hiperparámetros¶

La búsqueda de hiperparámetros y lags implica probar sistemáticamente diferentes valores de hiperparámetros (y/o lags) para encontrar la configuración óptima que ofrezca el mejor rendimiento. skforecast proporciona dos métodos diferentes para evaluar cada configuración candidata, implementados mediante las clases TimeSeriesFold y OneStepAheadFold:

• Backtesting (TimeSeriesFold): en este método, el modelo predice varios pasos a futuro en cada iteración, utilizando el mismo horizonte de predicción y la misma estrategia de reentrenamiento que se utilizarían si se desplegara el modelo (incluyendo la opción de no reentrenar en absoluto, es decir, refit=False, utilizado en los ejemplos siguientes por rapidez). De este modo, se simula un escenario de predicción real en el que el modelo se reentrena y actualiza a lo largo del tiempo.

• One-Step Ahead (OneStepAheadFold): Evalúa el modelo utilizando solo predicciones de un paso a futuro. Este método es más rápido porque requiere menos iteraciones, pero solo evalúa el rendimiento del modelo en el siguiente paso temporal (t+1).

Cada método utiliza una estrategia de evaluación diferente, por lo que pueden producir resultados distintos. Sin embargo, a largo plazo, se espera que ambos métodos converjan a selecciones similares de hiperparámetros óptimos. Se recomienda realizar un backtest del modelo final para obtener una estimación más precisa del rendimiento cuando se predicen varios pasos a futuro.

En este ejemplo, la búsqueda se realiza con la función bayesian_search_forecaster_multiseries, que utiliza Optuna para explorar el espacio de hiperparámetros de manera eficiente. Para mantener un tiempo de ejecución reducido en este ejemplo, solo se evalúan 10 combinaciones (n_trials=10); en una aplicación real, se recomienda un número mucho mayor de combinaciones para explorar adecuadamente el espacio de búsqueda.

# Crear forecaster
# ==============================================================================
window_features = RollingFeatures(stats=['mean', 'min', 'max'], window_sizes=7)
forecaster = ForecasterRecursiveMultiSeries(
                estimator            = LGBMRegressor(random_state=8520, verbose=-1),
                lags                 = 14,
                window_features      = window_features,
                categorical_features = 'auto',
                encoding             = 'ordinal'
            )

# Bayesian search con OneStepAheadFold
# ==============================================================================
def search_space(trial):
    search_space  = {
        'lags'             : trial.suggest_categorical('lags', [31, 62]),
        'n_estimators'     : trial.suggest_int('n_estimators', 200, 800, step=100),
        'learning_rate'    : trial.suggest_float('learning_rate', 0.01, 0.3, log=True),
        'num_leaves'       : trial.suggest_int('num_leaves', 15, 255, log=True),
        'max_depth'        : trial.suggest_int('max_depth', 3, 8),
        'min_child_samples': trial.suggest_int('min_child_samples', 5, 200, log=True),
        'colsample_bytree' : trial.suggest_float('colsample_bytree', 0.6, 1.0),
        'subsample'        : trial.suggest_float('subsample', 0.6, 1.0),
        'subsample_freq'   : trial.suggest_int('subsample_freq', 0, 5),
        'reg_alpha'        : trial.suggest_float('reg_alpha', 1e-3, 10.0, log=True),
        'reg_lambda'       : trial.suggest_float('reg_lambda', 1e-3, 10.0, log=True),
    }

    return search_space

cv = OneStepAheadFold(initial_train_size='2017-08-31')  # última fecha incluida en el conjunto de entrenamiento

results_search, best_trial = bayesian_search_forecaster_multiseries(
    forecaster        = forecaster,
    series            = {k: v.loc[:end_validation] for k, v in series_dict.items()},
    exog              = {k: v.loc[:end_validation, exog_features] for k, v in exog_dict.items()},
    cv                = cv,
    search_space      = search_space,
    n_trials          = 10,
    metric            = 'mean_absolute_error',
    random_state      = 123,
    suppress_warnings = True
)

best_params = results_search.at[0, 'params']
best_lags = results_search.at[0, 'lags']
print(f'Mejores lags: {best_lags}')
print(f'Mejores parámetros: {best_params}')
results_search.head(3)

Mejores 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]
Mejores parámetros: {'n_estimators': 300, 'learning_rate': 0.032160750764041963, 'num_leaves': 63, 'max_depth': 6, 'min_child_samples': 7, 'colsample_bytree': 0.6523579802656323, 'subsample': 0.7287922425873232, 'subsample_freq': 3, 'reg_alpha': 2.4323434680180993, 'reg_lambda': 0.16331624244561055}

Backtesting con datos de test¶

Una vez seleccionada la mejor combinación de lags e hiperparámetros, se evalúa la capacidad predictiva del modelo en el conjunto de test (noviembre y diciembre de 2017) utilizando la función backtesting_forecaster_multiseries. En cada fold, el forecaster predice 7 días cada vez sin reentrenar (refit=False). Dado que la búsqueda se ha ejecutado con el valor por defecto return_best=True, el forecaster se actualiza automáticamente con la mejor configuración encontrada, por lo que puede pasarse directamente a la función de backtesting.

# Backtesting
# ==============================================================================
cv = TimeSeriesFold(
    initial_train_size = '2017-10-31',  # última fecha de entrenamiento + validación
    steps              = 7,
    refit              = False
)
metrics, predictions = backtesting_forecaster_multiseries(
                          forecaster        = forecaster,
                          series            = series_dict,
                          exog              = exog_dict,
                          cv                = cv,
                          metric            = 'mean_absolute_error',
                          suppress_warnings = True
                      )

display(predictions.head())
display(metrics)

# Agregación de métricas para todos los edificios
# ==============================================================================
average_metric_all_buildings = metrics.query("levels == 'average'")['mean_absolute_error'].item()
errors_all_buildings = pd.merge(
    left     = data[['building_id', 'meter_reading']].reset_index(),
    right    = predictions.rename_axis('timestamp').reset_index(),
    left_on  = ['timestamp', 'building_id'],
    right_on = ['timestamp', 'level'],
    how      = 'inner',
    validate = '1:1'
).assign(error=lambda df: df['meter_reading'] - df['pred'])

sum_abs_errors_all_buildings = errors_all_buildings['error'].abs().sum()
sum_bias_all_buildings = errors_all_buildings['error'].sum()
print(f'Error absoluto medio promedio para todos los edificios: {average_metric_all_buildings:.0f}')
print(f'Suma de errores absolutos para todos los edificios (x 10,000): {sum_abs_errors_all_buildings/10000:.0f}')
print(f'Sesgo (x 10,000): {sum_bias_all_buildings/10000:.0f}')

Error absoluto medio promedio para todos los edificios: 311
Suma de errores absolutos para todos los edificios (x 10,000): 2990
Sesgo (x 10,000): -316

# Gráfico de predicciones vs reales para dos edificios aleatorios
# ==============================================================================
rng = np.random.default_rng(14793)
n_buildings = 2
selected_buildings = rng.choice(data['building_id'].unique(), size=n_buildings, replace=False)

fig, axs = plt.subplots(n_buildings, 1, figsize=(7, 4.5), sharex=True)
axs = axs.flatten()

for i, building in enumerate(selected_buildings):
    data.query('building_id == @building').loc[predictions.index.unique(), 'meter_reading'].plot(ax=axs[i], label='test')
    predictions.query('level == @building')['pred'].plot(ax=axs[i], label='predictions')
    axs[i].set_title(f'Building {building}', fontsize=10)
    axs[i].set_xlabel('')
    axs[i].legend()

fig.tight_layout()
plt.show();

Selección de predictores¶

La selección de predictores es el proceso de seleccionar un subconjunto de predictores relevantes (variables) para su uso en la construcción del modelo. Las técnicas de selección de predictores se utilizan por varias razones: para simplificar los modelos y hacerlos más fáciles de interpretar, para reducir el tiempo de entrenamiento, para evitar los problemas de dimensionalidad, para mejorar la generalización reduciendo el sobreajuste (formalmente, la reducción de la varianza), entre otros.

Skforecast es compatible con los métodos de selección implementados en scikit-learn. Existen varios métodos de selección de características, pero los más comunes son:

• Recursive feature elimination (RFE)

• Sequential Feature Selection (SFS)

• Feature selection based on threshold (SelectFromModel)

Skforecast proporciona la función select_features_multiseries para aplicar cualquiera de estos selectores de scikit-learn directamente sobre los lags, window features, variables exógenas y variables de calendario de un ForecasterRecursiveMultiSeries.

# Selección de predictores
# ==============================================================================
warnings.filterwarnings('ignore', message='X does not have valid feature names.*')
estimator = LGBMRegressor(n_estimators=50, num_leaves=8, max_bin=63, random_state=15926, verbose=-1, n_jobs=-1)
selector = RFECV(estimator=estimator, step=3, cv=3, n_jobs=1)
selected_lags, selected_window_features, selected_exog, _ = select_features_multiseries(
    forecaster      = forecaster,
    selector        = selector,
    series          = {k: v.loc[:end_validation] for k, v in series_dict.items()},
    exog            = {k: v.loc[:end_validation, exog_features] for k, v in exog_dict.items()},
    select_only     = None,
    force_inclusion = None,
    subsample       = 0.2,
    random_state    = 123,
    verbose         = True,
)

╭──────────────────────────────── MissingValuesWarning ────────────────────────────────╮
│ NaNs detected in `X_train`. Some estimators do not allow NaN values during training. │
│ If you want to drop them, set `forecaster.dropna_from_series = True`.                │
│                                                                                      │
│ Category : skforecast.exceptions.MissingValuesWarning                                │
│ Location :                                                                           │
│ C:\Users\Joaquin\miniconda3\envs\skforecast_24_py13\Lib\site-packages\skforecast\rec │
│ ursive\_forecaster_recursive_multiseries.py:1456                                     │
│ Suppress : warnings.simplefilter('ignore', category=MissingValuesWarning)            │
╰──────────────────────────────────────────────────────────────────────────────────────╯

Recursive feature elimination (RFECV)
-------------------------------------
Total number of records available: 959424
Total number of records used for feature selection: 191884
Number of features available: 87
    Lags            (n=62)
    Window features (n=3)
    Exog            (n=22)
    Calendar        (n=0)
Number of features selected: 45
    Lags            (n=28) : [1, 5, 7, 8, 9, 10, 13, 14, 15, 20, 21, 23, 26, 28, 34, 35, 37, 42, 48, 49, 54, 56, 57, 58, 59, 60, 61, 62]
    Window features (n=3) : ['roll_mean_7', 'roll_min_7', 'roll_max_7']
    Exog            (n=14) : ['primaryspaceusage', 'sub_primaryspaceusage', 'timezone', 'sqm', 'airTemperature', 'cloudCoverage', 'dewTemperature', 'windDirection', 'day_of_week_sin', 'day_of_week_cos', 'airTemperature_window_7D_mean', 'windSpeed_window_7D_mean', 'airTemperature_window_14D_mean', 'windSpeed_window_14D_mean']
    Calendar        (n=0) : []

# Backtesting del forecaster con predictores seleccionados
# ==============================================================================
selected_stats = [feature.split('_')[1] for feature in selected_window_features]
window_features = RollingFeatures(stats=selected_stats, window_sizes=7)
forecaster = ForecasterRecursiveMultiSeries(
                estimator            = LGBMRegressor(**best_params, random_state=8520, verbose=-1),
                lags                 = selected_lags,
                window_features      = window_features,
                categorical_features = 'auto',
                encoding             = 'ordinal'
            )
cv = TimeSeriesFold(
    initial_train_size = '2017-10-31',  # última fecha de entrenamiento + validación
    steps              = 7,
    refit              = False
)
metrics, predictions = backtesting_forecaster_multiseries(
                          forecaster        = forecaster,
                          series            = series_dict,
                          exog              = {k: v[selected_exog] for k, v in exog_dict.items()},
                          cv                = cv,
                          metric            = 'mean_absolute_error',
                          suppress_warnings = True
                      )

display(predictions.head())
display(metrics)

Se ha conseguido reducir el número de predictores sin que el rendimiento del modelo se vea comprometido. Esto permite simplificar el modelo y acelera el entrenamiento.

Las variables de calendario que aparecen como Calendar (n=0) en la salida de selección de predictores no se evalúan aquí como su propio grupo porque se pasaron como variables exógenas regulares en lugar de a través del parámetro calendar_features del forecaster; sus columnas cíclicas se evalúan como parte del conjunto de variables exógenas.

Predecir nuevas series (series desconocidas)¶

ForecasterRecursiveMultiSeries permite predecir series desconocidas (series no vistas durante el proceso de entrenamiento). Pueden darse dos escenarios:

• Existen datos históricos para la nueva serie:
  El usuario debe proporcionar un DataFrame con los datos históricos necesarios para crear los predictores (lags y rolling) utilizados durante el proceso de predicción. Este DataFrame debe contener una columna por cada serie que se desee predecir. Es posible incluir tanto series vistas durante el entrenamiento como series nuevas (no vistas).

• No existen datos históricos para la nueva serie:
  Este escenario es similar, pero la last_window para las nuevas series está compuesta completamente por valores NaN. Este enfoque solo es posible cuando se utilizan estimadores que soportan valores faltantes, ya que todos los lags serán NaN. Si tampoco se dispone de variables exógenas, todos los inputs del modelo son NaN, por lo que el modelo devuelve la misma predicción para cada paso y serie, como se muestra en el segundo ejemplo más abajo.

En ambos casos, las variables exógenas pueden estar disponibles o no. Si no están disponibles, se establecerán automáticamente como NaN.

# Predecir nuevas series con datos históricos disponibles
# ==============================================================================
# Crear un dataframe con las últimas observaciones de la nueva serie
forecaster.fit(
    series = series_dict_train,
    exog   = {k: v[selected_exog] for k, v in exog_dict_train.items()},
)
rng = np.random.default_rng(67584)
last_window_new = pd.DataFrame(
    {
        'series_new_1': rng.normal(loc=1000, scale=500, size=62),
        'series_new_2': rng.normal(loc=1000, scale=500, size=62),
    },
    index=pd.date_range('2017-07-01', periods=62, freq='D'),
)

# Si las variables exógenas no están disponibles para la nueva serie, se establecen automáticamente como NaN
forecaster.predict(
    steps       = 5,
    last_window = last_window_new,
    exog        = {k: v[selected_exog] for k, v in exog_dict_valid.items()}
)

╭──────────────────────────────── MissingValuesWarning ────────────────────────────────╮
│ NaNs detected in `X_train`. Some estimators do not allow NaN values during training. │
│ If you want to drop them, set `forecaster.dropna_from_series = True`.                │
│                                                                                      │
│ Category : skforecast.exceptions.MissingValuesWarning                                │
│ Location :                                                                           │
│ C:\Users\Joaquin\miniconda3\envs\skforecast_24_py13\Lib\site-packages\skforecast\rec │
│ ursive\_forecaster_recursive_multiseries.py:1456                                     │
│ Suppress : warnings.simplefilter('ignore', category=MissingValuesWarning)            │
╰──────────────────────────────────────────────────────────────────────────────────────╯

╭──────────────────────────────── UnknownLevelWarning ─────────────────────────────────╮
│ `levels` {'series_new_2', 'series_new_1'} were not included in training. Unknown     │
│ levels are encoded as NaN, which may cause the prediction to fail if the estimator   │
│ does not accept NaN values.                                                          │
│                                                                                      │
│ Category : skforecast.exceptions.UnknownLevelWarning                                 │
│ Location :                                                                           │
│ C:\Users\Joaquin\miniconda3\envs\skforecast_24_py13\Lib\site-packages\skforecast\uti │
│ ls\utils.py:1258                                                                     │
│ Suppress : warnings.simplefilter('ignore', category=UnknownLevelWarning)             │
╰──────────────────────────────────────────────────────────────────────────────────────╯

╭───────────────────────────────── MissingExogWarning ─────────────────────────────────╮
│ `exog` does not contain keys for levels {'series_new_2', 'series_new_1'}. Missing    │
│ levels are filled with NaN. Most of machine learning models do not allow missing     │
│ values. Prediction method may fail.                                                  │
│                                                                                      │
│ Category : skforecast.exceptions.MissingExogWarning                                  │
│ Location :                                                                           │
│ C:\Users\Joaquin\miniconda3\envs\skforecast_24_py13\Lib\site-packages\skforecast\uti │
│ ls\utils.py:1384                                                                     │
│ Suppress : warnings.simplefilter('ignore', category=MissingExogWarning)              │
╰──────────────────────────────────────────────────────────────────────────────────────╯

# Predecir nuevas series sin datos históricos disponibles
# ==============================================================================
# Crear un dataframe con la última ventana de la nueva serie con valores nan
last_window_new = pd.DataFrame(
    {
        'series_new_1': np.nan,
        'series_new_2': np.nan,
    },
    index=pd.date_range('2017-07-01', periods=62, freq='D'),
)

# Si las variables exógenas no están disponibles para la nueva serie, se establecen automáticamente como NaN
forecaster.predict(
    steps       = 5,
    last_window = last_window_new,
    exog        = {k: v[selected_exog] for k, v in exog_dict_valid.items()}
)

╭──────────────────────────────── UnknownLevelWarning ─────────────────────────────────╮
│ `levels` {'series_new_2', 'series_new_1'} were not included in training. Unknown     │
│ levels are encoded as NaN, which may cause the prediction to fail if the estimator   │
│ does not accept NaN values.                                                          │
│                                                                                      │
│ Category : skforecast.exceptions.UnknownLevelWarning                                 │
│ Location :                                                                           │
│ C:\Users\Joaquin\miniconda3\envs\skforecast_24_py13\Lib\site-packages\skforecast\uti │
│ ls\utils.py:1258                                                                     │
│ Suppress : warnings.simplefilter('ignore', category=UnknownLevelWarning)             │
╰──────────────────────────────────────────────────────────────────────────────────────╯

╭──────────────────────────────── MissingValuesWarning ────────────────────────────────╮
│ `last_window` has missing values. Most of machine learning models do not allow       │
│ missing values. Prediction method may either raise an error or return NaN            │
│ predictions.                                                                         │
│                                                                                      │
│ Category : skforecast.exceptions.MissingValuesWarning                                │
│ Location :                                                                           │
│ C:\Users\Joaquin\miniconda3\envs\skforecast_24_py13\Lib\site-packages\skforecast\uti │
│ ls\utils.py:1338                                                                     │
│ Suppress : warnings.simplefilter('ignore', category=MissingValuesWarning)            │
╰──────────────────────────────────────────────────────────────────────────────────────╯

╭───────────────────────────────── MissingExogWarning ─────────────────────────────────╮
│ `exog` does not contain keys for levels {'series_new_2', 'series_new_1'}. Missing    │
│ levels are filled with NaN. Most of machine learning models do not allow missing     │
│ values. Prediction method may fail.                                                  │
│                                                                                      │
│ Category : skforecast.exceptions.MissingExogWarning                                  │
│ Location :                                                                           │
│ C:\Users\Joaquin\miniconda3\envs\skforecast_24_py13\Lib\site-packages\skforecast\uti │
│ ls\utils.py:1384                                                                     │
│ Suppress : warnings.simplefilter('ignore', category=MissingExogWarning)              │
╰──────────────────────────────────────────────────────────────────────────────────────╯

Clustering series temporales¶

La idea que hay detrás de modelar varias series al mismo tiempo es poder capturar los patrones principales que rigen dichas series, reduciendo así el impacto del ruido que pueda haber en cada una de ellas. Esto significa que las series que se comportan de manera similar pueden beneficiarse de ser modelizadas juntas. Una forma de identificar posibles grupos de series es realizar un estudio de clustering antes de modelizarlas. Si como resultado del clustering se identifican grupos claros, es apropiado modelar cada uno de ellos por separado.

El clustering es una técnica de análisis no supervisado que agrupa un conjunto de observaciones en clústeres que contienen observaciones consideradas homogéneas, mientras que las observaciones en diferentes clústeres se consideran heterogéneas. Los algoritmos que agrupan series temporales se pueden dividir en dos grupos: aquellos que utilizan una transformación para crear variables antes de agrupar (clustering de series temporales basado en características) y aquellos que trabajan directamente en las series temporales (medidas de distancia elástica).

• Clustering basado en características de series temporales: se extraen variables que describen las características estructurales de cada serie temporal y luego se introducen en algoritmos de clustering. Estas variables se obtienen aplicando operaciones estadísticas que capturan mejor las características subyacentes: tendencia, estacionalidad, periodicidad, correlación serial, asimetría, curtosis, caos, no linealidad y auto-similitud.

• Medidas de distancia elástica: este enfoque trabaja directamente en las series temporales, ajustando o «reajustando» las series en comparación con otras. La medida más conocida de esta familia es el Dynamic Time Warping (DTW).

Para un ejemplo detallado de cómo el clustering de series temporales puede mejorar los modelos de forecasting, consulte Clustering de series temporales para mejorar modelos de forecasting.

Información de sesión¶

import session_info
session_info.show(html=False)

-----
feature_engine      1.9.4
lightgbm            4.7.0
matplotlib          3.10.9
numpy               2.4.6
optuna              4.9.0
pandas              2.3.3
session_info        v1.0.1
skforecast          0.25.0
sklearn             1.7.2
-----
IPython             9.15.0
jupyter_client      8.9.1
jupyter_core        5.9.1
-----
Python 3.13.14 | packaged by conda-forge | (main, Jun 12 2026, 09:44:26) [MSC v.1944 64 bit (AMD64)]
Windows-11-10.0.26200-SP0
-----
Session information updated at 2026-09-24 16:10

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 escalable: modelado de miles de series temporales con un único modelo global 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/py59-modelos-forecasting-escalables.html

¿Cómo citar skforecast?

Si utilizas skforecast, te agradeceríamos mucho que lo cites. ¡Muchas gracias!

Zenodo:

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

APA:

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

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

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.

• No-Comercial: No puedes utilizar el material para fines comerciales.

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