Pronóstico de ventas de autos en México

"""#10 Pronóstico de ventas de vehículos ligeros en México (y participación de Nissan).

Datos: INEGI – Registro Administrativo de la Industria Automotriz de Vehículos Ligeros (RAIAVL),
consultados directo de la API PX-Web de INEGI (venta al público por marca y modelo, mensual).
Modelos: ingenuo estacional, ETS, SARIMA, LightGBM con rezagos y ensamble, evaluados con walk-forward.
"""
import io
import json
import urllib.request
import warnings

import numpy as np
import pandas as pd
from lightgbm import LGBMRegressor
from scipy import stats
from statsmodels.stats.diagnostic import acorr_ljungbox
from statsmodels.tsa.holtwinters import ExponentialSmoothing
from statsmodels.tsa.seasonal import STL
from statsmodels.tsa.statespace.sarimax import SARIMAX

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

API = "https://www.inegi.org.mx/app/tabulados/pxwebv2/api/v1/es/RAIAVL/RAIAVL_8_9.px"
MESES = ["Enero", "Febrero", "Marzo", "Abril", "Mayo", "Junio", "Julio", "Agosto", "Septiembre",
         "Octubre", "Noviembre", "Diciembre"]
HORIZONTE = 12


def descargar_inegi() -> pd.DataFrame:
    """Descarga ventas mensuales por marca (año por año) y las guarda en caché."""
    cache = RAW / "inegi_raiavl_ventas_marca.csv"
    if cache.exists():
        return pd.read_csv(cache, parse_dates=["fecha"])
    RAW.mkdir(parents=True, exist_ok=True)
    meta = json.load(urllib.request.urlopen(API, timeout=120))
    anio, mes, modelo, unidades = [x["code"] for x in meta["variables"]]
    anios = meta["variables"][0]["values"]
    partes = []
    for a in anios:
        q = {"query": [{"code": anio, "selection": {"filter": "item", "values": [a]}},
                       {"code": mes, "selection": {"filter": "all", "values": ["*"]}},
                       {"code": modelo, "selection": {"filter": "all", "values": ["*"]}},
                       {"code": unidades, "selection": {"filter": "item", "values": ["0"]}}],  # 0 = Ventas
             "response": {"format": "csv"}}
        req = urllib.request.Request(API, data=json.dumps(q).encode(),
                                     headers={"Content-Type": "application/json"})
        txt = urllib.request.urlopen(req, timeout=300).read().decode("utf-8-sig")
        d = pd.read_csv(io.StringIO(txt), na_values=["-"])
        d.columns = ["anio", "mes", "modelo", "ventas"]
        d["marca"] = d["modelo"].str.split(" - ").str[0].str.strip()
        d["fecha"] = pd.to_datetime(dict(year=d.anio, month=d.mes.map(MESES.index) + 1, day=1))
        partes.append(d.groupby(["fecha", "marca"], as_index=False)["ventas"].sum(min_count=1).dropna())
        print(f"  ↓ INEGI {d.anio.iloc[0]}")
    df = pd.concat(partes, ignore_index=True)
    df.to_csv(cache, index=False)
    return df


def features_lag(s: pd.Series) -> pd.DataFrame:
    f = pd.DataFrame({"y": s})
    for l in [1, 2, 3, 6, 12, 24]:
        f[f"lag_{l}"] = s.shift(l)
    f["media_3"] = s.shift(1).rolling(3).mean()
    f["media_12"] = s.shift(1).rolling(12).mean()
    f["mes"] = s.index.month
    return f


def pronostico_lgbm(hist: pd.Series, h: int) -> np.ndarray:
    """Pronóstico recursivo: cada predicción alimenta los rezagos del siguiente mes."""
    s = hist.copy()
    f = features_lag(s).dropna()
    m = LGBMRegressor(n_estimators=400, learning_rate=0.03, num_leaves=8, min_child_samples=10,
                      random_state=SEED, verbose=-1).fit(f.drop(columns="y"), f["y"])
    for _ in range(h):
        nxt = s.index[-1] + pd.offsets.MonthBegin(1)
        s.loc[nxt] = np.nan
        x = features_lag(s).iloc[[-1]].drop(columns="y")
        s.loc[nxt] = m.predict(x)[0]
    return s.iloc[-h:].values


