Forecasting probabilístico: intervalos de predicción para forecasting multi-step

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

# Gráficos
# ==============================================================================
import matplotlib.pyplot as plt
from skforecast.plot import set_dark_theme, plot_prediction_intervals, plot_residuals
from statsmodels.graphics.tsaplots import plot_pacf
from statsmodels.tsa.stattools import pacf
import plotly.graph_objects as go
import plotly.io as pio
import plotly.offline as poff
pio.templates.default = "seaborn"
pio.renderers.default = 'notebook' 
poff.init_notebook_mode(connected=True)

# Modelado y forecasting
# ==============================================================================
import skforecast
from sklearn.linear_model import Ridge
from sklearn.feature_selection import RFECV
from sklearn.pipeline import make_pipeline
from feature_engine.datetime import DatetimeFeatures
from feature_engine.creation import CyclicalFeatures
from feature_engine.timeseries.forecasting import LagFeatures, WindowFeatures
from skforecast.recursive import ForecasterRecursive
from skforecast.model_selection import TimeSeriesFold, bayesian_search_forecaster, backtesting_forecaster
from skforecast.feature_selection import select_features
from skforecast.preprocessing import RollingFeatures
from skforecast.metrics import calculate_coverage

# Warnings
# ==============================================================================
import warnings
warnings.filterwarnings('once')
print('Versión skforecast:', skforecast.__version__)

# Descarga de datos
# ==============================================================================
datos = fetch_dataset(name="ett_m2")
datos = datos.resample(rule="1h", closed="left", label="right").mean()
datos

# Gráfico de la serie
# ==============================================================================
set_dark_theme()
plt.rcParams['lines.linewidth'] = 0.5
fig, ax = plt.subplots(figsize=(8, 3))
ax.plot(datos['OT'])
ax.set_title('OT');

# Gráfico de las series exógenas
# ==============================================================================
columnas_grafico = [
    'HUFL', 'HULL', 'MUFL', 'MULL', 'LUFL', 'LULL'
]
colores = plt.rcParams['axes.prop_cycle'].by_key()['color']
colores = colores[1:] + ["violet", "magenta"]    
fig, axs = plt.subplots(6, 1, figsize=(8, 8), sharex=True)
for i, columna in enumerate(columnas_grafico):
    axs[i].plot(datos[columna], label=columna, color=colores[i])
    axs[i].legend(loc='lower left', fontsize=8)
    axs[i].tick_params(axis='both', labelsize=8)

# Autocorrelación parcial
# ==============================================================================
n_lags = 24 * 7
columnas_grafico = ['OT', 'HUFL', 'HULL', 'MUFL', 'MULL', 'LUFL', 'LULL']
pacf_df = []
fig, axs = plt.subplots(4, 2, figsize=(8, 6))
axs = axs.ravel()
for i, columna in enumerate(columnas_grafico):
    valores_pacf = pacf(datos[columna], nlags=n_lags)
    valores_pacf = pd.DataFrame({
        'lag': range(1, len(valores_pacf)),
        'value': valores_pacf[1:],
        'variable': columna
    })
    pacf_df.append(valores_pacf)
    plot_pacf(datos[columna], lags=n_lags, ax=axs[i], method='burg')
    axs[i].set_title(columna, fontsize=10)
    axs[i].set_ylim(-0.5, 1.1)
plt.tight_layout()

# Top n lags con mayor autocorrelación parcial por variable
# ==============================================================================
n = 10
lags_principales = set()
for valores_pacf in pacf_df:
    variable = valores_pacf['variable'].iloc[0]
    lags = valores_pacf.nlargest(n, 'value')['lag'].tolist()
    lags_principales.update(lags)
    print(f"{variable}: {lags}")
lags_principales = list(lags_principales)
print(f"\nTodos los lags: {lags_principales}")

# Variables de calendario
# ==============================================================================
variables_a_extraer = ['year', 'month', 'week', 'day_of_week', 'hour']
calendar_transformer = DatetimeFeatures(
    variables           = 'index',
    features_to_extract = variables_a_extraer,
    drop_original       = False,
)

# Lags de las variables exógenas
# ==============================================================================
lag_transformer = LagFeatures(
    variables = ["HUFL", "MUFL", "MULL", "HULL", "LUFL", "LULL"],
    periods   = lags_principales,
)

