import json
import re
from pathlib import Path
import sys
import warnings

import numpy as np
import pandas as pd
import plotly.colors as pc
import plotly.express as px
import plotly.graph_objects as go
from IPython.display import Markdown, display
from plotly.subplots import make_subplots

warnings.filterwarnings("ignore")

pd.options.display.float_format = "{:,.2f}".format
px.defaults.template = "plotly_white"
px.defaults.height = 450


def resolve_project_root() -> Path:
    candidate_bases = [Path.cwd().resolve()]
    try:
        candidate_bases.insert(0, Path(__file__).resolve().parent)
    except NameError:
        pass

    for base in candidate_bases:
        for candidate in [base, *base.parents]:
            if (candidate / "src" / "preprocess" / "preprocess_config.py").exists():
                return candidate

    return Path.cwd().resolve()


def get_local_calendar_date(index: pd.DatetimeIndex) -> pd.Index:
    # Evitar problemas de medianoche por horario de verano al convertir timestamps con zona horaria a fechas calendario.
    naive_index = index.tz_localize(None) if index.tz is not None else index
    return pd.Index(naive_index.date, name="date")


PROJECT_ROOT = resolve_project_root()
if str(PROJECT_ROOT) not in sys.path:
    sys.path.append(str(PROJECT_ROOT))

from src.preprocess.preprocess_config import PATH_TO_SAVE_PROCESSED_DATA
from utils.report_utils import display_responsive_plot, setup_collapsible_code

setup_collapsible_code()


SERIES_NAME = "P.MONTT220"
REAL_DATA_PATH = Path(PATH_TO_SAVE_PROCESSED_DATA) / "lags" / "real_data.pkl"

DAY_LABELS = ["Lun", "Mar", "Mié", "Jue", "Vie", "Sáb", "Dom"]
MONTH_LABELS = ["Ene", "Feb", "Mar", "Abr", "May", "Jun", "Jul", "Ago", "Sep", "Oct", "Nov", "Dic"]

df_real_data = pd.read_pickle(REAL_DATA_PATH)
series = df_real_data[SERIES_NAME].sort_index().astype(float).rename("value")

df_series = (
    series.to_frame()
    .assign(
        date=lambda frame: get_local_calendar_date(frame.index),
        hour=lambda frame: frame.index.hour,
        day_of_week=lambda frame: frame.index.dayofweek,
        month=lambda frame: frame.index.month,
        year=lambda frame: frame.index.year,
    )
)

display(
    Markdown(
        f"""
### Serie cargada
- **Serie:** `{SERIES_NAME}`
- **Fuente:** `{REAL_DATA_PATH}`
- **Filas:** `{len(df_series):,}`
- **Rango temporal:** `{df_series.index.min()}` a `{df_series.index.max()}`
"""
    )
)

df_series.head()

def infer_series_frequency(ts: pd.Series):
    inferred = pd.infer_freq(ts.index[: min(len(ts), 500)])
    if inferred is not None:
        return inferred
    return ts.index.to_series().diff().dropna().mode().iloc[0]


def build_timestamp_summary(ts: pd.Series):
    freq = infer_series_frequency(ts)
    expected_index = pd.date_range(
        start=ts.index.min(),
        end=ts.index.max(),
        freq=freq,
        tz=ts.index.tz,
    )
    missing_timestamps = expected_index.difference(ts.index)

    summary = pd.DataFrame(
        {
            "metric": [
                "timezone",
                "inferred_frequency",
                "start",
                "end",
                "rows",
                "duplicate_timestamps",
                "missing_timestamps",
                "missing_values",
                "zero_share",
                "negative_share",
            ],
            "value": [
                str(ts.index.tz),
                str(freq),
                ts.index.min(),
                ts.index.max(),
                len(ts),
                int(ts.index.duplicated().sum()),
                int(len(missing_timestamps)),
                int(ts.isna().sum()),
                round(float((ts == 0).mean()), 4),
                round(float((ts < 0).mean()), 4),
            ],
        }
    )
    return summary, missing_timestamps


def compute_iqr_outliers(ts: pd.Series):
    q1, q3 = ts.quantile([0.25, 0.75])
    iqr = q3 - q1
    lower_bound = max(0.0, q1 - 1.5 * iqr)
    upper_bound = q3 + 1.5 * iqr
    outlier_mask = (ts < lower_bound) | (ts > upper_bound)

    return {
        "q1": float(q1),
        "q3": float(q3),
        "iqr": float(iqr),
        "lower_bound": float(lower_bound),
        "upper_bound": float(upper_bound),
        "outlier_mask": outlier_mask,
        "outlier_count": int(outlier_mask.sum()),
        "outlier_share": float(outlier_mask.mean()),
    }


def build_acf_frame(ts: pd.Series, max_lag: int = 24 * 14):
    return pd.DataFrame(
        {
            "lag": np.arange(1, max_lag + 1),
            "autocorr": [ts.autocorr(lag=lag) for lag in range(1, max_lag + 1)],
        }
    )


timestamp_summary, missing_timestamps = build_timestamp_summary(series)
distribution_summary = (
    series.describe(percentiles=[0.01, 0.05, 0.25, 0.50, 0.75, 0.95, 0.99]).to_frame("value")
)
outlier_stats = compute_iqr_outliers(series)

display(Markdown("### Resumen de timestamps y variable objetivo"))
display(timestamp_summary)
display(distribution_summary.T)

if len(missing_timestamps) > 0:
    display(Markdown("### Muestra de timestamps faltantes"))
    missing_preview = pd.Series(missing_timestamps[:10], name="timestamp_faltante")
    display(missing_preview.to_frame())

daily_mean = series.resample("D").mean().rename("daily_mean")
rolling_windows = [7, 14, 30]
rolling_stats_by_window = {
    window: {
        "mean": daily_mean.rolling(window, min_periods=7).mean().rename(f"rolling_mean_{window}d"),
        "std": daily_mean.rolling(window, min_periods=7).std().rename(f"rolling_std_{window}d"),
    }
    for window in rolling_windows
}
default_window = 30

fig = make_subplots(specs=[[{"secondary_y": True}]])
fig.add_trace(
    go.Scattergl(
        x=series.index,
        y=series,
        mode="lines",
        name="Serie horaria",
        line=dict(color="rgba(99, 110, 250, 0.28)", width=1),
    ),
    secondary_y=False,
)
fig.add_trace(
    go.Scatter(
        x=daily_mean.index,
        y=daily_mean,
        mode="lines",
        name="Media diaria",
        line=dict(color="#1f77b4", width=2),
    ),
    secondary_y=False,
)

for window in rolling_windows:
    rolling_mean = rolling_stats_by_window[window]["mean"]
    rolling_std = rolling_stats_by_window[window]["std"]
    is_default_window = window == default_window

    fig.add_trace(
        go.Scatter(
            x=rolling_mean.index,
            y=rolling_mean,
            mode="lines",
            name=f"Media móvil {window}D",
            line=dict(color="#d62728", width=3),
            visible=is_default_window,
        ),
        secondary_y=False,
    )
    fig.add_trace(
        go.Scatter(
            x=rolling_std.index,
            y=rolling_std,
            mode="lines",
            name=f"Desv. est. móvil {window}D",
            line=dict(color="#2ca02c", width=2, dash="dot"),
            visible=is_default_window,
        ),
        secondary_y=True,
    )

dropdown_buttons = []
base_trace_count = 2
traces_per_window = 2
total_trace_count = base_trace_count + len(rolling_windows) * traces_per_window

for selected_window in rolling_windows:
    visible = [True, True] + [False] * (total_trace_count - base_trace_count)
    window_index = rolling_windows.index(selected_window)
    start_idx = base_trace_count + window_index * traces_per_window
    visible[start_idx] = True
    visible[start_idx + 1] = True
    dropdown_buttons.append(
        dict(
            label=f"{selected_window} días",
            method="update",
            args=[
                {"visible": visible},
                {"title": f"{SERIES_NAME} - nivel, tendencia y volatilidad ({selected_window}D)"},
            ],
        )
    )

fig.update_layout(
    title=f"{SERIES_NAME} - nivel, tendencia y volatilidad ({default_window}D)",
    hovermode="x unified",
    legend=dict(orientation="h", y=1.08, x=0),
    updatemenus=[
        dict(
            buttons=dropdown_buttons,
            direction="down",
            x=1.01,
            xanchor="left",
            y=1.12,
            yanchor="top",
        )
    ],
)
fig.update_yaxes(title_text="Valor CMG", secondary_y=False)
fig.update_yaxes(title_text="Volatilidad diaria", secondary_y=True)
display_responsive_plot(fig)

