Pronóstico de la velocidad del viento a 24 horas en el Caribe colombiano

import json
import os
import re
import unicodedata
import warnings
from pathlib import Path

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from scipy import stats
from matplotlib.cm import ScalarMappable
from matplotlib.collections import PolyCollection
from matplotlib.colors import Normalize
from scipy.spatial import ConvexHull
from sklearn.cluster import DBSCAN, KMeans
from sklearn.compose import TransformedTargetRegressor
from sklearn.decomposition import PCA
from sklearn.dummy import DummyRegressor
from sklearn.ensemble import IsolationForest
from sklearn.feature_selection import mutual_info_regression
from sklearn.linear_model import LinearRegression
from sklearn.metrics import (mean_absolute_error, mean_absolute_percentage_error, mean_squared_error,
                             r2_score, silhouette_score)
from sklearn.model_selection import GroupKFold
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.svm import LinearSVR
from statsmodels.stats.diagnostic import het_breuschpagan
from statsmodels.stats.multitest import multipletests
from statsmodels.stats.outliers_influence import variance_inflation_factor
from statsmodels.tsa.seasonal import STL
from statsmodels.tsa.stattools import acf, adfuller, kpss, pacf
from IPython.display import display
warnings.filterwarnings("ignore")
SEED = 42
np.random.seed(SEED)
RNG = np.random.default_rng(SEED)
sns.set_theme(style="whitegrid", context="notebook")
plt.rcParams["figure.dpi"] = 110
pd.set_option("display.float_format", "{:.3f}".format)
pd.set_option("display.max_columns", 30)

# carpeta con el panel y las auditorías generados por scripts/procesar_todo.py; por defecto, la carpeta
# hermana ideam_viento/ (ver "Reproducibilidad"); se puede cambiar con la variable de entorno DATOS_PANEL
DATOS = Path(os.environ.get("DATOS_PANEL", "../ideam_viento"))
TEST_INI = pd.Timestamp("2025-01-01")  # conjunto de prueba reservado: todo 2025
H = 24                                 # horizonte de pronóstico (horas)
K_VECINOS = 4                          # vecinos para el rezago espacial y las matrices de pesos


def muestra(x, n):
    """Muestra aleatoria reproducible de hasta n filas (todas si hay menos)."""
    return x.sample(min(n, len(x)), random_state=SEED)


def haversine_matriz(lat, lon):
    la, lo = np.deg2rad(lat), np.deg2rad(lon)
    a = np.sin((la[:, None] - la) / 2) ** 2 + np.cos(la[:, None]) * np.cos(la) * np.sin((lo[:, None] - lo) / 2) ** 2
    return 2 * 6371 * np.arcsin(np.sqrt(a))


def pesos_knn(dist, k=K_VECINOS):
    d = dist.copy()
    np.fill_diagonal(d, np.inf)
    W = np.zeros_like(d)
    idx = np.argsort(d, axis=1)[:, :k]
    W[np.arange(len(d))[:, None], idx] = 1 / k
    return W


def moran(z, W):
    z = np.asarray(z, float) - np.mean(z)
    return len(z) / W.sum() * (z @ W @ z) / (z @ z)


# ---------- mapa base: límites departamentales (geoBoundaries, OpenStreetMap, licencia ODbL 1.0) ----------
def _normalizar(t):
    return unicodedata.normalize("NFKD", t).encode("ascii", "ignore").decode().upper()


def _anillos(feature):
    g = feature["geometry"]
    partes = [g["coordinates"]] if g["type"] == "Polygon" else g["coordinates"]
    return [np.array(parte[0]) for parte in partes]  # anillo exterior de cada polígono


GEO = json.load(open("mapas/colombia_departamentos.geojson", encoding="utf-8"))
POLIGONOS = {_normalizar(f["properties"]["shapeName"]): _anillos(f) for f in GEO["features"]}
CARIBE = ["ATLANTICO", "BOLIVAR", "CESAR", "CORDOBA", "LA GUAJIRA", "MAGDALENA", "SUCRE"]


def mapa_base(ax, valores=None, cmap="Blues", etiqueta=None, nombres=False):
    """Dibuja los departamentos; con `valores` (dict depto → número) dibuja una coropleta."""
    norma = Normalize(min(valores.values()), max(valores.values())) if valores else None
    for d, anillos in POLIGONOS.items():
        caribe = d in CARIBE
        if valores and d in valores:
            color = plt.get_cmap(cmap)(norma(valores[d]))
        else:
            color = "#ececec" if caribe else "#f8f8f8"
        ax.add_collection(PolyCollection(anillos, facecolor=color, edgecolor="gray",
                                         linewidth=.7 if caribe else .3, zorder=0))
        if nombres and caribe:
            x, y = np.concatenate(anillos).mean(axis=0)
            ax.annotate(d.title(), (x, y), fontsize=7, color="dimgray", ha="center", zorder=3)
    ax.set_xlim(-76.7, -70.9); ax.set_ylim(7.2, 12.7)
    ax.set_aspect(1 / np.cos(np.deg2rad(10)))  # corrige la escala de la longitud a 10 °N
    ax.set_facecolor("white")  # fuera de Colombia (mar y Venezuela) queda en blanco: el archivo solo trae Colombia
    ax.grid(False)
    if valores:
        plt.colorbar(ScalarMappable(norm=norma, cmap=cmap), ax=ax, label=etiqueta, shrink=.8)


def moran_perm(z, W, n=999, rng=RNG):
    I = moran(z, W)
    sim = np.array([moran(rng.permutation(z), W) for _ in range(n)])
    return I, (np.sum(sim >= I) + 1) / (n + 1)

P = pd.read_csv(DATOS / "panel_multivariado.csv", dtype={"CodigoEstacion": str}, parse_dates=["hora"])
P = P.rename(columns={"CodigoEstacion": "est"})
P["anio"] = P.hora.dt.year
print(f"Panel estación-hora: {len(P):,} filas | {P.est.nunique()} estaciones | "
      f"{P.hora.min():%Y-%m-%d} a {P.hora.max():%Y-%m-%d}")
P.head()

dicc = pd.DataFrame([
    ("est", "categórica (ID)", "—", "Código de la estación en el catálogo del IDEAM (10 dígitos)"),
    ("nombre, depto, municipio", "categórica", "—", "Nombre y ubicación administrativa de la estación"),
    ("lat, lon", "numérica continua", "grados (EPSG:4326)", "Coordenadas de la estación"),
    ("hora", "fecha-hora", "hora local", "Inicio de la hora de agregación"),
    ("vel", "numérica continua", "m/s", "Velocidad media del viento en la hora (objetivo a t+24 h)"),
    ("dir", "numérica circular", "grados (0–360)", "Dirección de la que sopla el viento, media vectorial de la hora"),
    ("dir_sin, dir_cos", "numérica continua", "adimensional",
     "Seno y coseno medios de la dirección en la hora: componentes del vector medio de dirección. "
     "No es unitario: su norma es dir_constancia"),
    ("dir_constancia", "numérica continua", "0–1", "Longitud del vector medio: 1 = dirección fija en la hora"),
    ("temp", "numérica continua", "°C", "Temperatura del aire a 2 m"),
    ("pres", "numérica continua", "hPa", "Presión atmosférica en la estación (no reducida al nivel del mar)"),
    ("categoria", "categórica", "—", "Categoría de la estación según el catálogo"),
    ("altitud", "numérica continua", "m s. n. m.", "Altitud de la estación según el catálogo"),
], columns=["variable", "tipo", "unidad", "significado"])
dicc["% de filas del panel con dato"] = [100.0, 100.0, 100.0, 100.0,
                                         *(P[["vel", "dir", "dir_sin", "dir_constancia", "temp", "pres",
                                              "categoria", "altitud"]].notna().mean() * 100).round(1)]
dicc

completas = P[["vel", "dir", "temp", "pres"]].notna().all(axis=1)
por_anio = pd.DataFrame({
    "estación-horas": P.groupby("anio").size(),
    "con las 4 variables": completas.groupby(P.anio).sum(),
    "estaciones": P.groupby("anio").est.nunique(),
    "estaciones con las 4": P[completas].groupby("anio").est.nunique(),
})
por_anio.loc["total"] = [len(P), completas.sum(), P.est.nunique(), P[completas].est.nunique()]
n_pred = 26  # predictoras del modelo final (sección 2.9)
print(f"Filas con las 4 variables: {completas.sum():,} | predictoras del modelo: {n_pred} | "
      f"n/p = {completas.sum() / n_pred:,.0f}")
print(f"Rango de la variable objetivo: {P.vel.min():.2f} a {P.vel.max():.2f} m/s")
por_anio

filas = []
for f in sorted((DATOS / "reportes").glob("*.txt")):
    if "_causal" in f.stem:  # los reportes de la limpieza causal (y su sensibilidad) no son de la depuración original
        continue
    var, anio = f.stem.split("_")
    t = f.read_text(encoding="utf-8")
    num = lambda patron: int(re.search(patron, t).group(1).replace(",", "")) if re.search(patron, t) else 0
    filas.append({"variable": var, "año": int(anio),
                  "filas crudas": num(r"Filas: ([\d,]+)"),
                  "fuera del Caribe continental": num(r"descartadas\): ([\d,]+)"),
                  "duplicadas": num(r"instante: ([\d,]+)"),
                  "duplicadas con valores distintos": num(r"valores distintos: ([\d,]+)"),
                  "fuera de rango físico": num(r"Lecturas fuera de \[[^\]]+\] \S+: ([\d,]+)"),
                  "pegadas o con picos": num(r"enmascaradas\): ([\d,]+)")})
calidad = pd.DataFrame(filas)
resumen_calidad = calidad.groupby("variable").sum(numeric_only=True).drop(columns="año")
resumen_calidad.loc["total"] = resumen_calidad.sum()
resumen_calidad

aud = pd.concat([pd.read_csv(DATOS / f"auditoria_{v}.csv", dtype={"CodigoEstacion": str}).assign(variable=v)
                 for v in ["velocidad", "direccion", "temperatura", "presion"]])
excl = aud[aud.excluida][["variable", "CodigoEstacion", "NombreEstacion", "anio", "frac_fuera", "max_v", "motivo"]]
print(f"Estación-año excluidas por sensor defectuoso: {len(excl)}")
lect_total = aud.lecturas.sum()
pegadas = aud.frac_pegado * aud.lecturas
retro = pd.DataFrame({
    "lecturas": [aud[aud.excluida].lecturas.sum(), pegadas.sum(), pegadas[aud.excluida].sum()],
}, index=["en estación-años excluidas (proporción anual de anomalías)",
          "en tramos pegados o picos de presión (tramo completo / mediana anual)",
          "contadas en ambos grupos (pegadas dentro de estación-años excluidos)"])
retro.loc["unión sin doble conteo"] = retro.lecturas.iloc[0] + retro.lecturas.iloc[1] - retro.lecturas.iloc[2]
retro["% de las lecturas auditadas"] = 100 * retro.lecturas / lect_total
print(f"Lecturas auditadas (después de quitar duplicados y filas fuera del Caribe): {lect_total:,.0f}")
display(retro.round(3))
altas = aud[aud.variable == "velocidad"].assign(n_fuera=lambda d: d.frac_fuera * d.lecturas)
top2 = altas.groupby("NombreEstacion").n_fuera.sum().nlargest(2)
print(f"Lecturas de velocidad fuera de `[0, 25] m/s`: {altas.n_fuera.sum():,.0f}; "
      f"{100 * top2.sum() / altas.n_fuera.sum():.0f} % vienen de {', '.join(top2.index.str.strip())}")
