Tasa de Banxico y demanda de crédito

"""#11 Impacto de la tasa de Banxico en la demanda de crédito automotriz y de consumo.

Datos: API SIE de Banco de México (requiere `BANXICO_TOKEN`; con caché en proyectos/data/raw/):
- SF61745  Tasa objetivo (diaria → promedio mensual)
- SP1      INPC general (para deflactar) · SP30578 inflación anual
- SF266105 Cartera vigente automotriz de la banca comercial (saldo nominal, millones de pesos)
- SF235712 Cartera vigente de consumo de la banca comercial
Métodos: empalme de saltos de nivel, proyecciones locales (Jordà) con errores HAC, VAR con impulso-respuesta y
causalidad de Granger, validación fuera de muestra (¿la tasa mejora el pronóstico?) y pronóstico a 12 meses.
"""
import os
import warnings

import httpx
import numpy as np
import pandas as pd
import statsmodels.api as sm
from scipy import stats
from statsmodels.tsa.api import VAR
from statsmodels.tsa.stattools import adfuller

from common import RAW, bars, export, line, table
from validacion import model_card

API = "https://www.banxico.org.mx/SieAPIRest/service/v1/series/{}/datos/{}/{}"
SERIES = {"SF61745": "tasa", "SP1": "inpc", "SP30578": "inflacion", "SF266105": "auto", "SF235712": "consumo"}
DESDE, HASTA = "2011-02-01", "2026-12-31"
H = 24                                     # horizonte de respuesta (meses)
CREDITOS = {"auto": "Crédito automotriz", "consumo": "Crédito al consumo"}


def descargar() -> pd.DataFrame | None:
    cache = RAW / "banxico_sie.csv"
    if cache.exists():
        return pd.read_csv(cache, parse_dates=["fecha"], index_col="fecha")
    token = os.environ.get("BANXICO_TOKEN")
    if not token:
        return None
    r = httpx.get(API.format(",".join(SERIES), DESDE, HASTA), headers={"Bmx-Token": token}, timeout=120)
    r.raise_for_status()
    partes = {}
    for s in r.json()["bmx"]["series"]:
        d = pd.DataFrame(s["datos"])
        d["fecha"] = pd.to_datetime(d["fecha"], format="%d/%m/%Y")
        d["dato"] = pd.to_numeric(d["dato"].str.replace(",", ""), errors="coerce")
        # La tasa es diaria: promedio del mes; las demás ya son mensuales
        partes[SERIES[s["idSerie"]]] = d.set_index("fecha")["dato"].resample("MS").mean()
    df = pd.DataFrame(partes).dropna()
    RAW.mkdir(parents=True, exist_ok=True)
    df.to_csv(cache, index_label="fecha")
    return df


def empalmar(s: pd.Series, umbral: float = 8.0) -> tuple[pd.Series, list[dict]]:
    """Corrige saltos de nivel (reclasificaciones, fusiones de entidades) que no son demanda: si un cambio
    mensual en log está a más de `umbral` desviaciones robustas (MAD) de la mediana, se reemplaza por la
    mediana y se reescala todo el tramo anterior. Devuelve la serie empalmada y los saltos detectados."""
    dl = np.log(s).diff()
    med, mad = dl.median(), 1.4826 * (dl - dl.median()).abs().median()
    z = (dl - med) / mad
    saltos = []
    ajustada = np.log(s).copy()
    for fecha in z[z.abs() > umbral].index:
        exceso = dl[fecha] - med
        ajustada[ajustada.index < fecha] += exceso
        saltos.append({"Mes": f"{fecha:%Y-%m}", "Cambio observado": f"{np.expm1(dl[fecha]):+.1%}",
                       "Desviaciones robustas": round(float(z[fecha]), 1)})
    return np.exp(ajustada), saltos


def proyeccion_local(dlog: pd.Series, dtasa: pd.Series, controles: pd.DataFrame, h_max: int = H) -> pd.DataFrame:
    """Jordà (2005): para cada h, log(y_{t+h}) − log(y_{t−1}) = a + b_h·Δtasa_t + controles_t + e.
    b_h es el efecto acumulado (en %) de subir 1 pp la tasa sobre el saldo h meses después. HAC (Newey-West)."""
    ly = np.log(dlog)
    filas = []
    for h in range(h_max + 1):
        y = (ly.shift(-h) - ly.shift(1)) * 100
        X = pd.concat([dtasa.rename("dtasa"), controles], axis=1)
        d = pd.concat([y.rename("y"), X], axis=1).dropna()
        m = sm.OLS(d.y, sm.add_constant(d.drop(columns="y"))).fit(cov_type="HAC", cov_kwds={"maxlags": h + 1})
        b, se = m.params["dtasa"], m.bse["dtasa"]
        filas.append({"h": h, "b": b, "lo": b - 1.96 * se, "hi": b + 1.96 * se, "p": m.pvalues["dtasa"]})
    return pd.DataFrame(filas)