def same_hour_rolling_mean(
    ts: pd.Series, window: int, min_periods: int | None = None
) -> list[go.Scatter]:
    """Media móvil sobre las `window` observaciones previas a la misma hora del día (shift(1) primero)."""
    if min_periods is None:
        min_periods = max(1, min(7, window))
    palette = pc.qualitative.Dark24
    traces = []
    for h in range(24):
        mask = ts.index.hour == h
        sub = ts[mask]
        rolled = sub.shift(1).rolling(window, min_periods=min_periods).mean()
        color = palette[h % len(palette)]
        traces.append(
            go.Scatter(
                x=rolled.index,
                y=rolled,
                mode="lines",
                name=f"{h:02d}:00",
                line=dict(color=color, width=1.8),
                legendgroup=str(h),
            )
        )
    return traces


SAME_HOUR_WINDOWS = (7, 14, 30)
fig_same_hour = go.Figure()
fig_same_hour.add_trace(
    go.Scattergl(
        x=series.index,
        y=series,
        mode="lines",
        name="Serie horaria",
        line=dict(color="rgba(99, 110, 250, 0.22)", width=1),
        legendgroup="bg",
        showlegend=True,
    )
)
for wi, w in enumerate(SAME_HOUR_WINDOWS):
    hour_traces = same_hour_rolling_mean(series, w)
    for h, tr in enumerate(hour_traces):
        tr.visible = wi == 0
        tr.showlegend = wi == 0
        fig_same_hour.add_trace(tr)

n_per_win = 24
buttons = []
for wi, w in enumerate(SAME_HOUR_WINDOWS):
    visible = [True]
    showlegend = [True]
    for wj in range(len(SAME_HOUR_WINDOWS)):
        on = wj == wi
        visible.extend([on] * n_per_win)
        showlegend.extend([on] * n_per_win)
    buttons.append(
        dict(
            label=f"{w} misma hora",
            method="update",
            args=[
                {"visible": visible, "showlegend": showlegend},
                {
                    "title": {
                        "text": (
                            f"{SERIES_NAME} — media móvil de {w} valores previos a la misma hora "
                            f"(una traza por hora del día)"
                        )
                    }
                },
            ],
        )
    )

fig_same_hour.update_layout(
    title=dict(
        text=(
            f"{SERIES_NAME} — media móvil de {SAME_HOUR_WINDOWS[0]} valores previos a la misma hora "
            f"(una traza por hora del día)"
        )
    ),
    hovermode="x unified",
    legend=dict(orientation="h", y=-0.22, x=0, font=dict(size=9)),
    yaxis_title="Valor CMG",
    xaxis_title="Tiempo",
    updatemenus=[
        dict(
            buttons=buttons,
            direction="down",
            showactive=True,
            x=0.0,
            xanchor="left",
            y=1.12,
            yanchor="top",
        )
    ],
)
#fig_same_hour.show()
display_responsive_plot(fig_same_hour)

outlier_mask = outlier_stats["outlier_mask"]
outlier_preview = series[outlier_mask].sort_values(ascending=False).head(15).to_frame("value")

fig = make_subplots(
    rows=1,
    cols=2,
    subplot_titles=("Distribución", "Diagrama de caja"),
    column_widths=[0.7, 0.3],
)
fig.add_trace(
    go.Histogram(
        x=series,
        nbinsx=70,
        name="Valores horarios",
        marker_color="#1f77b4",
        opacity=0.85,
    ),
    row=1,
    col=1,
)
fig.add_trace(
    go.Box(
        y=series,
        name=SERIES_NAME,
        marker_color="#d62728",
        boxpoints=False,
    ),
    row=1,
    col=2,
)
fig.update_layout(
    title=f"{SERIES_NAME} - distribución del objetivo y valores extremos",
    showlegend=False,
)
fig.update_xaxes(title_text="Valor CMG", row=1, col=1)
fig.update_yaxes(title_text="Cantidad", row=1, col=1)
fig.update_yaxes(title_text="Valor CMG", row=1, col=2)
display_responsive_plot(fig)

display(
    Markdown(
        f"""
### Detección de outliers
- **Método:** IQR sobre el objetivo horario
- **Límites:** `{outlier_stats["lower_bound"]:.2f}` a `{outlier_stats["upper_bound"]:.2f}`
- **Puntos marcados:** `{outlier_stats["outlier_count"]:,}` (`{outlier_stats["outlier_share"]:.1%}`)
- **Decisión actual:** mantenerlos por ahora y tratarlos como posibles regímenes de spike, no como datos erróneos
"""
    )
)
display(outlier_preview)

hourly_profile = df_series.groupby("hour")["value"].mean().reindex(range(24))
weekday_profile = df_series.groupby("day_of_week")["value"].mean().reindex(range(7))
monthly_profile = df_series.groupby("month")["value"].mean().reindex(range(1, 13))
weekday_hour_profile = (
    df_series.groupby(["hour", "day_of_week"])["value"]
    .mean()
    .unstack("day_of_week")
    .reindex(index=range(24), columns=range(7))
)

fig = make_subplots(
    rows=2,
    cols=2,
    subplot_titles=(
        "Promedio por hora",
        "Promedio por día de semana",
        "Promedio por mes",
        "Mapa de calor hora x día de semana",
    ),
    specs=[
        [{"type": "xy"}, {"type": "xy"}],
        [{"type": "xy"}, {"type": "heatmap"}],
    ],
)

fig.add_trace(
    go.Scatter(
        x=hourly_profile.index,
        y=hourly_profile.values,
        mode="lines+markers",
        name="Perfil horario",
        line=dict(color="#636efa", width=2),
    ),
    row=1,
    col=1,
)
fig.add_trace(
    go.Bar(
        x=list(range(7)),
        y=weekday_profile.values,
        name="Perfil semanal",
        marker_color="#00cc96",
    ),
    row=1,
    col=2,
)
fig.add_trace(
    go.Bar(
        x=list(range(1, 13)),
        y=monthly_profile.values,
        name="Perfil mensual",
        marker_color="#ef553b",
    ),
    row=2,
    col=1,
)
fig.add_trace(
    go.Heatmap(
        z=weekday_hour_profile.values,
        x=DAY_LABELS,
        y=weekday_hour_profile.index,
        colorscale="YlOrRd",
        colorbar=dict(title="CMG medio"),
        name="Mapa de calor",
    ),
    row=2,
    col=2,
)

fig.update_layout(
    title=f"{SERIES_NAME} - perfil de estacionalidad",
    showlegend=False,
    height=800,
)
fig.update_xaxes(title_text="Hora", row=1, col=1)
fig.update_yaxes(title_text="CMG medio", row=1, col=1)
fig.update_xaxes(
    title_text="Día de semana",
    tickmode="array",
    tickvals=list(range(7)),
    ticktext=DAY_LABELS,
    row=1,
    col=2,
)
fig.update_yaxes(title_text="CMG medio", row=1, col=2)
fig.update_xaxes(
    title_text="Mes",
    tickmode="array",
    tickvals=list(range(1, 13)),
    ticktext=MONTH_LABELS,
    row=2,
    col=1,
)
fig.update_yaxes(title_text="CMG medio", row=2, col=1)
fig.update_xaxes(title_text="Día de semana", row=2, col=2)
fig.update_yaxes(title_text="Hora", row=2, col=2)
display_responsive_plot(fig)

acf_frame = build_acf_frame(series, max_lag=24 * 14)
lag_frame = pd.DataFrame(
    {
        "value": series,
        "lag_1h": series.shift(1),
        "lag_24h": series.shift(24),
        "lag_168h": series.shift(168),
    }
).dropna()
lag_sample = lag_frame.sample(min(len(lag_frame), 4000), random_state=7)

acf_min = float(acf_frame["autocorr"].min())
acf_max = float(acf_frame["autocorr"].max())

fig = make_subplots(
    rows=2,
    cols=2,
    subplot_titles=(
        "Autocorrelación hasta 14 días",
        "Valor actual vs lag 1h",
        "Valor actual vs lag 24h",
        "Valor actual vs lag 168h",
    ),
)

fig.add_trace(
    go.Scatter(
        x=acf_frame["lag"],
        y=acf_frame["autocorr"],
        mode="lines",
        name="ACF",
        line=dict(color="#636efa", width=2),
    ),
    row=1,
    col=1,
)
for lag_marker, label, color in [(24, "24h", "#ef553b"), (168, "168h", "#00cc96")]:
    fig.add_trace(
        go.Scatter(
            x=[lag_marker, lag_marker],
            y=[acf_min, acf_max],
            mode="lines",
            line=dict(color=color, dash="dot"),
            name=label,
            showlegend=(lag_marker == 24),
        ),
        row=1,
        col=1,
    )