umbral = pd.read_csv("resultados/umbral_velocidad.csv")  # generado por scripts/revisar_umbral.py con los CSV crudos
print("Lecturas crudas de velocidad que eliminaría cada umbral (revisión sobre los archivos originales):")
display(umbral.set_index(["umbral (m/s)", "estaciones"]))
excl.sort_values(["variable", "anio"])

coord_est = P.groupby("est")[["lat", "lon"]].median().round(4)
grupos_coord = coord_est[coord_est.duplicated(["lat", "lon"], keep=False)].groupby(["lat", "lon"]).groups
for (la_, lo_), ests_ in grupos_coord.items():
    ests_ = list(ests_)
    nombres_ = P[P.est.isin(ests_)].groupby("est").nombre.first()
    ancho_ = P[P.est.isin(ests_)].pivot(index="hora", columns="est", values="vel")
    comunes = ancho_.dropna()
    r_ = comunes.corr().iloc[0, 1] if len(comunes) > 24 else np.nan
    print(f"Coordenadas ({la_}, {lo_}): {', '.join(nombres_.str.strip() + ' [' + nombres_.index + ']')}")
    print(f"  horas con dato: {ancho_.notna().sum().to_dict()} | horas en común: {len(comunes)} | "
          f"correlación en horas comunes: {r_:.2f}" if len(comunes) > 24 else
          f"  horas con dato: {ancho_.notna().sum().to_dict()} | no comparten horas: se registraron en periodos distintos")

P["mes"] = P.hora.dt.to_period("M")
fig, axs = plt.subplots(1, 3, figsize=(16, 6), sharey=True)
orden = P.groupby("est").vel.size().sort_values(ascending=False).index
for ax, v in zip(axs, ["dir", "temp", "pres"]):
    m = P.groupby(["est", "mes"])[v].apply(lambda s: s.notna().mean()).unstack().reindex(orden)
    sns.heatmap(m, ax=ax, cmap="viridis", vmin=0, vmax=1, cbar=ax is axs[-1],
                cbar_kws={"label": "fracción de horas con dato"}, xticklabels=12, yticklabels=False)
    ax.set(title=f"{v}: disponibilidad (estación × mes)", xlabel="mes", ylabel="estaciones" if ax is axs[0] else "")
fig.tight_layout(); plt.show()

mec = []
for v in ["dir", "temp", "pres"]:
    falta = P[v].isna()
    por_est = falta.groupby(P.est).mean()
    estructural = P.est.map(por_est).gt(0.95) & falta  # faltantes en estaciones sin ese sensor
    ks = stats.ks_2samp(muestra(P.vel[falta], 20000), muestra(P.vel[~falta], 20000))
    mec.append({"variable": v, "% faltante": 100 * falta.mean(),
                "% de faltantes en estaciones sin el sensor (>95 % faltante)": 100 * estructural.sum() / falta.sum(),
                "vel media si falta": P.vel[falta].mean(), "vel media si hay dato": P.vel[~falta].mean(),
                "KS (vel | falta vs. no falta)": ks.statistic, "p": ks.pvalue})
pd.DataFrame(mec).set_index("variable")

est = P.groupby("est").agg(nombre=("nombre", "first"), depto=("depto", "first"), lat=("lat", "median"),
                           lon=("lon", "median"), altitud=("altitud", "first"), categoria=("categoria", "first"),
                           horas=("vel", "size"), vel_media=("vel", "mean"))
fig, axs = plt.subplots(1, 2, figsize=(14, 4))
est.depto.value_counts().plot.bar(ax=axs[0], color="steelblue")
axs[0].set(title="Estaciones por departamento", ylabel="estaciones")
est.altitud.plot.hist(bins=30, ax=axs[1], color="steelblue")
axs[1].set(title="Altitud de las estaciones", xlabel="m s. n. m.")
fig.tight_layout(); plt.show()
print("Horas de datos por departamento (%):")
print((P.depto.value_counts(normalize=True) * 100).round(1).to_string())

PT = P[P.hora < TEST_INI].copy()  # EDA solo con entrenamiento
print(f"Entrenamiento (EDA): {len(PT):,} estación-horas 2020–2024 | "
      f"prueba reservada: {(P.hora >= TEST_INI).sum():,} estación-horas 2025")

coords = est[["lat", "lon"]]
DIST = haversine_matriz(coords.lat.to_numpy(), coords.lon.to_numpy())
d_ = DIST.copy(); np.fill_diagonal(d_, np.inf)
vecinos = {e: coords.index[np.argsort(d_[i])[:K_VECINOS]].tolist() for i, e in enumerate(coords.index)}
ancho_vel = P.pivot(index="hora", columns="est", values="vel")
vel_vecinos = pd.DataFrame({e: ancho_vel[v].mean(axis=1) for e, v in vecinos.items()})

# rejilla horaria completa por estación
G = (P.set_index("hora").groupby("est")[["vel", "dir_sin", "dir_cos", "dir_constancia", "temp", "pres", "altitud"]]
     .apply(lambda g: g.asfreq("h")))


def construir(g, vv=None):
    vv = vel_vecinos if vv is None else vv
    e = g.name
    g = g.droplevel(0)
    X = pd.DataFrame(index=g.index)
    v = g.vel
    for L in [0, 1, 2, 3, 6, 12, 23]:
        X[f"vel_l{L}"] = v.shift(L)
    X["vel_m24"] = v.rolling(24, min_periods=18).mean()
    # aproximación de las componentes este (u) y norte (w) del viento, convención "desde": velocidad media
    # × vector medio de dirección. Como ese vector tiene norma = constancia, equivale a ponderar por ella
    X["u"], X["w"] = -v * g.dir_sin, -v * g.dir_cos
    X["u_m24"] = X.u.rolling(24, min_periods=18).mean()
    X["w_m24"] = X.w.rolling(24, min_periods=18).mean()
    X["dir_constancia"] = g.dir_constancia
    X["temp"] = g.temp
    X["temp_d24"] = g.temp - g.temp.shift(24)
    X["temp_amp24"] = g.temp.rolling(24, min_periods=18).max() - g.temp.rolling(24, min_periods=18).min()
    X["pres"] = g.pres
    X["pres_d3"], X["pres_d24"] = g.pres - g.pres.shift(3), g.pres - g.pres.shift(24)
    X["vel_vecinos"] = vv[e].reindex(g.index)
    X["vel_vecinos_m24"] = X.vel_vecinos.rolling(24, min_periods=18).mean()
    t_y = g.index + pd.Timedelta(H, "h")
    X["hy_sin"], X["hy_cos"] = np.sin(2 * np.pi * t_y.hour / 24), np.cos(2 * np.pi * t_y.hour / 24)
    X["doy_sin"] = np.sin(2 * np.pi * t_y.dayofyear / 365.25)
    X["doy_cos"] = np.cos(2 * np.pi * t_y.dayofyear / 365.25)
    X["altitud"] = g.altitud.ffill().bfill()
    X["hora_y"], X["mes_y"] = t_y.hour, t_y.month
    X["y"] = v.shift(-H)
    X["t_y"] = t_y
    return X


D = G.groupby(level="est", group_keys=True).apply(construir).reset_index().rename(columns={"hora": "t"})
F_UNI = ["vel_l0", "vel_l1", "vel_l2", "vel_l3", "vel_l6", "vel_l12", "vel_l23", "vel_m24", "hy_sin", "hy_cos"]
F_MET = F_UNI + ["u", "w", "u_m24", "w_m24", "dir_constancia", "temp", "temp_d24", "temp_amp24",
                 "pres", "pres_d3", "pres_d24", "doy_sin", "doy_cos", "altitud"]
F_ESP = F_MET + ["vel_vecinos", "vel_vecinos_m24"]
D_todo = D                                   # antes de exigir casos completos (sección 3.2)
D = D.dropna(subset=F_ESP + ["y"]).reset_index(drop=True)
TR = D[D.t_y < TEST_INI]                      # el objetivo de entrenamiento siempre es anterior a 2025
TE = D[(D.t >= TEST_INI)]                      # prueba: pronósticos emitidos en 2025
frontera = len(D) - len(TR) - len(TE)
print(f"Conjunto de modelado (casos completos): {len(D):,} filas | {D.est.nunique()} estaciones")
print(f"  {frontera} filas de las últimas 24 h de 2024 quedan fuera de ambos conjuntos: su predictora es de 2024 "
      "pero su objetivo cae en 2025 (actúan como separación entre entrenamiento y prueba)")
print(f"  entrenamiento: {len(TR):,} filas ({TR.t.min():%Y-%m-%d} a {TR.t_y.max():%Y-%m-%d})")
print(f"  prueba:        {len(TE):,} filas ({TE.t.min():%Y-%m-%d} a {TE.t_y.max():%Y-%m-%d})")

y = PT.vel.dropna()
desc = y.describe(percentiles=[.01, .05, .25, .5, .75, .95, .99])
print(desc.round(3).to_string())
print(f"Asimetría: {stats.skew(y):.3f} | curtosis (exceso): {stats.kurtosis(y):.3f} | "
      f"horas con 0 m/s: {100 * (y == 0).mean():.2f} %")
_, lam = stats.yeojohnson(muestra(y, 50000))
print(f"λ de Yeo-Johnson: {lam:.3f} (1 = sin transformación, 0 ≈ logaritmo)")

fig, axs = plt.subplots(1, 3, figsize=(16, 4))
sns.histplot(y, bins=100, ax=axs[0], color="steelblue")
axs[0].set(title="Velocidad horaria del viento (train)", xlabel="m/s")
sns.histplot(np.log1p(y), bins=100, ax=axs[1], color="darkorange")
axs[1].set(title="log(1 + velocidad)", xlabel="log(1 + m/s)")
stats.probplot(muestra(y, 20000), dist="norm", plot=axs[2])
axs[2].set_title("QQ-plot (muestra de 20 000)")
fig.tight_layout(); plt.show()

fig, axs = plt.subplots(1, 3, figsize=(17, 4))
ciclo = PT.groupby(PT.hora.dt.hour).vel.quantile([.25, .5, .75]).unstack()
axs[0].fill_between(ciclo.index, ciclo[.25], ciclo[.75], alpha=.3, label="rango intercuartil")
axs[0].plot(ciclo.index, ciclo[.5], "o-", label="mediana")
axs[0].set(title="Ciclo diario", xlabel="hora local", ylabel="m/s"); axs[0].legend()
mes = PT.groupby(PT.hora.dt.month).vel.quantile([.25, .5, .75]).unstack()
axs[1].fill_between(mes.index, mes[.25], mes[.75], alpha=.3)
axs[1].plot(mes.index, mes[.5], "o-")
axs[1].set(title="Ciclo anual", xlabel="mes", ylabel="m/s", xticks=range(1, 13))
mapa_base(axs[2])
sc = axs[2].scatter(est.lon, est.lat, c=PT.groupby("est").vel.mean().reindex(est.index), s=40, cmap="viridis",
                    edgecolor="k", lw=.3)