def evaluacion_fuera_de_muestra(d: pd.DataFrame, credito: str, h: int = 12, inicio: int = 72) -> dict:
    """Pronóstico directo a 12 meses del crecimiento real, con origen móvil: AR (solo historia del crédito)
    vs. AR + tasa e inflación. Diebold-Mariano sobre errores cuadráticos (HAC por traslape de horizontes)."""
    ly = np.log(d[credito + "_real"])
    y = (ly.shift(-h) - ly) * 100
    base = pd.DataFrame({"g12": (ly - ly.shift(12)) * 100, "g1": (ly - ly.shift(1)) * 100})
    ampl = base.assign(tasa=d.tasa, dtasa12=d.tasa - d.tasa.shift(12), inflacion=d.inflacion)
    datos = pd.concat([y.rename("y"), ampl], axis=1).dropna()
    errores = {"AR (solo crédito)": [], "AR + tasa e inflación": []}
    fechas = []
    for t in range(inicio, len(datos)):
        tr = datos.iloc[: t - h + 1]            # solo observaciones cuyo objetivo ya se conocía en t
        te = datos.iloc[[t]]
        for nombre, cols in [("AR (solo crédito)", ["g12", "g1"]), ("AR + tasa e inflación", list(ampl.columns))]:
            m = sm.OLS(tr.y, sm.add_constant(tr[cols])).fit()
            errores[nombre].append(float(te.y.iloc[0] - m.predict(sm.add_constant(te[cols], has_constant="add")).iloc[0]))
        fechas.append(datos.index[t])
    e1, e2 = np.array(errores["AR (solo crédito)"]), np.array(errores["AR + tasa e inflación"])
    dif = e1 ** 2 - e2 ** 2
    nw = sm.OLS(dif, np.ones(len(dif))).fit(cov_type="HAC", cov_kwds={"maxlags": h})
    return {"rmse_ar": float(np.sqrt((e1 ** 2).mean())), "rmse_tasa": float(np.sqrt((e2 ** 2).mean())),
            "dm_t": float(nw.tvalues[0]), "dm_p": float(nw.pvalues[0]), "n": len(dif), "desde": fechas[0]}