fig.add_trace(
    go.Scattergl(
        x=lag_sample["lag_1h"],
        y=lag_sample["value"],
        mode="markers",
        marker=dict(size=4, color="#636efa", opacity=0.30),
        name="Lag 1h",
    ),
    row=1,
    col=2,
)
fig.add_trace(
    go.Scattergl(
        x=lag_sample["lag_24h"],
        y=lag_sample["value"],
        mode="markers",
        marker=dict(size=4, color="#ef553b", opacity=0.30),
        name="Lag 24h",
    ),
    row=2,
    col=1,
)
fig.add_trace(
    go.Scattergl(
        x=lag_sample["lag_168h"],
        y=lag_sample["value"],
        mode="markers",
        marker=dict(size=4, color="#00cc96", opacity=0.30),
        name="Lag 168h",
    ),
    row=2,
    col=2,
)

fig.update_layout(
    title=f"{SERIES_NAME} - estructura de dependencia útil para pronóstico",
    height=800,
    legend=dict(orientation="h", y=1.08, x=0),
)
fig.update_xaxes(title_text="Rezago", row=1, col=1)
fig.update_yaxes(title_text="Autocorrelación", row=1, col=1)
fig.update_xaxes(title_text="Lag 1h", row=1, col=2)
fig.update_yaxes(title_text="Valor actual", row=1, col=2)
fig.update_xaxes(title_text="Lag 24h", row=2, col=1)
fig.update_yaxes(title_text="Valor actual", row=2, col=1)
fig.update_xaxes(title_text="Lag 168h", row=2, col=2)
fig.update_yaxes(title_text="Valor actual", row=2, col=2)
display_responsive_plot(fig)

display(acf_frame[acf_frame["lag"].isin([1, 24, 168, 24 * 14])].reset_index(drop=True))

lag_summary = {
    lag: series.autocorr(lag=lag)
    for lag in [1, 24, 168]
}
peak_hour = int(hourly_profile.idxmax())
peak_weekday = DAY_LABELS[int(weekday_profile.idxmax())]
peak_month = MONTH_LABELS[int(monthly_profile.idxmax()) - 1]

display(
    Markdown(
        f"""
## Primeros hallazgos para `P.MONTT220`

- El objetivo es una serie **horaria** en `{series.index.tz}` **sin timestamps faltantes, sin duplicados y sin valores nulos** en la muestra actual.
- La distribución tiene **sesgo a la derecha** e incluye muchas observaciones cercanas a cero: **{(series == 0).mean():.1%}** de las filas son exactamente cero.
- Existe una clara **estacionalidad intradiaria**: el nivel promedio más alto aparece alrededor de las **{peak_hour}:00**.
- También hay una clara **estructura semanal y anual**: el día de semana con mayor media es **{peak_weekday}**, y el mes con mayor media es **{peak_month}**.
- El objetivo mantiene fuerte memoria del historial reciente: la autocorrelación es **{lag_summary[1]:.2f}** a 1 hora, **{lag_summary[24]:.2f}** a 24 horas, y **{lag_summary[168]:.2f}** a 168 horas.
- Las estadísticas diarias móviles muestran que el nivel **no es estacionario**, por lo que el modelado futuro debe usar **validación basada en tiempo** y permitir cambios de régimen.

### Qué parece útil para pronóstico
- Usar al menos **lag 1h, lag 24h y lag 168h** como features temporales base.
- Agregar features de calendario: **hora, día de semana y mes**.
- Mantener un flag para **regímenes de spike / valores extremos** en vez de descartarlos ciegamente.
- Revisitar la serie con drivers exógenos a continuación, porque el objetivo por sí solo ya muestra estacionalidad pero también fuertes cambios de nivel.
"""
    )
)

HOLDOUT_DAYS = 90
PLOT_WINDOW_DAYS = 14
HORIZONS = list(range(1, 25))
BASELINE_COLUMNS = ["current_value", "same_hour_previous_day", "same_hour_previous_week"]


def mean_absolute_error_np(y_true: np.ndarray, y_pred: np.ndarray) -> float:
    return float(np.mean(np.abs(y_true - y_pred)))


def rmse_np(y_true: np.ndarray, y_pred: np.ndarray) -> float:
    return float(np.sqrt(np.mean((y_true - y_pred) ** 2)))


def wape_np(y_true: np.ndarray, y_pred: np.ndarray) -> float:
    denominator = float(np.abs(y_true).sum())
    if denominator == 0:
        return np.nan
    return float(np.abs(y_true - y_pred).sum() / denominator)


def smape_np(y_true: np.ndarray, y_pred: np.ndarray) -> float:
    denominator = np.abs(y_true) + np.abs(y_pred)
    valid_mask = denominator > 0
    if not valid_mask.any():
        return np.nan
    return float(np.mean(2 * np.abs(y_true[valid_mask] - y_pred[valid_mask]) / denominator[valid_mask]) * 100)


holdout_start = series.index.max() - pd.Timedelta(days=HOLDOUT_DAYS)
horizon_eval_frames = {}
baseline_metrics = []

for horizon in HORIZONS:
    horizon_df = pd.DataFrame(
        {
            "actual": series.shift(-horizon),
            "current_value": series,
            "same_hour_previous_day": series.shift(24 - horizon),
            "same_hour_previous_week": series.shift(168 - horizon),
        }
    )
    horizon_df["forecast_origin"] = horizon_df.index
    horizon_df["target_timestamp"] = horizon_df.index + pd.to_timedelta(horizon, unit="h")
    horizon_df["origin_hour"] = horizon_df["forecast_origin"].dt.hour
    horizon_df["target_hour"] = horizon_df["target_timestamp"].dt.hour
    horizon_df["horizon"] = horizon
    horizon_df["split"] = np.where(horizon_df.index < holdout_start, "train", "test")
    horizon_df = horizon_df.dropna().copy()

    horizon_eval_frames[horizon] = horizon_df

    test_horizon_df = horizon_df.loc[horizon_df["split"] == "test"].copy()
    y_true = test_horizon_df["actual"].to_numpy()

    for baseline_name in BASELINE_COLUMNS:
        y_pred = test_horizon_df[baseline_name].to_numpy()
        baseline_metrics.append(
            {
                "horizon": horizon,
                "baseline": baseline_name,
                "mae": mean_absolute_error_np(y_true, y_pred),
                "rmse": rmse_np(y_true, y_pred),
                "wape_pct": wape_np(y_true, y_pred) * 100,
                "smape_pct": smape_np(y_true, y_pred),
            }
        )

baseline_metrics_df = pd.DataFrame(baseline_metrics)
best_by_horizon_df = (
    baseline_metrics_df.sort_values(["horizon", "wape_pct", "mae", "baseline"])
    .groupby("horizon", as_index=False)
    .first()
)
best_baseline_by_horizon = dict(zip(best_by_horizon_df["horizon"], best_by_horizon_df["baseline"]))

average_wape_df = (
    baseline_metrics_df.groupby("baseline", as_index=False)["wape_pct"]
    .mean()
    .sort_values("wape_pct")
    .rename(columns={"wape_pct": "mean_wape_pct"})
)
overall_best_baseline = average_wape_df.iloc[0]["baseline"]
baseline_win_counts = best_by_horizon_df["baseline"].value_counts().to_dict()
last_valid_origin = max(frame.index.max() for frame in horizon_eval_frames.values())

display(
    Markdown(
        f"""
### Diseño del holdout
- **Período de entrenamiento (origen de pronóstico):** `{series.index.min()}` a `{holdout_start - pd.Timedelta(hours=1)}`
- **Período de prueba (origen de pronóstico):** `{holdout_start}` a `{last_valid_origin}`
- **Horizontes evaluados:** `1h` a `24h`
- **Objetivo:** identificar el baseline univariado más fuerte antes de usar covariables externas
"""
    )
)
display(average_wape_df)
display(best_by_horizon_df[["horizon", "baseline", "wape_pct", "mae", "rmse"]])

display(
    Markdown(
        """
### Como leer WAPE
- **WAPE (Weighted Absolute Percentage Error)** se calcula como `sum(abs(y_real - y_pred)) / sum(abs(y_real))`.
- En este notebook se muestra como porcentaje, es decir, `WAPE * 100`.
- Si un baseline tiene **WAPE = 12%**, significa que el error absoluto acumulado equivale al `12%` del total real observado.
- **Mas bajo es mejor**, porque implica menos error relativo sobre la serie real.
"""
    )
)

