# -*- coding: utf-8 -*-
"""07 — Analyses complémentaires de la révision scientifique (editorial board, 2026-10-17).

Tourne sur le poste de travail à partir des sorties déjà rapatriées dans out/<version>/ (aucun accès à la base
Vrai Prix ni au cluster ; < 5 min) :
  targets.parquet, models_est.parquet (estimations par cible des modèles nommés M0–M13),
  tables/surface_E1.csv, tables/facteurs_E2.csv, index_full.csv.

Produit results/revision-<date>/*.csv|json et les fragments LaTeX correspondants dans paper/tables/ :
  (a) stabilité validation → test du classement des 2 256 configurations (corrélation de rang, rang en test de la
      configuration choisie en 2024, ensemble des configurations quasi optimales en validation et leur erreur en test) ;
  (b) bootstrap par grappes à 999 réplications des différences appariées (MdAPE et MAPE), statistique t robuste aux
      grappes de la différence de MAPE (type Diebold–Mariano), sensibilité au choix des grappes ;
  (c) intervalles de confiance bootstrap de la MdAPE de chaque modèle nommé ;
  (d) sensibilité aux bornes de la règle R8 (ratio prix/rôle) côté cibles ;
  (e) biais médian par trimestre civil (saisonnalité) pour les principaux modèles ;
  (f) croissance régionale de l'indice descriptif (vérification des chiffres du texte).
"""
import json, os, sys, time
from pathlib import Path
import numpy as np, pandas as pd
from scipy import stats

ROOT = Path(__file__).resolve().parents[1]
V = os.environ.get("CV_VERSION", "2026-10-10")
OUT = ROOT / "out" / V; TAB = OUT / "tables"
RES = ROOT / "results" / "revision-2026-10-17"; RES.mkdir(parents=True, exist_ok=True)
PT = ROOT / "paper" / "tables"
B = int(os.environ.get("CV_BOOT", 999)); rng = np.random.default_rng(20261017)
t0 = time.time()
def log(m): print(f"[{time.time() - t0:6.1f}s] {m}", flush=True)
def n_(x): return f"{int(round(x)):,}"
def esc(s): return str(s).replace("%", "\\%").replace("&", "\\&")
def tabular(rows, cols, align):
    L = ["\\begin{tabular}{" + align + "}", "\\toprule", " & ".join(cols) + " \\\\", "\\midrule"]
    L += [" & ".join(str(c) for c in r) + " \\\\" for r in rows]
    L += ["\\bottomrule", "\\end{tabular}"]
    return "\n".join(L) + "\n"

T = pd.read_parquet(OUT / "targets.parquet"); M = pd.read_parquet(OUT / "models_est.parquet")
assert (T["idx"].to_numpy() == M["idx"].to_numpy()).all()
TEST = (T["period"] == "test").to_numpy(); VAL = (T["period"] == "val").to_numpy()
lnP = T["ln_price"].to_numpy(float)
MODELS = ["M1", "M2", "M3", "M4", "M5", "M6", "M7", "M8a", "M8b", "M13", "M10", "M11a", "M11b", "M0", "M9", "M12"]
ERR = {m: np.exp(M[m].to_numpy(float) - lnP) - 1 for m in MODELS}
T["quarter"] = pd.to_datetime(T["date"]).dt.to_period("Q").astype(str)
T["cal_q"] = pd.to_datetime(T["date"]).dt.quarter
res = {}

# ---------------------------------------------------------------- (a) stabilité validation → test des 2 256 configurations
S = pd.read_csv(TAB / "surface_E1.csv").assign(src="E1"); F = pd.read_csv(TAB / "facteurs_E2.csv").assign(src="E2")
A = pd.concat([S, F], ignore_index=True)
ok = A["couverture_val"] >= 97
def spear(d): return float(stats.spearmanr(d["mdape_val"], d["mdape_test"]).statistic)
stab = {"n_configs": int(len(A)), "n_configs_cov97": int(ok.sum()),
        "spearman_E1": spear(S), "spearman_E2": spear(F), "spearman_all": spear(A), "spearman_cov97": spear(A[ok])}
A97 = A[ok].copy(); A97["rank_val"] = A97["mdape_val"].rank(method="min"); A97["rank_test"] = A97["mdape_test"].rank(method="min")
chosen = A97.sort_values("mdape_val").iloc[0]; oracle = A97.sort_values("mdape_test").iloc[0]
stab.update(chosen_val=float(chosen.mdape_val), chosen_test=float(chosen.mdape_test), chosen_rank_test=int(chosen.rank_test),
            oracle_test=float(oracle.mdape_test), oracle_desc=f"{oracle.sel}|k={int(oracle.k)}|d={oracle.d}|T={int(oracle['T'])}|{oracle.adj}|{oracle.w}|{oracle.recon}",
            regret_test=float(chosen.mdape_test - oracle.mdape_test))