# Estadísticas móviles de las variables exógenas
# ==============================================================================
wf_transformer = WindowFeatures(
    variables   = ["HUFL", "MUFL", "MULL", "HULL", "LUFL", "LULL"],
    window      = ["1D", "7D"],
    functions   = ["mean", "max", "min"],
    freq        = "1h",
)

# Codificación de variables cíclicas
# ==============================================================================
variables_a_codificar = ['year', 'month', 'week', 'day_of_week', 'hour']
valores_maximos = {
    "month": 12,
    "week": 52,
    "day_of_week": 7,
    "hour": 24,
}
codificador_ciclico = CyclicalFeatures(
                        variables     = variables_a_codificar,
                        max_values    = valores_maximos,
                        drop_original = True
                   )

exog_transformer = make_pipeline(
                        calendar_transformer,
                        lag_transformer,
                        wf_transformer,
                        codificador_ciclico
                   )
display(exog_transformer)

datos = exog_transformer.fit_transform(datos)
# Eliminación de valores nulos introducidos al crear los lags
datos = datos.dropna()
exog_features = datos.columns.difference(['OT']).tolist()
display(datos.head(3))

# División de los datos en train-validación-test
# ==============================================================================
fin_train = '2017-10-01 23:59:00'
fin_validacion = '2018-04-03 23:59:00'
datos_train = datos.loc[: fin_train, :]
datos_val   = datos.loc[fin_train:fin_validacion, :]
datos_test  = datos.loc[fin_validacion:, :]

print(f"Fechas train      : {datos_train.index.min()} --- {datos_train.index.max()}  (n={len(datos_train)})")
print(f"Fechas validación : {datos_val.index.min()} --- {datos_val.index.max()}  (n={len(datos_val)})")
print(f"Fechas test       : {datos_test.index.min()} --- {datos_test.index.max()}  (n={len(datos_test)})")

# Gráfico de las particiones
# ==============================================================================
fig, ax = plt.subplots(figsize=(8, 3))
ax.plot(datos_train['OT'], label='Train')
ax.plot(datos_val['OT'], label='Validación')
ax.plot(datos_test['OT'], label='Test')
ax.set_title('OT')
ax.legend();

# Gráfico de las particiones tras aplicar diferenciación
# ==============================================================================
fig, ax = plt.subplots(figsize=(8, 3))
ax.plot(datos_train['OT'].diff(1), label='Train')
ax.plot(datos_val['OT'].diff(1), label='Validación')
ax.plot(datos_test['OT'].diff(1), label='Test')
ax.set_title('OT diferenciada')
ax.legend();

# Crear forecaster
# ==============================================================================
window_features = RollingFeatures(stats=['mean', 'min', 'max'], window_sizes=24)
forecaster = ForecasterRecursive(
                 estimator       = Ridge(random_state=15926),
                 lags            = lags_principales,
                 differentiation = 1,
                 window_features = window_features,
                 binner_kwargs   = {'n_bins': 10}
             )

# Selección de predictores (autorregresivos y exógenas) con scikit-learn RFECV
# ==============================================================================
estimator = Ridge(random_state=15926)

selector = RFECV(
    estimator=estimator, step=1, cv=3, min_features_to_select=25, n_jobs=-1
)

lags_seleccionados, wf_seleccionadas, exog_seleccionadas, _ = select_features(
    forecaster      = forecaster,
    selector        = selector,
    y               = datos.loc[:fin_train, 'OT'],
    exog            = datos.loc[:fin_train, exog_features],
    select_only     = None,
    force_inclusion = None,
    subsample       = 0.5,
    random_state    = 123,
    verbose         = True,
)

# Establecer los lags seleccionados
# ==============================================================================
forecaster.set_lags(lags_seleccionados)

# Búsqueda de hiperparámetros
# ==============================================================================
cv = TimeSeriesFold(
        initial_train_size = len(datos.loc[:fin_train, :]),
        steps              = 24,  # todas las horas del día siguiente
     )

def search_space(trial):
    search_space  = {
        'alpha': trial.suggest_float('alpha', 1e-7, 1e+3, log=True)
    }
    return search_space