fig = make_subplots(rows=1, cols=2, subplot_titles=("División del origen de pronóstico", "Mejor baseline por horizonte"))
fig.add_trace(
    go.Scattergl(
        x=series.index,
        y=series,
        mode="lines",
        name=SERIES_NAME,
        line=dict(color="#1f77b4", width=1.2),
    ),
    row=1,
    col=1,
)
fig.add_vline(x=holdout_start, line_dash="dash", line_color="#d62728", row=1, col=1)
fig.add_annotation(
    x=holdout_start,
    y=float(series.max()),
    text="Inicio holdout",
    showarrow=True,
    arrowhead=2,
    ay=-40,
    xref="x1",
    yref="y1",
)
fig.add_trace(
    go.Bar(
        x=best_by_horizon_df["horizon"],
        y=best_by_horizon_df["wape_pct"],
        text=best_by_horizon_df["baseline"],
        name="Mejor baseline",
        marker_color="#00cc96",
        hovertemplate="Horizon=%{x}h<br>WAPE=%{y:.2f}%<br>Baseline=%{text}<extra></extra>",
    ),
    row=1,
    col=2,
)
fig.update_layout(
    title=f"{SERIES_NAME} - configuración holdout multi-horizonte",
    hovermode="x unified",
    height=500,
)
fig.update_xaxes(title_text="Timestamp", row=1, col=1)
fig.update_yaxes(title_text="Valor CMG", row=1, col=1)
fig.update_xaxes(title_text="Horizonte", row=1, col=2)
fig.update_yaxes(title_text="Mejor WAPE (%)", row=1, col=2)
display_responsive_plot(fig)

metric_options = [
    ("wape_pct", "WAPE (%)"),
    ("mae", "MAE"),
    ("rmse", "RMSE"),
    ("smape_pct", "sMAPE (%)"),
]

fig = go.Figure()
for metric_name, metric_label in metric_options:
    for baseline_name in BASELINE_COLUMNS:
        metric_slice = baseline_metrics_df.loc[baseline_metrics_df["baseline"] == baseline_name].sort_values("horizon")
        fig.add_trace(
            go.Scatter(
                x=metric_slice["horizon"],
                y=metric_slice[metric_name],
                mode="lines+markers",
                name=baseline_name,
                visible=(metric_name == "wape_pct"),
                hovertemplate="Horizon=%{x}h<br>" + metric_label + "=%{y:.2f}<br>Baseline=" + baseline_name + "<extra></extra>",
            )
        )

buttons = []
traces_per_metric = len(BASELINE_COLUMNS)
for metric_idx, (metric_name, metric_label) in enumerate(metric_options):
    visible = [False] * (len(metric_options) * traces_per_metric)
    start = metric_idx * traces_per_metric
    for offset in range(traces_per_metric):
        visible[start + offset] = True
    buttons.append(
        dict(
            label=metric_label,
            method="update",
            args=[
                {"visible": visible},
                {"title": f"{SERIES_NAME} - rendimiento de baselines por horizonte ({metric_label})", "yaxis": {"title": metric_label}},
            ],
        )
    )

fig.update_layout(
    title=f"{SERIES_NAME} - rendimiento de baselines por horizonte (WAPE %)",
    xaxis_title="Horizonte de pronóstico (horas adelante)",
    yaxis_title="WAPE (%)",
    hovermode="x unified",
    updatemenus=[
        dict(
            buttons=buttons,
            direction="down",
            x=1.02,
            xanchor="left",
            y=1.14,
            yanchor="top",
        )
    ],
    legend=dict(orientation="h", y=1.08, x=0),
)
display_responsive_plot(fig)

default_horizon = 24

diagnostic_frames = {}
for horizon in HORIZONS:
    winner_baseline = best_baseline_by_horizon[horizon]
    diagnostic_df = horizon_eval_frames[horizon].loc[horizon_eval_frames[horizon]["split"] == "test"].copy()
    diagnostic_df["prediction"] = diagnostic_df[winner_baseline]
    diagnostic_df["residual"] = diagnostic_df["actual"] - diagnostic_df["prediction"]
    diagnostic_df["abs_error"] = diagnostic_df["residual"].abs()

    focus_start = diagnostic_df["target_timestamp"].max() - pd.Timedelta(days=PLOT_WINDOW_DAYS)
    focus_df = diagnostic_df.loc[diagnostic_df["target_timestamp"] >= focus_start].copy()
    error_by_target_hour = diagnostic_df.groupby("target_hour")["abs_error"].mean().reindex(range(24))

    diagnostic_frames[horizon] = {
        "winner_baseline": winner_baseline,
        "focus_df": focus_df,
        "error_by_target_hour": error_by_target_hour,
    }

fig = make_subplots(
    rows=1,
    cols=2,
    subplot_titles=(
        f"Real vs mejor baseline (ventana {PLOT_WINDOW_DAYS}D)",
        "Error absoluto medio por hora de entrega",
    ),
    column_widths=[0.68, 0.32],
)

for horizon in HORIZONS:
    horizon_payload = diagnostic_frames[horizon]
    focus_df = horizon_payload["focus_df"]
    winner_baseline = horizon_payload["winner_baseline"]
    error_by_target_hour = horizon_payload["error_by_target_hour"]
    is_default = horizon == default_horizon

    fig.add_trace(
        go.Scatter(
            x=focus_df["target_timestamp"],
            y=focus_df["actual"],
            mode="lines",
            name=f"Real ({horizon}h)",
            line=dict(color="#1f77b4", width=2),
            visible=is_default,
        ),
        row=1,
        col=1,
    )
    fig.add_trace(
        go.Scatter(
            x=focus_df["target_timestamp"],
            y=focus_df["prediction"],
            mode="lines",
            name=f"{winner_baseline} ({horizon}h)",
            line=dict(color="#ef553b", width=2),
            visible=is_default,
        ),
        row=1,
        col=1,
    )
    fig.add_trace(
        go.Bar(
            x=list(range(24)),
            y=error_by_target_hour.values,
            name=f"MAE por hora de entrega ({horizon}h)",
            marker_color="#00cc96",
            visible=is_default,
        ),
        row=1,
        col=2,
    )
buttons = []
traces_per_horizon = 3
for horizon in HORIZONS:
    visible = [False] * (len(HORIZONS) * traces_per_horizon)
    start = (horizon - 1) * traces_per_horizon
    for offset in range(traces_per_horizon):
        visible[start + offset] = True
    buttons.append(
        dict(
            label=f"{horizon}h adelante",
            method="update",
            args=[
                {"visible": visible},
                {
                    "title": f"{SERIES_NAME} - diagnóstico por horizonte ({horizon}h adelante, mejor baseline: {best_baseline_by_horizon[horizon]})",
                },
            ],
        )
    )

fig.update_layout(
    title=f"{SERIES_NAME} - diagnóstico por horizonte ({default_horizon}h adelante, mejor baseline: {best_baseline_by_horizon[default_horizon]})",
    height=520,
    hovermode="x unified",
    updatemenus=[
        dict(
            buttons=buttons,
            direction="down",
            x=1.02,
            xanchor="left",
            y=1.14,
            yanchor="top",
        )
    ],
    legend=dict(orientation="h", y=1.08, x=0),
)
fig.update_xaxes(title_text="Timestamp objetivo", row=1, col=1)
fig.update_yaxes(title_text="Valor CMG", row=1, col=1)
fig.update_xaxes(title_text="Hora de entrega", row=1, col=2)
fig.update_yaxes(title_text="Error absoluto medio", row=1, col=2)
display_responsive_plot(fig)

best_1h_row = best_by_horizon_df.loc[best_by_horizon_df["horizon"] == 1].iloc[0]
best_24h_row = best_by_horizon_df.loc[best_by_horizon_df["horizon"] == 24].iloc[0]
win_summary = ", ".join(f"`{name}`: {count}" for name, count in baseline_win_counts.items())

display(
    Markdown(
        f"""
## Hallazgos del siguiente paso para `P.MONTT220`

- A través de los horizontes `1` a `24`, el mejor baseline promedio es **`{overall_best_baseline}`**.
- Los ganadores por horizonte se distribuyen como: {win_summary}.
- A **1h adelante**, el baseline más fuerte es **`{best_1h_row["baseline"]}`** con **WAPE {best_1h_row["wape_pct"]:.2f}%**.
- A **24h adelante**, el baseline más fuerte es **`{best_24h_row["baseline"]}`** con **WAPE {best_24h_row["wape_pct"]:.2f}%**.
- En esta serie, los horizontes muy cortos se comportan más como **persistencia**, mientras que la mayor parte de la ventana day-ahead se comporta más como **estacionalidad diaria**.

### Qué significa esto antes de agregar otras series
- El valor de los drivers externos debe evaluarse **horizonte por horizonte**, no con un puntaje agregado único.
- El siguiente paso útil es comparar variables externas candidatas contra la estructura residual para los horizontes donde estos baselines dejan de mejorar.
"""
    )
)