def run():
    df = descargar()
    if df is None:
        print("  · p11: sin BANXICO_TOKEN ni caché; se omite")
        return None
    d = df.copy()
    rupturas = []
    for c in CREDITOS:
        d[c], saltos = empalmar(d[c])
        rupturas += [{"Cartera": CREDITOS[c], **x} for x in saltos]
        d[c + "_real"] = d[c] / d.inpc * d.inpc.iloc[-1]          # pesos constantes del último mes
    d = d.dropna()
    meses = f"{d.index[0]:%Y-%m} a {d.index[-1]:%Y-%m}"
    g12 = {c: (np.log(d[c + "_real"]).diff(12) * 100) for c in CREDITOS}

    # Estacionariedad (ADF): niveles vs. diferencias, para justificar la especificación
    adf = []
    for nombre, s in [("log crédito automotriz real", np.log(d.auto_real)),
                      ("Δ log crédito automotriz real", np.log(d.auto_real).diff()),
                      ("Tasa objetivo", d.tasa), ("Δ tasa objetivo", d.tasa.diff()),
                      ("Inflación anual", d.inflacion), ("Δ inflación anual", d.inflacion.diff())]:
        p = adfuller(s.dropna(), autolag="AIC", result_object=False)[1]
        adf.append({"Serie": nombre, "p-valor ADF": round(p, 3), "¿Estacionaria (5 %)?": "Sí" if p < 0.05 else "No"})
    adf = pd.DataFrame(adf)

    # 1) Proyecciones locales: respuesta acumulada del saldo real a +1 pp de tasa
    lp = {}
    for c in CREDITOS:
        ly = np.log(d[c + "_real"])
        controles = pd.concat({f"g_l{k}": ly.diff().shift(k) * 100 for k in range(1, 4)} |
                              {f"dtasa_l{k}": d.tasa.diff().shift(k) for k in range(1, 4)} |
                              {"dinfl": d.inflacion.diff(), "dinfl_l1": d.inflacion.diff().shift(1)}, axis=1)
        lp[c] = proyeccion_local(d[c + "_real"], d.tasa.diff(), controles)

    # 2) VAR en diferencias (estacionarias): [Δ inflación, Δ tasa, Δ log crédito automotriz real]
    v = pd.DataFrame({"dinfl": d.inflacion.diff(), "dtasa": d.tasa.diff(),
                      "dauto": np.log(d.auto_real).diff() * 100}).dropna()
    with warnings.catch_warnings():
        warnings.simplefilter("ignore")
        var = VAR(v).fit(maxlags=12, ic="aic")
    p_var = var.k_ar
    irf = var.irf(H)
    # Choque ortogonal (Cholesky, orden: inflación → tasa → crédito) reescalado a +1 pp de tasa
    escala = 1 / irf.orth_irfs[0, 1, 1]
    acum = irf.orth_cum_effects[:, 2, 1] * escala
    lo_hi = var.irf_errband_mc(orth=True, repl=500, steps=H, signif=0.05, seed=42, cum=True)
    banda_var = [[h, round(lo_hi[0][h, 2, 1] * escala, 3), round(lo_hi[1][h, 2, 1] * escala, 3)] for h in range(H + 1)]
    granger = []
    for causa, efecto in [("dtasa", "dauto"), ("dinfl", "dauto"), ("dauto", "dtasa")]:
        t = var.test_causality(efecto, [causa], kind="f")
        granger.append({"Hipótesis nula": f"{causa} NO causa (Granger) a {efecto}", "F": round(t.test_statistic, 2),
                        "p-valor": round(t.pvalue, 4), "Conclusión": "Se rechaza" if t.pvalue < 0.05 else "No se rechaza"})
    granger = pd.DataFrame(granger)
    estable = bool(var.is_stable())
    resid_p = var.test_whiteness(nlags=max(p_var + 6, 12)).pvalue

    # 3) ¿La tasa mejora el pronóstico a 12 meses? (origen móvil, Diebold-Mariano)
    oos = {c: evaluacion_fuera_de_muestra(d, c) for c in CREDITOS}

    # 4) Pronóstico a 12 meses con el modelo que ganó fuera de muestra. No se publican escenarios de tasa:
    #    con un coeficiente endógeno (la tasa reacciona al ciclo) un "¿y si sube 100 pb?" sería engañoso.
    ult = d.index[-1]
    ly = np.log(d.auto_real)
    X = pd.DataFrame({"g12": ly.diff(12) * 100, "g1": ly.diff() * 100})
    y12 = (ly.shift(-12) - ly) * 100
    datos = pd.concat([y12.rename("y"), X], axis=1).dropna()
    m_ar = sm.OLS(datos.y, sm.add_constant(datos[["g12", "g1"]])).fit()
    g_hat = float(m_ar.predict(sm.add_constant(X.iloc[[-1]], has_constant="add")).iloc[0])
    q10, q90 = np.percentile(m_ar.resid, [10, 90])
    fut = pd.date_range(ult + pd.offsets.MonthBegin(), periods=12, freq="MS")
    saldo0 = d.auto_real.iloc[-1] / 1000                      # miles de millones de pesos constantes
    camino = lambda g: [round(saldo0 * np.exp(g / 100 * (i + 1) / 12), 2) for i in range(12)]
    hist = [[f"{f:%Y-%m}", round(v / 1000, 2)] for f, v in d.auto_real.iloc[-36:].items()]
    pron = [hist[-1]] + [[f"{f:%Y-%m}", v] for f, v in zip(fut, camino(g_hat))]
    banda_pron = [[hist[-1][0], hist[-1][1], hist[-1][1]]] + [
        [f"{f:%Y-%m}", lo, hi] for f, lo, hi in zip(fut, camino(g_hat + q10), camino(g_hat + q90))]
    esc = {"Histórico (empalmado)": hist, "Pronóstico AR a 12 meses": pron}

    # Gráficas y tablas
    serie_t = lambda s: [[f"{f:%Y-%m}", round(float(v), 2)] for f, v in s.dropna().items()]
    lp_series = {CREDITOS[c]: [[int(r.h), round(r.b, 3)] for r in lp[c].itertuples()] for c in CREDITOS}
    lp_bandas = {CREDITOS[c]: [[int(r.h), round(r.lo, 3), round(r.hi, 3)] for r in lp[c].itertuples()] for c in CREDITOS}
    b12 = {c: lp[c].set_index("h").loc[12] for c in CREDITOS}
    b24 = {c: lp[c].set_index("h").loc[24] for c in CREDITOS}
    tabla_lp = pd.DataFrame([{"Cartera": CREDITOS[c], "Horizonte (meses)": h,
                              "Efecto de +1 pp (%)": round(lp[c].set_index("h").loc[h, "b"], 2),
                              "IC 95 %": f"{lp[c].set_index('h').loc[h, 'lo']:.2f} a {lp[c].set_index('h').loc[h, 'hi']:.2f}",
                              "p-valor": round(lp[c].set_index("h").loc[h, "p"], 3)}
                             for c in CREDITOS for h in (6, 12, 18, 24)])
    tabla_oos = pd.DataFrame([{"Cartera": CREDITOS[c], "RMSE AR (pp)": round(o["rmse_ar"], 2),
                               "RMSE AR + tasa (pp)": round(o["rmse_tasa"], 2),
                               "Mejora": f"{1 - o['rmse_tasa'] / o['rmse_ar']:+.0%}",
                               "Diebold-Mariano p": round(o["dm_p"], 3), "Pronósticos": o["n"]}
                              for c, o in oos.items()])
    corr = d[["tasa", "inflacion"]].assign(**{CREDITOS[c]: g12[c] for c in CREDITOS}).dropna().corr().round(2)
    sig = lambda r: "significativo al 5 %" if r.p < 0.05 else "no significativo al 5 %"
    o_a = oos["auto"]

    result = {
        "title": "Tasa de Banxico y demanda de crédito",
        "summary": [
            f"Datos oficiales de la API SIE de Banxico ({meses}, {len(d)} meses): tasa objetivo, inflación y cartera "
            f"vigente de la banca comercial (automotriz y consumo), deflactada con el INPC. "
            + (f"Calidad de datos: {len(rupturas)} salto(s) de nivel que no es demanda, detectado(s) y empalmado(s) "
               f"({', '.join(r['Cartera'] + ' ' + r['Mes'] + ' ' + r['Cambio observado'] for r in rupturas)}), "
               "típicos de reclasificaciones o de entidades que se integran a la banca." if rupturas else ""),
            f"Proyecciones locales: tras subir la tasa 1 pp, el saldo automotriz real es {b12['auto'].b:+.1f} % a 12 meses "
            f"(IC {b12['auto'].lo:.1f} a {b12['auto'].hi:.1f}; {sig(b12['auto'])}); consumo {b12['consumo'].b:+.1f} % "
            f"({sig(b12['consumo'])}). "
            + ("El signo positivo NO significa que subir tasas genere crédito: Banxico sube la tasa cuando la economía y "
               "la demanda están fuertes (2022–2023 coincidió con la recuperación de ventas de autos tras la escasez de "
               "chips). Es endogeneidad, no causalidad." if b12["auto"].b > 0 else
               "El signo negativo es el esperado: tasas más altas encarecen el financiamiento."),
            f"VAR({p_var}) en diferencias ({'estable' if estable else 'NO estable'}): Granger tasa → crédito automotriz "
            f"p = {granger.iloc[0]['p-valor']:.3f} y crédito → tasa p = {granger.iloc[2]['p-valor']:.3f}"
            + (": la relación va en ambos sentidos, otra señal de política monetaria que reacciona al ciclo."
               if granger.iloc[2]["p-valor"] < 0.05 else "."),
            f"Prueba de fuego fuera de muestra ({o_a['n']} pronósticos a 12 meses desde {o_a['desde']:%Y-%m}): agregar "
            f"tasa e inflación {'reduce' if o_a['rmse_tasa'] < o_a['rmse_ar'] else 'NO reduce'} el error del crecimiento "
            f"automotriz ({o_a['rmse_ar']:.1f} → {o_a['rmse_tasa']:.1f} pp; Diebold-Mariano p = {o_a['dm_p']:.2f}). "
            f"Pronóstico AR del saldo real a 12 meses: {g_hat:+.1f} % (P10–P90 {g_hat + q10:+.1f} a {g_hat + q90:+.1f} %).",
            "Conclusión para la dirección: para planear originación conviene la inercia propia de la cartera, no la "
            "tasa. Medir el efecto causal de la tasa requiere choques de política identificados (sorpresas frente a "
            "las expectativas del mercado); por eso no se publican escenarios de ±100 pb.",
        ],
        "metrics": [
            {"label": "Asociación +1 pp a 12 m (auto)", "value": f"{b12['auto'].b:+.1f} %"},
            {"label": "Granger tasa → auto (p)", "value": f"{granger.iloc[0]['p-valor']:.3f}"},
            {"label": "Mejora RMSE fuera de muestra", "value": f"{1 - o_a['rmse_tasa'] / o_a['rmse_ar']:+.0%}"},
            {"label": "Meses analizados", "value": str(len(d))},
        ],
        "charts": [
            line("tasa", "Tasa objetivo e inflación anual", {"Tasa objetivo": serie_t(d.tasa),
                                                              "Inflación anual (INPC)": serie_t(d.inflacion)},
                 x_label="Mes", y_label="%", fmt="pct"),
            line("credito", "Crecimiento anual real de la cartera vigente", {CREDITOS[c]: serie_t(g12[c]) for c in CREDITOS},
                 x_label="Mes", y_label="% anual, pesos constantes", fmt="pct"),
            line("lp", "Respuesta acumulada del saldo real a +1 pp de tasa (proyecciones locales)", lp_series,
                 x_label="Meses después del cambio de tasa", y_label="% de cambio en el saldo real", fmt="dec",
                 bands=lp_bandas, subtitle="Banda: IC 95 % con errores HAC (Newey-West)"),
            line("irf", "VAR: impulso-respuesta acumulado del crédito automotriz", {
                "Crédito automotriz": [[h, round(float(acum[h]), 3)] for h in range(H + 1)]},
                 x_label="Meses", y_label="% acumulado ante +1 pp", fmt="dec",
                 bands={"Crédito automotriz": banda_var},
                 subtitle="Choque ortogonal (Cholesky: inflación → tasa → crédito); banda Monte Carlo 95 %"),
            line("pronostico", "Saldo automotriz real: pronóstico a 12 meses", esc, x_label="Mes",
                 y_label="Miles de millones de pesos constantes", fmt="num",
                 bands={"Pronóstico AR a 12 meses": banda_pron}, markers=[[f"{ult:%Y-%m}", "Último dato"]],
                 subtitle="Modelo AR directo (el mejor fuera de muestra); banda: percentiles 10–90 de sus errores"),
            bars("corr", "Correlación con la tasa objetivo (niveles)", corr.columns.drop("tasa"),
                 corr.loc["tasa"].drop("tasa").values, fmt="dec"),
        ],
        "tables": [table("Proyecciones locales: efecto acumulado de +1 pp de tasa", tabla_lp),
                   table("Causalidad de Granger (VAR en diferencias)", granger),
                   table("¿La tasa mejora el pronóstico a 12 meses? (origen móvil)", tabla_oos),
                   table("Raíz unitaria (ADF): por qué se modela en diferencias", adf)]
                  + ([table("Calidad de datos: saltos de nivel detectados y empalmados", pd.DataFrame(rupturas))]
                     if rupturas else []),
        "card": model_card(
            "Proyecciones locales de Jordà con errores HAC y VAR en diferencias (lags por AIC) con impulso-respuesta "
            "ortogonal; validación fuera de muestra con origen móvil y prueba de Diebold-Mariano.",
            ["ADF para decidir niveles vs. diferencias; residuales del VAR sin autocorrelación (Portmanteau "
             f"p = {resid_p:.2f}) y estabilidad {'verificada' if estable else 'NO verificada'}.",
             "Pronósticos directos a 12 meses, re-estimados mes a mes solo con información disponible en cada fecha."],
            ["El orden de Cholesky supone que la tasa no reacciona en el mismo mes al crédito.",
             "La cartera vigente (saldo) aproxima la demanda; también refleja oferta y castigos."],
            ["Es evidencia de asociación dinámica, no un experimento: otros choques (pandemia 2020, nearshoring) "
             "ocurren al mismo tiempo que los ciclos de tasa.",
             "Datos agregados de la banca comercial; una financiera automotriz cautiva puede responder distinto "
             "(promociones de tasa subsidiada por la armadora)."]),
    }
    sample = d[["tasa", "inflacion", "auto_real", "consumo_real"]].reset_index()
    sample["fecha"] = sample["fecha"].dt.strftime("%Y-%m")
    export("p11-banxico-credito", result, sample=sample.round(3))
    return result


if __name__ == "__main__":
    run()