resultados_busqueda, frozen_trial = bayesian_search_forecaster(
                                   forecaster   = forecaster,
                                   y            = datos.loc[:fin_validacion, 'OT'],
                                   exog         = datos.loc[:fin_validacion, exog_seleccionadas],
                                   search_space = search_space,
                                   cv           = cv,
                                   metric       = 'mean_absolute_error',
                                   n_trials     = 20
                               )
best_params = resultados_busqueda['params'].iloc[0]
best_lags = resultados_busqueda['lags'].iloc[0]
resultados_busqueda.head()

# Backtesting con datos de validación para obtener residuos out-sample
# ==============================================================================
cv = TimeSeriesFold(
        initial_train_size = len(datos.loc[:fin_train, :]),
        steps              = 24,  # todas las horas del día siguiente
     )

metrica_val, predicciones_val = backtesting_forecaster(
                                  forecaster = forecaster,
                                  y          = datos.loc[:fin_validacion, 'OT'],
                                  exog       = datos.loc[:fin_validacion, exog_seleccionadas],
                                  cv         = cv,
                                  metric     = 'mean_absolute_error'
                              )
display(metrica_val)
fig, ax = plt.subplots(figsize=(8, 3))
datos.loc[fin_train:fin_validacion, 'OT'].plot(ax=ax, label='Valor real')
predicciones_val['pred'].plot(ax=ax, label='Predicción')
ax.set_title("Backtesting en los datos de validación")
ax.legend();

# Distribución de los residuos out-sample
# ==============================================================================
residuos = datos.loc[predicciones_val.index, 'OT'] - predicciones_val['pred']
print(pd.Series(np.where(residuos < 0, 'negativo', 'positivo')).value_counts())
plt.rcParams.update({'font.size': 8})
_ = plot_residuals(residuals=residuos, figsize=(7, 4))

# Almacenar residuos out-sample en el forecaster
# ==============================================================================
forecaster.fit(y=datos.loc[:fin_train, 'OT'], exog=datos.loc[:fin_train, exog_seleccionadas])
forecaster.set_out_sample_residuals(
    y_true = datos.loc[predicciones_val.index, 'OT'], 
    y_pred = predicciones_val['pred']
)

# Intervalos de los bins de residuos (condicionados a los valores predichos) y número de residuos almacenados
# ==============================================================================
for k, v in forecaster.binner_intervals_.items():
    print(f"{k}: {v}  n={len(forecaster.out_sample_residuals_by_bin_[k])}")

# Distribución de los residuos out-sample por bin
# ==============================================================================
residuos_out_sample_por_bin_df = pd.DataFrame(
    dict(
        [(k, pd.Series(v))
         for k, v in forecaster.out_sample_residuals_by_bin_.items()]
    )
)

fig, ax = plt.subplots(figsize=(6, 3))
residuos_out_sample_por_bin_df.boxplot(
    ax=ax,
    flierprops=dict(marker='o', markerfacecolor='white', markeredgecolor='black', markersize=4)
)
ax.set_title("Distribución de los residuos por bin", fontsize=12)
ax.set_xlabel("Bin", fontsize=10)
ax.set_ylabel("Residuos", fontsize=10)
plt.show();

# Backtesting con intervalos de predicción en datos de test usando residuos out-sample
# ==============================================================================
cv = TimeSeriesFold(
        initial_train_size = len(datos.loc[:fin_validacion, :]),
        steps              = 24,  # todas las horas del día siguiente
     )

metrica, predicciones = backtesting_forecaster(
                          forecaster              = forecaster,
                          y                       = datos['OT'],
                          exog                    = datos[exog_seleccionadas],
                          cv                      = cv,
                          metric                  = 'mean_absolute_error',
                          interval                = [0.1, 0.9],  # 80% intervalo de predicción
                          n_boot                  = 150,
                          interval_method         = 'bootstrapping',
                          use_in_sample_residuals = False,  # Usar residuos out-sample 
                          use_binned_residuals    = True,   # Intervalos condicionados a los valores predichos
                      )
predicciones