OTHER_DATA_PATH = Path(PATH_TO_SAVE_PROCESSED_DATA) / "lags" / "other_data_w_lags.pkl"
COLUMN_TYPE_NAMES_PATH = Path(PATH_TO_SAVE_PROCESSED_DATA) / "metadata" / "column_type_names.json"
KNOWN_FEATURE_SELECTION_PATH = Path(PATH_TO_SAVE_PROCESSED_DATA) / "metadata" / "feature_selection_v6_PMONTT220_all_24h.json"
MANUAL_KNOWN_FEATURE_SELECTION = {
    "PMontt220": {
        "known_values": [
            'PMontt220_PRG_latest',
            'Charrua220_PRG_latest',
            "CIRUELOS220_rad",
            "CIRUELOS220_temp",
            "GDT_Dem",
            "GD2_Dem",
            "GDT_GRat_EO",
            "GDT_GRat_HE",
            "GDT_GRat_HP",
            "GDT_GRat_PV",
            "GDT_GRat_TER",
            "GDT_GRat_PV",
            "GD2_GRat_EO",
            "GD2_GRat_HE",
            "GD2_GRat_HP",
            "GDT_GP_TER",
            "GDT_GP_HE",
            "GDT_GP_HP",
            "GDT_GP_PV",
            "GDT_GP_EO",
            "GDT_GP_TER",
            "GD2_GP_TER",
            "GD2_GP_HE",
            "GD2_GP_HP",
            "GD2_GP_PV",
            "GD2_GP_EO",
            "GD2_GP_TER",
            "TCO_HE",
            "TCO_TER",
            "GD2_GCurt_EO",
            "GDT_GCurt_PV",
            "GDT_GCurt_EO",
            "GDT_D-GP_PV+EO",
            "GDT_D-GP_+HP",
            "GDT_D-GP_+HE",
            "GD2_D-GP_PV+EO",
            "GD2_D-GP_+HP",
            "GD2_D-GP_+HE",
            # 'GDT_SSCC_CPF+_EO',
            # 'GDT_SSCC_CPF+_HE',
            # 'GDT_SSCC_CPF+_PV',
            # 'GDT_SSCC_CPF+_TER',
            # 'GDT_SSCC_CPF-_EO',
            # 'GDT_SSCC_CPF-_HE',
            # 'GDT_SSCC_CPF-_PV',
            # 'GDT_SSCC_CPF-_TER',
            # 'GDT_SSCC_CSF+_EO',
            # 'GDT_SSCC_CSF+_HE',
            # 'GDT_SSCC_CSF+_PV',
            # 'GDT_SSCC_CSF+_TER',
            # 'GDT_SSCC_CSF-_EO',
            # 'GDT_SSCC_CSF-_HE',
            # 'GDT_SSCC_CSF-_PV',
            # 'GDT_SSCC_CSF-_TER',
            # 'GDT_SSCC_CTF+_EO',
            # 'GDT_SSCC_CTF+_HE',
            # 'GDT_SSCC_CTF+_HP',
            # 'GDT_SSCC_CTF+_PV',
            # 'GDT_SSCC_CTF+_TER',
            # 'GDT_SSCC_CTF-_EO',
            # 'GDT_SSCC_CTF-_HE',
            # 'GDT_SSCC_CTF-_HP',
            # 'GDT_SSCC_CTF-_PV',
            # 'GDT_SSCC_CTF-_TER',
            # 'GD2_SSCC_CPF+_EO',
            # 'GD2_SSCC_CPF+_HE',
            # 'GD2_SSCC_CPF+_PV',
            # 'GD2_SSCC_CPF+_TER',
            # 'GD2_SSCC_CPF-_EO',
            # 'GD2_SSCC_CPF-_HE',
            # 'GD2_SSCC_CPF-_PV',
            # 'GD2_SSCC_CPF-_TER',
            # 'GD2_SSCC_CSF+_EO',
            # 'GD2_SSCC_CSF+_HE',
            # 'GD2_SSCC_CSF+_PV',
            # 'GD2_SSCC_CSF+_TER',
            # 'GD2_SSCC_CSF-_EO',
            # 'GD2_SSCC_CSF-_HE',
            # 'GD2_SSCC_CSF-_PV',
            # 'GD2_SSCC_CSF-_TER',
            # 'GD2_SSCC_CTF+_EO',
            # 'GD2_SSCC_CTF+_HE',
            # 'GD2_SSCC_CTF+_HP',
            # 'GD2_SSCC_CTF+_PV',
            # 'GD2_SSCC_CTF+_TER',
            # 'GD2_SSCC_CTF-_EO',
            # 'GD2_SSCC_CTF-_HE',
            # 'GD2_SSCC_CTF-_HP',
            # 'GD2_SSCC_CTF-_PV',
            # 'GD2_SSCC_CTF-_TER',
            # 'GDT_D-GRat_PV+EO_RAMP',
            # 'GDT_D-GRat_+HP_RAMP',
            # 'GDT_D-GRat_+HE_RAMP',
            # 'GDT_D-GRat_+TER_RAMP',
            # 'GD2_D-GRat_PV+EO_RAMP',
            # 'GD2_D-GRat_+HP_RAMP',
            # 'GD2_D-GRat_+HE_RAMP',
            # 'GD2_D-GRat_+TER_RAMP',
            'GDT_D-GRat_PV+EO',
            'GDT_D-GRat_+HP',
            'GDT_D-GRat_+HE',
            'GDT_D-GRat_+TER',
            'GD2_D-GRat_PV+EO',
            'GD2_D-GRat_+HP',
            'GD2_D-GRat_+HE',
            'GD2_D-GRat_+TER',
            # 'GDT_D-GP_PV+EO_RAMP',
            # 'GDT_D-GP_+HP_RAMP',
            # 'GDT_D-GP_+HE_RAMP',
            # 'GDT_D-GP_+TER_RAMP',
            # 'GD2_D-GP_PV+EO_RAMP',
            # 'GD2_D-GP_+HP_RAMP',
            # 'GD2_D-GP_+HE_RAMP',
            # 'GD2_D-GP_+TER_RAMP',
            'GDT_D-GP_PV+EO',
            'GDT_D-GP_+HP',
            'GDT_D-GP_+HE',
            'GDT_D-GP_+TER',
            'GD2_D-GP_PV+EO',
            'GD2_D-GP_+HP',
            'GD2_D-GP_+HE',
            'GD2_D-GP_+TER',
        ]
    }
}


def normalize_feature_name(feature_name: str) -> str:
    return feature_name.replace(".", "")


def infer_feature_group(feature_name: str, base_group_lookup: dict) -> str:
    lag_match = re.match(r"^(.*)_Lag\d+h$", feature_name)
    if lag_match:
        base_feature = lag_match.group(1)
        if normalize_feature_name(base_feature) == normalize_feature_name(SERIES_NAME):
            return "target_lag"
        base_group = base_group_lookup.get(base_feature)
        return f"{base_group}_lag" if base_group else "lagged_known"
    return base_group_lookup.get(feature_name, "known_other")


df_other_data = pd.read_pickle(OTHER_DATA_PATH).sort_index()
column_type_names = json.loads(COLUMN_TYPE_NAMES_PATH.read_text(encoding="utf-8"))

base_group_lookup = {}
for group_name, group_columns in column_type_names.items():
    for group_column in group_columns:
        base_group_lookup[group_column] = group_name

normalized_feature_lookup = {normalize_feature_name(column): column for column in df_other_data.columns}

manual_feature_aliases = MANUAL_KNOWN_FEATURE_SELECTION.get("PMontt220", {}).get("known_values", [])
selected_feature_aliases = list(dict.fromkeys(manual_feature_aliases))
selection_source_label = "manual known_values list for PMontt220"

if not selected_feature_aliases and KNOWN_FEATURE_SELECTION_PATH.exists():
    known_feature_selection_payload = json.loads(KNOWN_FEATURE_SELECTION_PATH.read_text(encoding="utf-8"))
    selected_feature_aliases = known_feature_selection_payload.get("known_reals", [])
    selection_source_label = KNOWN_FEATURE_SELECTION_PATH.name
elif not selected_feature_aliases:
    selection_source_label = "all columns from other_data_w_lags.pkl"

feature_rows = []
missing_selected_features = []

if selected_feature_aliases:
    for feature_alias in selected_feature_aliases:
        actual_feature = normalized_feature_lookup.get(feature_alias)
        if actual_feature is None:
            missing_selected_features.append(feature_alias)
            continue
        feature_rows.append(
            {
                "feature_alias": feature_alias,
                "feature_actual": actual_feature,
                "group": infer_feature_group(actual_feature, base_group_lookup),
            }
        )
else:
    for actual_feature in df_other_data.columns:
        feature_rows.append(
            {
                "feature_alias": normalize_feature_name(actual_feature),
                "feature_actual": actual_feature,
                "group": infer_feature_group(actual_feature, base_group_lookup),
            }
        )

feature_catalog_df = pd.DataFrame(feature_rows).drop_duplicates(subset=["feature_actual"]).reset_index(drop=True)
df_other_data_to_test = df_other_data[feature_catalog_df["feature_actual"].tolist()].copy()

group_counts_df = (
    feature_catalog_df.groupby("group", as_index=False)
    .size()
    .rename(columns={"size": "feature_count"})
    .sort_values(["feature_count", "group"], ascending=[False, True])
)