near = []
for tol in [0.1, 0.25, 0.5, 1.0]:
    g = A97[A97["mdape_val"] <= chosen.mdape_val + tol]
    near.append(dict(tol=tol, n=int(len(g)), test_min=float(g.mdape_test.min()), test_max=float(g.mdape_test.max()), test_median=float(g.mdape_test.median())))
stab["near_optimal"] = near
top = A97.sort_values("mdape_val").head(10)
stab["top10_val"] = top[["src", "sel", "k", "d", "T", "adj", "w", "recon", "mdape_val", "mdape_test", "couverture_test", "rank_test"]].to_dict("records")
res["selection_stability"] = stab
LAB = {"geo": "nearest", "sim": "most similar", "mix": "compromise", "adj": "least adjusted", "none": "time only", "hedo": "hedonic grid",
       "role": "assessment ratio", "gbm": "GBM adjustments", "conv": "conventional grid", "hedo_reg": "regional grid", "local": "local grid", "raw": "none",
       "equal": "equal", "dist": "1/distance", "rec": "recency", "all": "product", "median": "median", "mean": "mean"}
rows = [[LAB[r.sel], int(r.k), f"{r.d:g}", int(r["T"]), LAB[r.adj], f"{LAB[r.w]}/{LAB[r.recon]}", f"{r.mdape_val:.2f}", f"{r.mdape_test:.2f}", int(r.rank_test)] for _, r in top.iterrows()]
(PT / "selection_top10.tex").write_text(tabular(rows, ["Selection", "$k$", "$d$ (km)", "$t$ (mo)", "Adjustment", "Weights/recon.", "MdAPE 2024", "MdAPE 2025--26", "Rank 2025--26"], "lrrrllrrr"))
rows = [[f"$\\le {d['tol']:.2f}$", d["n"], f"{d['test_min']:.2f}", f"{d['test_median']:.2f}", f"{d['test_max']:.2f}"] for d in near]
(PT / "selection_near.tex").write_text(tabular(rows, ["Validation MdAPE within (pts of best)", "Configurations", "Test MdAPE min", "median", "max"], "lrrrr"))
log(f"(a) stabilité : Spearman all {stab['spearman_all']:.4f}, rang test de M0 {stab['chosen_rank_test']}/{stab['n_configs_cov97']}, regret {stab['regret_test']:.2f}")

# ---------------------------------------------------------------- bootstrap par grappes : outils
def cluster_ids(key):
    return pd.factorize(key)[0]
CLUSTERS = {"municipality x quarter": cluster_ids(T["code_mun"].astype(str) + "|" + T["quarter"]),
            "MRC x quarter": cluster_ids(T["mrc"].astype(str) + "|" + T["quarter"]),
            "municipality": cluster_ids(T["code_mun"].astype(str)),
            "region x quarter": cluster_ids(T["region"].astype(str) + "|" + T["quarter"])}
def boot_weights(cl, mask, B):
    """Multiplicités bootstrap (grappes tirées avec remise) pour les observations de `mask` ; matrice B x n."""
    c = cl[mask]; ids, inv = np.unique(c, return_inverse=True); nc = len(ids)
    W = np.empty((B, mask.sum()), np.int32)
    for b in range(B):
        pick = rng.integers(0, nc, nc); mult = np.bincount(pick, minlength=nc); W[b] = mult[inv]
    return W
def wmedian_sorted(v_sorted, w_sorted_rows):
    """Médiane pondérée pour chaque ligne de poids (observations déjà triées par v)."""
    cw = np.cumsum(w_sorted_rows, axis=1); half = cw[:, -1] / 2.0
    pos = np.array([np.searchsorted(cw[b], half[b]) for b in range(cw.shape[0])]); return v_sorted[np.minimum(pos, len(v_sorted) - 1)]
def ci(x): return [float(np.percentile(x, 2.5)), float(np.percentile(x, 97.5))]