plt.colorbar(sc, ax=axs[2], label="velocidad media (m/s)")
axs[2].set(title="Velocidad media por estación", xlabel="longitud", ylabel="latitud")
fig.tight_layout(); plt.show()

num = PT[["vel", "temp", "pres", "dir_constancia"]]


def resumen_num(s):
    s = s.dropna()
    q1, q3 = s.quantile([.25, .75])
    iqr = q3 - q1
    return pd.Series({"n": len(s), "media": s.mean(), "mediana": s.median(), "desv. est.": s.std(),
                      "p1": s.quantile(.01), "p99": s.quantile(.99), "asimetría": stats.skew(s),
                      "curtosis": stats.kurtosis(s),
                      "% atípicos IQR": 100 * ((s < q1 - 1.5 * iqr) | (s > q3 + 1.5 * iqr)).mean(),
                      "p normalidad (D'Agostino, n=5000)": stats.normaltest(muestra(s, 5000)).pvalue})


tabla_uni = num.apply(resumen_num).T
tabla_uni.loc["altitud (por estación)"] = resumen_num(est.altitud)
tabla_uni

fig, axs = plt.subplots(2, 3, figsize=(16, 8))
for ax, v, u in zip(axs.flat, ["vel", "temp", "pres", "dir_constancia"], ["m/s", "°C", "hPa", "0–1"]):
    sns.histplot(muestra(PT[v].dropna(), 200000), bins=80, ax=ax, color="steelblue")
    ax.set(title=v, xlabel=u)
sns.boxplot(data=muestra(PT, 200000), x="depto", y="vel", ax=axs[1, 2], showfliers=False,
            color="lightsteelblue")
axs[1, 2].set(title="Velocidad por departamento", xlabel="", ylabel="m/s")
axs[1, 2].tick_params(axis="x", rotation=45)
# rosa de vientos
axs[1, 1].remove()
ax = fig.add_subplot(2, 3, 5, projection="polar")
d = PT[["dir", "vel"]].dropna()
sect = (np.floor(((d.dir + 11.25) % 360) / 22.5)).astype(int)
cuenta = sect.value_counts(normalize=True).sort_index()
vmed = d.vel.groupby(sect).mean().sort_index()
barras = ax.bar(np.deg2rad(cuenta.index * 22.5), cuenta.values, width=np.deg2rad(22.5), bottom=0,
                color=plt.cm.viridis(vmed.values / vmed.max()), edgecolor="k", lw=.3)
ax.set_theta_zero_location("N"); ax.set_theta_direction(-1)
ax.set_title("Rosa de vientos (color = velocidad media)", pad=15)
fig.tight_layout(); plt.show()

ang = np.arctan2(PT.dir_sin.mean(), PT.dir_cos.mean())
print(f"Dirección media vectorial: {np.rad2deg(ang) % 360:.0f}° | longitud del vector medio: "
      f"{np.hypot(PT.dir_sin.mean(), PT.dir_cos.mean()):.2f}")

cat = pd.DataFrame({"estaciones": est.categoria.value_counts(),
                    "% de horas": (PT.categoria.value_counts(normalize=True) * 100).round(1)})
cat

m = PT[["vel", "temp", "pres", "dir_sin", "dir_cos", "dir_constancia", "altitud", "lat", "lon"]].dropna()
m = muestra(m, 200000)
fig, axs = plt.subplots(1, 2, figsize=(15, 6))
sns.heatmap(m.corr("pearson"), annot=True, fmt=".2f", cmap="RdBu_r", center=0, ax=axs[0], cbar=False)
axs[0].set_title("Pearson")
sns.heatmap(m.corr("spearman"), annot=True, fmt=".2f", cmap="RdBu_r", center=0, ax=axs[1], cbar=False)
axs[1].set_title("Spearman")
fig.tight_layout(); plt.show()

fig, axs = plt.subplots(1, 3, figsize=(17, 4.5))
for ax, v in zip(axs[:2], ["temp", "pres"]):
    hb = ax.hexbin(m[v], m.vel, gridsize=60, bins="log", cmap="viridis", mincnt=1)
    ax.set(xlabel=v, ylabel="vel (m/s)", title=f"vel vs. {v}")
axs[2].scatter(est.altitud, PT.groupby("est").vel.mean().reindex(est.index), c="steelblue", edgecolor="k", lw=.3)
axs[2].set(xlabel="altitud (m)", ylabel="velocidad media (m/s)", title="Velocidad media vs. altitud (por estación)")
fig.tight_layout(); plt.show()

est["vel_media"] = PT.groupby("est").vel.mean().reindex(est.index)
est["temp_media"] = PT.groupby("est").temp.mean().reindex(est.index)
pruebas = []
for cat_v in ["depto", "categoria"]:
    for num_v in ["vel_media", "temp_media", "altitud"]:
        e = est.dropna(subset=[num_v])
        grupos = [g[num_v].to_numpy() for _, g in e.groupby(cat_v) if len(g) >= 2]
        H_, p = stats.kruskal(*grupos)
        pruebas.append({"categórica": cat_v, "numérica": num_v, "H": H_, "p": p,
                        "ε²": H_ / (len(e) - 1), "grupos": len(grupos)})
kw = pd.DataFrame(pruebas)
kw["p (Holm)"] = multipletests(kw.p, method="holm")[1]
kw

tabla = pd.crosstab(est.depto, est.categoria)
chi2, p, gl, _ = stats.chi2_contingency(tabla)
v_cramer = np.sqrt(chi2 / (tabla.values.sum() * (min(tabla.shape) - 1)))
print(f"Departamento × categoría (por estación): χ² = {chi2:.1f}, gl = {gl}, p = {p:.3g}, V de Cramér = {v_cramer:.2f}")
tabla

mu = muestra(TR, 20000)
mi = mutual_info_regression(mu[F_ESP], mu.y, random_state=SEED)
asoc = pd.DataFrame({"Spearman con y": mu[F_ESP].corrwith(mu.y, method="spearman"),
                     "información mutua": mi}, index=F_ESP).sort_values("información mutua", ascending=False)
Xv = StandardScaler().fit_transform(mu[F_ESP])
Xv = np.column_stack([np.ones(len(Xv)), Xv])
asoc["VIF"] = pd.Series([variance_inflation_factor(Xv, i + 1) for i in range(len(F_ESP))], index=F_ESP)
asoc

fig, ax = plt.subplots(figsize=(12, 10))
sns.heatmap(mu[F_ESP + ["y"]].corr(), cmap="RdBu_r", center=0, ax=ax, square=True, cbar_kws={"shrink": .6})
ax.set_title("Correlación de Pearson entre predictoras y objetivo (train)")
fig.tight_layout(); plt.show()

Xs = StandardScaler().fit_transform(muestra(TR[F_ESP], 100000))
pca = PCA(random_state=SEED).fit(Xs)
var_acum = np.cumsum(pca.explained_variance_ratio_)
fig, axs = plt.subplots(1, 2, figsize=(15, 4.5))
axs[0].bar(range(1, len(var_acum) + 1), pca.explained_variance_ratio_, color="steelblue", label="individual")
axs[0].plot(range(1, len(var_acum) + 1), var_acum, "o-", color="darkorange", label="acumulada")
axs[0].axhline(.9, ls="--", c="gray"); axs[0].legend()
axs[0].set(title="PCA de las predictoras (estandarizadas)", xlabel="componente", ylabel="varianza explicada")
cargas = pd.DataFrame(pca.components_[:3].T, index=F_ESP, columns=["PC1", "PC2", "PC3"])
sns.heatmap(cargas, cmap="RdBu_r", center=0, ax=axs[1], annot=True, fmt=".2f", annot_kws={"size": 7})
axs[1].set_title("Cargas de las 3 primeras componentes")
fig.tight_layout(); plt.show()
print(f"Componentes para explicar el 90 % de la varianza: {np.argmax(var_acum >= .9) + 1} de {len(F_ESP)}")

# atípicos multivariados relativos a cada estación: variables estandarizadas dentro de su estación
cols4 = ["vel", "temp", "pres", "dir_constancia"]
z = PT.dropna(subset=cols4).groupby("est")[cols4].transform(lambda s: (s - s.mean()) / s.std()).dropna()
zs = muestra(z, 100000)
inv = np.linalg.inv(np.cov(zs.T))
d2 = np.einsum("ij,jk,ik->i", zs.values, inv, zs.values)
umbral = stats.chi2.ppf(.999, df=len(cols4))
iso = IsolationForest(random_state=SEED, contamination="auto").fit(zs)
atip_iso = iso.predict(zs) == -1
print(f"Mahalanobis (χ² 99.9 %): {100 * (d2 > umbral).mean():.2f} % de atípicos multivariados")
print(f"Isolation Forest: {100 * atip_iso.mean():.2f} % | coincidencia con Mahalanobis entre los atípicos: "
      f"{100 * (atip_iso & (d2 > umbral)).sum() / max(1, (d2 > umbral).sum()):.0f} %")
print("Perfil medio (z) de los atípicos de Mahalanobis:")
print(zs[d2 > umbral].mean().round(2).to_string())

perfil = pd.DataFrame({
    "vel_media": PT.groupby("est").vel.mean(),
    "vel_desv": PT.groupby("est").vel.std(),
    "amplitud_diaria": PT.groupby(["est", PT.hora.dt.hour]).vel.mean().groupby("est").agg(lambda s: s.max() - s.min()),
    "constancia_dir": PT.groupby("est").dir_constancia.mean(),
    "temp_media": PT.groupby("est").temp.mean(),
}).join(est[["altitud", "lat", "lon"]]).dropna()
Zp = StandardScaler().fit_transform(perfil.drop(columns=["lat", "lon"]))
sil = {k: silhouette_score(Zp, KMeans(k, n_init=20, random_state=SEED).fit_predict(Zp)) for k in range(2, 7)}
k_opt = max(sil, key=sil.get)
perfil["grupo"] = KMeans(k_opt, n_init=20, random_state=SEED).fit_predict(Zp)
print("Silueta por k:", {k: round(s, 3) for k, s in sil.items()}, "→ k =", k_opt)
fig, ax = plt.subplots(figsize=(8, 7))
mapa_base(ax)
sns.scatterplot(data=perfil, x="lon", y="lat", hue="grupo", size="vel_media", palette="Set1", ax=ax, sizes=(30, 250))
ax.set_title(f"Grupos de estaciones por régimen de viento (k-means, k = {k_opt})")
plt.show()
perfil.groupby("grupo").mean().round(2).assign(estaciones=perfil.grupo.value_counts().sort_index())

tr_a, va_a = TR[TR.t_y < "2024-01-01"], TR[(TR.t >= "2024-01-01")]
uni = {f: r2_score(va_a.y, LinearRegression().fit(tr_a[[f]], tr_a.y).predict(va_a[[f]])) for f in F_ESP}
tabla_fuga = pd.Series(uni, name="R² univariado (2020–23 → 2024)").sort_values(ascending=False).to_frame()