def covid(idx: pd.DatetimeIndex) -> pd.DataFrame:
    """Variables de intervención: cierre (abr–jun 2020) y recuperación con escasez de chips (jul-2020 a dic-2021)."""
    return pd.DataFrame({"cierre": ((idx >= "2020-04-01") & (idx <= "2020-06-01")).astype(float),
                         "escasez": ((idx >= "2020-07-01") & (idx <= "2021-12-01")).astype(float)}, index=idx)


def pronostico_sarima(hist: pd.Series, h: int, alpha=0.05):
    """SARIMA(1,1,1)(0,1,1)12 en logaritmo con intervención COVID; devuelve media e intervalo de predicción."""
    with warnings.catch_warnings():   # avisos de convergencia del optimizador en orígenes cortos
        warnings.simplefilter("ignore")
        m = SARIMAX(np.log(hist), exog=covid(hist.index), order=(1, 1, 1), seasonal_order=(0, 1, 1, 12)).fit(disp=False)
    fi = pd.date_range(hist.index[-1] + pd.offsets.MonthBegin(1), periods=h, freq="MS")
    f = m.get_forecast(h, exog=covid(fi))
    ci = np.exp(f.conf_int(alpha=alpha).values)
    return np.exp(f.predicted_mean.values), ci[:, 0], ci[:, 1], m


def pronostico_ets(hist: pd.Series, h: int) -> np.ndarray:
    m = ExponentialSmoothing(np.log(hist), trend="add", damped_trend=True, seasonal="add",
                             seasonal_periods=12).fit()
    return np.exp(m.forecast(h).values)


def mape(y, p):
    return float(np.mean(np.abs((y - p) / y)) * 100)


def mase(y, p, hist):
    """MASE: error relativo al ingenuo estacional dentro de la historia (< 1 = mejor que el ingenuo)."""
    escala = np.mean(np.abs(hist.values[12:] - hist.values[:-12]))
    return float(np.mean(np.abs(y - p)) / escala)