# ---------------------------------------------------------------- (c) IC bootstrap de la MdAPE de chaque modèle (grappes municipalité x trimestre)
Wmain = boot_weights(CLUSTERS["municipality x quarter"], TEST, B)
log(f"poids bootstrap principaux : {Wmain.shape}")
mod_ci = {}
for m in MODELS:
    a = np.abs(ERR[m][TEST]); okm = ~np.isnan(a)
    v = a[okm]; o = np.argsort(v); vs = v[o]; W = Wmain[:, okm][:, o]
    med = wmedian_sorted(vs, W) * 100; mean = (W @ vs) / W.sum(axis=1) * 100
    mod_ci[m] = dict(mdape=float(np.median(v) * 100), ci_mdape=ci(med), mape=float(v.mean() * 100), ci_mape=ci(mean), n=int(okm.sum()))
res["model_ci"] = mod_ci
order = sorted(MODELS, key=lambda m: mod_ci[m]["mdape"])
rows = [[m, f"{mod_ci[m]['mdape']:.2f}", f"[{mod_ci[m]['ci_mdape'][0]:.2f}; {mod_ci[m]['ci_mdape'][1]:.2f}]", f"{mod_ci[m]['mape']:.2f}", f"[{mod_ci[m]['ci_mape'][0]:.2f}; {mod_ci[m]['ci_mape'][1]:.2f}]", n_(mod_ci[m]["n"])] for m in order]
(PT / "competition_ci.tex").write_text(tabular(rows, ["Model", "MdAPE", "95\\% CI", "MAPE", "95\\% CI", "$n$"], "lrrrrr"))
log("(c) IC par modèle écrits")

# ---------------------------------------------------------------- (b) différences appariées : B = 999, t robuste aux grappes, autres grappes
PAIRS = [("M7", "M2"), ("M7", "M6"), ("M8a", "M7"), ("M8b", "M7"), ("M11b", "M0"), ("M11a", "M0"), ("M12", "M11b"), ("M13", "M7"), ("M3", "M2"), ("M9", "M7"), ("M1", "M2"),
         ("M0", "M13"), ("M0", "M9"), ("M0", "M7"), ("M13", "M11b")]
def paired(a_, b_, mask, cl_name, B_):
    cl = CLUSTERS[cl_name]
    a = np.abs(a_[mask]); b = np.abs(b_[mask]); okp = ~np.isnan(a) & ~np.isnan(b)
    a, b = a[okp], b[okp]; d = a - b
    sub = np.zeros(mask.sum(), bool); sub[okp] = True
    W = Wmain[:, okp] if (cl_name == "municipality x quarter" and B_ == B) else boot_weights(cl, _mask_sub(mask, okp), B_)
    oa = np.argsort(a); ob = np.argsort(b)
    med_d = (wmedian_sorted(a[oa], W[:, oa]) - wmedian_sorted(b[ob], W[:, ob])) * 100
    mean_d = (W @ d) / W.sum(axis=1) * 100
    # statistique t robuste aux grappes (Diebold–Mariano sur la perte absolue en pourcentage)
    c = cl[_mask_sub(mask, okp)]; ids, inv = np.unique(c, return_inverse=True)
    dbar = d.mean(); S_c = np.bincount(inv, d - dbar); se = np.sqrt((S_c ** 2).sum()) / len(d)
    return dict(n=int(len(d)), n_clusters=int(len(ids)), diff_mdape=float((np.median(a) - np.median(b)) * 100), ci_mdape=ci(med_d),
                diff_mape=float(dbar * 100), ci_mape=ci(mean_d), t_cluster=float(dbar / se), share_A_better=float((a < b).mean() * 100))
def _mask_sub(mask, okp):
    m = np.zeros(len(mask), bool); m[np.flatnonzero(mask)[okp]] = True; return m
pairs_main = {f"{a}-{b}": paired(ERR[a], ERR[b], TEST, "municipality x quarter", B) for a, b in PAIRS}
res["pairs_B999"] = pairs_main
rows = [[f"{a} $-$ {b}", f"{p['diff_mdape']:.2f} [{p['ci_mdape'][0]:.2f}; {p['ci_mdape'][1]:.2f}]", f"{p['diff_mape']:.2f} [{p['ci_mape'][0]:.2f}; {p['ci_mape'][1]:.2f}]", f"{p['t_cluster']:.1f}", f"{p['share_A_better']:.1f}", n_(p["n"]), n_(p["n_clusters"])]
        for (a, b), p in zip(PAIRS, pairs_main.values())]
(PT / "pairs_b999.tex").write_text(tabular(rows, ["Comparison", "$\\Delta$ MdAPE [95\\% CI]", "$\\Delta$ MAPE [95\\% CI]", "$t$", "$A$ better (\\%)", "$n$", "Clusters"], "lrrrrrr"))
log("(b) paires B=999 écrites")
# sensibilité au choix des grappes (B = 399) pour trois comparaisons clés
SENS = [("M7", "M2"), ("M11b", "M0"), ("M12", "M11b"), ("M7", "M6")]
sens = {}
for cl_name in CLUSTERS:
    for a, b in SENS:
        sens[f"{cl_name}|{a}-{b}"] = paired(ERR[a], ERR[b], TEST, cl_name, 399)
