Riesgo de crédito: Deep Learning vs. gradient boosting

"""#6 Riesgo de crédito con Deep Learning vs. gradient boosting (comparación honesta).

Datos: Home Credit Default Risk (Kaggle; requiere aceptar las reglas de la competencia y KAGGLE_API_TOKEN).
307,511 solicitudes con 8 % de incumplimiento y tablas relacionales: buró de crédito, solicitudes previas
y pagos de cuotas. Se construyen agregados por cliente con DuckDB (SQL) y se comparan:
regresión logística (línea base), LightGBM, una red neuronal tabular en PyTorch con embeddings para
variables categóricas, y su ensamble. Métricas de discriminación, calibración, estabilidad entre semillas
y costo (tiempo de entrenamiento e inferencia), para decidir qué conviene en producción.
"""
import os
import subprocess
import sys
import time
from pathlib import Path
import zipfile

import duckdb
import numpy as np
import pandas as pd
import torch
from lightgbm import LGBMClassifier, early_stopping
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import brier_score_loss, roc_auc_score
from sklearn.model_selection import train_test_split
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import QuantileTransformer, StandardScaler
from torch import nn

from common import RAW, SEED, bars, cumulative_gain_curve, export, ks_stat, line, table
from validacion import bootstrap_ci, calibracion_lineal, graficas_validacion, model_card

DIR = RAW / "kaggle_hc"
ARCHIVOS = ["application_train.csv", "bureau.csv", "previous_application.csv", "installments_payments.csv"]
SEMILLAS = [SEED, SEED + 1, SEED + 2]


def descargar():
    DIR.mkdir(parents=True, exist_ok=True)
    for f in ARCHIVOS:
        if not (DIR / f"{f}.zip").exists() and not (DIR / f).exists():
            subprocess.run([str(Path(sys.executable).parent / "kaggle"), "competitions", "download",
                            "home-credit-default-risk", "-f", f, "-p", str(DIR)], check=True)
    for f in ARCHIVOS:
        if not (DIR / f).exists():
            with zipfile.ZipFile(DIR / f"{f}.zip") as z:
                z.extractall(DIR)