def run():
    df = descargar_inegi()
    total = df.groupby("fecha")["ventas"].sum().asfreq("MS")
    total = total[total > 0].asfreq("MS")
    nissan = df[df.marca.str.upper() == "NISSAN"].groupby("fecha")["ventas"].sum().reindex(total.index)
    share = (nissan / total * 100)
    y = total.loc["2010":]

    # --- 1) Validación walk-forward: 12 orígenes mensuales consecutivos, horizonte de 12 meses ---
    #     (ventana expansiva; cada origen solo ve el pasado)
    N_ORIG = 12
    origenes = [len(y) - HORIZONTE - k for k in range(N_ORIG - 1, -1, -1)]
    modelos = ["Estacional ingenuo", "ETS amortiguado", "SARIMA", "LightGBM (rezagos)", "Ensamble SARIMA+ETS"]
    filas, por_h, cobertura = [], {m: np.zeros(HORIZONTE) for m in modelos}, []
    for o in origenes:
        hist, real = y.iloc[:o], y.iloc[o:o + HORIZONTE].values
        sar, lo, hi, _ = pronostico_sarima(hist, HORIZONTE)
        ets = pronostico_ets(hist, HORIZONTE)
        preds = {"Estacional ingenuo": hist.iloc[-12:].values, "ETS amortiguado": ets, "SARIMA": sar,
                 "LightGBM (rezagos)": pronostico_lgbm(hist, HORIZONTE), "Ensamble SARIMA+ETS": (sar + ets) / 2}
        for mname, pr in preds.items():
            filas.append({"origen": y.index[o], "modelo": mname, "MAPE": mape(real, pr), "MASE": mase(real, pr, hist)})
            por_h[mname] += np.abs((real - pr) / real) * 100 / N_ORIG
        cobertura.append(np.mean((real >= lo) & (real <= hi)))
    ev = pd.DataFrame(filas)
    res = ev.groupby("modelo", sort=False)[["MAPE", "MASE"]].agg(["mean", "std"])
    mejor = res[("MAPE", "mean")].idxmin()
    cob = float(np.mean(cobertura) * 100)
    # Prueba pareada del mejor modelo contra el ingenuo estacional (errores por origen)
    retador = res.drop(index="Estacional ingenuo")[("MAPE", "mean")].idxmin()
    p_dm = float(stats.wilcoxon(ev[ev.modelo == retador].MAPE.values,
                                ev[ev.modelo == "Estacional ingenuo"].MAPE.values).pvalue)

    # --- 2) Diagnóstico de residuos del SARIMA y descomposición estacional ---
    _, _, _, m_full = pronostico_sarima(y, 1)
    lb = acorr_ljungbox(m_full.resid.iloc[13:], lags=[12, 24], return_df=True)
    stl = STL(np.log(y), period=12, robust=True).fit()
    est_mes = pd.Series(np.exp(stl.seasonal.values) * 100 - 100, index=y.index.month).groupby(level=0).mean()

    # --- 3) Pronóstico a futuro con intervalos (80 % y 95 %) ---
    fut_idx = pd.date_range(y.index[-1] + pd.offsets.MonthBegin(1), periods=HORIZONTE, freq="MS")
    fut, lo95, hi95, _ = pronostico_sarima(y, HORIZONTE, 0.05)
    _, lo80, hi80, _ = pronostico_sarima(y, HORIZONTE, 0.20)
    ets_f = pronostico_ets(y, HORIZONTE)
    hist, real = y.iloc[:-HORIZONTE], y.iloc[-HORIZONTE:]
    p_sar, b_lo, b_hi, _ = pronostico_sarima(hist, HORIZONTE)
    p_ets = pronostico_ets(hist, HORIZONTE)

    def pts(idx, vals):
        return [[d.strftime("%Y-%m"), float(round(v))] for d, v in zip(idx, vals)]
    def banda(idx, lo, hi):
        return [[d.strftime("%Y-%m"), float(round(a)), float(round(b))] for d, a, b in zip(idx, lo, hi)]

    ult = y.index[-1]
    yoy = y.iloc[-12:].sum() / y.iloc[-24:-12].sum() - 1
    anual = total.groupby(total.index.year).sum()
    anual = anual[anual.index < ult.year]
    share_anual = (nissan.groupby(nissan.index.year).sum() / total.groupby(total.index.year).sum() * 100)
    share_anual = share_anual[share_anual.index < ult.year]
    ult12 = df[df.fecha > ult - pd.DateOffset(months=12)].groupby("marca")["ventas"].sum()
    marcas = (ult12 / ult12.sum() * 100).sort_values(ascending=False).head(10)
    res_tbl = pd.DataFrame({"Modelo": res.index,
                            "MAPE % (media ± desv.)": [f"{m:.2f} ± {s:.2f}" for m, s in res["MAPE"].values],
                            "MASE (media ± desv.)": [f"{m:.2f} ± {s:.2f}" for m, s in res["MASE"].values]})
    lb_tbl = pd.DataFrame({"Rezagos": lb.index, "Estadístico Q": lb.lb_stat.round(2), "p-valor": lb.lb_pvalue.round(4)})
    meses_c = ["Ene", "Feb", "Mar", "Abr", "May", "Jun", "Jul", "Ago", "Sep", "Oct", "Nov", "Dic"]

    result = {
        "title": "Pronóstico de ventas de autos en México",
        "summary": [
            f"Serie oficial INEGI RAIAVL de {total.index[0]:%Y} a {ult:%m/%Y}: últimos 12 meses "
            f"{y.iloc[-12:].sum():,.0f} unidades ({yoy:+.1%} vs. año previo).",
            f"Walk-forward con {N_ORIG} orígenes × 12 meses y 5 modelos: el mejor es {mejor} con MAPE "
            f"{res.loc[mejor, ('MAPE', 'mean')]:.1f} % y MASE {res.loc[mejor, ('MASE', 'mean')]:.2f}. Mejor modelo "
            f"estadístico: {retador} ({res.loc[retador, ('MAPE', 'mean')]:.1f} %); diferencia vs. ingenuo "
            f"Wilcoxon p = {p_dm:.3f}.",
            "En el régimen estable reciente el ingenuo estacional es muy difícil de vencer; la intervención COVID "
            "mejora al SARIMA y le da intervalos útiles. "
            f"Cobertura empírica del IC 95 %: {cob:.0f} %"
            + (" (intervalos conservadores)" if cob > 97 else "") +
            f"; Ljung-Box (12 rezagos) p = {lb.lb_pvalue.iloc[0]:.3f}"
            + (" → residuos sin autocorrelación." if lb.lb_pvalue.iloc[0] > 0.05 else " → queda estructura sin modelar."),
            f"Pronóstico próximos 12 meses: {fut.sum():,.0f} unidades (IC 95 % mensual en la gráfica). Diciembre vende "
            f"{est_mes.loc[12]:+.0f} % vs. un mes típico. Nissan lidera con {share.iloc[-12:].mean():.1f} % del mercado.",
        ],
        "metrics": [
            {"label": "Ventas 12 m", "value": f"{y.iloc[-12:].sum() / 1e6:.2f} M"},
            {"label": f"MAPE {mejor.split(' (')[0]}", "value": f"{res.loc[mejor, ('MAPE', 'mean')]:.1f} %"},
            {"label": "Cobertura IC 95 %", "value": f"{cob:.0f} %"},
            {"label": "Participación Nissan", "value": f"{share.iloc[-12:].mean():.1f} %"},
        ],
        "charts": [
            line("serie", "Ventas mensuales y pronóstico SARIMA con intervalos",
                 {"Real": pts(y.index[-60:], y.iloc[-60:].values), "Pronóstico": pts(fut_idx, fut),
                  "ETS": pts(fut_idx, ets_f)},
                 x_label="Mes", y_label="Unidades", fmt="int", bands={"Pronóstico": banda(fut_idx, lo95, hi95)},
                 subtitle="Banda: intervalo de predicción al 95 %"),
            line("backtest", "Backtest de los últimos 12 meses (con IC 95 % del SARIMA)",
                 {"Real": pts(real.index, real.values), "SARIMA": pts(real.index, p_sar), "ETS": pts(real.index, p_ets)},
                 x_label="Mes", y_label="Unidades", fmt="int", bands={"SARIMA": banda(real.index, b_lo, b_hi)}),
            line("horiz", "Error (MAPE %) según el horizonte de pronóstico",
                 {m: [[h + 1, round(float(v), 2)] for h, v in enumerate(por_h[m])] for m in
                  ["Estacional ingenuo", "SARIMA", "ETS amortiguado", "LightGBM (rezagos)"]},
                 x_label="Meses hacia adelante", y_label="MAPE %", fmt="pct"),
            bars("mape", "MAPE promedio por modelo (walk-forward)", res.index, res[("MAPE", "mean")].round(2).values, fmt="pct"),
            bars("estac", "Estacionalidad: % sobre / bajo un mes típico (STL)", meses_c, est_mes.round(1).values, fmt="pct"),
            line("share", "Participación de mercado de Nissan por año",
                 {"Nissan": [[str(a), round(float(v), 2)] for a, v in share_anual.items()]},
                 x_label="Año", y_label="% de las ventas totales", fmt="pct"),
            bars("marcas", "Participación por marca, últimos 12 meses", marcas.index, marcas.values.round(2), fmt="pct"),
            bars("anual", "Ventas anuales del mercado", anual.index.astype(str), anual.values, fmt="int"),
        ],
        "tables": [table("Walk-forward: error por modelo (12 orígenes)", res_tbl),
                   table("Diagnóstico de residuos SARIMA (Ljung-Box)", lb_tbl)],
        "card": model_card(
            "Ingenuo estacional, ETS amortiguado, SARIMA(1,1,1)(0,1,1)12 en logaritmo con intervención COVID, "
            "LightGBM con rezagos (recursivo) y ensamble.",
            [f"Validación walk-forward con ventana expansiva: {N_ORIG} orígenes consecutivos × 12 meses.",
             "MAPE y MASE (escala libre y comparable con el ingenuo); error por horizonte.",
             "Cobertura empírica de los intervalos de predicción; Ljung-Box sobre residuos; STL para estacionalidad.",
             "Prueba de Wilcoxon pareada contra el ingenuo estacional."],
            ["La estructura estacional de los últimos años se mantiene.", "No se incluyen variables exógenas (tasas, tipo de cambio)."],
            ["2020 (COVID) es un choque atípico dentro de la historia de entrenamiento.",
             "Las cifras del último mes son preliminares en INEGI y pueden revisarse."]),
    }
    sample = pd.DataFrame({"fecha": total.index.strftime("%Y-%m"), "ventas": total.values,
                           "nissan": nissan.values})
    export("p10-autos-mexico", result, sample)
    return result


if __name__ == "__main__":
    run()