res["pairs_cluster_sensitivity"] = sens
rows = []
for cl_name in CLUSTERS:
    r = [cl_name.replace(" x ", " $\\times$ "), n_(sens[f"{cl_name}|M7-M2"]["n_clusters"])]
    for a, b in SENS:
        p = sens[f"{cl_name}|{a}-{b}"]; r.append(f"[{p['ci_mdape'][0]:.2f}; {p['ci_mdape'][1]:.2f}]")
    rows.append(r)
(PT / "pairs_cluster_sens.tex").write_text(tabular(rows, ["Cluster definition", "Clusters"] + [f"{a} $-$ {b}" for a, b in SENS], "lr" + "r" * len(SENS)))
log("(b') sensibilité aux grappes écrite")

# ---------------------------------------------------------------- (d) sensibilité aux bornes de R8 (cibles)
ratio = (T["amount"] / T["role_vente"]).to_numpy()
KEY = ["M2", "M7", "M13", "M0", "M10", "M11a", "M11b", "M12"]
r8 = {}
for lab, lo, hi in [("$[0.33, 3]$ (baseline)", 0.33, 3.0), ("$[0.40, 2.5]$", 0.40, 2.5), ("$[0.50, 2.0]$", 0.50, 2.0), ("$[0.67, 1.5]$", 0.67, 1.5)]:
    m = TEST & (ratio >= lo) & (ratio <= hi)
    r8[lab] = dict(n=int(m.sum()), share=float(m.sum() / TEST.sum() * 100), **{k: float(np.nanmedian(np.abs(ERR[k][m])) * 100) for k in KEY})
res["r8_sensitivity"] = r8
rows = [[lab, n_(d["n"]), f"{d['share']:.1f}"] + [f"{d[k]:.1f}" for k in KEY] for lab, d in r8.items()]
(PT / "r8_sensitivity.tex").write_text(tabular(rows, ["Target price/assessment ratio", "$n$", "Share (\\%)"] + KEY, "lrr" + "r" * len(KEY)))
log("(d) sensibilité R8 écrite")

# ---------------------------------------------------------------- (e) biais médian par trimestre civil
seas = {}
for q in [1, 2, 3, 4]:
    m = TEST & (T["cal_q"].to_numpy() == q)
    seas[f"Q{q}"] = dict(n=int(m.sum()), **{k: float(np.nanmedian(ERR[k][m]) * 100) for k in KEY})
res["seasonal_bias"] = seas
rows = [[q, n_(d["n"])] + [f"{d[k]:+.1f}" for k in KEY] for q, d in seas.items()]
(PT / "seasonal_bias.tex").write_text(tabular(rows, ["Calendar quarter", "$n$"] + KEY, "lr" + "r" * len(KEY)))
log("(e) biais saisonnier écrit")

# ---------------------------------------------------------------- (f) indice descriptif : croissance régionale
I = pd.read_csv(OUT / "index_full.csv", dtype={"region": str})
NAMES = {"QC": "Quebec (province)", "06": "Montréal", "03": "Capitale-Nationale", "16": "Montérégie", "13": "Laval", "07": "Outaouais", "15": "Laurentides", "14": "Lanaudière", "05": "Estrie"}
idx = {}
for r, nm in NAMES.items():
    g = I[I.region == r].set_index("ym")["log_idx"]
    if len(g):
        peak = g.loc["2021-01":"2022-12"].idxmax(); trough = g.loc[peak:"2023-12"].idxmin()
        idx[nm] = dict(g_2021_2026=float((np.exp(g.loc["2026-09"] - g.loc["2021-01"]) - 1) * 100), g_2025_2026=float((np.exp(g.loc["2026-09"] - g.loc["2024-12"]) - 1) * 100),
                       peak_2022=peak, trough_2023=trough, drawdown=float((np.exp(g.loc[trough] - g.loc[peak]) - 1) * 100))
res["index_growth"] = idx
rows = [[nm, f"{d['g_2021_2026']:.1f}", f"{d['peak_2022']}", f"{d['trough_2023']}", f"{d['drawdown']:.1f}", f"{d['g_2025_2026']:.1f}"] for nm, d in idx.items()]
(PT / "index_growth.tex").write_text(tabular(rows, ["Region", "2021-01 to 2026-09 (\\%)", "Peak", "Trough", "Peak to trough (\\%)", "2024-12 to 2026-09 (\\%)"], "lrrrrr"))
log("(f) indice écrit")