def features() -> pd.DataFrame:
    """Una fila por solicitud: variables de la solicitud + agregados de buró, historial y pagos (SQL en DuckDB)."""
    cache = DIR / "features.parquet"
    if cache.exists():
        return pd.read_parquet(cache)
    descargar()
    con = duckdb.connect()
    for f in ARCHIVOS:
        con.execute(f"CREATE VIEW {f.split('.')[0]} AS SELECT * FROM read_csv_auto('{DIR / f}', sample_size=-1)")
    df = con.execute("""
        WITH buro AS (
            SELECT SK_ID_CURR,
                   COUNT(*)                                             AS buro_creditos,
                   SUM(CASE WHEN CREDIT_ACTIVE = 'Active' THEN 1 ELSE 0 END) AS buro_activos,
                   SUM(AMT_CREDIT_SUM_DEBT)                             AS buro_deuda,
                   MAX(CREDIT_DAY_OVERDUE)                              AS buro_max_dias_vencido,
                   AVG(DAYS_CREDIT)                                     AS buro_antig_media,
                   SUM(AMT_CREDIT_SUM_OVERDUE)                          AS buro_monto_vencido
            FROM bureau GROUP BY 1),
        previas AS (
            SELECT SK_ID_CURR,
                   COUNT(*)                                                  AS prev_solicitudes,
                   AVG(CASE WHEN NAME_CONTRACT_STATUS = 'Approved' THEN 1.0 ELSE 0 END) AS prev_tasa_aprobacion,
                   AVG(CASE WHEN NAME_CONTRACT_STATUS = 'Refused' THEN 1.0 ELSE 0 END)  AS prev_tasa_rechazo,
                   AVG(AMT_APPLICATION)                                      AS prev_monto_medio,
                   AVG(AMT_CREDIT / NULLIF(AMT_APPLICATION, 0))              AS prev_otorgado_vs_pedido
            FROM previous_application GROUP BY 1),
        pagos AS (
            SELECT SK_ID_CURR,
                   COUNT(*)                                                         AS pagos_n,
                   AVG(CASE WHEN DAYS_ENTRY_PAYMENT > DAYS_INSTALMENT THEN 1.0 ELSE 0 END) AS pagos_tasa_tardios,
                   AVG(GREATEST(DAYS_ENTRY_PAYMENT - DAYS_INSTALMENT, 0))            AS pagos_dias_atraso_medio,
                   MAX(DAYS_ENTRY_PAYMENT - DAYS_INSTALMENT)                         AS pagos_max_atraso,
                   AVG(CASE WHEN AMT_PAYMENT < AMT_INSTALMENT THEN 1.0 ELSE 0 END)   AS pagos_tasa_incompletos
            FROM installments_payments GROUP BY 1)
        SELECT a.*, buro.* EXCLUDE (SK_ID_CURR), previas.* EXCLUDE (SK_ID_CURR), pagos.* EXCLUDE (SK_ID_CURR)
        FROM application_train a
        LEFT JOIN buro USING (SK_ID_CURR) LEFT JOIN previas USING (SK_ID_CURR) LEFT JOIN pagos USING (SK_ID_CURR)
    """).df()
    # Limpieza de dominio: 365243 días de empleo es un código de "pensionado/sin dato", no 1,000 años
    df["EMPLEO_SIN_DATO"] = (df.DAYS_EMPLOYED == 365243).astype(int)
    df.loc[df.DAYS_EMPLOYED == 365243, "DAYS_EMPLOYED"] = np.nan
    df["credito_ingreso"] = df.AMT_CREDIT / df.AMT_INCOME_TOTAL
    df["cuota_ingreso"] = df.AMT_ANNUITY / df.AMT_INCOME_TOTAL
    df["credito_bien"] = df.AMT_CREDIT / df.AMT_GOODS_PRICE
    df["plazo_implicito"] = df.AMT_CREDIT / df.AMT_ANNUITY
    df["empleo_edad"] = df.DAYS_EMPLOYED / df.DAYS_BIRTH
    ext = df[["EXT_SOURCE_1", "EXT_SOURCE_2", "EXT_SOURCE_3"]]
    df["ext_media"], df["ext_min"], df["ext_nulos"] = ext.mean(axis=1), ext.min(axis=1), ext.isna().sum(axis=1)
    df.to_parquet(cache)
    return df


# --------------------------------------------------------------------------- red neuronal tabular
class RedTabular(nn.Module):
    """Embeddings para cada categórica (dimensión ≈ 1.6·n^0.56, regla de fast.ai) + numéricas → MLP."""

    def __init__(self, cardinalidades: list[int], n_num: int, ocultas=(256, 128), dropout=0.3):
        super().__init__()
        self.emb = nn.ModuleList([nn.Embedding(c + 1, min(16, round(1.6 * c ** 0.56))) for c in cardinalidades])
        d_emb = sum(e.embedding_dim for e in self.emb)
        capas, d = [], d_emb + n_num
        for h in ocultas:
            capas += [nn.Linear(d, h), nn.BatchNorm1d(h), nn.ReLU(), nn.Dropout(dropout)]
            d = h
        self.mlp = nn.Sequential(*capas, nn.Linear(d, 1))
        self.drop_emb = nn.Dropout(0.1)

    def forward(self, x_cat, x_num):
        e = torch.cat([m(x_cat[:, i]) for i, m in enumerate(self.emb)], dim=1)
        return self.mlp(torch.cat([self.drop_emb(e), x_num], dim=1)).squeeze(1)


