"""Phase 6 — robustness suite (analysis_plan §6; identification threats and sensitivity).

Every variant is re-derived in memory from data/interim (price_log, vehicles,
listings_full) and data/processed (listings, risk_set); no existing module or
processed file is modified. Blocks:

  1. Hedonic H1/H4 (Model A2: splines + make-model FE, dealer-clustered):
     sample variants (sampled portals, Kijiji, 1 000-500 000 $ price filter),
     region joint Wald under each, heterogeneity of the linear age/km slopes,
     composition reweighting to the otogo make-model-year mix (pyfixest weights).
  2. Hazards H2/H5 (discrete-time cloglog, GLM, dealer-clustered): markdown
     thresholds 2 %/200 $ and 0.5 %/50 $, transient rule 24 h / 72 h, exit
     confirmation 3 / 14 d, S2L vs S2. Listing-days are subsampled to <= MAX_ROWS
     with a fixed seed (reported).
  3. AMLI H6: rank standardisation, MIN_LISTING_DAYS = 100, 2-week aggregation,
     province vs region, leave-one-component-out predictive validation.
  4. Holm correction over the six primary preregistered tests.
  5. paper/tables/tab09_robustness.tex (booktabs fragment).
"""
from __future__ import annotations

import json
import time
import warnings

import numpy as np
import pandas as pd
import polars as pl
import pyfixest as pf
import statsmodels.api as sm
from scipy.stats import chi2, norm

from . import amli, config as C, features, hazards, hedonic, pricepaths

OUT = C.RESULTS / "robustness"
warnings.filterwarnings("ignore", message=".*multicollinearity.*")   # na6/nk6 (redundant orthogonalised spline columns) as in hedonic.py
SEED = 36
MAX_ROWS = 800_000
HAZ_TERMS = ["ldens", "resid_A2", "lprice", "lkm"]
NUM_B = ["ldens", "resid_A2", "pricepos_entry", "nocomp_entry", "lprice", "lkm", "age_years_entry", "own_stock"]
NUM_C = ["ldens", "resid_A2", "pricepos_25", "nocomp", "prior_cuts", "lprice", "lkm", "age_years_entry", "own_stock"]
CATS = ["age_bin", "powertrain", "body_filled", "source_class", "region", "cal_week"]
_T0 = time.time()


def _log(msg: str) -> None:
    print(f"[{time.time() - _T0:7.1f}s] {msg}", flush=True)


# =============================================================================== 1. hedonic

def _wald(m, cols: list[str]) -> dict:
    names = list(m._coefnames)
    idx = [names.index(c) for c in cols if c in names]
    if not idx:
        return {"stat": None, "df": 0, "p": None}
    b = np.asarray(m.coef())[idx]; V = np.asarray(m._vcov)[np.ix_(idx, idx)]
    W = float(b @ np.linalg.solve(V, b))
    return {"stat": W, "df": len(idx), "p": float(chi2.sf(W, len(idx)))}


def _term(m, term: str) -> dict:
    t = m.tidy().reset_index(); t = t.rename(columns={t.columns[0]: "term"})
    r = t[t["term"] == term]
    if r.empty:
        return {"b": None, "se": None, "lo": None, "hi": None, "p": None}
    r = r.iloc[0]
    return {"b": float(r["Estimate"]), "se": float(r["Std. Error"]), "lo": float(r["2.5%"]), "hi": float(r["97.5%"]), "p": float(r["Pr(>|t|)"])}


def _region_premia(m) -> pd.DataFrame:
    t = m.tidy().reset_index(); t = t.rename(columns={t.columns[0]: "term"})
    rp = t[t["term"].str.startswith("C(region)")].copy()
    rp["region"] = rp["term"].str.extract(r"\[T\.(.*)\]")
    return rp[["region", "Estimate", "Std. Error", "2.5%", "97.5%"]].rename(columns={"Estimate": "est", "Std. Error": "se", "2.5%": "lo", "97.5%": "hi"})