with open(RES / "revision_results.json", "w") as f: json.dump(res, f, indent=2, default=lambda o: float(o) if isinstance(o, (np.floating,)) else (int(o) if isinstance(o, np.integer) else str(o)))
pd.DataFrame(top).to_csv(RES / "selection_top10.csv", index=False)
pd.DataFrame([dict(pair=k, **v) for k, v in pairs_main.items()]).to_csv(RES / "pairs_B999.csv", index=False)
pd.DataFrame(mod_ci).T.to_csv(RES / "model_ci.csv")
pd.DataFrame(r8).T.to_csv(RES / "r8_sensitivity.csv"); pd.DataFrame(seas).T.to_csv(RES / "seasonal_bias.csv"); pd.DataFrame(idx).T.to_csv(RES / "index_growth.csv")
# macros pour le texte
mac = {"spearmanAll": f"{stab['spearman_all']:.3f}", "spearmanCov": f"{stab['spearman_cov97']:.3f}", "mZeroRankTest": str(stab["chosen_rank_test"]), "nConfigsCov": n_(stab["n_configs_cov97"]),
       "oracleTest": f"{stab['oracle_test']:.1f}\\%", "regretTest": f"{stab['regret_test']:.1f}", "nearQuarterN": str(near[1]["n"]), "nearQuarterMin": f"{near[1]['test_min']:.1f}\\%", "nearQuarterMax": f"{near[1]['test_max']:.1f}\\%",
       "ciMzeroLo": f"{mod_ci['M0']['ci_mdape'][0]:.1f}", "ciMzeroHi": f"{mod_ci['M0']['ci_mdape'][1]:.1f}", "ciMelevenBLo": f"{mod_ci['M11b']['ci_mdape'][0]:.1f}", "ciMelevenBHi": f"{mod_ci['M11b']['ci_mdape'][1]:.1f}",
       "tMsevenMtwo": f"{pairs_main['M7-M2']['t_cluster']:.0f}", "tMelevenBMzero": f"{pairs_main['M11b-M0']['t_cluster']:.0f}", "tMtwelveMelevenB": f"{pairs_main['M12-M11b']['t_cluster']:.1f}", "tMsevenMsix": f"{pairs_main['M7-M6']['t_cluster']:.0f}",
       "tMzeroMthirteen": f"{pairs_main['M0-M13']['t_cluster']:.1f}", "dMzeroMthirteen": f"{pairs_main['M0-M13']['diff_mdape']:.2f}",
       "rEightTightShare": f"{r8['$[0.50, 2.0]$']['share']:.1f}\\%", "rEightTightMzero": f"{r8['$[0.50, 2.0]$']['M0']:.1f}\\%", "rEightTightMelevenB": f"{r8['$[0.50, 2.0]$']['M11b']:.1f}\\%",
       "idxMtl": f"{idx['Montréal']['g_2025_2026']:.0f}\\%", "idxQc": f"{idx['Capitale-Nationale']['g_2025_2026']:.0f}\\%", "idxMonteregie": f"{idx['Montérégie']['g_2025_2026']:.0f}\\%",
       "idxProvAll": f"{idx['Quebec (province)']['g_2021_2026']:.0f}\\%", "idxMtlAll": f"{idx['Montréal']['g_2021_2026']:.0f}\\%", "idxEstrieAll": f"{idx['Estrie']['g_2021_2026']:.0f}\\%", "idxDrawdown": f"{abs(idx['Quebec (province)']['drawdown']):.0f}\\%",
       "seasQoneMseven": f"{seas['Q1']['M7']:+.1f}\\%", "seasQthreeMseven": f"{seas['Q3']['M7']:+.1f}\\%", "seasQoneMelevenB": f"{seas['Q1']['M11b']:+.1f}\\%", "seasQthreeMelevenB": f"{seas['Q3']['M11b']:+.1f}\\%",
       "seasQoneMzero": f"{seas['Q1']['M0']:+.1f}\\%", "seasQthreeMzero": f"{seas['Q3']['M0']:+.1f}\\%", "seasQoneMten": f"{seas['Q1']['M10']:+.1f}\\%", "seasQthreeMten": f"{seas['Q3']['M10']:+.1f}\\%"}
(PT / "numbers_revision.tex").write_text("% généré par pipeline/07_revision_analyses.py\n" + "\n".join(f"\\newcommand{{\\{k}}}{{{v}}}" for k, v in mac.items()) + "\n")
log("terminé")