def entrenar_red(Xc_tr, Xn_tr, y_tr, Xc_va, Xn_va, y_va, cards, semilla, epocas=30, paciencia=4):
    torch.manual_seed(semilla)
    red = RedTabular(cards, Xn_tr.shape[1])
    opt = torch.optim.AdamW(red.parameters(), lr=2e-3, weight_decay=1e-5)
    sched = torch.optim.lr_scheduler.ReduceLROnPlateau(opt, mode="max", factor=0.5, patience=1)
    perdida = nn.BCEWithLogitsLoss()
    tc, tn, ty = torch.tensor(Xc_tr), torch.tensor(Xn_tr), torch.tensor(y_tr, dtype=torch.float32)
    mejor, estado, sin_mejora, historia = -1, None, 0, []
    gen = torch.Generator().manual_seed(semilla)
    for ep in range(epocas):
        red.train()
        orden = torch.randperm(len(ty), generator=gen)
        total = 0.0
        for i in range(0, len(orden), 1024):
            idx = orden[i:i + 1024]
            opt.zero_grad()
            l = perdida(red(tc[idx], tn[idx]), ty[idx])
            l.backward()
            opt.step()
            total += l.item() * len(idx)
        auc = roc_auc_score(y_va, predecir_red(red, Xc_va, Xn_va))
        historia.append([ep + 1, round(total / len(ty), 4), round(auc, 4)])
        sched.step(auc)
        if auc > mejor + 1e-4:
            mejor, estado, sin_mejora = auc, {k: v.clone() for k, v in red.state_dict().items()}, 0
        else:
            sin_mejora += 1
            if sin_mejora >= paciencia:          # early stopping sobre AUC de validación
                break
    red.load_state_dict(estado)
    return red, historia


@torch.no_grad()
def predecir_red(red, Xc, Xn):
    red.eval()
    out = [torch.sigmoid(red(torch.tensor(Xc[i:i + 8192]), torch.tensor(Xn[i:i + 8192]))) for i in range(0, len(Xc), 8192)]
    return torch.cat(out).numpy()


def hay_datos() -> bool:
    if (DIR / "features.parquet").exists():
        return True
    if all((DIR / f).exists() or (DIR / f"{f}.zip").exists() for f in ARCHIVOS):
        return True
    return bool(os.environ.get("KAGGLE_API_TOKEN")) or (Path.home() / ".kaggle").exists()


