TeX 52.9%
Python 45.7%
Shell 1.4%
1# -*- coding: utf-8 -*-2"""07 — Analyses complémentaires de la révision scientifique (editorial board, 2026-10-17).34Tourne sur le poste de travail à partir des sorties déjà rapatriées dans out/<version>/ (aucun accès à la base5Vrai Prix ni au cluster ; < 5 min) :6 targets.parquet, models_est.parquet (estimations par cible des modèles nommés M0–M13),7 tables/surface_E1.csv, tables/facteurs_E2.csv, index_full.csv.89Produit results/revision-<date>/*.csv|json et les fragments LaTeX correspondants dans paper/tables/ :10 (a) stabilité validation → test du classement des 2 256 configurations (corrélation de rang, rang en test de la11 configuration choisie en 2024, ensemble des configurations quasi optimales en validation et leur erreur en test) ;12 (b) bootstrap par grappes à 999 réplications des différences appariées (MdAPE et MAPE), statistique t robuste aux13 grappes de la différence de MAPE (type Diebold–Mariano), sensibilité au choix des grappes ;14 (c) intervalles de confiance bootstrap de la MdAPE de chaque modèle nommé ;15 (d) sensibilité aux bornes de la règle R8 (ratio prix/rôle) côté cibles ;16 (e) biais médian par trimestre civil (saisonnalité) pour les principaux modèles ;17 (f) croissance régionale de l'indice descriptif (vérification des chiffres du texte).18"""19import json, os, sys, time20from pathlib import Path21import numpy as np, pandas as pd22from scipy import stats2324ROOT = Path(__file__).resolve().parents[1]25V = os.environ.get("CV_VERSION", "2026-10-10")26OUT = ROOT / "out" / V; TAB = OUT / "tables"27RES = ROOT / "results" / "revision-2026-10-17"; RES.mkdir(parents=True, exist_ok=True)28PT = ROOT / "paper" / "tables"29B = int(os.environ.get("CV_BOOT", 999)); rng = np.random.default_rng(20261017)30t0 = time.time()31def log(m): print(f"[{time.time() - t0:6.1f}s] {m}", flush=True)32def n_(x): return f"{int(round(x)):,}"33def esc(s): return str(s).replace("%", "\\%").replace("&", "\\&")34def tabular(rows, cols, align):35 L = ["\\begin{tabular}{" + align + "}", "\\toprule", " & ".join(cols) + " \\\\", "\\midrule"]36 L += [" & ".join(str(c) for c in r) + " \\\\" for r in rows]37 L += ["\\bottomrule", "\\end{tabular}"]38 return "\n".join(L) + "\n"3940T = pd.read_parquet(OUT / "targets.parquet"); M = pd.read_parquet(OUT / "models_est.parquet")41assert (T["idx"].to_numpy() == M["idx"].to_numpy()).all()42TEST = (T["period"] == "test").to_numpy(); VAL = (T["period"] == "val").to_numpy()43lnP = T["ln_price"].to_numpy(float)44MODELS = ["M1", "M2", "M3", "M4", "M5", "M6", "M7", "M8a", "M8b", "M13", "M10", "M11a", "M11b", "M0", "M9", "M12"]45ERR = {m: np.exp(M[m].to_numpy(float) - lnP) - 1 for m in MODELS}46T["quarter"] = pd.to_datetime(T["date"]).dt.to_period("Q").astype(str)47T["cal_q"] = pd.to_datetime(T["date"]).dt.quarter48res = {}4950# ---------------------------------------------------------------- (a) stabilité validation → test des 2 256 configurations51S = pd.read_csv(TAB / "surface_E1.csv").assign(src="E1"); F = pd.read_csv(TAB / "facteurs_E2.csv").assign(src="E2")52A = pd.concat([S, F], ignore_index=True)53ok = A["couverture_val"] >= 9754def spear(d): return float(stats.spearmanr(d["mdape_val"], d["mdape_test"]).statistic)55stab = {"n_configs": int(len(A)), "n_configs_cov97": int(ok.sum()),56 "spearman_E1": spear(S), "spearman_E2": spear(F), "spearman_all": spear(A), "spearman_cov97": spear(A[ok])}57A97 = A[ok].copy(); A97["rank_val"] = A97["mdape_val"].rank(method="min"); A97["rank_test"] = A97["mdape_test"].rank(method="min")58chosen = A97.sort_values("mdape_val").iloc[0]; oracle = A97.sort_values("mdape_test").iloc[0]59stab.update(chosen_val=float(chosen.mdape_val), chosen_test=float(chosen.mdape_test), chosen_rank_test=int(chosen.rank_test),60 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}",61 regret_test=float(chosen.mdape_test - oracle.mdape_test))62near = []63for tol in [0.1, 0.25, 0.5, 1.0]:64 g = A97[A97["mdape_val"] <= chosen.mdape_val + tol]65 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())))66stab["near_optimal"] = near67top = A97.sort_values("mdape_val").head(10)68stab["top10_val"] = top[["src", "sel", "k", "d", "T", "adj", "w", "recon", "mdape_val", "mdape_test", "couverture_test", "rank_test"]].to_dict("records")69res["selection_stability"] = stab70LAB = {"geo": "nearest", "sim": "most similar", "mix": "compromise", "adj": "least adjusted", "none": "time only", "hedo": "hedonic grid",71 "role": "assessment ratio", "gbm": "GBM adjustments", "conv": "conventional grid", "hedo_reg": "regional grid", "local": "local grid", "raw": "none",72 "equal": "equal", "dist": "1/distance", "rec": "recency", "all": "product", "median": "median", "mean": "mean"}73rows = [[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()]74(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"))75rows = [[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]76(PT / "selection_near.tex").write_text(tabular(rows, ["Validation MdAPE within (pts of best)", "Configurations", "Test MdAPE min", "median", "max"], "lrrrr"))77log(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}")7879# ---------------------------------------------------------------- bootstrap par grappes : outils80def cluster_ids(key):81 return pd.factorize(key)[0]82CLUSTERS = {"municipality x quarter": cluster_ids(T["code_mun"].astype(str) + "|" + T["quarter"]),83 "MRC x quarter": cluster_ids(T["mrc"].astype(str) + "|" + T["quarter"]),84 "municipality": cluster_ids(T["code_mun"].astype(str)),85 "region x quarter": cluster_ids(T["region"].astype(str) + "|" + T["quarter"])}86def boot_weights(cl, mask, B):87 """Multiplicités bootstrap (grappes tirées avec remise) pour les observations de `mask` ; matrice B x n."""88 c = cl[mask]; ids, inv = np.unique(c, return_inverse=True); nc = len(ids)89 W = np.empty((B, mask.sum()), np.int32)90 for b in range(B):91 pick = rng.integers(0, nc, nc); mult = np.bincount(pick, minlength=nc); W[b] = mult[inv]92 return W93def wmedian_sorted(v_sorted, w_sorted_rows):94 """Médiane pondérée pour chaque ligne de poids (observations déjà triées par v)."""95 cw = np.cumsum(w_sorted_rows, axis=1); half = cw[:, -1] / 2.096 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)]97def ci(x): return [float(np.percentile(x, 2.5)), float(np.percentile(x, 97.5))]9899# ---------------------------------------------------------------- (c) IC bootstrap de la MdAPE de chaque modèle (grappes municipalité x trimestre)100Wmain = boot_weights(CLUSTERS["municipality x quarter"], TEST, B)101log(f"poids bootstrap principaux : {Wmain.shape}")102mod_ci = {}103for m in MODELS:104 a = np.abs(ERR[m][TEST]); okm = ~np.isnan(a)105 v = a[okm]; o = np.argsort(v); vs = v[o]; W = Wmain[:, okm][:, o]106 med = wmedian_sorted(vs, W) * 100; mean = (W @ vs) / W.sum(axis=1) * 100107 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()))108res["model_ci"] = mod_ci109order = sorted(MODELS, key=lambda m: mod_ci[m]["mdape"])110rows = [[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]111(PT / "competition_ci.tex").write_text(tabular(rows, ["Model", "MdAPE", "95\\% CI", "MAPE", "95\\% CI", "$n$"], "lrrrrr"))112log("(c) IC par modèle écrits")113114# ---------------------------------------------------------------- (b) différences appariées : B = 999, t robuste aux grappes, autres grappes115PAIRS = [("M7", "M2"), ("M7", "M6"), ("M8a", "M7"), ("M8b", "M7"), ("M11b", "M0"), ("M11a", "M0"), ("M12", "M11b"), ("M13", "M7"), ("M3", "M2"), ("M9", "M7"), ("M1", "M2"),116 ("M0", "M13"), ("M0", "M9"), ("M0", "M7"), ("M13", "M11b")]117def paired(a_, b_, mask, cl_name, B_):118 cl = CLUSTERS[cl_name]119 a = np.abs(a_[mask]); b = np.abs(b_[mask]); okp = ~np.isnan(a) & ~np.isnan(b)120 a, b = a[okp], b[okp]; d = a - b121 sub = np.zeros(mask.sum(), bool); sub[okp] = True122 W = Wmain[:, okp] if (cl_name == "municipality x quarter" and B_ == B) else boot_weights(cl, _mask_sub(mask, okp), B_)123 oa = np.argsort(a); ob = np.argsort(b)124 med_d = (wmedian_sorted(a[oa], W[:, oa]) - wmedian_sorted(b[ob], W[:, ob])) * 100125 mean_d = (W @ d) / W.sum(axis=1) * 100126 # statistique t robuste aux grappes (Diebold–Mariano sur la perte absolue en pourcentage)127 c = cl[_mask_sub(mask, okp)]; ids, inv = np.unique(c, return_inverse=True)128 dbar = d.mean(); S_c = np.bincount(inv, d - dbar); se = np.sqrt((S_c ** 2).sum()) / len(d)129 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),130 diff_mape=float(dbar * 100), ci_mape=ci(mean_d), t_cluster=float(dbar / se), share_A_better=float((a < b).mean() * 100))131def _mask_sub(mask, okp):132 m = np.zeros(len(mask), bool); m[np.flatnonzero(mask)[okp]] = True; return m133pairs_main = {f"{a}-{b}": paired(ERR[a], ERR[b], TEST, "municipality x quarter", B) for a, b in PAIRS}134res["pairs_B999"] = pairs_main135rows = [[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"])]136 for (a, b), p in zip(PAIRS, pairs_main.values())]137(PT / "pairs_b999.tex").write_text(tabular(rows, ["Comparison", "$\\Delta$ MdAPE [95\\% CI]", "$\\Delta$ MAPE [95\\% CI]", "$t$", "$A$ better (\\%)", "$n$", "Clusters"], "lrrrrrr"))138log("(b) paires B=999 écrites")139# sensibilité au choix des grappes (B = 399) pour trois comparaisons clés140SENS = [("M7", "M2"), ("M11b", "M0"), ("M12", "M11b"), ("M7", "M6")]141sens = {}142for cl_name in CLUSTERS:143 for a, b in SENS:144 sens[f"{cl_name}|{a}-{b}"] = paired(ERR[a], ERR[b], TEST, cl_name, 399)145res["pairs_cluster_sensitivity"] = sens146rows = []147for cl_name in CLUSTERS:148 r = [cl_name.replace(" x ", " $\\times$ "), n_(sens[f"{cl_name}|M7-M2"]["n_clusters"])]149 for a, b in SENS:150 p = sens[f"{cl_name}|{a}-{b}"]; r.append(f"[{p['ci_mdape'][0]:.2f}; {p['ci_mdape'][1]:.2f}]")151 rows.append(r)152(PT / "pairs_cluster_sens.tex").write_text(tabular(rows, ["Cluster definition", "Clusters"] + [f"{a} $-$ {b}" for a, b in SENS], "lr" + "r" * len(SENS)))153log("(b') sensibilité aux grappes écrite")154155# ---------------------------------------------------------------- (d) sensibilité aux bornes de R8 (cibles)156ratio = (T["amount"] / T["role_vente"]).to_numpy()157KEY = ["M2", "M7", "M13", "M0", "M10", "M11a", "M11b", "M12"]158r8 = {}159for 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)]:160 m = TEST & (ratio >= lo) & (ratio <= hi)161 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})162res["r8_sensitivity"] = r8163rows = [[lab, n_(d["n"]), f"{d['share']:.1f}"] + [f"{d[k]:.1f}" for k in KEY] for lab, d in r8.items()]164(PT / "r8_sensitivity.tex").write_text(tabular(rows, ["Target price/assessment ratio", "$n$", "Share (\\%)"] + KEY, "lrr" + "r" * len(KEY)))165log("(d) sensibilité R8 écrite")166167# ---------------------------------------------------------------- (e) biais médian par trimestre civil168seas = {}169for q in [1, 2, 3, 4]:170 m = TEST & (T["cal_q"].to_numpy() == q)171 seas[f"Q{q}"] = dict(n=int(m.sum()), **{k: float(np.nanmedian(ERR[k][m]) * 100) for k in KEY})172res["seasonal_bias"] = seas173rows = [[q, n_(d["n"])] + [f"{d[k]:+.1f}" for k in KEY] for q, d in seas.items()]174(PT / "seasonal_bias.tex").write_text(tabular(rows, ["Calendar quarter", "$n$"] + KEY, "lr" + "r" * len(KEY)))175log("(e) biais saisonnier écrit")176177# ---------------------------------------------------------------- (f) indice descriptif : croissance régionale178I = pd.read_csv(OUT / "index_full.csv", dtype={"region": str})179NAMES = {"QC": "Quebec (province)", "06": "Montréal", "03": "Capitale-Nationale", "16": "Montérégie", "13": "Laval", "07": "Outaouais", "15": "Laurentides", "14": "Lanaudière", "05": "Estrie"}180idx = {}181for r, nm in NAMES.items():182 g = I[I.region == r].set_index("ym")["log_idx"]183 if len(g):184 peak = g.loc["2021-01":"2022-12"].idxmax(); trough = g.loc[peak:"2023-12"].idxmin()185 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),186 peak_2022=peak, trough_2023=trough, drawdown=float((np.exp(g.loc[trough] - g.loc[peak]) - 1) * 100))187res["index_growth"] = idx188rows = [[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()]189(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"))190log("(f) indice écrit")191192with 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)))193pd.DataFrame(top).to_csv(RES / "selection_top10.csv", index=False)194pd.DataFrame([dict(pair=k, **v) for k, v in pairs_main.items()]).to_csv(RES / "pairs_B999.csv", index=False)195pd.DataFrame(mod_ci).T.to_csv(RES / "model_ci.csv")196pd.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")197# macros pour le texte198mac = {"spearmanAll": f"{stab['spearman_all']:.3f}", "spearmanCov": f"{stab['spearman_cov97']:.3f}", "mZeroRankTest": str(stab["chosen_rank_test"]), "nConfigsCov": n_(stab["n_configs_cov97"]),199 "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}\\%",200 "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}",201 "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}",202 "tMzeroMthirteen": f"{pairs_main['M0-M13']['t_cluster']:.1f}", "dMzeroMthirteen": f"{pairs_main['M0-M13']['diff_mdape']:.2f}",203 "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}\\%",204 "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}\\%",205 "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}\\%",206 "seasQoneMseven": f"{seas['Q1']['M7']:+.1f}\\%", "seasQthreeMseven": f"{seas['Q3']['M7']:+.1f}\\%", "seasQoneMelevenB": f"{seas['Q1']['M11b']:+.1f}\\%", "seasQthreeMelevenB": f"{seas['Q3']['M11b']:+.1f}\\%",207 "seasQoneMzero": f"{seas['Q1']['M0']:+.1f}\\%", "seasQthreeMzero": f"{seas['Q3']['M0']:+.1f}\\%", "seasQoneMten": f"{seas['Q1']['M10']:+.1f}\\%", "seasQthreeMten": f"{seas['Q3']['M10']:+.1f}\\%"}208(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")209log("terminé")210