display(
    Markdown(
        f"""
### Covariables conocidas candidatas
- **Dataframe fuente:** `{OTHER_DATA_PATH.name}`
- **Fuente de selección:** `{selection_source_label}`
- **Features conocidos candidatos retenidos:** `{len(df_other_data_to_test.columns)}`
- **Rango temporal:** `{df_other_data_to_test.index.min()}` a `{df_other_data_to_test.index.max()}`
- **Supuesto importante:** cada columna en `df_other_data_to_test` se trata como conocida al origen del pronóstico
"""
    )
)
display(group_counts_df)
display(feature_catalog_df.sort_values(["group", "feature_alias"]).reset_index(drop=True))

if missing_selected_features:
    display(
        Markdown(
            "### Features seleccionados no encontrados en el dataframe actual\n"
            + "\n".join(f"- `{feature}`" for feature in missing_selected_features)
        )
    )

fig = px.bar(
    group_counts_df,
    x="group",
    y="feature_count",
    title=f"{SERIES_NAME} - conteo de covariables conocidas por grupo fuente",
    text="feature_count",
)
fig.update_layout(xaxis_title="Grupo de features", yaxis_title="Cantidad")
display_responsive_plot(fig)

if df_other_data_to_test.empty:
    raise RuntimeError("No hay covariables conocidas candidatas disponibles para el screening.")

train_feature_frame = df_other_data_to_test.loc[df_other_data_to_test.index < holdout_start].copy()
feature_catalog_by_actual = feature_catalog_df.set_index("feature_actual")

screening_frame = (
    pd.concat([series.rename("target"), train_feature_frame], axis=1)
    .dropna()
    .copy()
)

if screening_frame.empty:
    raise RuntimeError("No hay observaciones completas entre la serie objetivo y las covariables conocidas antes del holdout.")

target_corr = screening_frame[df_other_data_to_test.columns].corrwith(screening_frame["target"])

contemporaneous_screen_df = (
    feature_catalog_df.assign(target_corr=lambda frame: frame["feature_actual"].map(target_corr))
    .dropna(subset=["target_corr"])
    .assign(target_corr_abs=lambda frame: frame["target_corr"].abs())
    .sort_values(["target_corr_abs", "feature_alias"], ascending=[False, True])
    .reset_index(drop=True)
)

if contemporaneous_screen_df.empty:
    raise RuntimeError("No se pudieron calcular correlaciones validas para las covariables conocidas.")

top_target_features_df = contemporaneous_screen_df.head(12)[
    ["feature_alias", "group", "target_corr", "target_corr_abs"]
].copy()

top_feature_by_group_df = (
    contemporaneous_screen_df.sort_values(["group", "target_corr_abs", "feature_alias"], ascending=[True, False, True])
    .groupby("group", as_index=False)
    .first()[["group", "feature_alias", "target_corr", "target_corr_abs"]]
    .sort_values(["target_corr_abs", "group"], ascending=[False, True])
    .reset_index(drop=True)
)

group_target_summary_df = (
    contemporaneous_screen_df.groupby("group", as_index=False)
    .agg(
        feature_count=("feature_actual", "size"),
        max_target_corr_abs=("target_corr_abs", "max"),
        mean_target_corr_abs=("target_corr_abs", "mean"),
    )
    .sort_values(["max_target_corr_abs", "mean_target_corr_abs", "group"], ascending=[False, False, True])
    .reset_index(drop=True)
)

display(
    Markdown(
        f"""
### Screening contemporaneo de covariables conocidas
- **Muestras de entrenamiento:** `{len(screening_frame):,}`
- **Features con correlacion valida:** `{len(contemporaneous_screen_df):,}`
- **Cobertura temporal:** `{screening_frame.index.min()}` a `{screening_frame.index.max()}`
- **Criterio:** correlacion lineal contemporanea de cada covariable conocida con la serie objetivo.
- **Motivo para no separar por horizonte:** estas covariables se conocen al origen del pronostico, por lo que se prioriza un solo pool comun de candidatos.
- **Lectura correcta:** este ranking sirve para priorizar features; la validacion final sigue dependiendo del backtest del modelo.
"""
    )
)
display(Markdown("### Covariables conocidas con mayor correlacion absoluta al objetivo"))
display(top_target_features_df)
display(Markdown("### Mejor covariable por grupo"))
display(top_feature_by_group_df)

plot_feature_df = top_target_features_df.sort_values("target_corr_abs", ascending=True).copy()

fig = px.bar(
    plot_feature_df,
    x="target_corr",
    y="feature_alias",
    color="group",
    orientation="h",
    title=f"{SERIES_NAME} - screening contemporaneo de covariables conocidas",
    labels={"target_corr": "Correlacion con objetivo", "feature_alias": "Feature", "group": "Grupo"},
)
fig.add_vline(x=0, line_dash="dash", line_color="grey")
fig.update_layout(yaxis_title="")
display_responsive_plot(fig)

# fig = px.bar(
#     group_target_summary_df,
#     x="group",
#     y="max_target_corr_abs",
#     color="feature_count",
#     title=f"{SERIES_NAME} - senal por grupo en covariables conocidas",
#     text="max_target_corr_abs",
#     labels={
#         "group": "Grupo de features",
#         "max_target_corr_abs": "Max |corr| con objetivo",
#         "feature_count": "Cantidad de features",
#     },
# )
# fig.update_traces(texttemplate="%{text:.2f}", textposition="outside")
# fig.update_layout(yaxis_range=[0, min(1.05, group_target_summary_df["max_target_corr_abs"].max() * 1.15)])
# #display_responsive_plot(fig)

top_target_row = contemporaneous_screen_df.iloc[0]
strongest_group_row = group_target_summary_df.iloc[0]

display(
    Markdown(
        f"""
## Hallazgos multivariados para `{SERIES_NAME}`

- El screening multivariado uso **{len(df_other_data_to_test.columns)} covariables conocidas** de `df_other_data_to_test`, evaluadas solo en la **ventana de entrenamiento** para evitar filtracion del holdout.
- Como estas covariables ya estan disponibles al origen del pronostico, el analisis se hizo de forma **contemporanea** y con un **ranking unico**, en lugar de separar por horizonte.
- La covariable con señal lineal mas fuerte sobre el objetivo es **`{top_target_row["feature_alias"]}`** del grupo **`{top_target_row["group"]}`** con correlacion **{top_target_row["target_corr"]:.2f}**.
- A nivel de grupos, la familia con mayor correlacion absoluta maxima es **`{strongest_group_row["group"]}`** con **|corr| {strongest_group_row["max_target_corr_abs"]:.2f}**.
- Este screening define un pool base de exogenas conocidas; el peso final de cada variable todavia debe confirmarse con el modelo y su backtest.

### Que significa esto para el siguiente paso de modelado
- Mantener un **pool comun de covariables conocidas** para todos los horizontes, en vez de construir un ranking distinto por cada uno.
- Dejar que el modelo y el backtest determinen si el peso de cada covariable cambia por horizonte.
- Priorizar en el siguiente paso las variables y grupos mejor rankeados aqui, y contrastarlos con la seleccion de Lasso que sigue.
"""
    )
)

from sklearn.linear_model import LassoCV
from sklearn.preprocessing import StandardScaler

_lasso_frame = (
    pd.concat([series.rename("target"), df_other_data_to_test], axis=1)
    .loc[lambda df: df.index < holdout_start]
    .dropna()
)

_X = _lasso_frame[df_other_data_to_test.columns]
_y = _lasso_frame["target"]

_scaler = StandardScaler()
_X_scaled = _scaler.fit_transform(_X)

with warnings.catch_warnings():
    warnings.simplefilter("ignore")
    _lasso = LassoCV(cv=5, max_iter=10_000, random_state=42)
    _lasso.fit(_X_scaled, _y)

lasso_df = pd.DataFrame(
    {
        "feature_alias": [feature_catalog_by_actual.loc[f, "feature_alias"] for f in df_other_data_to_test.columns],
        "feature_actual": df_other_data_to_test.columns,
        "group": [feature_catalog_by_actual.loc[f, "group"] for f in df_other_data_to_test.columns],
        "coefficient": _lasso.coef_,
        "abs_coefficient": np.abs(_lasso.coef_),
    }
)
lasso_selected_df = lasso_df[lasso_df["abs_coefficient"] > 0].sort_values("abs_coefficient", ascending=False).reset_index(drop=True)

display(
    Markdown(
        f"### Selección de features con Lasso para `{SERIES_NAME}`\n"
        f"- **Muestras de entrenamiento:** {len(_y)}\n"
        f"- **Features candidatos:** {len(df_other_data_to_test.columns)}\n"
        f"- **Seleccionados (coef no-cero):** {len(lasso_selected_df)}\n"
        f"- **Alpha óptimo (CV):** {_lasso.alpha_:.4f}\n"
        f"- **R² (entrenamiento):** {_lasso.score(_X_scaled, _y):.4f}"
    )
)
display(lasso_selected_df[["feature_alias", "group", "coefficient", "abs_coefficient"]])