# demostración: una predictora con fuga, la velocidad de la hora t + 23 (la hora justo anterior al
# objetivo, que aún no se conoce al emitir el pronóstico). Se busca por fecha exacta en la rejilla horaria
clave = pd.MultiIndex.from_arrays([D.est, D.t + pd.Timedelta(H - 1, "h")])
D_f = D.assign(fuga=G.vel.reindex(clave).to_numpy()).dropna(subset=["fuga"])
trf, vaf = D_f[D_f.t_y < "2024-01-01"], D_f[(D_f.t >= "2024-01-01") & (D_f.t_y < TEST_INI)]
r2_fuga = r2_score(vaf.y, LinearRegression().fit(trf[["fuga"]], trf.y).predict(vaf[["fuga"]]))
print(f"Con una variable filtrada del futuro (velocidad de la hora t+23, anterior al objetivo): R² = {r2_fuga:.3f}")

chequeos = pd.Series({
    "Objetivos de train con fecha ≥ 2025-01-01": int((TR.t_y >= TEST_INI).sum()),
    "Filas duplicadas (estación, t)": int(D.duplicated(["est", "t"]).sum()),
    "Estaciones de prueba también presentes en train (%)": round(100 * TE.est.isin(TR.est).mean(), 1),
    "Escalado ajustado dentro del Pipeline": "sí",
}, name="resultado").to_frame()
display(tabla_fuga.head(10)); chequeos

dups = PT.duplicated(["est", "hora"]).sum()
huecos = PT.sort_values(["est", "hora"]).groupby("est").hora.diff().dt.total_seconds().div(3600).dropna()
huecos = huecos[huecos > 1]
h_vel = PT.groupby(PT.hora.dt.hour).vel.median().idxmax()
h_tmp = PT.groupby(PT.hora.dt.hour).temp.median().idxmax()
print(f"Zona horaria: la temperatura mediana es máxima a las {h_tmp} h y el viento a las {h_vel} h. En hora local "
      "ambos máximos se esperan a primera hora de la tarde; si las marcas estuvieran en UTC aparecerían 5 h más "
      "tarde (18–20 h). Suponemos hora local de Colombia (UTC−5, sin horario de verano); la fuente no especifica "
      "la zona horaria.")
print("Frecuencia: horaria regular")
print(f"Marcas de tiempo duplicadas (estación, hora): {dups}")
print(f"Huecos (> 1 h) dentro de cada estación: {len(huecos):,} | mediana {huecos.median():.0f} h | "
      f"p90 {huecos.quantile(.9):.0f} h | máx {huecos.max() / 24:.0f} días")
cob = pd.read_csv(DATOS / "cobertura_velocidad.csv", dtype={"CodigoEstacion": str}).set_index(["CodigoEstacion", "nombre", "depto"])
fig, ax = plt.subplots(figsize=(10, 12))
horas_anio = pd.Series({c: 8784 if int(c) % 4 == 0 else 8760 for c in cob.columns})  # 2020 y 2024 son bisiestos
sns.heatmap(cob.sort_values("2024", ascending=False).div(horas_anio), cmap="viridis", vmin=0, vmax=1, ax=ax,
            yticklabels=[n for _, n, _ in cob.sort_values("2024", ascending=False).index],
            cbar_kws={"label": "fracción del año con dato"})
ax.set(title="Cobertura de velocidad por estación y año", xlabel="año", ylabel="")
ax.tick_params(axis="y", labelsize=6)
fig.tight_layout(); plt.show()

fig, axs = plt.subplots(3, 1, figsize=(15, 10))
diaria = PT.groupby(PT.hora.dt.floor("D")).vel.mean()
semanal = diaria.resample("W").mean()
axs[0].plot(diaria.index, diaria, lw=.5, color="lightsteelblue", label="media diaria de la red")
axs[0].plot(semanal.index, semanal, lw=1.2, color="steelblue", label="media semanal")
axs[0].plot(diaria.rolling(30, center=True).mean(), color="darkorange", lw=2, label="media móvil 30 días")
axs[0].set(title="Velocidad media de la red, 2020–2024 (diaria, semanal y móvil de 30 días)", ylabel="m/s")
axs[0].legend(ncol=3)
mensual = PT.groupby([PT.hora.dt.to_period("M").dt.to_timestamp(), "depto"]).vel.mean().unstack()
mensual.plot(ax=axs[1], lw=1.5)
axs[1].set(title="Media mensual por departamento", ylabel="m/s", xlabel=""); axs[1].legend(ncol=4, fontsize=8)
sub = muestra(PT, 300000)
sns.boxplot(x=sub.hora.dt.dayofweek, y=sub.vel, ax=axs[2], showfliers=False, color="lightsteelblue")
axs[2].set(title="Velocidad por día de la semana (0 = lunes)", xlabel="", ylabel="m/s")
fig.tight_layout(); plt.show()

def pettitt(x):
    x = np.asarray(x, float)
    n = len(x)
    U = np.array([np.sign(x[t + 1:, None] - x[None, :t + 1]).sum() for t in range(n - 1)])
    k = np.argmax(np.abs(U))
    K = np.abs(U[k])
    return k, K, min(1.0, 2 * np.exp(-6 * K ** 2 / (n ** 3 + n ** 2)))


# media mensual por estación (solo meses con al menos 15 días con dato) y anomalía respecto a su climatología
dias_mes = PT.groupby([PT.hora.dt.to_period("M"), "est"]).hora.apply(lambda h: h.dt.day.nunique()).unstack()
men = PT.groupby([PT.hora.dt.to_period("M"), "est"]).vel.mean().unstack().where(dias_mes >= 15)
clim_est = men.groupby(men.index.month).transform("mean")
anom_est = men - clim_est
n_est_mes = anom_est.notna().sum(axis=1)
anom_red = anom_est.mean(axis=1)[n_est_mes >= 15]
presentes_siempre = anom_est.loc[anom_red.index].notna().all().sum()
k_cp, K_cp, p_cp = pettitt(anom_red.to_numpy())
fecha_cp = anom_red.index[k_cp + 1].to_timestamp()
fig, (ax, ax2) = plt.subplots(2, 1, figsize=(15, 5.5), sharex=True, gridspec_kw={"height_ratios": [3, 1]})
x_ = anom_red.index.to_timestamp()
ax.bar(x_, anom_red, width=20, color=np.where(anom_red > 0, "steelblue", "indianred"))
ax.axvspan(pd.Timestamp("2020-08-01"), pd.Timestamp("2023-02-28"), color="blue", alpha=.07, label="La Niña (ONI)")
ax.axvspan(pd.Timestamp("2023-06-01"), pd.Timestamp("2024-05-31"), color="red", alpha=.07, label="El Niño (ONI)")
ax.axvline(fecha_cp, color="k", ls="--", label=f"punto de cambio de Pettitt (p = {p_cp:.3f})")
ax.set(title="Anomalía mensual media de las estaciones (cada una respecto a su propia climatología)", ylabel="m/s")
ax.legend(loc="upper left", fontsize=8)
ax2.bar(x_, n_est_mes.loc[anom_red.index], width=20, color="gray")
ax2.set(ylabel="estaciones", title="Estaciones que aportan a cada mes")
fig.tight_layout(); plt.show()
print(f"Meses analizados: {len(anom_red)} | estaciones por mes: {n_est_mes.loc[anom_red.index].min()} a "
      f"{n_est_mes.loc[anom_red.index].max()} | estaciones presentes en todos esos meses: {presentes_siempre}")
media_antes, media_despues = anom_red[:k_cp + 1].mean(), anom_red[k_cp + 1:].mean()
print(f"Pettitt: cambio más probable en {fecha_cp:%Y-%m} (p = {p_cp:.3f}); "
      f"anomalía media antes {media_antes:+.2f} m/s y después {media_despues:+.2f} m/s")
enso = pd.Series(np.select([(x_ >= "2020-08-01") & (x_ <= "2023-02-28"), (x_ >= "2023-06-01") & (x_ <= "2024-05-31")],
                           ["La Niña", "El Niño"], "neutral"), index=anom_red.index)
print("Anomalía media por fase ENSO (m/s):", anom_red.groupby(enso).mean().round(3).to_dict())

e_rep = PT.groupby("est").vel.count().idxmax()
nombre_rep = est.loc[e_rep, "nombre"]
s_h = G.loc[e_rep].vel
s_h = s_h[s_h.index < TEST_INI]
s_d = s_h.resample("D").mean().where(s_h.resample("D").count() >= 18)


def tramo_continuo(s, limite):
    s = s.interpolate(limit=limite, limit_area="inside")
    grupos = s.isna().cumsum()
    largo = s.notna().groupby(grupos).sum()
    return s[grupos == largo.idxmax()].dropna()


sd = tramo_continuo(s_d, 7)
stl = STL(sd, period=365, robust=True).fit()
fig = stl.plot(); fig.set_size_inches(14, 8)
fig.suptitle(f"STL de la velocidad media diaria — {nombre_rep} ({sd.index.min():%Y-%m} a {sd.index.max():%Y-%m})", y=1.01)
plt.show()
fuerza_est = max(0, 1 - stl.resid.var() / (stl.seasonal + stl.resid).var())
fuerza_ten = max(0, 1 - stl.resid.var() / (stl.trend + stl.resid).var())
print(f"Fuerza de la estacionalidad anual: {fuerza_est:.2f} | fuerza de la tendencia: {fuerza_ten:.2f}")

sh = tramo_continuo(s_h, 3).iloc[-24 * 60:]
stl_h = STL(sh, period=24, robust=True).fit()
print(f"Fuerza del ciclo diario (60 días horarios): "
      f"{max(0, 1 - stl_h.resid.var() / (stl_h.seasonal + stl_h.resid).var()):.2f}")

fig, axs = plt.subplots(1, 2, figsize=(15, 4))
axs[0].plot(sd.rolling(30).mean(), label="media móvil 30 d")
axs[0].plot(sd.rolling(30).std(), label="desviación móvil 30 d")
axs[0].set(title=f"Media y varianza móviles — {nombre_rep}", ylabel="m/s"); axs[0].legend()
res = []
for e_, g in PT.groupby("est"):
    s = g.set_index("hora").vel.resample("D").mean()
    s = s.where(g.set_index("hora").vel.resample("D").count() >= 18)
    s = tramo_continuo(s, 7) if s.notna().sum() > 30 else pd.Series(dtype=float)
    if len(s) < 365:
        continue
    anom_s = s - s.groupby(s.index.month).transform("mean")  # sin el ciclo anual (media por mes)
    res.append({"est": e_, "días": len(s), "p ADF": adfuller(s, autolag="AIC")[1],
                "p KPSS": kpss(s, regression="c", nlags="auto")[1],
                "p ADF (anomalía)": adfuller(anom_s, autolag="AIC")[1],
                "p KPSS (anomalía)": kpss(anom_s, regression="c", nlags="auto")[1]})