def hedonic_frame(flt: pl.Expr) -> pd.DataFrame:
    """Same preparation as hedonic.load() but with an arbitrary sample filter."""
    L = pl.read_parquet(C.PROCESSED / "listings.parquet").filter(flt)
    df = L.select(["listing_id", "dealer_id", "source", "source_class", "make", "model", "trim", "year", "entry_price", "mileage_km",
                   "transmission", "drivetrain", "powertrain", "body_filled", "region", "entry_ts", "age_years_entry"]).to_pandas()
    df = df.rename(columns={"body_filled": "body"})
    df["body"] = df["body"].fillna("unknown"); df["trim"] = df["trim"].fillna("(none)"); df["region"] = df["region"].fillna("(missing)")
    df["logp"] = np.log(df["entry_price"]); df["age"] = df["age_years_entry"]; df["km"] = df["mileage_km"] / 1e4
    df["week"] = ((df["entry_ts"] - C.WINDOW_START_TS) // (7 * 86400)).astype(int)
    df["trim"] = df["trim"].str.replace(r"\s+", " ", regex=True).str.strip().str.slice(0, 60)
    mm_s = df["make"] + "|" + df["model"]; mmt_s = mm_s + "|" + df["trim"]
    df["mm"] = pd.factorize(mm_s)[0]; df["mm_label"] = mm_s; df["mmt"] = pd.factorize(mmt_s)[0]
    df["mmy_label"] = mm_s + "|" + df["year"].astype(str)
    df["region"] = pd.Categorical(df["region"], categories=["Montréal"] + sorted(r for r in df["region"].unique() if r != "Montréal"))
    for col, ref in (("powertrain", "ICE"), ("transmission", "Automatic"), ("drivetrain", "AWD/4x4"), ("source_class", "direct_dealer")):
        levels = ([ref] if ref in set(df[col]) else []) + sorted(x for x in df[col].unique() if x != ref)
        df[col] = pd.Categorical(df[col], categories=levels)
    return df.dropna(subset=["age", "km", "logp"])


def _rhs(acols, kcols, controls) -> str:
    return " + ".join(["age", "km"] + acols + kcols + [f"C({c})" for c in controls])


def fit_A2(df: pd.DataFrame, controls=hedonic.CONTROLS, weights: str | None = None, region: bool = True, fe: str = "mm") -> tuple[dict, pd.DataFrame]:
    """A2 (fe='mm': make-model FE) or A4 (fe='mmt': make-model-trim FE) with the same spline RHS."""
    X, acols, kcols, _ = hedonic.spline_basis(df)
    fml = f"logp ~ {_rhs(acols, kcols, controls)}" + (" + C(region)" if region else "") + f" | {fe} + week"
    kw = {"weights": weights} if weights else {}
    m = pf.feols(fml, data=X, vcov={"CRV1": "dealer_id"}, fixef_rm="none", **kw)
    res = {"n": int(m._N), "n_dealers": int(X["dealer_id"].nunique()), "r2": float(m._r2), "rmse": float(np.sqrt(np.mean(m.resid() ** 2))),
           "age": _term(m, "age"), "km": _term(m, "km"),
           "H1_nonlinear_joint": _wald(m, acols + kcols), "H1_nonlinear_age": _wald(m, acols), "H1_nonlinear_km": _wald(m, kcols)}
    names = list(m._coefnames)
    res["H4_region_joint"] = _wald(m, [n for n in names if n.startswith("C(region)")]) if region else None
    return res, (_region_premia(m) if region else pd.DataFrame())


def otogo_weights(df: pd.DataFrame, cap: float = 10.0) -> tuple[pd.Series, dict]:
    """Post-stratification weights: direct-dealer + banner listings reweighted to the otogo
    make-model-year distribution (fallback to the make-model cell; cells absent from otogo
    get weight 0 and are dropped); other sources keep weight 1."""
    tgt = df[df["source"] == "otogo"]; smp = df[df["source_class"].isin(["direct_dealer", "banner"])]
    t_mmy = tgt["mmy_label"].value_counts(normalize=True); s_mmy = smp["mmy_label"].value_counts(normalize=True)
    t_mm = tgt["mm_label"].value_counts(normalize=True); s_mm = smp["mm_label"].value_counts(normalize=True)
    w = pd.Series(1.0, index=df.index)
    r_mmy = (t_mmy / s_mmy).reindex(smp["mmy_label"]).values
    r_mm = (t_mm / s_mm).reindex(smp["mm_label"]).values
    ws = np.where(np.isfinite(r_mmy), r_mmy, np.where(np.isfinite(r_mm), r_mm, 0.0))
    n_fallback_mm = int((~np.isfinite(r_mmy) & np.isfinite(r_mm)).sum()); n_zero = int((ws == 0).sum()); n_capped = int((ws > cap).sum())
    ws = np.minimum(ws, cap)
    pos = ws > 0
    ws[pos] = ws[pos] / ws[pos].mean()          # mean 1 among retained dealer/banner listings
    w.loc[smp.index] = ws
    info = {"n_target_otogo": int(len(tgt)), "n_reweighted": int(len(smp)), "n_fallback_mm_cell": n_fallback_mm, "n_zero_weight_dropped": n_zero,
            "n_capped_at_10": n_capped, "weight_max": float(ws.max()), "weight_cv": float(ws[pos].std() / ws[pos].mean()), "cap": cap}
    return w, info


def run_hedonic() -> dict:
    OUT.mkdir(parents=True, exist_ok=True)
    out, rows, prem = {}, [], []
    base_flt = pl.col("s1_hedonic")
    s1core = (pl.col("f_auto") & pl.col("f_canonical") & pl.col("f_year") & pl.col("f_km") & pl.col("f_mm") & pl.col("f_nocollision"))
    variants = {
        "baseline_S1": base_flt,
        "plus_sampled_portals": pl.col("s1_all_sources") & (pl.col("source_class") != "private_marketplace"),
        "plus_kijiji": pl.col("s1_all_sources") & (pl.col("source_class") != "sampled_portal"),
        "plus_sampled_and_kijiji": pl.col("s1_all_sources"),
        "price_1000_500000": s1core & pl.col("entry_price").is_between(1000.0, 500000.0) & pl.col("f_not_sampled_or_private"),
    }
    base_prem = None
    for name, flt in variants.items():
        df = hedonic_frame(flt)
        res, rp = fit_A2(df)
        r4, rp4 = fit_A2(df, fe="mmt"); res["A4_mmt"] = {k: r4[k] for k in ("n", "r2", "age", "km", "H4_region_joint")}; res["H4_region_joint_A4"] = r4["H4_region_joint"]
        res["source_classes"] = sorted(df["source_class"].unique().tolist())
        if base_prem is None:
            base_prem = rp.set_index("region")["est"]
        res["region_premia_corr_with_baseline"] = float(rp.set_index("region")["est"].reindex(base_prem.index).corr(base_prem))
        res["region_premia_max_abs"] = float(rp["est"].abs().max())
        out[name] = res
        rows.append({"variant": name, **_flat_hed(res)}); rp.insert(0, "variant", name); prem.append(rp)
        _log(f"hedonic {name}: n={res['n']} age={res['age']['b']:.4f} km={res['km']['b']:.4f} H1p={res['H1_nonlinear_joint']['p']:.2e} H4p(A2)={res['H4_region_joint']['p']:.3f} H4p(A4)={res['H4_region_joint_A4']['p']:.3f}")
    # --- (f) composition reweighting --------------------------------------------
    df = hedonic_frame(base_flt)
    w, winfo = otogo_weights(df)
    df["w_otogo"] = w.values
    dfw = df[df["w_otogo"] > 0].copy()
    res_w, rp_w = fit_A2(dfw, weights="w_otogo")
    r4w, _ = fit_A2(dfw, weights="w_otogo", fe="mmt"); res_w["A4_mmt"] = {k: r4w[k] for k in ("n", "r2", "age", "km", "H4_region_joint")}; res_w["H4_region_joint_A4"] = r4w["H4_region_joint"]
    res_w["weights"] = winfo
    res_w["region_premia_corr_with_baseline"] = float(rp_w.set_index("region")["est"].reindex(base_prem.index).corr(base_prem)); res_w["region_premia_max_abs"] = float(rp_w["est"].abs().max())
    out["reweighted_to_otogo_mmy"] = res_w; rows.append({"variant": "reweighted_to_otogo_mmy", **_flat_hed(res_w)}); rp_w.insert(0, "variant", "reweighted_to_otogo_mmy"); prem.append(rp_w)
    _log(f"hedonic reweighted: n={res_w['n']} age={res_w['age']['b']:.4f} H4p={res_w['H4_region_joint']['p']:.3f} {winfo}")
    # direct+banner only: unweighted vs reweighted
    dd = df[df["source_class"].isin(["direct_dealer", "banner"])].copy()
    res_dd, rp_dd = fit_A2(dd); out["direct_banner_only_unweighted"] = res_dd; rows.append({"variant": "direct_banner_only_unweighted", **_flat_hed(res_dd)})
    ddw = dd[dd["w_otogo"] > 0].copy(); res_ddw, rp_ddw = fit_A2(ddw, weights="w_otogo"); out["direct_banner_only_reweighted"] = res_ddw; rows.append({"variant": "direct_banner_only_reweighted", **_flat_hed(res_ddw)})
    for nm, rp_ in (("direct_banner_only_unweighted", rp_dd), ("direct_banner_only_reweighted", rp_ddw)):
        rp_.insert(0, "variant", nm); prem.append(rp_)
        out[nm]["region_premia_corr_with_baseline"] = float(rp_.set_index("region")["est"].reindex(base_prem.index).corr(base_prem))
    _log(f"hedonic direct+banner: unweighted age={res_dd['age']['b']:.4f} H4p={res_dd['H4_region_joint']['p']:.3f}; reweighted age={res_ddw['age']['b']:.4f} H4p={res_ddw['H4_region_joint']['p']:.3f}")
    pd.DataFrame(rows).to_csv(OUT / "hedonic_variants.csv", index=False)
    pd.concat(prem).to_csv(OUT / "hedonic_region_premia_variants.csv", index=False)
    # --- (e) heterogeneity -------------------------------------------------------
    het_rows, het = [], {}
    df["yband"] = pd.cut(df["year"], bins=[0, 2015, 2020, 9999], labels=["<=2015", "2016-2020", "2021+"]).astype(str)
    groups = {"powertrain": ["ICE", "BEV", "HEV"], "yband": ["<=2015", "2016-2020", "2021+"], "body": ["SUV", "Sedan", "Pickup"],
              "source_class": ["direct_dealer", "full_portal", "banner"]}
    for gvar, levels in groups.items():
        controls = [c for c in hedonic.CONTROLS if c != gvar]
        for lv in levels:
            sub = df[df[gvar] == lv]
            if len(sub) < 1000:
                continue
            r, _ = fit_A2(sub, controls=controls)
            het[f"{gvar}={lv}"] = {k: r[k] for k in ("n", "n_dealers", "r2", "age", "km", "H1_nonlinear_joint", "H4_region_joint")}
            het_rows.append({"group_var": gvar, "group": lv, "n": r["n"], "age_b": r["age"]["b"], "age_se": r["age"]["se"], "age_lo": r["age"]["lo"], "age_hi": r["age"]["hi"],
                             "km_b": r["km"]["b"], "km_se": r["km"]["se"], "km_lo": r["km"]["lo"], "km_hi": r["km"]["hi"], "H1_joint_p": r["H1_nonlinear_joint"]["p"], "H4_region_p": r["H4_region_joint"]["p"]})
            _log(f"hedonic het {gvar}={lv}: n={r['n']} age={r['age']['b']:.4f} ({r['age']['se']:.4f}) km={r['km']['b']:.4f} ({r['km']['se']:.4f})")
    out["heterogeneity"] = het
    pd.DataFrame(het_rows).to_csv(OUT / "hedonic_heterogeneity.csv", index=False)
    pub = json.load(open(C.RESULTS / "hedonic" / "summary.json"))
    out["published"] = {"H1_nonlinear_joint_A2": pub["tests"]["H1_nonlinear_joint"], "H4_region_joint_A4": pub["tests"]["H4_region_joint"], "n_S1": pub["n_S1"]}
    return out


def _flat_hed(r: dict) -> dict:
    return {"n": r["n"], "n_dealers": r["n_dealers"], "r2": r["r2"], "rmse": r["rmse"], "age_b": r["age"]["b"], "age_se": r["age"]["se"], "age_lo": r["age"]["lo"], "age_hi": r["age"]["hi"],
            "km_b": r["km"]["b"], "km_se": r["km"]["se"], "km_lo": r["km"]["lo"], "km_hi": r["km"]["hi"],
            "H1_joint_stat": r["H1_nonlinear_joint"]["stat"], "H1_joint_p": r["H1_nonlinear_joint"]["p"],
            "H4_region_stat": r["H4_region_joint"]["stat"], "H4_region_df": r["H4_region_joint"]["df"], "H4_region_p": r["H4_region_joint"]["p"],
            "H4_region_p_A4_mmt": (r.get("H4_region_joint_A4") or {}).get("p"), "age_b_A4": ((r.get("A4_mmt") or {}).get("age") or {}).get("b"),
            "region_premia_max_abs": r.get("region_premia_max_abs"), "region_premia_corr_with_baseline": r.get("region_premia_corr_with_baseline")}


# =============================================================================== 2. hazards

def rebuild_events(transient_seconds: float, min_pct: float, min_abs: float) -> tuple[pd.DataFrame, pd.DataFrame]:
    """Re-run the sustained-path algorithm under alternative rules for the S2L universe.
    Returns (per-listing: first_cut_day, n_cuts, oscillating, entry_price) and (listing_id, t) cut days."""
    LF = pl.read_parquet(C.INTERIM / "listings_full.parquet").filter(pl.col("s2l_dynamic_with_stock")).select(["uid", "listing_id"])
    V = pl.read_parquet(C.INTERIM / "vehicles.parquet").select(["uid", "last_seen"])
    P = pl.read_parquet(C.INTERIM / "price_log.parquet").filter(pl.col("uid").is_in(LF["uid"].to_list()))
    ls = dict(zip(V["uid"].to_list(), V["last_seen"].to_list()))
    lid = dict(zip(LF["uid"].to_list(), LF["listing_id"].to_list()))
    grouped = P.sort(["uid", "ts"]).group_by("uid", maintain_order=True).agg([pl.col("ts"), pl.col("price")])
    rows, cut_rows = [], []
    for uid, ts, pr in grouped.iter_rows():
        ps = pricepaths.sustain(ts, pr, ls.get(uid, max(ts)), transient_seconds=transient_seconds)
        ev = pricepaths.classify(ps.levels, min_pct=min_pct, min_abs=min_abs)
        cuts = [e for e in ev if e["kind"] == "cut"]
        l = lid[uid]
        for e in cuts:
            cut_rows.append((l, int(features.day_of(np.array([e["ts"]]))[0])))
        rows.append({"listing_id": l, "first_cut_day_v": int(features.day_of(np.array([cuts[0]["ts"]]))[0]) if cuts else np.nan,
                     "n_cuts_v": len(cuts), "oscillating_v": ps.n_transients >= C.OSCILLATION_MAX_TRANSIENTS,
                     "entry_price_v": ev[0]["price"] if ev else np.nan, "first_cut_pct_v": cuts[0]["change_pct"] if cuts else np.nan})
    per = pd.DataFrame(rows)
    cd = pd.DataFrame(cut_rows, columns=["listing_id", "t"]).drop_duplicates(); cd["cut_v"] = 1.0
    return per, cd


def _cloglog(d: pd.DataFrame, y: str, numeric: list[str], cats: list[str]) -> dict:
    """Same estimator as hazards.cloglog (GLM cloglog, dealer-clustered) without writing to results/hazards."""
    t0 = time.time()
    d = d.dropna(subset=numeric + [y]).copy()
    X = hazards.design(d, numeric, cats); X = X.loc[:, (X != 0).any(axis=0)]
    model = sm.GLM(d[y].values.astype(float), X, family=sm.families.Binomial(link=sm.families.links.CLogLog()))
    res = model.fit(cov_type="cluster", cov_kwds={"groups": pd.factorize(d["dealer_id"])[0]}, maxiter=100)
    idx = [i for i, c in enumerate(X.columns) if c.startswith("age_bin_")]
    b = res.params.values[idx]; V = res.cov_params().values[np.ix_(idx, idx)]; W = float(b @ np.linalg.solve(V, b))
    ci = res.conf_int()
    return {"n_rows": int(len(d)), "n_events": int(d[y].sum()), "n_listings": int(d["listing_id"].nunique()), "n_dealers": int(d["dealer_id"].nunique()),
            "llf": float(res.llf), "joint_age_bins": {"stat": W, "df": len(idx), "p": float(chi2.sf(W, len(idx)))}, "seconds": round(time.time() - t0, 1),
            "coef": {c: {"b": float(res.params[c]), "se": float(res.bse[c]), "p": float(res.pvalues[c]), "lo": float(ci.loc[c, 0]), "hi": float(ci.loc[c, 1])} for c in HAZ_TERMS + ["prior_cuts", "pricepos_25", "pricepos_entry"] if c in X.columns}}


def _subsample(d: pd.DataFrame, max_rows: int = MAX_ROWS) -> pd.DataFrame:
    if len(d) <= max_rows:
        return d
    return d.sample(n=max_rows, random_state=SEED)


def _apply_events(d: pd.DataFrame, per: pd.DataFrame, cd: pd.DataFrame, drop_new_oscillating: bool) -> pd.DataFrame:
    """Rebuild first-cut risk set, first_cut, cut and prior_cuts from re-derived events."""
    d = d.drop(columns=[c for c in ("first_cut_day_v", "n_cuts_v", "oscillating_v", "entry_price_v", "first_cut_pct_v", "cut_v") if c in d.columns])
    d = d.merge(per, on="listing_id", how="left").merge(cd, on=["listing_id", "t"], how="left")
    if drop_new_oscillating:
        d = d[~d["oscillating_v"].fillna(False).astype(bool)]
    d["first_cut_day"] = d["first_cut_day_v"]
    d["in_first_cut_risk"] = d["first_cut_day"].isna() | (d["t"] <= d["first_cut_day"])
    d["first_cut"] = (d["t"] == d["first_cut_day"]).astype(float)
    d["cut"] = d["cut_v"].fillna(0.0)
    d = d.sort_values(["listing_id", "t"])
    d["prior_cuts"] = d.groupby("listing_id")["cut"].cumsum() - d["cut"]
    return d


def run_hazards() -> dict:
    OUT.mkdir(parents=True, exist_ok=True)
    out = {"subsample": {"max_rows": MAX_ROWS, "seed": SEED, "unit": "listing-days, simple random sample without replacement"}}
    _log("hazards: loading S2L risk set")
    full = hazards.load("s2l")                      # all rows of the risk set (S2L universe)
    # NOTE: in risk_set.parquet, `first_cut` is NULL for listings that never get a cut (t == NULL first_cut_day).
    # hazards.cloglog drops NaN outcomes, so the published Model B was fitted on listings that eventually got a cut only
    # (n_listings == n_events in results/hazards/summary.json). Here NULL first_cut is coded 0 (listing-day without a cut).
    out["note_first_cut_null"] = {"risk_set_first_cut_null_rows_S2": int(full.loc[full["s2_dynamic"], "first_cut"].isna().sum()),
                                  "published_B_n_rows": None, "published_B_n_listings": None, "published_B_n_events": None}
    pub = json.load(open(C.RESULTS / "hazards" / "summary.json"))
    out["note_first_cut_null"].update({"published_B_n_rows": pub["B_first_markdown_cloglog"]["n_rows"], "published_B_n_listings": pub["B_first_markdown_cloglog"]["n_listings"],
                                       "published_B_n_events": pub["B_first_markdown_cloglog"]["n_events"]})
    full["first_cut"] = full["first_cut"].fillna(0).astype(float)
    full["exit"] = full["exit"].astype(float); full["cut"] = full["cut"].astype(float)
    L = pl.read_parquet(C.PROCESSED / "listings.parquet").select(["listing_id", "last_seen_ts", "active"]).to_pandas()
    full = full.merge(L, on="listing_id", how="left")
    full["last_t"] = full.groupby("listing_id")["t"].transform("max")
    base_events = pl.read_parquet(C.PROCESSED / "price_events.parquet")
    S2 = full[full["s2_dynamic"]].copy()
    # reproduction check of exit coding (7-day rule): must equal the stored `exit`
    ex7 = ((S2["t"] == S2["last_t"]) & (S2["active"] == 0) & (S2["last_seen_ts"] <= C.SNAPSHOT_TS - 7 * 86400)).astype(float)
    out["exit_recode_check"] = {"stored_exits_S2": int(S2["exit"].sum()), "recoded_7d_exits_S2": int(ex7.sum()), "disagreements": int((ex7 != S2["exit"]).sum())}
    _log(f"hazards: S2 rows={len(S2)} exits={int(S2['exit'].sum())} recode check={out['exit_recode_check']}")

    results = {"B": {}, "C": {}}
    def fitBC(name: str, d: pd.DataFrame, note: str):
        dB = _subsample(d[d["in_first_cut_risk"]]); dC = _subsample(d)
        rB = _cloglog(dB, "first_cut", NUM_B, CATS); rB["note"] = note; rB["n_rows_before_subsample"] = int(d["in_first_cut_risk"].sum())
        rB["n_events_before_subsample"] = int(d.loc[d["in_first_cut_risk"], "first_cut"].sum())
        results["B"][name] = rB
        _log(f"B {name}: n={rB['n_rows']} ev={rB['n_events']} ageP={rB['joint_age_bins']['p']:.2e} " + " ".join(f"{k}={rB['coef'][k]['b']:+.3f}({rB['coef'][k]['se']:.3f})" for k in HAZ_TERMS))
        rC = _cloglog(dC, "exit", NUM_C, CATS); rC["note"] = note; rC["n_rows_before_subsample"] = int(len(d)); rC["n_events_before_subsample"] = int(d["exit"].sum())
        results["C"][name] = rC
        _log(f"C {name}: n={rC['n_rows']} ev={rC['n_events']} ageP={rC['joint_age_bins']['p']:.2e} " + " ".join(f"{k}={rC['coef'][k]['b']:+.3f}({rC['coef'][k]['se']:.3f})" for k in HAZ_TERMS))

    # published baselines (full sample)
    results["B"]["published_full"] = {k: pub["B_first_markdown_cloglog"][k] for k in ("n_rows", "n_events", "n_listings", "joint_age_bins")}
    results["B"]["published_full"]["coef"] = {k: pub["B_first_markdown_cloglog"]["coef"][k] for k in HAZ_TERMS}
    results["B"]["published_full"]["note"] = "results/hazards/summary.json as published (fitted on listings with a NULL-free first_cut, i.e. eventually-cut listings only)"
    results["C"]["published_full"] = {k: pub["C_exit_cloglog"][k] for k in ("n_rows", "n_events", "n_listings", "joint_age_bins")}
    results["C"]["published_full"]["coef"] = {k: pub["C_exit_cloglog"]["coef"][k] for k in HAZ_TERMS}
    results["C"]["published_full"]["note"] = "results/hazards/summary.json as published (full S2 risk set)"
    # corrected B on the full first-cut risk set (no subsample) — reference for the NULL-coding issue
    rBfull = _cloglog(S2[S2["in_first_cut_risk"]], "first_cut", NUM_B, CATS); rBfull["note"] = "full S2 first-cut risk set, NULL first_cut coded 0, no subsample"
    results["B"]["baseline_full_nullfix"] = rBfull
    _log(f"B baseline_full_nullfix: n={rBfull['n_rows']} ev={rBfull['n_events']} ageP={rBfull['joint_age_bins']['p']:.2e} " + " ".join(f"{k}={rBfull['coef'][k]['b']:+.3f}({rBfull['coef'][k]['se']:.3f})" for k in HAZ_TERMS))

    # --- baseline rules on the subsample ---------------------------------------------
    fitBC("baseline_sub", S2, "S2; 1 %/100 $; 48 h; exit 7 d; subsample")
    # reproduction check: re-deriving events under the baseline rules must give the stored first_cut_day
    per0, _ = rebuild_events(C.TRANSIENT_SECONDS, C.MARKDOWN_MIN_PCT, C.MARKDOWN_MIN_ABS)
    stored = S2.drop_duplicates("listing_id").set_index("listing_id")["first_cut_day"]
    chk = per0.set_index("listing_id")["first_cut_day_v"].reindex(stored.index)
    out["events_baseline_reproduction"] = {"listings_S2": int(len(stored)), "first_cut_day_disagreements": int(((chk.fillna(-1) != stored.fillna(-1))).sum()),
                                           "listings_with_cut_S2L": int(per0["first_cut_day_v"].notna().sum())}
    _log(f"events baseline reproduction: {out['events_baseline_reproduction']}")
    # --- (a) markdown thresholds, (b) transient rule -------------------------------------
    for name, (tr, pct, ab), drop_osc in [("threshold_2pct_200", (C.TRANSIENT_SECONDS, 0.02, 200.0), False), ("threshold_0.5pct_50", (C.TRANSIENT_SECONDS, 0.005, 50.0), False),
                                           ("transient_24h", (24 * 3600, C.MARKDOWN_MIN_PCT, C.MARKDOWN_MIN_ABS), True), ("transient_72h", (72 * 3600, C.MARKDOWN_MIN_PCT, C.MARKDOWN_MIN_ABS), True)]:
        per, cd = rebuild_events(tr, pct, ab)
        d = _apply_events(S2, per, cd, drop_new_oscillating=drop_osc)
        info = {"transient_seconds": tr, "min_pct": pct, "min_abs": ab, "listings_with_cut": int(per["first_cut_day_v"].notna().sum()), "oscillating_under_rule": int(per["oscillating_v"].sum()),
                "listings_dropped_as_newly_oscillating": int(S2["listing_id"].nunique() - d["listing_id"].nunique()) if drop_osc else 0,
                "entry_price_changed_share": float((np.abs(per["entry_price_v"] - per["listing_id"].map(S2.drop_duplicates("listing_id").set_index("listing_id")["entry_price"])) > 0.5).mean())}
        out[f"events_{name}"] = info
        _log(f"events {name}: {info}")
        fitBC(name, d, f"S2; rebuilt events ({tr/3600:.0f} h, {pct*100:g} %/{ab:g} $); lprice/resid_A2 from baseline entry price; subsample")
    # --- (c) exit confirmation 3 / 14 days ----------------------------------------------
    for k in (3, 14):
        d = S2.copy()
        d["exit"] = ((d["t"] == d["last_t"]) & (d["active"] == 0) & (d["last_seen_ts"] <= C.SNAPSHOT_TS - k * 86400)).astype(float)
        out[f"exit_confirm_{k}d"] = {"confirmed_exits_S2": int(d["exit"].sum()), "baseline_7d": int(S2["exit"].sum())}
        fitBC(f"exit_confirm_{k}d", d, f"S2; exit = inactive and last_seen <= snapshot - {k} d (exit on last observed listing-day); subsample")
    # --- (d) S2L ---------------------------------------------------------------------
    out["S2L"] = {"rows": int(len(full)), "listings": int(full["listing_id"].nunique()), "left_censored_listings": int(full.loc[full["left_censored"], "listing_id"].nunique()),
                  "note": "left-censored listings enter with listing_age measured from first_seen (lower bound)"}
    fitBC("S2L_with_stock", full, "S2L (S2 + left-censored stock, age from first_seen); subsample")
    out["results"] = results
    # flat CSV
    rows = []
    for mdl in ("B", "C"):
        for nm, r in results[mdl].items():
            row = {"model": mdl, "variant": nm, "n_rows": r["n_rows"], "n_events": r["n_events"], "n_listings": r["n_listings"], "joint_age_stat": r["joint_age_bins"]["stat"], "joint_age_p": r["joint_age_bins"]["p"]}
            for k in HAZ_TERMS:
                c = r["coef"].get(k, {})
                row.update({f"{k}_b": c.get("b"), f"{k}_se": c.get("se"), f"{k}_p": c.get("p")})
            row["note"] = r.get("note"); rows.append(row)
    pd.DataFrame(rows).to_csv(OUT / "hazards_variants.csv", index=False)
    return out


# =============================================================================== 3. AMLI

def _validate(Z: pd.DataFrame, key: list[str], zcols: list[str], eval_weeks) -> dict:
    """Leave-exit-out (equal-weight mean of zcols) vs next-period exit, market FE, cluster by market."""
    Z = Z.copy(); Z["idx"] = Z[zcols].mean(axis=1)
    Z["mkt"] = Z[key].astype(str).agg("|".join, axis=1)
    Z = Z.sort_values(["mkt", "week"]); Z["exit_next"] = Z.groupby("mkt")["exit"].shift(-1)
    d = Z[Z["week"].isin(list(eval_weeks))].dropna(subset=["idx", "exit_next", "exit"])
    nm = d["mkt"].nunique()
    if nm > 1:
        m = pf.feols("exit_next ~ idx | mkt", data=d, vcov={"CRV1": "mkt"})
    else:
        m = pf.feols("exit_next ~ idx", data=d, vcov="hetero")
    t = m.tidy().reset_index(); t = t.rename(columns={t.columns[0]: "term"}); r = t[t["term"] == "idx"].iloc[0]
    b, se = float(r["Estimate"]), float(r["Std. Error"])
    return {"b": b, "se": se, "lo": float(r["2.5%"]), "hi": float(r["97.5%"]), "p_normal": float(2 * norm.sf(abs(b / se))) if se > 0 else None, "n": int(m._N), "n_markets": int(nm), "r2": float(m._r2)}


def _rank_index(comp: pd.DataFrame, cols: list[str]) -> pd.DataFrame:
    """Rank standardisation: signed normal scores of the within-panel rank (replaces z-scores)."""
    Z = comp.copy()
    for c in cols:
        r = Z[c].rank(); n = Z[c].notna().sum()
        Z[f"z_{c}"] = amli.SIGN[c] * norm.ppf((r - 0.5) / n)
    zc = [f"z_{c}" for c in cols]
    ok = Z[zc].notna().sum(axis=1) >= min(3, len(cols))
    Z["amli_equal"] = Z.loc[ok, zc].mean(axis=1)
    return Z


def run_amli() -> dict:
    OUT.mkdir(parents=True, exist_ok=True)
    out, rows = {}, []
    d = amli.listing_days()
    others = [f"z_{c}" for c in amli.PRIMARY if c != "exit"]
    ev_w = range(6, amli.WEEKS)
    comp = amli.components(d, ["region"])
    Z, meta = amli.build_index(comp, amli.PRIMARY)
    base = _validate(Z, ["region"], others, ev_w)
    out["baseline"] = {"validation_noexit_fe": base, "pca_explained": meta["weights"]["pca_explained"], "n_market_weeks": int(Z["amli_equal"].notna().sum()),
                       "min_listing_days": amli.MIN_LISTING_DAYS, "weeks": [int(Z["week"].min()), int(Z["week"].max())]}
    rows.append({"variant": "baseline_z_region_weekly", **base, "corr_equal_with_baseline": 1.0})
    _log(f"AMLI baseline: slope={base['b']:.5f} ({base['se']:.5f}) n={base['n']}")
    # (a) rank standardisation
    Zr = _rank_index(comp, amli.PRIMARY)
    vr = _validate(Zr, ["region"], others, ev_w)
    out["rank_standardisation"] = {"validation_noexit_fe": vr, "corr_equal_with_baseline": float(Zr["amli_equal"].corr(Z["amli_equal"])),
                                   "spearman_equal_with_baseline": float(Zr["amli_equal"].corr(Z["amli_equal"], method="spearman"))}
    rows.append({"variant": "rank_standardisation", **vr, "corr_equal_with_baseline": out["rank_standardisation"]["corr_equal_with_baseline"]})
    _log(f"AMLI rank: slope={vr['b']:.5f} ({vr['se']:.5f}) corr={out['rank_standardisation']['corr_equal_with_baseline']:.3f}")
    # (b) MIN_LISTING_DAYS = 100
    comp100 = comp.copy(); comp100.loc[comp100["ld"] < 100, amli.PRIMARY + ["stale"]] = np.nan
    Z100, m100 = amli.build_index(comp100, amli.PRIMARY)
    v100 = _validate(Z100, ["region"], others, ev_w)
    out["min_listing_days_100"] = {"validation_noexit_fe": v100, "n_market_weeks": int(Z100["amli_equal"].notna().sum()), "n_market_weeks_baseline": int(Z["amli_equal"].notna().sum()),
                                   "cells_between_30_and_100_listing_days": int(((comp["ld"] >= 30) & (comp["ld"] < 100)).sum()), "min_listing_days_in_panel": float(comp["ld"].min()),
                                   "ld_quantiles": {q: float(comp["ld"].quantile(q)) for q in (0.0, 0.05, 0.25, 0.5)},
                                   "corr_equal_with_baseline": float(Z100["amli_equal"].corr(Z["amli_equal"])), "pca_explained": m100["weights"]["pca_explained"]}
    rows.append({"variant": "min_listing_days_100", **v100, "corr_equal_with_baseline": out["min_listing_days_100"]["corr_equal_with_baseline"]})
    _log(f"AMLI minld100: slope={v100['b']:.5f} ({v100['se']:.5f}) mw={out['min_listing_days_100']['n_market_weeks']}")
    # (c) 2-week aggregation (periods: weeks 2-3, 4-5, 6-7, 8-9, 10)
    d2 = d.copy(); d2["week"] = (d2["week"] - amli.WEEK_MIN) // 2
    comp2 = amli.components(d2, ["region"])
    Z2, m2 = amli.build_index(comp2, amli.PRIMARY, fit_weeks=range(0, 2))
    v2 = _validate(Z2, ["region"], others, range(2, 5))
    # comparable weekly index averaged over the same 2-week periods
    Zw = Z.copy(); Zw["period"] = (Zw["week"] - amli.WEEK_MIN) // 2
    avg = Zw.groupby(["region", "period"])["amli_equal"].mean().rename("amli_weekly_avg").reset_index().rename(columns={"period": "week"})
    cmp2 = Z2.merge(avg, on=["region", "week"], how="left")
    out["two_week_aggregation"] = {"validation_noexit_fe": v2, "periods": "weeks 2-3, 4-5, 6-7, 8-9, 10 (last period = 1 week)", "fit_periods": [0, 1], "eval_periods": [2, 3, 4],
                                   "corr_with_weekly_index_averaged": float(cmp2["amli_equal"].corr(cmp2["amli_weekly_avg"])), "pca_explained": m2["weights"]["pca_explained"],
                                   "n_market_periods": int(Z2["amli_equal"].notna().sum())}
    rows.append({"variant": "two_week_aggregation", **v2, "corr_equal_with_baseline": out["two_week_aggregation"]["corr_with_weekly_index_averaged"]})
    _log(f"AMLI 2-week: slope={v2['b']:.5f} ({v2['se']:.5f}) n={v2['n']}")
    # (d) province vs region
    dP = d.copy(); dP["market"] = "Quebec"
    compP = amli.components(dP, ["market"]); ZP, mP = amli.build_index(compP, amli.PRIMARY)
    wavg = Z.dropna(subset=["amli_equal"]).assign(wx=lambda x: x["amli_equal"] * x["ld"]).groupby("week").agg(wx=("wx", "sum"), ld=("ld", "sum"))
    wavg["region_weighted_equal"] = wavg["wx"] / wavg["ld"]
    cmpP = ZP.set_index("week")[["amli_equal", "amli_pca"]].join(wavg["region_weighted_equal"])
    vP = _validate(ZP, ["market"], others, ev_w)
    vP_all = _validate(ZP, ["market"], others, range(amli.WEEK_MIN, amli.WEEKS))
    out["province_vs_region"] = {"corr_province_equal_vs_region_weighted_equal": float(cmpP["amli_equal"].corr(cmpP["region_weighted_equal"])),
                                 "spearman": float(cmpP["amli_equal"].corr(cmpP["region_weighted_equal"], method="spearman")), "n_weeks": int(cmpP["amli_equal"].notna().sum()),
                                 "province_validation_pooled_eval_weeks": vP, "province_validation_pooled_all_weeks": vP_all, "region_validation": base,
                                 "note": "province index has one market: pooled OLS with HC SEs on 4 (eval) / 8 (all) weekly points, descriptive only",
                                 "series": cmpP.reset_index().to_dict(orient="list")}
    rows.append({"variant": "province_level_pooled", **vP, "corr_equal_with_baseline": out["province_vs_region"]["corr_province_equal_vs_region_weighted_equal"]})
    _log(f"AMLI province: corr={out['province_vs_region']['corr_province_equal_vs_region_weighted_equal']:.3f} slope(pooled)={vP['b']:.5f} n={vP['n']}")
    # (e) leave-one-component-out predictive validation (on top of leave-exit-out)
    loo = {}
    for c in [x for x in amli.PRIMARY if x != "exit"]:
        cols = [z for z in others if z != f"z_{c}"]
        loo[c] = {"baseline_z": _validate(Z, ["region"], cols, ev_w), "rank": _validate(Zr, ["region"], cols, ev_w), "min_ld_100": _validate(Z100, ["region"], cols, ev_w),
                  "two_week": _validate(Z2, ["region"], cols, range(2, 5))}
        rows.append({"variant": f"loo_drop_{c}", **loo[c]["baseline_z"], "corr_equal_with_baseline": None})
        _log(f"AMLI LOO drop {c}: slope={loo[c]['baseline_z']['b']:.5f} ({loo[c]['baseline_z']['se']:.5f})")
    out["leave_one_component_out"] = loo
    pubA = json.load(open(C.RESULTS / "amli" / "summary.json"))
    out["published"] = {"region_validation_noexit_fe": pubA["region_validation"]["noexit_fe"], "listing_days": pubA["listing_days"]}
    pd.DataFrame(rows).to_csv(OUT / "amli_variants.csv", index=False)
    return out


# =============================================================================== 4. Holm

def holm_primary_tests() -> pd.DataFrame:
    hed = json.load(open(C.RESULTS / "hedonic" / "summary.json")); haz = json.load(open(C.RESULTS / "hazards" / "summary.json"))
    comp = json.load(open(C.RESULTS / "competition" / "summary.json")); am = json.load(open(C.RESULTS / "amli" / "summary.json"))
    v = am["region_validation"]["noexit_fe"]; p6 = float(2 * norm.sf(abs(v["b"] / v["se"])))
    tests = [
        ("H1", "hedonic A2: joint Wald on nonlinear age/km spline terms", "results/hedonic/summary.json tests.H1_nonlinear_joint", hed["tests"]["H1_nonlinear_joint"]["p"]),
        ("H2", "Model B cloglog: joint test of listing-age bins", "results/hazards/summary.json B_first_markdown_cloglog.joint_age_bins", haz["B_first_markdown_cloglog"]["joint_age_bins"]["p"]),
        ("H3", "Model D LPM 25 km: joint test of competitor-cut lags", "results/competition/summary.json main_25km.lags_joint", comp["main_25km"]["lags_joint"]["p"]),
        ("H4", "hedonic A4: joint Wald on region effects", "results/hedonic/summary.json tests.H4_region_joint", hed["tests"]["H4_region_joint"]["p"]),
        ("H5", "Model C cloglog: coefficient on entry hedonic residual (resid_A2)", "results/hazards/summary.json C_exit_cloglog.coef.resid_A2", haz["C_exit_cloglog"]["coef"]["resid_A2"]["p"]),
        ("H6", "AMLI leave-exit-out index vs next-week exit (region FE), normal p from b/se", "results/amli/summary.json region_validation.noexit_fe", p6),
    ]
    df = pd.DataFrame(tests, columns=["hypothesis", "test", "source", "p_raw"])
    m = len(df); df = df.sort_values("p_raw").reset_index(drop=True); df["rank"] = np.arange(1, m + 1)
    # Holm step-down: adj_(i) = max_{j<=i} min(1, (m-j+1) p_(j))
    adj = np.maximum.accumulate(np.minimum(1.0, (m - df["rank"].values + 1) * df["p_raw"].values))
    df["p_holm"] = adj; df["reject_0.05_raw"] = df["p_raw"] < 0.05; df["reject_0.05_holm"] = df["p_holm"] < 0.05
    df = df.sort_values("hypothesis").reset_index(drop=True)
    df.to_csv(OUT / "holm_primary_tests.csv", index=False)
    return df


# =============================================================================== 5. table

def _fp(p) -> str:
    if p is None or not np.isfinite(p):
        return "--"
    return "$<$0.001" if p < 0.001 else f"{p:.3f}"


def _fe(b, se=None, d=3) -> str:
    if b is None or not np.isfinite(b):
        return "--"
    s = f"{b:.{d}f}"
    return s + (f" ({se:.{d}f})" if se is not None and np.isfinite(se) else "")


def _delta(a, b, d=3) -> str:
    if a is None or b is None or not (np.isfinite(a) and np.isfinite(b)):
        return "--"
    return f"{b - a:+.{d}f}"


def _verdict(pa, pb) -> str:
    if pa is None or pb is None:
        return "--"
    return "same" if (pa < 0.05) == (pb < 0.05) else "flips"


def _esc(s: str) -> str:
    return s.replace("%", r"\%").replace("$", r"\$").replace("_", r"\_")


def write_table(S: dict | None = None) -> None:
    if S is None:
        S = json.load(open(OUT / "summary.json"))
    H, Z, A, holm = S["hedonic"], S["hazards"]["results"], S["amli"], S["holm"]
    holm = pd.DataFrame(holm) if not isinstance(holm, pd.DataFrame) else holm
    b = H["baseline_S1"]
    rows = []
    def hed_row(label, key, what):
        v = H[key]
        if what == "age":
            rows.append((label, "age slope (log \\$/yr)", _fe(b["age"]["b"], b["age"]["se"], 4), _fe(v["age"]["b"], v["age"]["se"], 4), _delta(b["age"]["b"], v["age"]["b"], 4)))
        elif what == "km":
            rows.append((label, "km slope (log \\$/10k km)", _fe(b["km"]["b"], b["km"]["se"], 4), _fe(v["km"]["b"], v["km"]["se"], 4), _delta(b["km"]["b"], v["km"]["b"], 4)))
        elif what == "H4":
            rows.append((label, "H4 region joint $p$ (A4 trim FE / A2)", f"{_fp(b['H4_region_joint_A4']['p'])} / {_fp(b['H4_region_joint']['p'])}", f"{_fp(v['H4_region_joint_A4']['p'])} / {_fp(v['H4_region_joint']['p'])}", _verdict(b["H4_region_joint_A4"]["p"], v["H4_region_joint_A4"]["p"])))
        elif what == "H1":
            rows.append((label, "H1 nonlinearity Wald $p$", _fp(b["H1_nonlinear_joint"]["p"]), _fp(v["H1_nonlinear_joint"]["p"]), _verdict(b["H1_nonlinear_joint"]["p"], v["H1_nonlinear_joint"]["p"])))
    rows.append(("\\textit{Hedonic (Model A2 splines, make--model FE, dealer-clustered)}",))
    hed_row("+ sampled portals (autotrader, cargurus)", "plus_sampled_portals", "age"); hed_row("+ sampled portals", "plus_sampled_portals", "H4")
    hed_row("+ Kijiji private sellers", "plus_kijiji", "age"); hed_row("+ Kijiji private sellers", "plus_kijiji", "H4")
    hed_row("Price filter 1\\,000--500\\,000\\,\\$", "price_1000_500000", "age"); hed_row("Price filter 1\\,000--500\\,000\\,\\$", "price_1000_500000", "H4")
    hed_row("Reweighted to otogo make--model--year mix", "reweighted_to_otogo_mmy", "age"); hed_row("Reweighted to otogo make--model--year mix", "reweighted_to_otogo_mmy", "H4")
    het = H["heterogeneity"]
    def het_row(label, k1, k2, what):
        a, c = het.get(k1), het.get(k2)
        if a and c:
            l1, l2 = _esc(k1.split("=", 1)[1]).replace("<=", "$\\leq$"), _esc(k2.split("=", 1)[1]).replace("<=", "$\\leq$")
            rows.append((label, f"{what} slope: {l1} vs {l2}", _fe(a[what]["b"], a[what]["se"], 4), _fe(c[what]["b"], c[what]["se"], 4), _delta(a[what]["b"], c[what]["b"], 4)))
    het_row("Heterogeneity: powertrain", "powertrain=ICE", "powertrain=BEV", "age")
    het_row("Heterogeneity: model-year band", "yband=<=2015", "yband=2021+", "age")
    het_row("Heterogeneity: body", "body=SUV", "body=Pickup", "km")
    het_row("Heterogeneity: source class", "source_class=direct_dealer", "source_class=full_portal", "age")
    n_hed = len(rows)
    # hazards
    rows.append(("\\textit{Hazards (cloglog, dealer-clustered; $\\leq$800k listing-day subsample)}",))
    bB, bC = Z["B"]["baseline_sub"], Z["C"]["baseline_sub"]
    def haz_row(label, model, key, term):
        base = bB if model == "B" else bC; v = Z[model][key]
        if term == "age":
            rows.append((label, f"{model}: listing-age bins joint $p$", _fp(base["joint_age_bins"]["p"]), _fp(v["joint_age_bins"]["p"]), _verdict(base["joint_age_bins"]["p"], v["joint_age_bins"]["p"])))
        else:
            rows.append((label, f"{model}: coef.\\ {_esc(term)}", _fe(base["coef"][term]["b"], base["coef"][term]["se"]), _fe(v["coef"][term]["b"], v["coef"][term]["se"]), _delta(base["coef"][term]["b"], v["coef"][term]["b"])))
    # (rows comparing the pre-E-016 survivor-only Model B with the full risk set were removed once the
    #  pipeline bug was fixed: the two estimates are now identical by construction)
    haz_row("Markdown $\\geq$2\\,\\%/200\\,\\$", "B", "threshold_2pct_200", "ldens")
    haz_row("Markdown $\\geq$0.5\\,\\%/50\\,\\$", "B", "threshold_0.5pct_50", "ldens")
    haz_row("Transient rule 24\\,h", "B", "transient_24h", "ldens"); haz_row("Transient rule 72\\,h", "B", "transient_72h", "ldens")
    haz_row("Exit confirmation 3\\,d", "C", "exit_confirm_3d", "resid_A2")
    haz_row("Exit confirmation 14\\,d", "C", "exit_confirm_14d", "resid_A2")
    haz_row("S2L (left-censored stock included)", "C", "S2L_with_stock", "resid_A2"); haz_row("S2L (left-censored stock included)", "C", "S2L_with_stock", "ldens")
    # AMLI
    rows.append(("\\textit{AMLI (leave-exit-out index $\\to$ next-week exit intensity, region FE)}",))
    bA = A["baseline"]["validation_noexit_fe"]
    def am_row(label, v):
        rows.append((label, "slope $\\hat\\beta$ (s.e.)", _fe(bA["b"], bA["se"], 4), _fe(v["b"], v["se"], 4), _delta(bA["b"], v["b"], 4)))
    am_row("Rank standardisation", A["rank_standardisation"]["validation_noexit_fe"])
    am_row("Min.\\ 100 listing-days per cell", A["min_listing_days_100"]["validation_noexit_fe"])
    am_row("Two-week aggregation", A["two_week_aggregation"]["validation_noexit_fe"])
    loo = A["leave_one_component_out"]; worst = min(loo, key=lambda c: loo[c]["baseline_z"]["b"]); best = max(loo, key=lambda c: loo[c]["baseline_z"]["b"])
    am_row(f"Leave-one-component-out (min: drop {_esc(worst)})", loo[worst]["baseline_z"]); am_row(f"Leave-one-component-out (max: drop {_esc(best)})", loo[best]["baseline_z"])
    pv = A["province_vs_region"]
    rows.append(("Province vs region aggregation", "corr.\\ equal-weight index (weekly)", "1.000", f"{pv['corr_province_equal_vs_region_weighted_equal']:.3f}", "--"))
    # Holm
    rows.append(("\\textit{Multiple testing (six primary tests)}",))
    hm = holm
    rows.append(("Holm correction", "tests rejected at 5\\,\\% (raw $\\to$ Holm)", str(int(hm["reject_0.05_raw"].sum())), str(int(hm["reject_0.05_holm"].sum())), f"{int(hm['reject_0.05_raw'].sum() - hm['reject_0.05_holm'].sum())} lost"))
    assert sum(1 for r in rows if len(r) == 5) <= 30, "table exceeds 30 data rows"
    lines = ["\\begin{tabular}{llrrr}", "\\toprule", "Check & Quantity & Baseline & Variant & Change \\\\", "\\midrule"]
    for i, r in enumerate(rows):
        if len(r) == 1:
            if i > 0:
                lines.append("\\addlinespace")
            lines.append(f"\\multicolumn{{5}}{{l}}{{{r[0]}}} \\\\")
        else:
            lines.append(" & ".join(r) + " \\\\")
    lines += ["\\bottomrule", "\\end{tabular}"]
    C.TABLES.mkdir(parents=True, exist_ok=True)
    (C.TABLES / "tab09_robustness.tex").write_text("\n".join(lines) + "\n")


# =============================================================================== driver

def run() -> dict:
    OUT.mkdir(parents=True, exist_ok=True)
    t0 = time.time(); S = {"generated": time.strftime("%Y-%m-%d %H:%M:%S"), "seed": SEED}
    _log("=== 1. hedonic"); t1 = time.time(); S["hedonic"] = run_hedonic(); S["hedonic"]["seconds"] = round(time.time() - t1, 1)
    _log("=== 2. hazards"); t1 = time.time(); S["hazards"] = run_hazards(); S["hazards"]["seconds"] = round(time.time() - t1, 1)
    _log("=== 3. AMLI"); t1 = time.time(); S["amli"] = run_amli(); S["amli"]["seconds"] = round(time.time() - t1, 1)
    _log("=== 4. Holm"); S["holm"] = holm_primary_tests().to_dict(orient="records")
    S["seconds_total"] = round(time.time() - t0, 1)
    json.dump(S, open(OUT / "summary.json", "w"), indent=2, default=float)
    _log("=== 5. table"); write_table(S)
    _log(f"done in {S['seconds_total']} s")
    return S


if __name__ == "__main__":
    run()