def run():
    if not hay_datos():
        print("  · p06: sin datos de Kaggle ni KAGGLE_API_TOKEN; se omite")
        return None
    torch.set_num_threads(4)
    df = features()
    y = df.pop("TARGET").values
    df = df.drop(columns="SK_ID_CURR")
    # Atributos sensibles fuera del modelo (no discriminación)
    df = df.drop(columns=["CODE_GENDER"])
    cats = [c for c in df.columns if df[c].dtype == object]
    nums = [c for c in df.columns if c not in cats]

    idx = np.arange(len(df))
    i_tr, i_te = train_test_split(idx, test_size=0.2, stratify=y, random_state=SEED)
    i_tr, i_va = train_test_split(i_tr, test_size=0.2, stratify=y[i_tr], random_state=SEED)
    y_tr, y_va, y_te = y[i_tr], y[i_va], y[i_te]

    # Codificación común: categorías → enteros (0 = desconocida/nula), aprendida solo en entrenamiento
    mapas = {c: {v: k + 1 for k, v in enumerate(df[c].iloc[i_tr].dropna().unique())} for c in cats}
    Xcat = np.stack([df[c].map(mapas[c]).fillna(0).astype(np.int64).values for c in cats], axis=1)
    cards = [len(mapas[c]) for c in cats]
    Xnum_raw = df[nums].astype(float).replace([np.inf, -np.inf], np.nan)

    resultados, scores, tiempos, historia_red = [], {}, [], None

    def registrar(nombre, p_te, t_fit, t_pred, extra=""):
        auc, lo, hi = bootstrap_ci(y_te, p_te, n=300)
        pend, _ = calibracion_lineal(y_te, p_te)
        resultados.append({"Modelo": nombre, "AUC": round(auc, 4), "IC 95 %": f"{lo:.3f}–{hi:.3f}",
                           "Gini": round(2 * auc - 1, 3), "KS": round(ks_stat(y_te, p_te), 3),
                           "Brier": round(brier_score_loss(y_te, p_te), 4), "Pendiente calib.": round(pend, 2),
                           "Entrenamiento (s)": round(t_fit, 1), "Inferencia por 1,000 (ms)": round(t_pred, 2),
                           "Notas": extra})
        scores[nombre] = p_te

    # 1) Línea base: logística sobre numéricas imputadas (mediana) + indicadores de nulo + one-hot de categóricas
    med = Xnum_raw.iloc[i_tr].median()
    Xlr = pd.concat([Xnum_raw.fillna(med), Xnum_raw.isna().astype(int).add_suffix("_nulo"),
                     pd.get_dummies(df[cats], dummy_na=True, dtype=int)], axis=1).values
    t = time.time()
    lr = make_pipeline(StandardScaler(), LogisticRegression(max_iter=2000, C=0.1)).fit(Xlr[i_tr], y_tr)
    t_fit = time.time() - t
    t = time.time()
    p = lr.predict_proba(Xlr[i_te])[:, 1]
    registrar("Regresión logística", p, t_fit, (time.time() - t) / len(i_te) * 1000 * 1000)

    # 2) LightGBM: varias semillas para medir estabilidad; early stopping en validación
    Xgb = pd.concat([Xnum_raw, df[cats].astype("category")], axis=1)
    p_gbm, auc_gbm_sem, t_gbm = [], [], []
    for s in SEMILLAS:
        t = time.time()
        gbm = LGBMClassifier(n_estimators=3000, learning_rate=0.03, num_leaves=48, min_child_samples=80,
                             subsample=0.8, subsample_freq=1, colsample_bytree=0.6, reg_lambda=2,
                             random_state=s, verbose=-1)
        gbm.fit(Xgb.iloc[i_tr], y_tr, eval_set=[(Xgb.iloc[i_va], y_va)], eval_metric="auc",
                callbacks=[early_stopping(100, verbose=False)])
        t_gbm.append(time.time() - t)
        t = time.time()
        pp = gbm.predict_proba(Xgb.iloc[i_te])[:, 1]
        t_pred_gbm = (time.time() - t) / len(i_te) * 1000 * 1000
        p_gbm.append(pp)
        auc_gbm_sem.append(roc_auc_score(y_te, pp))
    registrar("LightGBM", p_gbm[0], np.mean(t_gbm), t_pred_gbm, f"{gbm.best_iteration_} árboles")
    imp = pd.Series(gbm.booster_.feature_importance("gain"), index=Xgb.columns).sort_values(ascending=False)

    # 3) Red neuronal tabular (PyTorch): numéricas con transformación cuantil→normal + indicadores de nulo
    qt = QuantileTransformer(output_distribution="normal", n_quantiles=1000, subsample=200_000, random_state=SEED)
    Xn = np.nan_to_num(qt.fit(Xnum_raw.iloc[i_tr]).transform(Xnum_raw), nan=0.0)
    Xn = np.hstack([Xn, Xnum_raw.isna().values.astype(float)]).astype(np.float32)
    p_red, auc_red_sem, t_red = [], [], []
    for s in SEMILLAS:
        t = time.time()
        red, hist = entrenar_red(Xcat[i_tr], Xn[i_tr], y_tr, Xcat[i_va], Xn[i_va], y_va, cards, s)
        t_red.append(time.time() - t)
        t = time.time()
        pp = predecir_red(red, Xcat[i_te], Xn[i_te])
        t_pred_red = (time.time() - t) / len(i_te) * 1000 * 1000
        p_red.append(pp)
        auc_red_sem.append(roc_auc_score(y_te, pp))
        historia_red = historia_red or hist
    n_param = sum(p.numel() for p in red.parameters())
    registrar("Red neuronal (PyTorch, embeddings)", p_red[0], np.mean(t_red), t_pred_red, f"{n_param:,} parámetros")

    # 4) Ensamble: promedio en escala logit (no requiere ajustar nada sobre el set de prueba)
    logit = lambda p: np.log(np.clip(p, 1e-6, 1 - 1e-6) / (1 - np.clip(p, 1e-6, 1 - 1e-6)))
    ens = 1 / (1 + np.exp(-(logit(p_gbm[0]) + logit(p_red[0])) / 2))
    registrar("Ensamble LightGBM + red", ens, np.mean(t_gbm) + np.mean(t_red), t_pred_gbm + t_pred_red,
              "promedio de logits")

    res = pd.DataFrame(resultados)
    # Diferencias pareadas por bootstrap (mismas solicitudes remuestreadas)
    rng = np.random.default_rng(SEED)
    boots = [rng.integers(0, len(y_te), len(y_te)) for _ in range(300)]

    def pareada(a, b):
        d = np.array([roc_auc_score(y_te[i], scores[a][i]) - roc_auc_score(y_te[i], scores[b][i]) for i in boots])
        return {"Comparación": f"{a} − {b}", "Δ AUC": round(roc_auc_score(y_te, scores[a]) - roc_auc_score(y_te, scores[b]), 4),
                "IC 95 %": f"{np.percentile(d, 2.5):+.4f} a {np.percentile(d, 97.5):+.4f}",
                "P(Δ ≤ 0)": round(float((d <= 0).mean()), 3)}
    pareadas = pd.DataFrame([pareada("LightGBM", "Red neuronal (PyTorch, embeddings)"),
                             pareada("LightGBM", "Regresión logística"),
                             pareada("Ensamble LightGBM + red", "LightGBM")])
    estabilidad = pd.DataFrame({"Semilla": SEMILLAS, "AUC LightGBM": np.round(auc_gbm_sem, 4),
                                "AUC red neuronal": np.round(auc_red_sem, 4)})
    r = res.set_index("Modelo")
    g, n_ = r.loc["LightGBM"], r.loc["Red neuronal (PyTorch, embeddings)"]
    d_gr = pareadas.iloc[0]
    d_en = pareadas.iloc[2]
    gana_ens = d_en["P(Δ ≤ 0)"] < 0.025
    pierde_ens = d_en["P(Δ ≤ 0)"] > 0.975

    result = {
        "title": "Riesgo de crédito: Deep Learning vs. gradient boosting",
        "summary": [
            f"{len(df):,} solicitudes reales ({y.mean():.1%} de incumplimiento) y 4 tablas relacionales unidas con SQL "
            f"(DuckDB): buró, solicitudes previas y {len(nums) + len(cats)} variables por cliente, incluyendo historial "
            "de pagos tardíos e incompletos.",
            f"LightGBM gana: AUC {g['AUC']:.3f} vs. red neuronal {n_['AUC']:.3f} (Δ pareado {d_gr['Δ AUC']:+.3f}, "
            f"IC {d_gr['IC 95 %']}) y logística {r.loc['Regresión logística', 'AUC']:.3f}. Entre semillas, la red varía "
            f"±{np.std(auc_red_sem):.4f} y LightGBM ±{np.std(auc_gbm_sem):.4f}.",
            f"Costo en CPU (4 núcleos): entrenar toma {g['Entrenamiento (s)']:.0f} s con LightGBM y {n_['Entrenamiento (s)']:.0f} s "
            f"con la red; inferir 1,000 solicitudes, {g['Inferencia por 1,000 (ms)']:.1f} ms vs. "
            f"{n_['Inferencia por 1,000 (ms)']:.1f} ms. El costo no decide: la red pide escalado cuantil, embeddings, "
            "early stopping y más ajuste para quedar atrás por poco.",
            (f"El ensamble sí agrega: {d_en['Δ AUC']:+.4f} de AUC sobre LightGBM (IC {d_en['IC 95 %']}); vale la pena solo "
             "si ese punto extra paga mantener dos modelos." if gana_ens else
             f"El ensamble {'empeora ligeramente' if pierde_ens else 'no mejora de forma concluyente'} a LightGBM "
             f"({d_en['Δ AUC']:+.4f}, IC {d_en['IC 95 %']}): la red no aporta información que el boosting no tenga."),
            "Decisión para producción: LightGBM. Mejor discriminación con una diferencia que excluye el cero, "
            "explicable con SHAP y más fácil de monitorear; la red tabular se justifica cuando hay texto, secuencias o "
            "imágenes, no en tablas.",
        ],
        "metrics": [
            {"label": "AUC LightGBM", "value": f"{g['AUC']:.3f}"},
            {"label": "AUC red neuronal", "value": f"{n_['AUC']:.3f}"},
            {"label": "Δ AUC (boosting − red)", "value": f"{d_gr['Δ AUC']:+.3f}"},
            {"label": "Solicitudes", "value": f"{len(df):,}"},
        ],
        "charts": [
            bars("auc", "AUC en prueba por modelo", res.Modelo, res.AUC.values, fmt="dec"),
            *graficas_validacion(y_te, {k: scores[k] for k in ["LightGBM", "Red neuronal (PyTorch, embeddings)",
                                                                "Regresión logística"]}),
            line("curva_red", "Red neuronal: pérdida de entrenamiento y AUC de validación por época",
                 {"AUC validación": [[e, a] for e, _, a in historia_red]}, x_label="Época", y_label="AUC", fmt="dec",
                 subtitle=f"Early stopping: se conserva la mejor época ({max(historia_red, key=lambda h: h[2])[0]})"),
            line("ganancia", "Ganancia acumulada: % de incumplimientos capturados", {
                k: cumulative_gain_curve(y_te, scores[k]) for k in ["LightGBM", "Red neuronal (PyTorch, embeddings)"]},
                 x_label="% de solicitudes revisadas (mayor riesgo primero)", y_label="% de incumplimientos", fmt="pct"),
            bars("tiempo", "Tiempo de entrenamiento (s, CPU)", res.Modelo, res["Entrenamiento (s)"].values, fmt="num"),
            bars("imp", "LightGBM: importancia por ganancia (top 12)", imp.index[:12],
                 (imp.values[:12] / imp.sum() * 100).round(2), fmt="pct"),
        ],
        "tables": [table("Comparación en el set de prueba (61,503 solicitudes)", res),
                   table("Diferencias pareadas de AUC (bootstrap)", pareadas),
                   table("Estabilidad entre semillas", estabilidad)],
        "card": model_card(
            "Agregados relacionales en SQL (DuckDB) → logística, LightGBM con early stopping y red tabular en PyTorch "
            "(embeddings categóricos + MLP con BatchNorm y dropout); ensamble por promedio de logits.",
            ["Partición estratificada 64/16/20: entrenamiento, validación (early stopping) y prueba intocable.",
             "AUC con IC bootstrap, diferencias pareadas, KS, Brier y pendiente de calibración; 3 semillas por modelo."],
            ["Sin dimensión temporal en los datos públicos: la partición es aleatoria, no fuera de tiempo.",
             "Género excluido del modelo (no discriminación)."],
            ["Búsqueda de hiperparámetros acotada para ambos modelos: mismo presupuesto, comparación justa.",
             "Cartera de un prestamista de consumo internacional, no mexicana: la conclusión es metodológica."]),
    }
    muestra = pd.DataFrame({"ext_media": df.ext_media, "credito_ingreso": df.credito_ingreso,
                            "cuota_ingreso": df.cuota_ingreso, "edad": -df.DAYS_BIRTH / 365.25,
                            "pagos_tasa_tardios": df.pagos_tasa_tardios, "buro_activos": df.buro_activos,
                            "incumple": y}).iloc[i_te].sample(4000, random_state=SEED).round(4)
    export("p06-riesgo-deep-learning", result, sample=muestra)
    return result


if __name__ == "__main__":
    run()