# ── Gráfico de barras: coeficientes Lasso seleccionados ──────────────────────
if not lasso_selected_df.empty:
    _plot_df = lasso_selected_df.head(min(25, len(lasso_selected_df))).copy()
    _plot_df = _plot_df.sort_values("abs_coefficient", ascending=True)

    fig = px.bar(
        _plot_df,
        x="coefficient",
        y="feature_alias",
        color="group",
        orientation="h",
        title=f"{SERIES_NAME} - features seleccionados por Lasso (coeficientes estandarizados)",
        labels={"coefficient": "Coeficiente Lasso (estandarizado)", "feature_alias": "Feature"},
    )
    fig.add_vline(x=0, line_dash="dash", line_color="grey")
    fig.update_layout(yaxis_title="")
    display_responsive_plot(fig)
else:
    display(Markdown("_Lasso no seleccionó ningún feature — todos los coeficientes se redujeron a cero._"))

from statistics import NormalDist

_lasso_pred = pd.Series(_lasso.predict(_X_scaled), index=_lasso_frame.index, name="prediction")
lasso_residuals_df = pd.DataFrame(
    {
        "target": _y,
        "prediction": _lasso_pred,
    }
)
lasso_residuals_df["residual"] = lasso_residuals_df["target"] - lasso_residuals_df["prediction"]
lasso_residuals_df["abs_residual"] = lasso_residuals_df["residual"].abs()

_residual_values = lasso_residuals_df["residual"].to_numpy()
_residual_mean = float(np.mean(_residual_values))
_residual_std = float(np.std(_residual_values, ddof=1)) if len(_residual_values) > 1 else 0.0
_residual_rmse = rmse_np(lasso_residuals_df["target"].to_numpy(), lasso_residuals_df["prediction"].to_numpy())
_residual_mae = mean_absolute_error_np(lasso_residuals_df["target"].to_numpy(), lasso_residuals_df["prediction"].to_numpy())
_residual_autocorr_1 = lasso_residuals_df["residual"].autocorr(lag=1)

_sorted_residuals = np.sort(_residual_values)
_qq_probabilities = (np.arange(1, len(_sorted_residuals) + 1) - 0.5) / len(_sorted_residuals)
_theoretical_quantiles = np.array([NormalDist().inv_cdf(prob) for prob in _qq_probabilities])
if _residual_std > 0:
    _sample_quantiles = (_sorted_residuals - _residual_mean) / _residual_std
else:
    _sample_quantiles = np.zeros_like(_sorted_residuals)

_qq_min = float(min(_theoretical_quantiles.min(), _sample_quantiles.min()))
_qq_max = float(max(_theoretical_quantiles.max(), _sample_quantiles.max()))

display(
    Markdown(
        f"""
### Diagnóstico residual del Lasso
- **Muestras analizadas:** `{len(lasso_residuals_df):,}` (entrenamiento)
- **MAE:** `{_residual_mae:.2f}`
- **RMSE:** `{_residual_rmse:.2f}`
- **Media del residual:** `{_residual_mean:.2f}`
- **Desv. estándar del residual:** `{_residual_std:.2f}`
- **Autocorrelación lag 1:** `{_residual_autocorr_1:.2f}`
"""
    )
)

fig = make_subplots(
    rows=2,
    cols=2,
    subplot_titles=(
        "Residuales en el tiempo",
        "Residuales vs predicción",
        "Distribución de residuales",
        "Q-Q plot de residuales estandarizados",
    ),
)

fig.add_trace(
    go.Scattergl(
        x=lasso_residuals_df.index,
        y=lasso_residuals_df["residual"],
        mode="lines",
        name="Residual",
        line=dict(color="#1f77b4", width=1.2),
        showlegend=False,
    ),
    row=1,
    col=1,
)
fig.add_hline(y=0, line_dash="dash", line_color="grey", row=1, col=1)

fig.add_trace(
    go.Scattergl(
        x=lasso_residuals_df["prediction"],
        y=lasso_residuals_df["residual"],
        mode="markers",
        name="Residual vs pred",
        marker=dict(color="#ff7f0e", size=5, opacity=0.55),
        showlegend=False,
    ),
    row=1,
    col=2,
)
fig.add_hline(y=0, line_dash="dash", line_color="grey", row=1, col=2)

fig.add_trace(
    go.Histogram(
        x=lasso_residuals_df["residual"],
        nbinsx=40,
        marker_color="#00cc96",
        opacity=0.85,
        showlegend=False,
        name="Hist residual",
    ),
    row=2,
    col=1,
)

fig.add_trace(
    go.Scattergl(
        x=_theoretical_quantiles,
        y=_sample_quantiles,
        mode="markers",
        marker=dict(color="#ab63fa", size=5, opacity=0.6),
        name="Q-Q",
        showlegend=False,
    ),
    row=2,
    col=2,
)
fig.add_trace(
    go.Scattergl(
        x=[_qq_min, _qq_max],
        y=[_qq_min, _qq_max],
        mode="lines",
        line=dict(color="grey", dash="dash"),
        name="45°",
        showlegend=False,
    ),
    row=2,
    col=2,
)

fig.update_xaxes(title_text="Fecha", row=1, col=1)
fig.update_yaxes(title_text="Residual", row=1, col=1)
fig.update_xaxes(title_text="Predicción Lasso", row=1, col=2)
fig.update_yaxes(title_text="Residual", row=1, col=2)
fig.update_xaxes(title_text="Residual", row=2, col=1)
fig.update_yaxes(title_text="Frecuencia", row=2, col=1)
fig.update_xaxes(title_text="Cuantil teórico normal", row=2, col=2)
fig.update_yaxes(title_text="Cuantil residual estandarizado", row=2, col=2)
fig.update_layout(
    title=f"{SERIES_NAME} - diagnostico residual del modelo Lasso",
    bargap=0.05,
    height=850,
)
display_responsive_plot(fig)

	value	date	hour	day_of_week	month	year
2023-05-06 17:00:00-04:00	212.41	2023-05-06	17	5	5	2023
2023-05-06 18:00:00-04:00	247.31	2023-05-06	18	5	5	2023
2023-05-06 19:00:00-04:00	230.83	2023-05-06	19	5	5	2023
2023-05-06 20:00:00-04:00	192.14	2023-05-06	20	5	5	2023
2023-05-06 21:00:00-04:00	192.14	2023-05-06	21	5	5	2023

	value
2024-02-05 15:00:00-03:00	371.66
2024-02-05 14:00:00-03:00	371.66
2024-02-05 16:00:00-03:00	371.66
2024-02-01 15:00:00-03:00	349.73
2024-01-25 15:00:00-03:00	338.87
2025-06-17 12:00:00-04:00	338.65
2024-02-14 15:00:00-03:00	335.49
2024-02-14 13:00:00-03:00	335.49
2024-02-14 14:00:00-03:00	335.49
2024-02-14 12:00:00-03:00	335.49
2024-02-14 16:00:00-03:00	335.49
2024-02-15 15:00:00-03:00	335.42
2024-02-15 16:00:00-03:00	335.42
2024-02-15 14:00:00-03:00	334.85
2025-06-17 13:00:00-04:00	332.31

	lag	autocorr
0	1	0.90
1	24	0.52
2	168	0.42
3	336	0.33

	horizon	baseline	wape_pct	mae	rmse
0	1	current_value	20.38	18.45	33.82
1	2	current_value	33.07	29.95	48.53
2	3	current_value	41.22	37.34	56.18
3	4	same_hour_previous_day	46.79	42.39	66.89
4	5	same_hour_previous_day	46.80	42.41	66.91
5	6	same_hour_previous_day	46.80	42.42	66.92
6	7	same_hour_previous_day	46.81	42.44	66.94
7	8	same_hour_previous_day	46.81	42.45	66.95
8	9	same_hour_previous_day	46.81	42.47	66.97
9	10	same_hour_previous_day	46.79	42.48	66.98
10	11	same_hour_previous_day	46.76	42.47	66.99
11	12	same_hour_previous_day	46.74	42.47	66.99
12	13	same_hour_previous_day	46.71	42.46	66.99
13	14	same_hour_previous_day	46.68	42.45	67.00
14	15	same_hour_previous_day	46.65	42.45	67.00
15	16	same_hour_previous_day	46.62	42.44	67.01
16	17	same_hour_previous_day	46.59	42.43	67.01
17	18	same_hour_previous_day	46.56	42.43	67.02
18	19	same_hour_previous_day	46.54	42.43	67.03
19	20	same_hour_previous_day	46.55	42.44	67.04
20	21	same_hour_previous_day	46.56	42.46	67.06
21	22	same_hour_previous_day	46.58	42.48	67.07
22	23	same_hour_previous_day	46.58	42.49	67.09
23	24	current_value	46.60	42.51	67.10

	group	feature_count