est_test = pd.DataFrame(res)
resumen_est = pd.DataFrame({
    "serie original": [100 * (est_test["p ADF"] < .05).mean(), 100 * (est_test["p KPSS"] > .05).mean()],
    "sin ciclo anual": [100 * (est_test["p ADF (anomalía)"] < .05).mean(), 100 * (est_test["p KPSS (anomalía)"] > .05).mean()],
}, index=["% estaciones estacionarias según ADF (p < 0.05)", "% estacionarias según KPSS (p > 0.05)"]).round(1)
print(f"Estaciones con ≥ 365 días continuos: {len(est_test)}")
display(resumen_est)
axs[1].hist(est_test["p KPSS"], bins=20, color="steelblue")
axs[1].axvline(.05, c="r", ls="--"); axs[1].set(title="Valores p de KPSS por estación (serie diaria)", xlabel="p")
fig.tight_layout(); plt.show()

shc = tramo_continuo(s_h, 3)
fig, axs = plt.subplots(1, 3, figsize=(17, 4))
axs[0].stem(acf(shc, nlags=72), markerfmt=" ")
axs[0].set(title=f"ACF horaria — {nombre_rep}", xlabel="rezago (h)")
axs[1].stem(pacf(shc, nlags=72, method="ywm"), markerfmt=" ")
axs[1].set(title="PACF horaria", xlabel="rezago (h)")
# correlación de X(t) con vel(t + k), buscando vel(t + k) por fecha exacta en la rejilla horaria
# (no por posición de fila: con huecos, k filas no son k horas)
lags = list(range(0, 49, 4))
base_cc = muestra(TR, 200000)
cc = {v: [] for v in ["vel_l0", "temp", "pres", "u", "w", "vel_vecinos"]}
for k in lags:
    b_k = base_cc[base_cc.t + pd.Timedelta(k, "h") < TEST_INI]  # nunca consultar datos de 2025
    yk = pd.Series(G.vel.reindex(pd.MultiIndex.from_arrays([b_k.est, b_k.t + pd.Timedelta(k, "h")])).to_numpy())
    for v in cc:
        cc[v].append(pd.Series(b_k[v].to_numpy()).corr(yk))
for v, c in cc.items():
    axs[2].plot(lags, c, "o-", ms=3, label=v)
axs[2].set(title="Correlación de X(t) con vel(t+k)", xlabel="k (h)"); axs[2].legend(fontsize=8)
fig.tight_layout(); plt.show()

deriva = []
for v in ["vel", "temp", "pres", "dir_constancia"]:
    base = muestra(PT[PT.anio == 2020][v].dropna(), 20000)
    for a in range(2021, 2025):
        otro = muestra(PT[PT.anio == a][v].dropna(), 20000)
        deriva.append({"variable": v, "año": a, "KS vs. 2020": stats.ks_2samp(base, otro).statistic})
deriva = pd.DataFrame(deriva).pivot(index="variable", columns="año", values="KS vs. 2020")
concepto = TR.groupby(TR.t.dt.year)[["vel_l0", "y"]].apply(lambda g: g.corr().iloc[0, 1]).rename("corr(vel_t, vel_t+24)")
display(deriva.round(3)); print(concepto.round(3).to_string())

fig, axs = plt.subplots(1, 2, figsize=(14, 4))
sns.histplot(perfil.vel_media, bins=20, ax=axs[0], color="steelblue")
axs[0].set(title="Velocidad media por estación", xlabel="m/s")
sns.histplot(perfil.amplitud_diaria, bins=20, ax=axs[1], color="darkorange")
axs[1].set(title="Amplitud del ciclo diario por estación", xlabel="m/s")
fig.tight_layout(); plt.show()

print("Sistema de referencia: EPSG:4326 (WGS84, grados), el del catálogo del IDEAM")
print(f"Latitud: {est.lat.min():.3f} a {est.lat.max():.3f} | longitud: {est.lon.min():.3f} a {est.lon.max():.3f}")
print(f"Coordenadas fuera de rango válido: {(~est.lat.between(-90, 90) | ~est.lon.between(-180, 180)).sum()} | "
      f"(0, 0): {((est.lat == 0) & (est.lon == 0)).sum()} | posibles lat/lon invertidas: {(est.lat.abs() > est.lon.abs()).sum()}")
print(f"Estaciones con coordenadas duplicadas: {est.duplicated(['lat', 'lon']).sum()}")
fuera_deptos = (~est.depto.isin(["ATLANTICO", "BOLIVAR", "CESAR", "CORDOBA", "LA GUAJIRA", "MAGDALENA", "SUCRE"])).sum()
print(f"Caja envolvente de las estaciones: lat {est.lat.min():.2f}–{est.lat.max():.2f} °N, "
      f"lon {est.lon.min():.2f}–{est.lon.max():.2f} °. Estaciones con departamento fuera de los 7 seleccionados: {fuera_deptos}.")
sur = est[est.lat < 8]
print(f"Estaciones al sur de 8 °N (sur de Bolívar y Córdoba, dentro de los departamentos): {len(sur)} — "
      f"{', '.join(sur.nombre.str.strip() + ' (' + sur.depto + ')')}")

lat0 = np.deg2rad(est.lat.mean())
xk = 6371 * np.deg2rad(est.lon) * np.cos(lat0)
yk = 6371 * np.deg2rad(est.lat)
area = ConvexHull(np.column_stack([xk, yk])).volume  # en 2D, 'volume' es el área (km²)
nn = np.sort(np.where(np.eye(len(DIST)) == 1, np.inf, DIST), axis=1)[:, 0]
R_ce = nn.mean() / (0.5 / np.sqrt(len(est) / area))
print(f"Distancia al vecino más cercano: mediana {np.median(nn):.1f} km, máx {nn.max():.1f} km")
print(f"Índice de Clark-Evans R = {R_ce:.2f} (R < 1 agrupado, R ≈ 1 aleatorio, R > 1 regular) | área de estudio ≈ {area:,.0f} km²")

# función K de Ripley (estimador sin corrección de borde) y su versión L(r) − r; bajo CSR, L(r) − r ≈ 0
xy = np.column_stack([xk, yk])
dxy = np.sqrt(((xy[:, None] - xy[None]) ** 2).sum(-1))
np.fill_diagonal(dxy, np.inf)
radios = np.arange(5, 155, 5)
n_e = len(xy)
K_r = np.array([area / (n_e * (n_e - 1)) * (dxy < r).sum() for r in radios])
L_r = np.sqrt(K_r / np.pi) - radios
# envolvente de CSR con la MISMA ventana y el mismo estimador que los datos: puntos uniformes dentro del
# polígono convexo de las estaciones (muestreo por rechazo) y su misma área
from matplotlib.path import Path as Trazado
casco = ConvexHull(xy)
ventana = Trazado(xy[casco.vertices])
sims = []
for _ in range(99):
    xs = np.empty((0, 2))
    while len(xs) < n_e:
        cand = np.column_stack([RNG.uniform(xk.min(), xk.max(), 4 * n_e), RNG.uniform(yk.min(), yk.max(), 4 * n_e)])
        xs = np.vstack([xs, cand[ventana.contains_points(cand)]])
    xs = xs[:n_e]
    ds = np.sqrt(((xs[:, None] - xs[None]) ** 2).sum(-1)); np.fill_diagonal(ds, np.inf)
    sims.append(np.sqrt(np.array([area / (n_e * (n_e - 1)) * (ds < r).sum() for r in radios]) / np.pi) - radios)
sims = np.array(sims)
fuera = radios[(L_r < sims.min(0)) | (L_r > sims.max(0))]
print(f"Ripley: radios donde la curva observada sale de la envolvente CSR: "
      f"{', '.join(map(str, fuera)) + ' km' if len(fuera) else 'ninguno'}")

# clusters espaciales con DBSCAN (distancia haversine)
eps_km = 40
db = DBSCAN(eps=eps_km / 6371, min_samples=3, metric="haversine").fit(np.deg2rad(est[["lat", "lon"]].to_numpy()))
est["cluster_dbscan"] = db.labels_
print(f"DBSCAN (eps = {eps_km} km, mínimo 3 estaciones): {len(set(db.labels_) - {-1})} clusters; "
      f"{(db.labels_ == -1).sum()} estaciones aisladas (ruido)")

fig, axs = plt.subplots(1, 2, figsize=(16, 5))
axs[0].fill_between(radios, sims.min(0), sims.max(0), color="gray", alpha=.3, label="envolvente CSR (99 simulaciones)")
axs[0].plot(radios, L_r, "o-", color="steelblue", label="estaciones")
axs[0].axhline(0, c="k", lw=.6)
axs[0].set(title="Función L de Ripley: L(r) − r", xlabel="r (km)", ylabel="L(r) − r (km)"); axs[0].legend()
mapa_base(axs[1])
for c_, g in est.groupby("cluster_dbscan"):
    axs[1].scatter(g.lon, g.lat, s=45, edgecolor="k", lw=.3, label="aislada" if c_ == -1 else f"cluster {c_}",
                   color="lightgray" if c_ == -1 else None)
axs[1].set(title=f"Clusters espaciales de estaciones (DBSCAN haversine, eps = {eps_km} km)", xlabel="longitud", ylabel="latitud")
axs[1].legend(fontsize=7, ncol=2)
fig.tight_layout(); plt.show()

fig, axs = plt.subplots(1, 3, figsize=(20, 6.5))
mapa_base(axs[0], nombres=True)
sc = axs[0].scatter(est.lon, est.lat, c=est.vel_media, s=est.horas / est.horas.max() * 250 + 20, cmap="viridis",
                    edgecolor="k", lw=.4, zorder=2)
plt.colorbar(sc, ax=axs[0], label="velocidad media (m/s)", shrink=.8)
axs[0].set(title="Estaciones sobre el mapa base\n(tamaño = horas con dato; color = velocidad media)",
           xlabel="longitud", ylabel="latitud")
mapa_base(axs[1])
hb = axs[1].hexbin(est.lon, est.lat, gridsize=12, cmap="Oranges", mincnt=1, alpha=.85, zorder=2)
plt.colorbar(hb, ax=axs[1], label="estaciones por celda", shrink=.8)
axs[1].set(title="Densidad de estaciones (hexbin)", xlabel="longitud", ylabel="latitud")
# coropleta: velocidad media por departamento (media de sus estaciones, 2020–2024)
por_dep = est.dropna(subset=["vel_media"]).groupby("depto").vel_media.mean().to_dict()
mapa_base(axs[2], valores=por_dep, cmap="viridis", etiqueta="velocidad media de sus estaciones (m/s)", nombres=True)
axs[2].scatter(est.lon, est.lat, s=6, c="k", zorder=2)
axs[2].set(title="Coropleta por departamento", xlabel="longitud", ylabel="latitud")
fig.tight_layout(); plt.show()
print("Mapa base: límites departamentales de geoBoundaries (fuente OpenStreetMap, licencia ODbL 1.0).")
print("Velocidad media por departamento (m/s):", {k: round(v, 2) for k, v in sorted(por_dep.items(), key=lambda x: -x[1])})

e_ok = est.dropna(subset=["vel_media"])
W = pesos_knn(haversine_matriz(e_ok.lat.to_numpy(), e_ok.lon.to_numpy()))
I, p_I = moran_perm(e_ok.vel_media.to_numpy(), W)
print(f"I de Moran global (velocidad media por estación): I = {I:.3f}, p = {p_I:.3f} (esperado sin estructura: {-1 / (len(e_ok) - 1):.3f})")