# Gráfico de los intervalos de predicción
# ==============================================================================
plt.rcParams['lines.linewidth'] = 1
fig, ax = plt.subplots(figsize=(9, 4))
plot_prediction_intervals(
    predictions     = predicciones,
    y_true          = datos_test,
    target_variable = "OT",
    initial_x_zoom  = None,
    title           = "Valor real vs predicción en los datos de test",
    xaxis_title     = "Fecha",
    yaxis_title     = "OT",
    ax              = ax
)
fill_between_obj = ax.collections[0]
fill_between_obj.set_facecolor('white')
fill_between_obj.set_alpha(0.3)

# Cobertura del intervalo predicho
# ==============================================================================
cobertura = calculate_coverage(
    y_true      = datos.loc[fin_validacion:, 'OT'],
    lower_bound = predicciones["lower_bound"], 
    upper_bound = predicciones["upper_bound"]
)
print(f"Cobertura del intervalo predicho: {round(100 * cobertura, 2)} %")

# Área del intervalo
# ==============================================================================
area = (predicciones["upper_bound"] - predicciones["lower_bound"]).sum()
print(f"Área del intervalo: {round(area, 2)}")

display(metrica)

# Gráfico interactivo
# ==============================================================================
fig = go.Figure([
        go.Scatter(name='Predicción', x=predicciones.index, y=predicciones['pred'], mode='lines'),
        go.Scatter(name='Valor real', x=datos_test.index, y=datos_test['OT'], mode='lines'),
        go.Scatter(
            name='Límite superior', x=predicciones.index, y=predicciones['upper_bound'],
            mode='lines', marker=dict(color="#444"), line=dict(width=0), showlegend=False
        ),
        go.Scatter(
            name='Límite inferior', x=predicciones.index, y=predicciones['lower_bound'],
            marker=dict(color="#444"), line=dict(width=0), mode='lines',
            fillcolor='rgba(68, 68, 68, 0.3)', fill='tonexty', showlegend=False
        )
    ])
fig.update_layout(
    title="Valor real vs predicción en los datos de test",
    xaxis_title= "Fecha",
    yaxis_title= "OT",
    width=800,
    height=400,
    margin=dict(l=10, r=10, t=35, b=10),
    hovermode="x",
    legend=dict(orientation="h", yanchor="top", y=1.1, xanchor="left", x=0.001)
)
fig.show()

# Intervalos de predicción para diferentes coberturas nominales
# ==============================================================================
cuantiles = [
    0.025, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6,
    0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 0.975
]
intervalos = [
    [0.025, 0.975], [0.05, 0.95], [0.1, 0.9], [0.15, 0.85], [0.2, 0.8],
    [0.25, 0.75], [0.3, 0.7], [0.35, 0.65], [0.4, 0.6], [0.45, 0.55]
]
coberturas_observadas = []
areas_observadas = []

metrica, predicciones = backtesting_forecaster(
                          forecaster              = forecaster,
                          y                       = datos['OT'],
                          exog                    = datos[exog_seleccionadas],
                          cv                      = cv,
                          metric                  = 'mean_absolute_error',
                          interval                = cuantiles,
                          n_boot                  = 150,
                          interval_method         = 'bootstrapping',
                          use_in_sample_residuals = False,  # Usar residuos out-sample
                          use_binned_residuals    = True,   # Intervalos condicionados a los valores predichos
                      )
predicciones.head()

# Calcular cobertura y área observada para cada intervalo
# ==============================================================================
for intervalo in intervalos:
    cobertura_observada = calculate_coverage(
                            y_true      = datos.loc[fin_validacion:, 'OT'],
                            lower_bound = predicciones[f"q_{intervalo[0]}"], 
                            upper_bound = predicciones[f"q_{intervalo[1]}"]
                        )
    area_observada = (predicciones[f"q_{intervalo[1]}"] - predicciones[f"q_{intervalo[0]}"]).sum()
    coberturas_observadas.append(100 * cobertura_observada)
    areas_observadas.append(area_observada)

resultados = pd.DataFrame({
              'Intervalo': intervalos,
              'Cobertura nominal': [100 * (intervalo[1] - intervalo[0]) for intervalo in intervalos],
              'Cobertura observada': coberturas_observadas,
              'Área': areas_observadas
          })
resultados.round(2)

import session_info
session_info.show(html=False)