4	generacion_prg	10
3	disponibilidad	8
5	residuals_prg	8
6	residuals_ratings	8
1	curtailments_prg	3
0	cmg_prg	2
2	demanda_prg	2
7	tco	2
8	weather	2

CMG EDA - P.MONTT220¶

Serie cargada¶

Resumen de timestamps y variable objetivo¶

Detección de outliers¶

Primeros hallazgos para `P.MONTT220`¶

Qué parece útil para pronóstico¶

Siguiente paso - pronosticabilidad multi-horizonte (t+1 a t+24)¶

Diseño del holdout¶

Como leer WAPE¶

Hallazgos del siguiente paso para `P.MONTT220`¶

Qué significa esto antes de agregar otras series¶

EDA multivariado - covariables conocidas al momento de la predicción¶

Covariables conocidas candidatas¶

Screening contemporaneo de covariables conocidas¶

Covariables conocidas con mayor correlacion absoluta al objetivo¶

Mejor covariable por grupo¶

Hallazgos multivariados para `P.MONTT220`¶

Que significa esto para el siguiente paso de modelado¶

Regresión lineal Lasso: serie de interés vs todas las series de `df_other_data`¶

Selección de features con Lasso para `P.MONTT220`¶

Diagnóstico residual del Lasso¶

Conclusiones¶

1. Análisis univariado¶

2. Análisis multivariado¶

3. Implicancias para modelado¶

	metric	value
0	timezone	America/Santiago
1	inferred_frequency	h
2	start	2023-05-06 17:00:00-04:00
3	end	2026-02-06 00:00:00-03:00
4	rows	24151
5	duplicate_timestamps	0
6	missing_timestamps	0
7	missing_values	0
8	zero_share	0.25
9	negative_share	0.00

	baseline	mean_wape_pct
1	same_hour_previous_day	46.69
2	same_hour_previous_week	56.31
0	current_value	59.55

	feature_alias	feature_actual	group
0	Charrua220_PRG_latest	Charrua220_PRG_latest	cmg_prg
1	PMontt220_PRG_latest	PMontt220_PRG_latest	cmg_prg
2	GD2_GCurt_EO	GD2_GCurt_EO	curtailments_prg
3	GDT_GCurt_EO	GDT_GCurt_EO	curtailments_prg
4	GDT_GCurt_PV	GDT_GCurt_PV	curtailments_prg
5	GD2_Dem	GD2_Dem	demanda_prg
6	GDT_Dem	GDT_Dem	demanda_prg
7	GD2_GRat_EO	GD2_GRat_EO	disponibilidad
8	GD2_GRat_HE	GD2_GRat_HE	disponibilidad
9	GD2_GRat_HP	GD2_GRat_HP	disponibilidad
10	GDT_GRat_EO	GDT_GRat_EO	disponibilidad
11	GDT_GRat_HE	GDT_GRat_HE	disponibilidad
12	GDT_GRat_HP	GDT_GRat_HP	disponibilidad
13	GDT_GRat_PV	GDT_GRat_PV	disponibilidad
14	GDT_GRat_TER	GDT_GRat_TER	disponibilidad
15	GD2_GP_EO	GD2_GP_EO	generacion_prg
16	GD2_GP_HE	GD2_GP_HE	generacion_prg
17	GD2_GP_HP	GD2_GP_HP	generacion_prg
18	GD2_GP_PV	GD2_GP_PV	generacion_prg
19	GD2_GP_TER	GD2_GP_TER	generacion_prg
20	GDT_GP_EO	GDT_GP_EO	generacion_prg
21	GDT_GP_HE	GDT_GP_HE	generacion_prg
22	GDT_GP_HP	GDT_GP_HP	generacion_prg
23	GDT_GP_PV	GDT_GP_PV	generacion_prg
24	GDT_GP_TER	GDT_GP_TER	generacion_prg
25	GD2_D-GP_+HE	GD2_D-GP_+HE	residuals_prg
26	GD2_D-GP_+HP	GD2_D-GP_+HP	residuals_prg
27	GD2_D-GP_+TER	GD2_D-GP_+TER	residuals_prg
28	GD2_D-GP_PV+EO	GD2_D-GP_PV+EO	residuals_prg
29	GDT_D-GP_+HE	GDT_D-GP_+HE	residuals_prg
30	GDT_D-GP_+HP	GDT_D-GP_+HP	residuals_prg
31	GDT_D-GP_+TER	GDT_D-GP_+TER	residuals_prg
32	GDT_D-GP_PV+EO	GDT_D-GP_PV+EO	residuals_prg
33	GD2_D-GRat_+HE	GD2_D-GRat_+HE	residuals_ratings
34	GD2_D-GRat_+HP	GD2_D-GRat_+HP	residuals_ratings
35	GD2_D-GRat_+TER	GD2_D-GRat_+TER	residuals_ratings
36	GD2_D-GRat_PV+EO	GD2_D-GRat_PV+EO	residuals_ratings
37	GDT_D-GRat_+HE	GDT_D-GRat_+HE	residuals_ratings
38	GDT_D-GRat_+HP	GDT_D-GRat_+HP	residuals_ratings
39	GDT_D-GRat_+TER	GDT_D-GRat_+TER	residuals_ratings
40	GDT_D-GRat_PV+EO	GDT_D-GRat_PV+EO	residuals_ratings
41	TCO_HE	TCO_HE	tco
42	TCO_TER	TCO_TER	tco
43	CIRUELOS220_rad	CIRUELOS220_rad	weather
44	CIRUELOS220_temp	CIRUELOS220_temp	weather

	feature_alias	group	target_corr	target_corr_abs
0	PMontt220_PRG_latest	cmg_prg	0.60	0.60
1	TCO_HE	tco	0.40	0.40
2	Charrua220_PRG_latest	cmg_prg	0.39	0.39
3	GD2_D-GP_+HE	residuals_prg	0.37	0.37
4	GD2_Dem	demanda_prg	0.31	0.31
5	GDT_Dem	demanda_prg	0.29	0.29
6	GD2_GP_TER	generacion_prg	0.29	0.29
7	GDT_GP_TER	generacion_prg	0.28	0.28
8	GDT_D-GP_+HE	residuals_prg	0.28	0.28
9	GD2_D-GP_+HP	residuals_prg	0.21	0.21
10	GD2_D-GP_PV+EO	residuals_prg	0.20	0.20
11	GD2_GRat_EO	disponibilidad	-0.20	0.20

	feature_alias	group	coefficient	abs_coefficient
0	PMontt220_PRG_latest	cmg_prg	28.20	28.20
1	GD2_D-GP_+HE	residuals_prg	11.95	11.95
2	GD2_Dem	demanda_prg	5.02	5.02
3	CIRUELOS220_temp	weather	4.77	4.77
4	GD2_GP_EO	generacion_prg	-2.07	2.07
5	TCO_HE	tco	2.04	2.04
6	GDT_Dem	demanda_prg	0.95	0.95
7	GD2_GRat_EO	disponibilidad	-0.72	0.72
8	GDT_GP_HP	generacion_prg	0.00	0.00

CMG EDA - P.MONTT220¶

Serie cargada¶

Resumen de timestamps y variable objetivo¶

Detección de outliers¶

Primeros hallazgos para P.MONTT220¶

Qué parece útil para pronóstico¶

Siguiente paso - pronosticabilidad multi-horizonte (t+1 a t+24)¶

Diseño del holdout¶

Como leer WAPE¶

Hallazgos del siguiente paso para P.MONTT220¶

Qué significa esto antes de agregar otras series¶

EDA multivariado - covariables conocidas al momento de la predicción¶

Covariables conocidas candidatas¶

Screening contemporaneo de covariables conocidas¶

Covariables conocidas con mayor correlacion absoluta al objetivo¶

Mejor covariable por grupo¶

Hallazgos multivariados para P.MONTT220¶

Que significa esto para el siguiente paso de modelado¶

Regresión lineal Lasso: serie de interés vs todas las series de df_other_data¶

Selección de features con Lasso para P.MONTT220¶

Diagnóstico residual del Lasso¶

Conclusiones¶

1. Análisis univariado¶

2. Análisis multivariado¶

3. Implicancias para modelado¶

Primeros hallazgos para `P.MONTT220`¶

Hallazgos del siguiente paso para `P.MONTT220`¶

Hallazgos multivariados para `P.MONTT220`¶

Regresión lineal Lasso: serie de interés vs todas las series de `df_other_data`¶

Selección de features con Lasso para `P.MONTT220`¶