z = e_ok.vel_media.to_numpy() - e_ok.vel_media.mean()
m2 = (z ** 2).mean()
lag = W @ z
Ii = z / m2 * lag
p_loc = []
for i in range(len(z)):
    otros = np.delete(z, i)
    sim = np.array([z[i] / m2 * (RNG.choice(otros, K_VECINOS, replace=False).mean()) for _ in range(999)])
    p_loc.append((np.sum(np.abs(sim) >= abs(Ii[i])) + 1) / 1000)
e_ok = e_ok.assign(I_local=Ii, p_local=p_loc,
                   cuadrante=np.select([(z > 0) & (lag > 0), (z < 0) & (lag < 0), (z > 0) & (lag < 0)],
                                       ["alto-alto", "bajo-bajo", "alto-bajo"], "bajo-alto"))
# 70+ pruebas locales simultáneas: se corrige con Benjamini-Hochberg (FDR)
e_ok["p_fdr"] = multipletests(e_ok.p_local, method="fdr_bh")[1]
print(f"Estaciones con p local < 0.05 sin corrección: {(e_ok.p_local < .05).sum()} | con FDR < 0.05: {(e_ok.p_fdr < .05).sum()}")
e_ok["LISA"] = np.where(e_ok.p_fdr < .05, e_ok.cuadrante, "no significativo")
print(e_ok.LISA.value_counts().to_string())

fig, axs = plt.subplots(1, 2, figsize=(16, 6))
colores = {"alto-alto": "red", "bajo-bajo": "blue", "alto-bajo": "pink", "bajo-alto": "lightblue", "no significativo": "lightgray"}
mapa_base(axs[0])
for k_, g in e_ok.groupby("LISA"):
    axs[0].scatter(g.lon, g.lat, c=colores[k_], label=k_, s=70, edgecolor="k", lw=.4)
axs[0].legend(); axs[0].set(title="LISA (FDR < 0.05): agrupamientos locales de velocidad media", xlabel="longitud", ylabel="latitud")

# semivariograma empírico
Dm = haversine_matriz(e_ok.lat.to_numpy(), e_ok.lon.to_numpy())
iu = np.triu_indices(len(e_ok), 1)
gam = 0.5 * (e_ok.vel_media.to_numpy()[:, None] - e_ok.vel_media.to_numpy())[iu] ** 2
dist_p = Dm[iu]
bins = np.arange(0, 400, 25)
vg = pd.DataFrame({"d": dist_p, "g": gam}).assign(b=pd.cut(dist_p, bins)).groupby("b", observed=True).agg(
    d=("d", "mean"), g=("g", "mean"), n=("g", "size"))
axs[1].plot(vg.d, vg.g, "o-")
axs[1].axhline(e_ok.vel_media.var(), ls="--", c="gray", label="varianza total")
axs[1].set(title="Semivariograma empírico de la velocidad media", xlabel="distancia (km)", ylabel="semivarianza (m²/s²)")
axs[1].legend()
fig.tight_layout(); plt.show()

# hotspots de Getis-Ord Gi* (k vecinos incluyendo la propia estación), con corrección FDR
x_g = e_ok.vel_media.to_numpy()
W_g = (W > 0).astype(float)
np.fill_diagonal(W_g, 1)
n_g, xbar, s_g = len(x_g), x_g.mean(), x_g.std()
wsum = W_g.sum(1)
gi = (W_g @ x_g - xbar * wsum) / (s_g * np.sqrt((n_g * (W_g ** 2).sum(1) - wsum ** 2) / (n_g - 1)))
p_gi = multipletests(2 * stats.norm.sf(np.abs(gi)), method="fdr_bh")[1]
e_ok["Gi*"] = np.where(p_gi < .05, np.where(gi > 0, "hotspot", "coldspot"), "no significativo")
print("Getis-Ord Gi* (FDR < 0.05):", e_ok["Gi*"].value_counts().to_dict())
print("Hotspots:", ", ".join(e_ok[e_ok["Gi*"] == "hotspot"].nombre.str.strip()))

# efecto de la escala de agregación (MAUP): la misma relación a distintas escalas espaciales
e_m = est.dropna(subset=["vel_media", "temp_media"])
por_depto = e_m.groupby("depto")[["vel_media", "temp_media", "altitud"]].mean()
mau = pd.DataFrame({
    "estación (n = %d)" % len(e_m): [stats.spearmanr(e_m.vel_media, e_m.temp_media)[0],
                                     stats.spearmanr(e_m.vel_media, e_m.altitud)[0]],
    "departamento (n = %d)" % len(por_depto): [stats.spearmanr(por_depto.vel_media, por_depto.temp_media)[0],
                                               stats.spearmanr(por_depto.vel_media, por_depto.altitud)[0]],
}, index=["Spearman(velocidad, temperatura)", "Spearman(velocidad, altitud)"])
print("Efecto de la escala de agregación (MAUP): correlaciones con estaciones vs. con medias por departamento")
mau

seca = PT.hora.dt.month.isin([12, 1, 2, 3, 4])
temporada = pd.DataFrame({"seca-ventosa (dic–abr)": PT[seca].groupby("est").vel.mean(),
                          "lluviosa (may–nov)": PT[~seca].groupby("est").vel.mean()}).join(est[["lat", "lon"]])
fig, axs = plt.subplots(1, 3, figsize=(19, 5.5))
vmax = temporada.iloc[:, :2].max().max()
for ax, c in zip(axs[:2], temporada.columns[:2]):
    mapa_base(ax)
    sc = ax.scatter(temporada.lon, temporada.lat, c=temporada[c], cmap="viridis", vmin=0, vmax=vmax, s=60, edgecolor="k", lw=.3)
    ax.set(title=f"Velocidad media — temporada {c}", xlabel="longitud", ylabel="latitud")
plt.colorbar(sc, ax=axs[1], label="m/s")
cambio = (temporada.iloc[:, 0] - temporada.iloc[:, 1])
t_ok = temporada.dropna()  # estaciones con datos en ambas temporadas
W_t = pesos_knn(haversine_matriz(t_ok.lat.to_numpy(), t_ok.lon.to_numpy()))
I_s, p_s = moran_perm(t_ok.iloc[:, 0].to_numpy(), W_t)
I_l, p_l = moran_perm(t_ok.iloc[:, 1].to_numpy(), W_t)
print(f"I de Moran: temporada seca {I_s:.3f} (p = {p_s:.3f}) | temporada lluviosa {I_l:.3f} (p = {p_l:.3f})")
print(f"Diferencia seca − lluviosa por estación: mediana {cambio.median():.2f} m/s; "
      f"{100 * (cambio > 0).mean():.0f} % de las estaciones es más ventosa en la temporada seca")

# decaimiento de la correlación temporal entre estaciones con la distancia
dia = PT.groupby(["est", PT.hora.dt.floor("D")]).vel.mean().unstack(0)
anom = dia - dia.groupby(dia.index.month).transform("mean")
C = anom.corr(min_periods=200)
ests = [e for e in C.index if e in est.index]
Dc = haversine_matriz(est.loc[ests, "lat"].to_numpy(), est.loc[ests, "lon"].to_numpy())
iu = np.triu_indices(len(ests), 1)
pares = pd.DataFrame({"d": Dc[iu], "r": C.loc[ests, ests].to_numpy()[iu]}).dropna()
pares["b"] = pd.cut(pares.d, np.arange(0, 450, 25))
dec = pares.groupby("b", observed=True).agg(d=("d", "mean"), r=("r", "median"))
axs[2].scatter(pares.d, pares.r, s=4, alpha=.3, color="gray")
axs[2].plot(dec.d, dec.r, "o-", color="darkorange", label="mediana por bloque de 25 km")
axs[2].set(title="Correlación de anomalías diarias entre pares de estaciones", xlabel="distancia (km)", ylabel="r")
axs[2].legend()
fig.tight_layout(); plt.show()

def preparar_panel(ruta):
    # lee un panel estación-hora y construye el conjunto de modelado (mismas predictoras que en 2.0)
    P_ = pd.read_csv(ruta, dtype={"CodigoEstacion": str}, parse_dates=["hora"]).rename(columns={"CodigoEstacion": "est"})
    co = P_.groupby("est")[["lat", "lon"]].first()
    di = haversine_matriz(co.lat.to_numpy(), co.lon.to_numpy())
    np.fill_diagonal(di, np.inf)
    vec = {e: co.index[np.argsort(di[i])[:K_VECINOS]].tolist() for i, e in enumerate(co.index)}
    an = P_.pivot(index="hora", columns="est", values="vel")
    vv = pd.DataFrame({e: an[v].mean(axis=1) for e, v in vec.items()})
    G_ = (P_.set_index("hora").groupby("est")[["vel", "dir_sin", "dir_cos", "dir_constancia", "temp", "pres", "altitud"]]
          .apply(lambda g: g.asfreq("h")))
    todo = G_.groupby(level="est", group_keys=True).apply(construir, vv=vv).reset_index().rename(columns={"hora": "t"})
    return P_, todo, todo.dropna(subset=F_ESP + ["y"]).reset_index(drop=True)


PC, D_todo_c, D_c = preparar_panel(DATOS / "panel_multivariado_causal.csv")
D_o, TR_o, TE_o, D_todo_o = D, TR, TE, D_todo            # procesamiento original (EDA y comparación en 3.6)
D, D_todo = D_c, D_todo_c                                # desde aquí, todo el modelado usa la limpieza causal
TR, TE = D[D.t_y < TEST_INI], D[D.t >= TEST_INI]
completas_c = PC[["vel", "dir", "temp", "pres"]].notna().all(axis=1)
print(f"Panel causal: {len(PC):,} estación-horas con velocidad ({PC.est.nunique()} estaciones); "
      f"{completas_c.sum():,} con las 4 variables")
print(f"Conjunto de modelado causal: {len(D):,} filas | entrenamiento {len(TR):,} | prueba {len(TE):,} "
      f"({TE.est.nunique()} estaciones) | {len(D) - len(TR) - len(TE)} filas de separación en el cambio de año")

def svr(C=0.1, log=False):
    base = Pipeline([("escalado", StandardScaler()),
                     ("svr", LinearSVR(C=C, epsilon=0.0, loss="squared_epsilon_insensitive", dual=False,
                                       max_iter=20000, random_state=SEED))])
    return TransformedTargetRegressor(regressor=base, func=np.log1p, inverse_func=np.expm1) if log else base


def pliegues(datos):
    for a in range(2021, 2025):
        yield a, datos[datos.t_y < f"{a}-01-01"], datos[(datos.t >= f"{a}-01-01") & (datos.t_y < f"{a + 1}-01-01")]


def pliegues_2024(datos):
    """4 pliegues cronológicos (trimestres de 2024) con ventana de entrenamiento creciente."""
    cortes = pd.to_datetime(["2024-01-01", "2024-04-01", "2024-07-01", "2024-10-01", "2025-01-01"])
    for k in range(4):
        ini, fin = cortes[k], cortes[k + 1]
        yield f"2024-T{k + 1}", datos[datos.t_y < ini], datos[(datos.t >= ini) & (datos.t_y < fin)]


def validar(generador):
    filas = []
    for log in [False, True]:
        for C in [0.01, 0.1, 1.0]:
            for nombre, tr, va in generador(TR):
                rm = np.sqrt(mean_squared_error(va.y, svr(C, log).fit(tr[F_ESP], tr.y).predict(va[F_ESP])))
                filas.append({"log(1+y)": log, "C": C, "pliegue": nombre, "RMSE": rm,
                              "mejora": 1 - rm / np.sqrt(mean_squared_error(va.y, va.vel_l0))})
    r = pd.DataFrame(filas)
    t = r.pivot_table(index=["log(1+y)", "C"], columns="pliegue", values="RMSE")
    t["RMSE medio"] = t.mean(axis=1)
    t["mejora media vs. persistencia"] = r.groupby(["log(1+y)", "C"]).mejora.mean()
    return t


tabla_cv = validar(pliegues_2024)
LOG_OPT, C_OPT = tabla_cv["RMSE medio"].idxmin()
print(f"Configuración elegida (validación cruzada, 4 trimestres de 2024): log(1+y) = {LOG_OPT}, C = {C_OPT}")
display(tabla_cv)
tabla_cv_anual = validar(pliegues)
print("Referencia de estabilidad, pliegues por años completos (no se usan para elegir): "
      f"la configuración de menor RMSE medio sería log(1+y) = {tabla_cv_anual['RMSE medio'].idxmin()[0]}, "
      f"C = {tabla_cv_anual['RMSE medio'].idxmin()[1]}")
tabla_cv_anual

def bootstrap_dias(y, pred, pers, t, B=1000, bloque=7, rng=np.random.default_rng(SEED)):
    """IC 95 % remuestreando bloques de 7 días (respeta la autocorrelación)."""
    df = pd.DataFrame({"d": pd.to_datetime(t).dt.floor("D").to_numpy(), "y": y, "e": y - pred, "ep": y - pers})
    agg = df.groupby("d").agg(n=("y", "size"), sy=("y", "sum"), sy2=("y", lambda v: (v ** 2).sum()),
                              se=("e", lambda v: (v ** 2).sum()), sa=("e", lambda v: v.abs().sum()),
                              sep=("ep", lambda v: (v ** 2).sum())).to_numpy()
    nb = int(np.ceil(len(agg) / bloque))
    out = []
    for _ in range(B):
        ini = rng.integers(0, len(agg) - bloque + 1, nb)
        n, sy, sy2, se, sa, sep = agg[(ini[:, None] + np.arange(bloque)).ravel()].sum(axis=0)
        out.append([1 - se / (sy2 - sy ** 2 / n), np.sqrt(se / n), sa / n, 1 - np.sqrt(se / sep)])
    return pd.DataFrame(out, columns=["R2", "RMSE", "MAE", "mejora"]).quantile([.025, .975]).T


def metricas(y, pred, pers, grupos):
    rm = np.sqrt(mean_squared_error(y, pred))
    m = y >= 0.5  # MAPE solo con viento ≥ 0.5 m/s: con calma el error porcentual no está definido
    r2e = pd.DataFrame({"g": grupos, "y": y, "p": pred}).groupby("g").apply(lambda d: r2_score(d.y, d.p)).median()
    return {"R2": r2_score(y, pred), "R2 por estación (mediana)": r2e, "RMSE": rm,
            "MAE": mean_absolute_error(y, pred), "MAPE (y ≥ 0.5 m/s, %)": 100 * mean_absolute_percentage_error(y[m], pred[m]),
            "Mejora vs. persistencia": 1 - rm / np.sqrt(mean_squared_error(y, pers))}


y_te, pers_te = TE.y.to_numpy(), TE.vel_l0.to_numpy()
clim = TR.groupby(["est", "mes_y", "hora_y"]).y.mean()
clim_te = clim.reindex(pd.MultiIndex.from_arrays([TE.est, TE.mes_y, TE.hora_y])).to_numpy()
clim_te = np.where(np.isnan(clim_te), TR.y.mean(), clim_te)
modelos = {
    "Dummy (media de train)": DummyRegressor().fit(TR[F_ESP], TR.y).predict(TE[F_ESP]),
    "Persistencia 24 h": pers_te,
    "Climatología estación × mes × hora": clim_te,
}
final = {}
for nombre, F in [("SVR solo velocidad", F_UNI), ("SVR + meteorología", F_MET), ("SVR + meteorología + espacial", F_ESP)]:
    final[nombre] = svr(C_OPT, LOG_OPT).fit(TR[F], TR.y)
    modelos[nombre] = final[nombre].predict(TE[F])
resultados = pd.DataFrame({k: metricas(y_te, v, pers_te, TE.est.to_numpy()) for k, v in modelos.items()}).T
resultados

pred_te = modelos["SVR + meteorología + espacial"]
ic = bootstrap_dias(y_te, pred_te, pers_te, TE.t)
ic.columns = ["IC 2.5 %", "IC 97.5 %"]
ic.insert(0, "estimación", [r2_score(y_te, pred_te), np.sqrt(mean_squared_error(y_te, pred_te)),
                            mean_absolute_error(y_te, pred_te), resultados.loc["SVR + meteorología + espacial", "Mejora vs. persistencia"]])
ic_pers = bootstrap_dias(y_te, pers_te, pers_te, TE.t)
print(f"R² de la persistencia en prueba: IC 95 % [{ic_pers.loc['R2'].iloc[0]:.3f}, {ic_pers.loc['R2'].iloc[1]:.3f}]")
display(ic)

# sensibilidad a la longitud del bloque
sens = pd.DataFrame({f"bloque {b} días": bootstrap_dias(y_te, pred_te, pers_te, TE.t, bloque=b).loc["mejora"]
                     for b in [3, 7, 14]}).T
sens.columns = ["mejora IC 2.5 %", "mejora IC 97.5 %"]
display(sens)


def diferencia_pareada(y, p_a, p_b, t, B=1000, bloque=7, rng=np.random.default_rng(SEED)):
    """IC 95 % de RMSE(a) − RMSE(b) sobre los mismos bloques de días (comparación pareada)."""
    df = pd.DataFrame({"d": pd.to_datetime(t).dt.floor("D").to_numpy(), "ea": (y - p_a) ** 2, "eb": (y - p_b) ** 2})
    agg = df.groupby("d").agg(n=("ea", "size"), sa=("ea", "sum"), sb=("eb", "sum")).to_numpy()
    nb = int(np.ceil(len(agg) / bloque))
    dif = []
    for _ in range(B):
        ini = rng.integers(0, len(agg) - bloque + 1, nb)
        n, sa, sb = agg[(ini[:, None] + np.arange(bloque)).ravel()].sum(axis=0)
        dif.append(np.sqrt(sa / n) - np.sqrt(sb / n))
    return np.sqrt(np.mean((y - p_a) ** 2)) - np.sqrt(np.mean((y - p_b) ** 2)), *np.percentile(dif, [2.5, 97.5])


pares = [("SVR + meteorología + espacial", "SVR solo velocidad"), ("SVR + meteorología + espacial", "SVR + meteorología"),
         ("SVR + meteorología", "SVR solo velocidad")]
tabla_pares = pd.DataFrame([dict(zip(["comparación", "ΔRMSE (m/s)", "IC 2.5 %", "IC 97.5 %"],
                                     [f"{a} − {b}", *diferencia_pareada(y_te, modelos[a], modelos[b], TE.t)]))
                            for a, b in pares]).set_index("comparación")
tabla_pares["¿significativa?"] = np.where((tabla_pares["IC 2.5 %"] > 0) | (tabla_pares["IC 97.5 %"] < 0), "sí", "no")
tabla_pares

r2_est = pd.DataFrame({"est": TE.est, "y": y_te, "p": pred_te, "pers": pers_te}).groupby("est").apply(
    lambda d: pd.Series({"R2 SVR": r2_score(d.y, d.p), "R2 persistencia": r2_score(d.y, d.pers), "n": len(d)}))
fig, axs = plt.subplots(1, 2, figsize=(15, 4.5))
axs[0].scatter(r2_est["R2 persistencia"], r2_est["R2 SVR"], s=r2_est.n / r2_est.n.max() * 150 + 10, alpha=.7)
lim = [min(r2_est.min(numeric_only=True)[["R2 SVR", "R2 persistencia"]].min(), -.5), 1]
axs[0].plot(lim, lim, "k--", lw=.8)
axs[0].set(xlabel="R² persistencia", ylabel="R² SVR", title="R² por estación en prueba (2025)", xlim=lim, ylim=lim)
print(f"El SVR supera a la persistencia en {100 * (r2_est['R2 SVR'] > r2_est['R2 persistencia']).mean():.0f} % de las estaciones")

grupos = []
for k, (_, b) in enumerate(GroupKFold(n_splits=5).split(D, groups=D.est)):
    ests_k = set(D.est.iloc[b])
    tr_k, te_k = TR[~TR.est.isin(ests_k)], TE[TE.est.isin(ests_k)]
    if te_k.empty:
        continue
    p_k = svr(C_OPT, LOG_OPT).fit(tr_k[F_ESP], tr_k.y).predict(te_k[F_ESP])
    grupos.append({"pliegue": k + 1, "estaciones de prueba": te_k.est.nunique(), "filas": len(te_k),
                   **metricas(te_k.y.to_numpy(), p_k, te_k.vel_l0.to_numpy(), te_k.est.to_numpy())})
grupos = pd.DataFrame(grupos).set_index("pliegue")
grupos.loc["media"] = grupos.mean()
axs[1].bar(grupos.index[:-1].astype(str), grupos["Mejora vs. persistencia"].iloc[:-1], color="steelblue")
axs[1].set(title="Estaciones no usadas para ajustar (GroupKFold): mejora vs. persistencia", xlabel="pliegue")
fig.tight_layout(); plt.show()
grupos

BUFFER_KM = 30
est_mod = est.loc[est.index.isin(D.est.unique())]
xy_mod = np.column_stack([6371 * np.deg2rad(est_mod.lon) * np.cos(lat0), 6371 * np.deg2rad(est_mod.lat)])
bloque = pd.Series(KMeans(6, n_init=20, random_state=SEED).fit_predict(xy_mod), index=est_mod.index)
D_mod = haversine_matriz(est_mod.lat.to_numpy(), est_mod.lon.to_numpy())
filas_b = []
for b in sorted(bloque.unique()):
    en_b = bloque.index[bloque == b]
    cerca = est_mod.index[(D_mod[:, bloque.to_numpy() == b] < BUFFER_KM).any(axis=1)]
    tr_b = TR[~TR.est.isin(set(en_b) | set(cerca))]
    te_b = TE[TE.est.isin(en_b)]
    if te_b.est.nunique() < 2:
        continue
    p_b = svr(C_OPT, LOG_OPT).fit(tr_b[F_MET], tr_b.y).predict(te_b[F_MET])
    filas_b.append({"bloque": b, "estaciones de prueba": te_b.est.nunique(),
                    "excluidas por el buffer": len(set(cerca) - set(en_b)), "filas de prueba": len(te_b),
                    **metricas(te_b.y.to_numpy(), p_b, te_b.vel_l0.to_numpy(), te_b.est.to_numpy())})
bloques = pd.DataFrame(filas_b).set_index("bloque")
bloques.loc["media"] = bloques.mean()
fig, ax = plt.subplots(figsize=(8, 7))
mapa_base(ax)
for b, g in est_mod.groupby(bloque):
    ax.scatter(g.lon, g.lat, s=50, edgecolor="k", lw=.3, label=f"bloque {b}")
ax.set(title="Bloques geográficos para la validación espacial", xlabel="longitud", ylabel="latitud"); ax.legend(fontsize=8)
plt.show()
bloques

print(f"Predicciones negativas en prueba: {(pred_te < 0).sum()} de {len(pred_te):,} "
      f"({100 * (pred_te < 0).mean():.2f} %); mínimo {pred_te.min():.2f} m/s")
tramos = pd.cut(y_te, [0, 1, 2, 4, 6, np.inf], right=False, labels=["< 1", "1–2", "2–4", "4–6", "≥ 6"])
por_tramo = pd.DataFrame({"y": y_te, "svr": pred_te, "pers": pers_te, "tramo": tramos}).groupby(
    "tramo", observed=True).apply(lambda d: pd.Series({
        "filas": len(d), "RMSE SVR": np.sqrt(np.mean((d.y - d.svr) ** 2)),
        "RMSE persistencia": np.sqrt(np.mean((d.y - d.pers) ** 2)), "sesgo SVR (ŷ − y)": np.mean(d.svr - d.y)}))
por_tramo

D_vel = D_todo.dropna(subset=F_UNI + ["y"])
TR_v, TE_v = D_vel[D_vel.t_y < TEST_INI], D_vel[D_vel.t >= TEST_INI]
p_v = svr(C_OPT, LOG_OPT).fit(TR_v[F_UNI], TR_v.y).predict(TE_v[F_UNI])
cobertura = pd.DataFrame({
    "filas de prueba": [len(TE), len(TE_v)], "estaciones de prueba": [TE.est.nunique(), TE_v.est.nunique()],
    **{k: [resultados.loc["SVR solo velocidad", k], v] for k, v in
       metricas(TE_v.y.to_numpy(), p_v, TE_v.vel_l0.to_numpy(), TE_v.est.to_numpy()).items()}},
    index=["SVR solo velocidad, casos completos (4 variables)", "SVR solo velocidad, todas las filas con velocidad"])
cobertura

res = pd.DataFrame({"est": TE.est.to_numpy(), "t": TE.t.to_numpy(), "pred": pred_te, "res": y_te - pred_te})
fig, axs = plt.subplots(1, 3, figsize=(17, 4.5))
sns.histplot(res.res, bins=100, ax=axs[0], color="steelblue")
axs[0].set(title="Distribución de residuos", xlabel="y − ŷ (m/s)")
stats.probplot(muestra(res.res, 20000), dist="norm", plot=axs[1])
axs[1].set_title("QQ-plot de residuos")
axs[2].hexbin(res.pred, res.res, gridsize=60, bins="log", cmap="viridis", mincnt=1)
axs[2].axhline(0, c="r", lw=.8)
axs[2].set(title="Residuos vs. predicción", xlabel="ŷ (m/s)", ylabel="residuo")
fig.tight_layout(); plt.show()
mb = muestra(res, 20000)
bp = het_breuschpagan(mb.res, np.column_stack([np.ones(len(mb)), mb.pred]))
print(f"Asimetría: {stats.skew(res.res):.2f} | curtosis: {stats.kurtosis(res.res):.2f} | "
      f"Breusch-Pagan: LM = {bp[0]:.0f}, p = {bp[1]:.2g}")

acf_est = []
for e_, g in res.groupby("est"):
    s = g.set_index("t").res.asfreq("h")
    if s.notna().sum() > 500:
        acf_est.append(acf(s, nlags=48, missing="conservative"))
acf_m = np.mean(acf_est, axis=0)

dias = res.assign(d=res.t.dt.floor("D")).groupby(["d", "est"]).res.mean().reset_index()
I_dia = []
for d_, g in dias.groupby("d"):
    if len(g) < 10:
        continue
    c = est.loc[g.est, ["lat", "lon"]]
    Wd = pesos_knn(haversine_matriz(c.lat.to_numpy(), c.lon.to_numpy()))
    I_dia.append(moran_perm(g.res.to_numpy(), Wd, n=99))
I_dia = np.array(I_dia)
em = res.groupby("est").res.mean()
I_est, p_est = moran_perm(em.to_numpy(), pesos_knn(haversine_matriz(est.loc[em.index, "lat"].to_numpy(),
                                                                     est.loc[em.index, "lon"].to_numpy())))
fig, axs = plt.subplots(1, 2, figsize=(15, 4))
axs[0].bar(range(49), acf_m, color="steelblue")
axs[0].axhspan(-1.96 / np.sqrt(len(res) / len(acf_est)), 1.96 / np.sqrt(len(res) / len(acf_est)), color="gray", alpha=.3)
axs[0].set(title="ACF media de los residuos por estación", xlabel="rezago (h)")
axs[1].hist(I_dia[:, 0], bins=30, color="steelblue")
axs[1].set(title="I de Moran diario de los residuos (2025)", xlabel="I")
fig.tight_layout(); plt.show()
print(f"ACF de residuos: rezago 1 = {acf_m[1]:.2f}, rezago 24 = {acf_m[24]:.2f}")
print(f"I de Moran diario: media {I_dia[:, 0].mean():.3f}; significativo (p < 0.05) en {100 * (I_dia[:, 1] < .05).mean():.0f} % de los días")
print(f"I de Moran del residuo medio por estación: {I_est:.3f} (p = {p_est:.3f})")

base_lc, val_lc = TR[TR.t_y < "2024-01-01"], TR[TR.t >= "2024-01-01"]
fechas = np.sort(base_lc.t.unique())
lc = []
for frac in [.05, .1, .2, .4, .6, .8, 1.0]:
    sub = base_lc[base_lc.t >= fechas[int(len(fechas) * (1 - frac))]]
    mod = svr(C_OPT, LOG_OPT).fit(sub[F_ESP], sub.y)
    lc.append({"filas de entrenamiento": len(sub),
               "RMSE train": np.sqrt(mean_squared_error(sub.y, mod.predict(sub[F_ESP]))),
               "RMSE validación 2024": np.sqrt(mean_squared_error(val_lc.y, mod.predict(val_lc[F_ESP])))})
lc = pd.DataFrame(lc)
fig, ax = plt.subplots(figsize=(8, 4))
ax.plot(lc["filas de entrenamiento"], lc["RMSE train"], "o-", label="entrenamiento")
ax.plot(lc["filas de entrenamiento"], lc["RMSE validación 2024"], "s-", label="validación 2024")
ax.set(xscale="log", xlabel="filas de entrenamiento (escala log)", ylabel="RMSE (m/s)", title="Curva de aprendizaje")
ax.legend(); plt.show()
lc

reg = final["SVR + meteorología + espacial"]
svr_fit = reg.regressor_.named_steps["svr"] if LOG_OPT else reg.named_steps["svr"]
coef = pd.Series(svr_fit.coef_, index=F_ESP).sort_values(key=np.abs, ascending=False)
fig, ax = plt.subplots(figsize=(9, 7))
coef.iloc[::-1].plot.barh(ax=ax, color=np.where(coef.iloc[::-1] > 0, "steelblue", "indianred"))
ax.set(title="Coeficientes estandarizados del SVR final" + (" (escala log(1+y))" if LOG_OPT else ""),
       xlabel="coeficiente")
plt.show()
coef.to_frame("coeficiente").T

pred_o = svr(C_OPT, LOG_OPT).fit(TR_o[F_ESP], TR_o.y).predict(TE_o[F_ESP])
comp = pd.DataFrame({
    "original (depuración retrospectiva)": {"filas de prueba": len(TE_o), "estaciones de prueba": TE_o.est.nunique(),
                                            **metricas(TE_o.y.to_numpy(), pred_o, TE_o.vel_l0.to_numpy(), TE_o.est.to_numpy())},
    "causal (principal)": {"filas de prueba": len(TE), "estaciones de prueba": TE.est.nunique(),
                           **metricas(y_te, pred_te, pers_te, TE.est.to_numpy())},
}).T
display(comp)


def comparar_comunes(TE_a, pred_a, TE_b, pred_b, nombre_a, nombre_b):
    a = TE_a[["est", "t", "y", "vel_l0"]].assign(pred=pred_a)
    b = TE_b[["est", "t", "y", "vel_l0"]].assign(pred=pred_b)
    com = a.merge(b, on=["est", "t"], suffixes=("_a", "_b"))
    igual = np.isclose(com.y_a, com.y_b, atol=1e-6)
    print(f"Filas comunes: {len(com):,} ({100 * len(com) / len(TE_a):.1f} % de {nombre_a} y "
          f"{100 * len(com) / len(TE_b):.1f} % de {nombre_b}); con el mismo objetivo: {100 * igual.mean():.1f} %")
    cm = com[igual]
    dif, lo_, hi_ = diferencia_pareada(cm.y_b.to_numpy(), cm.pred_b.to_numpy(), cm.pred_a.to_numpy(), cm.t)
    print(f"Sobre esas {len(cm):,} filas: RMSE({nombre_b}) − RMSE({nombre_a}) = {dif:+.4f} m/s "
          f"(IC 95 %: {lo_:+.4f} a {hi_:+.4f})")
    return pd.DataFrame({
        nombre_a: metricas(cm.y_a.to_numpy(), cm.pred_a.to_numpy(), cm.vel_l0_a.to_numpy(), cm.est.to_numpy()),
        nombre_b: metricas(cm.y_b.to_numpy(), cm.pred_b.to_numpy(), cm.vel_l0_b.to_numpy(), cm.est.to_numpy())}).T


comparar_comunes(TE_o, pred_o, TE, pred_te, "original", "causal")

_, _, D_alt = preparar_panel(DATOS / "panel_multivariado_causal_temp_alt.csv")
TR_alt, TE_alt = D_alt[D_alt.t_y < TEST_INI], D_alt[D_alt.t >= TEST_INI]
pred_alt = svr(C_OPT, LOG_OPT).fit(TR_alt[F_ESP], TR_alt.y).predict(TE_alt[F_ESP])
display(pd.DataFrame({
    "causal, temperatura: fuera de [3, 45] °C, fijado antes de 2025 (principal)": {
        "filas de prueba": len(TE), "estaciones de prueba": TE.est.nunique(),
        **metricas(y_te, pred_te, pers_te, TE.est.to_numpy())},
    "causal, temperatura: solo > 45 °C, fijado viendo 2025": {
        "filas de prueba": len(TE_alt), "estaciones de prueba": TE_alt.est.nunique(),
        **metricas(TE_alt.y.to_numpy(), pred_alt, TE_alt.vel_l0.to_numpy(), TE_alt.est.to_numpy())},
}).T)
comparar_comunes(TE_alt, pred_alt, TE, pred_te, "criterio fijado viendo 2025", "criterio principal")

import sklearn, scipy, statsmodels, matplotlib
print("Versiones:", {m.__name__: m.__version__ for m in [np, pd, sklearn, scipy, statsmodels, matplotlib, sns]},
      "| semilla:", SEED)