TeX 52.9%
Python 45.7%
Shell 1.4%
1# -*- coding: utf-8 -*-2"""05 — Analyse : métriques par configuration (validation 2024 → choix ; test 2025-26 → verdict),3surface Erreur(k, d, t), effets des facteurs (sélection, ajustement, pondération, réconciliation),4compétition de modèles M1-M13, hypothèses H1-H5, strates (région, densité, prix), intervalles,5règles conditionnelles. Sorties : tables/*.csv, results.json, fig/*.pdf"""6from scipy.optimize import nnls7import matplotlib; matplotlib.use("Agg")8import matplotlib.pyplot as plt9from matplotlib import rcParams10from common import *1112# ---------------------------------------------------------------- style (palette de référence, marques fines)13PAL = ["#2a78d6", "#eb6834", "#1baf7a", "#eda100", "#e87ba4", "#008300", "#6250d6", "#e34948"]14SEQ = ["#cde2fb", "#9ec5f4", "#6da7ec", "#3987e5", "#256abf", "#184f95", "#0d366b"]15rcParams.update({"font.family": "Helvetica Neue", "font.size": 9, "axes.edgecolor": "#c3c2b7", "axes.linewidth": 0.6,16 "axes.grid": True, "grid.color": "#e8e7e3", "grid.linewidth": 0.5, "axes.spines.top": False,17 "axes.spines.right": False, "legend.frameon": False, "figure.dpi": 150, "savefig.bbox": "tight",18 "axes.titleweight": "bold", "axes.titlesize": 10, "axes.labelcolor": "#52514e", "xtick.color": "#52514e",19 "ytick.color": "#52514e", "figure.facecolor": "#fcfcfb", "axes.facecolor": "#fcfcfb"})20from matplotlib.colors import LinearSegmentedColormap21CMAP = LinearSegmentedColormap.from_list("seqblue", SEQ)22TAB = OUT / "tables"; TAB.mkdir(exist_ok=True); FIG = OUT / "fig"23def savefig(fig, name): fig.savefig(FIG / f"{name}.pdf"); fig.savefig(FIG / f"{name}.png"); plt.close(fig)2425# ---------------------------------------------------------------- chargement26T = load_parquet("targets.parquet"); S = load_parquet("sales.parquet").set_index("idx")27E1 = pd.read_csv(OUT / "configs_E1.csv"); E2 = pd.read_csv(OUT / "configs_E2.csv")28est1 = np.load(OUT / "est_E1.npy"); est2 = np.load(OUT / "est_E2.npy"); sd2 = np.load(OUT / "sd_E2.npy"); nav = np.load(OUT / "n_avail.npy")29meta = load_json("index_meta.json"); regions = meta["regions"]; M_EVAL0 = meta["M_EVAL0"]; IDX = np.load(OUT / "index_asof.npy")30rid = {r: i for i, r in enumerate(regions)}31lnP = T["ln_price"].to_numpy(np.float32)32VAL = (T["period"] == "val").to_numpy(); TEST = (T["period"] == "test").to_numpy()33log(f"{len(T):,} cibles (val {VAL.sum():,}, test {TEST.sum():,}) ; E1 {est1.shape[1]} ; E2 {est2.shape[1]}")3435def err_of(est): return (np.exp(est - lnP[:, None]) - 1).astype(np.float32)36def metrics(err, sel):37 e = err[sel]; a = np.abs(e)38 ok = ~np.isnan(a); n = ok.sum(axis=0)39 with np.errstate(all="ignore"):40 return pd.DataFrame({"mdape": np.nanmedian(a, axis=0) * 100, "mape": np.nanmean(a, axis=0) * 100,41 "biais": np.nanmedian(e, axis=0) * 100, "biais_moy": np.nanmean(e, axis=0) * 100,42 "p5": np.nanmean(a <= 0.05, axis=0) * 100, "p10": np.nanmean(a <= 0.10, axis=0) * 100,43 "p20": np.nanmean(a <= 0.20, axis=0) * 100, "couverture": n / sel.sum() * 100, "n": n})44err1 = err_of(est1); err2 = err_of(est2)45del est146res = {}4748# ---------------------------------------------------------------- densité & strates49T["dens"] = pd.cut(T["n_2km_12m"], [-1, 4, 19, 49, 199, 10 ** 9], labels=["< 5", "5-19", "20-49", "50-199", "≥ 200"])50T["prix_q"] = pd.qcut(T["amount"], 5, labels=["Q1 (bas)", "Q2", "Q3", "Q4", "Q5 (haut)"])51T["trimestre"] = pd.to_datetime(T["date"]).dt.to_period("Q").astype(str)52T["age_cl"] = pd.cut(T["age"], [-1, 3, 15, 40, 70, 1000], labels=["≤ 3 ans", "4-15", "16-40", "41-70", "> 70"])53T["cluster"] = pd.factorize(T["code_mun"].astype(str) + "|" + T["trimestre"])[0]5455# ---------------------------------------------------------------- E1 : surface Erreur(k, d, t)56mv = metrics(err1, VAL).add_suffix("_val"); mt = metrics(err1, TEST).add_suffix("_test")57surf = pd.concat([E1, mv, mt], axis=1)58surf.to_csv(TAB / "surface_E1.csv", index=False)59ok = surf["couverture_val"] >= 9760best_e1 = surf[ok].sort_values("mdape_val").iloc[0]61res["E1_best_val"] = best_e1.to_dict()62log(f"E1 meilleure (val, couverture ≥ 97 %) : {best_e1[['sel','d','T','k','adj']].to_dict()} → val {best_e1.mdape_val:.2f} %, test {best_e1.mdape_test:.2f} %")63for (s, a), g in surf[ok].groupby(["sel", "adj"]):64 b = g.sort_values("mdape_val").iloc[0]65 res.setdefault("E1_best_by_sel_adj", {})[f"{s}|{a}"] = dict(k=int(b.k), d=float(b.d), T=int(b["T"]), mdape_val=float(b.mdape_val), mdape_test=float(b.mdape_test), couverture_test=float(b.couverture_test))6667# figure : MdAPE vs k selon d (T = 12, geo, hedo, test)68fig, axes = plt.subplots(1, 2, figsize=(9, 3.4), sharey=True)69for ax, a in zip(axes, ["none", "hedo"]):70 for j, d in enumerate([0.5, 1, 2, 5, 10, 25]):71 g = surf[(surf.sel == "geo") & (surf["T"] == 12) & (surf.adj == a) & (surf.d == d)].sort_values("k")72 ax.plot(g.k, g.mdape_test, color=PAL[j], lw=2, marker="o", ms=4, label=f"d ≤ {d} km")73 ax.set_xscale("log"); ax.set_xticks([1, 3, 5, 10, 20, 50]); ax.set_xticklabels([1, 3, 5, 10, 20, 50])74 ax.set_xlabel("nombre de comparables k"); ax.set_title({"none": "sans ajustement (temps seulement)", "hedo": "grille hédonique"}[a])75axes[0].set_ylabel("MdAPE, test 2025-26 (%)"); axes[1].legend(ncol=2, fontsize=8)76savefig(fig, "f1_mdape_vs_k")77# figure : carte d × T au meilleur k (geo, hedo)78fig, axes = plt.subplots(1, 2, figsize=(8, 3.2))79for ax, a in zip(axes, ["none", "hedo"]):80 g = surf[(surf.sel == "geo") & (surf.adj == a)].groupby(["d", "T"])["mdape_test"].min().unstack("T")81 im = ax.imshow(g.to_numpy(), cmap=CMAP, aspect="auto")82 ax.set_xticks(range(len(g.columns))); ax.set_xticklabels([f"{t} mois" for t in g.columns]); ax.set_yticks(range(len(g.index))); ax.set_yticklabels([f"{d} km" for d in g.index])83 for i in range(g.shape[0]):84 for j in range(g.shape[1]):85 v = g.iloc[i, j]; ax.text(j, i, f"{v:.1f}", ha="center", va="center", fontsize=7.5, color="white" if v < np.nanpercentile(g.to_numpy(), 40) else "#0b0b0b")86 ax.grid(False); ax.set_title({"none": "sans ajustement", "hedo": "grille hédonique"}[a] + " — MdAPE au meilleur k")87savefig(fig, "f2_surface_d_T")8889# ---------------------------------------------------------------- E2 : facteurs90mv2 = metrics(err2, VAL).add_suffix("_val"); mt2 = metrics(err2, TEST).add_suffix("_test")91fac = pd.concat([E2, mv2, mt2], axis=1); fac.to_csv(TAB / "facteurs_E2.csv", index=False)92anchor = (fac.k == 5) & (fac.d == 5) & (fac["T"] == 12)93base = anchor & (fac.sel == "geo") & (fac.w == "equal") & (fac.recon == "median")94res["E2_adjust_anchor"] = fac[base].set_index("adj")[["mdape_val", "mdape_test", "mape_test", "biais_test", "p10_test", "couverture_test"]].to_dict("index")95res["E2_sel_anchor"] = fac[anchor & (fac.adj == "hedo") & (fac.w == "equal") & (fac.recon == "median")].set_index("sel")[["mdape_val", "mdape_test", "p10_test"]].to_dict("index")96res["E2_wr_anchor"] = {f"{r.w}|{r.recon}": dict(mdape_val=r.mdape_val, mdape_test=r.mdape_test, p10_test=r.p10_test)97 for r in fac[anchor & (fac.adj == "hedo") & (fac.sel == "geo")].itertuples()}98# effets marginaux moyens (moyenne de MdAPE test sur toutes les autres modalités)99marg = {}100for f in ["sel", "adj", "w", "recon"]:101 marg[f] = fac.groupby(f)["mdape_test"].mean().round(2).to_dict()102res["E2_marginal"] = marg103anc_tab = fac.groupby(["k", "d", "T"])["mdape_test"].min().round(2).to_dict()104res["E2_anchor_best"] = {f"{k}|{d}|{t}": v for (k, d, t), v in anc_tab.items()}105# figure : ajustements à l'ancre (val et test)106fig, ax = plt.subplots(figsize=(7, 3.2))107g = fac[base].set_index("adj").loc[["raw", "none", "conv", "hedo", "hedo_reg", "local", "gbm", "role"]]108lab = {"raw": "aucun", "none": "temps seul", "conv": "grille conv.", "hedo": "hédonique", "hedo_reg": "hédo. régionale", "local": "locale", "gbm": "GBM", "role": "ratio au rôle"}109x = np.arange(len(g)); ax.bar(x - 0.2, g.mdape_val, 0.38, color=PAL[0], label="validation 2024"); ax.bar(x + 0.2, g.mdape_test, 0.38, color=PAL[1], label="test 2025-26")110for i, (v1, v2) in enumerate(zip(g.mdape_val, g.mdape_test)):111 ax.text(i - 0.2, v1 + 0.2, f"{v1:.1f}", ha="center", fontsize=7); ax.text(i + 0.2, v2 + 0.2, f"{v2:.1f}", ha="center", fontsize=7)112ax.set_xticks(x); ax.set_xticklabels([lab[a] for a in g.index]); ax.set_ylabel("MdAPE (%)"); ax.legend(); ax.set_title("Type d'ajustement — 5 comparables les plus proches, ≤ 5 km, ≤ 12 mois, médiane")113savefig(fig, "f3_ajustements")114115# ---------------------------------------------------------------- modèles de référence (M10, M11) et compétition116def level_shift(m_asof):117 mt_ = T["month"].to_numpy(); r = T["region"].map(rid).to_numpy()118 return IDX[mt_ - M_EVAL0, r, mt_ - 1] - IDX[mt_ - M_EVAL0, r, m_asof - 1]119Sx = S.loc[T["idx"]]120XF = np.column_stack([Sx.ln_floor, Sx.ln_land, Sx.age, Sx.age ** 2 / 100, Sx.lien == 1, Sx.lien >= 2, Sx.genre == 1, Sx.genre == 2, Sx.genre == 3, Sx.etages, Sx.neuf]).astype(float)121M10 = np.full(len(T), np.nan); M11 = np.full(len(T), np.nan); M11b = np.full(len(T), np.nan)122for a in load_json("asof.json")["asof"]:123 mask = (T["asof"] == a).to_numpy(); period = str(a)124 if not mask.any(): continue125 mod = load_json(f"models_{period}.json"); beta = np.array([mod["beta_global"][f] for f in mod["feats"]])126 shift = level_shift(mod["m_asof"])127 fe = (T["region"].map(mod["region_fe"]).fillna(0) + T["mrc"].map(mod["mrc_eff"]).fillna(0) + T["code_mun"].map(mod["mun_eff"]).fillna(0)128 + (T["code_mun"].astype(str) + "|" + T["unite_voisinage"].astype(str)).map(mod["uv_eff"]).fillna(0)).to_numpy()129 M10[mask] = (XF @ beta + fe + shift)[mask]130 M11[mask] = (np.load(OUT / f"gbm_role_pred_{period}.npy")[T["idx"].to_numpy()] + shift)[mask]131 M11b[mask] = (np.load(OUT / f"gbm_pred_{period}.npy")[T["idx"].to_numpy()] + shift)[mask]132133# M10 : recentrage médian calibré sur la validation 2024 (la MCO vise la moyenne du log, distribution asymétrique)134M10 = M10 + np.nanmedian((lnP - M10)[VAL])135def e2col(sel, k, d, Tm, adj, w, recon):136 m = (E2.sel == sel) & (E2.k == k) & (E2.d == d) & (E2["T"] == Tm) & (E2.adj == adj) & (E2.w == w) & (E2.recon == recon)137 return int(E2.index[m][0])138def e1col(sel, k, d, Tm, adj):139 m = (E1.sel == sel) & (E1.k == k) & (E1.d == d) & (E1["T"] == Tm) & (E1.adj == adj)140 return int(E1.index[m][0])141MODELS = {142 "M1 médiane des 3 plus proches (temps seul)": ("E2", e2col("geo", 3, 2, 6, "none", "equal", "median")),143 "M2 moyenne des 5 plus proches (temps seul)": ("E2", e2col("geo", 5, 5, 12, "none", "equal", "mean")),144 "M3 sélection par similarité physique": ("E2", e2col("sim", 5, 5, 12, "none", "equal", "median")),145 "M4 pondération par distance": ("E2", e2col("geo", 5, 5, 12, "none", "dist", "mean")),146 "M5 pondération par récence": ("E2", e2col("geo", 5, 5, 12, "none", "rec", "mean")),147 "M6 grille d'ajustement conventionnelle": ("E2", e2col("geo", 5, 5, 12, "conv", "equal", "median")),148 "M7 grille d'ajustement hédonique": ("E2", e2col("geo", 5, 5, 12, "hedo", "equal", "median")),149 "M8a grille locale (régression sur le voisinage)": ("E2", e2col("geo", 5, 5, 12, "local", "equal", "median")),150 "M8b ajustements non linéaires (GBM)": ("E2", e2col("geo", 5, 5, 12, "gbm", "equal", "median")),151 "M13 ratio au rôle (5 plus proches, même rôle)": ("E2", e2col("geo", 5, 5, 12, "role", "equal", "median")),152 "M10 hédonique globale directe": ("X", M10),153 "M11a gradient boosting (sans rôle)": ("X", M11b),154 "M11b gradient boosting (avec rôle)": ("X", M11),155}156def est_of(spec):157 kind, c = spec158 return est2[:, c].astype(float) if kind == "E2" else c159# M0 : meilleure configuration empirique choisie sur la validation (E1 ∪ E2, couverture ≥ 90 %)160cand_all = pd.concat([surf.assign(src="E1"), fac.assign(src="E2")], ignore_index=True)161cand_ok = cand_all[cand_all.couverture_val >= 97].sort_values("mdape_val")162b0 = cand_ok.iloc[0]163M0_est = (est2[:, int(b0.cfg)] if b0.src == "E2" else np.load(OUT / "est_E1.npy", mmap_mode="r")[:, int(b0.cfg)]).astype(float)164res["M0"] = b0[["src", "sel", "d", "T", "k", "adj", "w", "recon", "mdape_val", "mdape_test", "couverture_test"]].to_dict()165res["M0_sans_contrainte"] = cand_all[cand_all.couverture_val >= 50].sort_values("mdape_val").iloc[0][["src", "sel", "d", "T", "k", "adj", "w", "recon", "mdape_val", "mdape_test", "couverture_val", "couverture_test"]].to_dict()166MODELS["M0 meilleure configuration empirique (choisie en 2024, couverture ≥ 97 %)"] = ("X", M0_est)167# M9 : réconciliation optimisée = combinaison convexe (NNLS sur la validation) d'estimations de comparables168stack_cols = {a: e2col("geo", 5, 5, 12, a, "equal", "median") for a in ["none", "conv", "hedo", "hedo_reg", "local", "gbm", "role"]}169def stack(names, mats, fit_mask):170 """Combinaison convexe apprise par NNLS sur la validation ; si une composante manque pour une cible,171 les poids des composantes disponibles sont renormalisés."""172 A = np.column_stack(mats); okr = ~np.isnan(A).any(axis=1)173 w, _ = nnls(A[okr & fit_mask], lnP[okr & fit_mask].astype(float)); w = w / w.sum()174 W = np.where(np.isnan(A), 0.0, w[None, :]); den = W.sum(axis=1)175 out = np.nansum(A * W, axis=1) / np.where(den > 0, den, np.nan)176 return out, dict(zip(names, w.round(3)))177M9_est, w9 = stack(list(stack_cols), [est2[:, c].astype(float) for c in stack_cols.values()], VAL)178res["M9_weights"] = w9; MODELS["M9 réconciliation optimisée (poids appris en 2024)"] = ("X", M9_est)179M12_est, w12 = stack(["gbm_role", "comps_gbm", "comps_hedo"], [M11, est2[:, stack_cols["gbm"]].astype(float), est2[:, stack_cols["hedo"]].astype(float)], VAL)180res["M12_weights"] = w12; MODELS["M12 hybride comparables + GBM (poids appris en 2024)"] = ("X", M12_est)181182rows = []183ERR = {}184for name, spec in MODELS.items():185 e = (np.exp(est_of(spec) - lnP) - 1).astype(np.float32); ERR[name] = e186 for per, mask in [("val", VAL), ("test", TEST)]:187 m = metrics(e[:, None], mask).iloc[0].to_dict(); m.update(modele=name, periode=per); rows.append(m)188comp = pd.DataFrame(rows); comp.to_csv(TAB / "competition_modeles.csv", index=False)189pd.DataFrame({n.split(" ")[0]: est_of(s_) for n, s_ in MODELS.items()}).assign(idx=T["idx"].to_numpy()).to_parquet(OUT / "models_est.parquet", index=False)190res["competition"] = comp[comp.periode == "test"].set_index("modele")[["mdape", "mape", "biais", "p5", "p10", "p20", "couverture"]].round(2).to_dict("index")191res["competition_val"] = comp[comp.periode == "val"].set_index("modele")[["mdape", "p10", "couverture"]].round(2).to_dict("index")192log("compétition (test) :\n" + comp[comp.periode == "test"][["modele", "mdape", "mape", "biais", "p10", "couverture"]].round(2).to_string(index=False))193fig, ax = plt.subplots(figsize=(7.5, 5))194g = comp[comp.periode == "test"].sort_values("mdape", ascending=False)195cols = [PAL[0] if ("M1" in n and "M1 " not in n) or n.startswith(("M10", "M11", "M12")) else PAL[2] if n.startswith(("M0", "M9")) else PAL[1] for n in g.modele]196cols = [PAL[0] if n.startswith(("M10", "M11", "M12")) else (PAL[2] if n.startswith(("M0", "M9")) else PAL[1]) for n in g.modele]197ax.barh(g.modele, g.mdape, color=cols, height=0.65)198for i, v in enumerate(g.mdape): ax.text(v + 0.15, i, f"{v:.1f} %", va="center", fontsize=8)199ax.set_xlabel("MdAPE, test 2025-26 (%)"); ax.set_title("Compétition de modèles (orange : comparables ; vert : optimisés ; bleu : références externes)")200savefig(fig, "f4_competition")201202# ---------------------------------------------------------------- bootstrap par grappes (municipalité × trimestre) des différences appariées203rng = np.random.default_rng(7)204def paired(eA, eB, mask, B=400):205 a = np.abs(eA[mask]); b = np.abs(eB[mask]); cl = T["cluster"].to_numpy()[mask]206 okp = ~np.isnan(a) & ~np.isnan(b); a, b, cl = a[okp], b[okp], cl[okp]207 d = a - b208 cl_ids, inv = np.unique(cl, return_inverse=True); nc = len(cl_ids)209 sums = np.bincount(inv, d, nc); cnts = np.bincount(inv, minlength=nc)210 boots_mean = []; boots_med = []211 order = np.argsort(inv); d_sorted = d[order]; starts = np.searchsorted(inv[order], np.arange(nc)); ends = np.append(starts[1:], len(d))212 for _ in range(B):213 pick = rng.integers(0, nc, nc)214 boots_mean.append(sums[pick].sum() / cnts[pick].sum())215 if _ < 150:216 idxs = np.concatenate([np.arange(starts[p], ends[p]) for p in pick])217 boots_med.append(np.median(a[order][idxs]) - np.median(b[order][idxs]))218 return dict(n=int(len(d)), diff_mape=float(d.mean() * 100), ci_mape=[float(np.percentile(boots_mean, 2.5) * 100), float(np.percentile(boots_mean, 97.5) * 100)],219 diff_mdape=float((np.median(a) - np.median(b)) * 100), ci_mdape=[float(np.percentile(boots_med, 2.5) * 100), float(np.percentile(boots_med, 97.5) * 100)],220 part_A_meilleur=float((a < b).mean() * 100))221PAIRS = [222 ("M7 grille d'ajustement hédonique", "M2 moyenne des 5 plus proches (temps seul)"),223 ("M7 grille d'ajustement hédonique", "M6 grille d'ajustement conventionnelle"),224 ("M8a grille locale (régression sur le voisinage)", "M7 grille d'ajustement hédonique"),225 ("M8b ajustements non linéaires (GBM)", "M7 grille d'ajustement hédonique"),226 ("M11b gradient boosting (avec rôle)", "M0 meilleure configuration empirique (choisie en 2024, couverture ≥ 97 %)"),227 ("M11a gradient boosting (sans rôle)", "M0 meilleure configuration empirique (choisie en 2024, couverture ≥ 97 %)"),228 ("M12 hybride comparables + GBM (poids appris en 2024)", "M11b gradient boosting (avec rôle)"),229 ("M13 ratio au rôle (5 plus proches, même rôle)", "M7 grille d'ajustement hédonique"),230 ("M3 sélection par similarité physique", "M2 moyenne des 5 plus proches (temps seul)"),231 ("M9 réconciliation optimisée (poids appris en 2024)", "M7 grille d'ajustement hédonique"),232 ("M1 médiane des 3 plus proches (temps seul)", "M2 moyenne des 5 plus proches (temps seul)"),233]234res["paires"] = {f"{a} vs {b}": paired(ERR[a], ERR[b], TEST) for a, b in PAIRS}235pd.DataFrame([dict(A=a, B=b, **res["paires"][f"{a} vs {b}"]) for a, b in PAIRS]).to_csv(TAB / "paires_bootstrap.csv", index=False)236237# ---------------------------------------------------------------- strates238def strata(err, by):239 out = []240 for lev, g in T[TEST].groupby(by, observed=True):241 m = np.zeros(len(T), bool); m[g.index] = True242 r = metrics(err[:, None], m).iloc[0].to_dict(); r[by] = str(lev); out.append(r)243 return pd.DataFrame(out)244KEY = ["M2 moyenne des 5 plus proches (temps seul)", "M7 grille d'ajustement hédonique", "M8a grille locale (régression sur le voisinage)",245 "M13 ratio au rôle (5 plus proches, même rôle)", "M0 meilleure configuration empirique (choisie en 2024, couverture ≥ 97 %)", "M10 hédonique globale directe", "M11b gradient boosting (avec rôle)", "M12 hybride comparables + GBM (poids appris en 2024)"]246for by in ["region_nom", "dens", "prix_q", "trimestre", "age_cl", "lien"]:247 tabs = []248 for name in KEY:249 s = strata(ERR[name], by); s["modele"] = name.split(" ")[0]; tabs.append(s)250 st = pd.concat(tabs); st.to_csv(TAB / f"strates_{by}.csv", index=False)251 res[f"strates_{by}"] = st.pivot(index=by, columns="modele", values="mdape").round(2).to_dict("index")252 res[f"strates_{by}_n"] = st[st.modele == "M7"].set_index(by)["n"].to_dict()253# figure : densité254fig, ax = plt.subplots(figsize=(6.5, 3.4))255st = pd.read_csv(TAB / "strates_dens.csv")256order = ["< 5", "5-19", "20-49", "50-199", "≥ 200"]257for j, (name, lab_) in enumerate([("M2", "M2 5 plus proches"), ("M7", "M7 grille hédonique"), ("M0", "M0 meilleure config."), ("M11b", "M11b GBM"), ("M12", "M12 hybride")]):258 g = st[st.modele == name].set_index("dens").loc[order]259 ax.plot(order, g.mdape, color=PAL[j], lw=2, marker="o", ms=4, label=lab_)260ax.set_xlabel("ventes disponibles dans un rayon de 2 km au cours des 12 mois précédents"); ax.set_ylabel("MdAPE, test (%)"); ax.legend(fontsize=8); ax.set_title("Précision selon la liquidité du sous-marché")261savefig(fig, "f5_densite")262263# ---------------------------------------------------------------- hypothèses264H = {}265# H1 : proche mais dissemblable vs loin mais semblable (k = 5, T = 12, sans ajustement)266H["H1"] = {"geo_d1": surf.loc[e1col("geo", 5, 1, 12, "none"), ["mdape_test", "couverture_test"]].to_dict(),267 "sim_d10": surf.loc[e1col("sim", 5, 10, 12, "none"), ["mdape_test", "couverture_test"]].to_dict(),268 "geo_d10": surf.loc[e1col("geo", 5, 10, 12, "none"), ["mdape_test", "couverture_test"]].to_dict(),269 "mix_d10": surf.loc[e1col("mix", 5, 10, 12, "none"), ["mdape_test", "couverture_test"]].to_dict(),270 "sim_d10_hedo": surf.loc[e1col("sim", 5, 10, 12, "hedo"), ["mdape_test", "couverture_test"]].to_dict(),271 "geo_d10_hedo": surf.loc[e1col("geo", 5, 10, 12, "hedo"), ["mdape_test", "couverture_test"]].to_dict(),272 "mix_d10_hedo": surf.loc[e1col("mix", 5, 10, 12, "hedo"), ["mdape_test", "couverture_test"]].to_dict()}273# H2 : k optimal selon la densité (geo, d = 10, T = 12, hedo et none)274h2 = {}275for a in ["none", "hedo"]:276 for lev in order:277 m = TEST & (T["dens"] == lev).to_numpy()278 cols_ = [e1col("geo", k, 10, 12, a) for k in [1, 3, 5, 7, 10, 15, 20, 30, 50]]279 md = np.nanmedian(np.abs(err1[m][:, cols_]), axis=0) * 100280 h2[f"{a}|{lev}"] = dict(zip([1, 3, 5, 7, 10, 15, 20, 30, 50], md.round(2)))281H["H2"] = h2282fig, axes = plt.subplots(1, 2, figsize=(9, 3.3), sharey=True)283for ax, a in zip(axes, ["none", "hedo"]):284 for j, lev in enumerate(order):285 v = h2[f"{a}|{lev}"]; ax.plot(list(v), list(v.values()), color=PAL[j], lw=2, marker="o", ms=3.5, label=f"densité {lev}")286 ax.set_xscale("log"); ax.set_xticks([1, 3, 5, 10, 20, 50]); ax.set_xticklabels([1, 3, 5, 10, 20, 50]); ax.set_xlabel("k"); ax.set_title({"none": "temps seul", "hedo": "grille hédonique"}[a] + " — d ≤ 10 km, 12 mois")287axes[0].set_ylabel("MdAPE, test (%)"); axes[1].legend(fontsize=8)288savefig(fig, "f6_k_par_densite")289# H3 : local vs hédonique fixe selon la densité, aux trois ancres290h3 = {}291for (k, d, Tm) in [(3, 2, 6), (5, 5, 12), (10, 10, 24)]:292 for lev in order:293 m = TEST & (T["dens"] == lev).to_numpy()294 h3[f"{k}|{d}|{Tm}|{lev}"] = {a: float(np.nanmedian(np.abs(err2[m][:, e2col("geo", k, d, Tm, a, "equal", "median")])) * 100) for a in ["conv", "hedo", "hedo_reg", "local", "gbm"]}295H["H3"] = h3296# H4 : pondération selon la densité (ancre 5/5/12, hedo)297h4 = {}298for lev in order:299 m = TEST & (T["dens"] == lev).to_numpy()300 h4[lev] = {f"{w}|{r}": float(np.nanmedian(np.abs(err2[m][:, e2col("geo", 5, 5, 12, "hedo", w, r)])) * 100) for (w, r) in [("equal", "median"), ("equal", "mean"), ("dist", "mean"), ("rec", "mean"), ("sim", "mean"), ("all", "mean"), ("all", "median")]}301H["H4"] = h4302# H5 : méthodes simples vs ML selon l'homogénéité (dispersion des prix ajustés des 10 voisins à 10 km / 24 mois)303het = sd2[:, e2col("geo", 10, 10, 24, "hedo", "equal", "median")].astype(float)304T["heterog"] = pd.qcut(pd.Series(np.where(np.isnan(het), np.nanmax(het), het)), 3, labels=["homogène", "moyen", "hétérogène"])305h5 = {}306for lev in ["homogène", "moyen", "hétérogène"]:307 m = TEST & (T["heterog"] == lev).to_numpy()308 h5[lev] = {n.split(" ")[0]: float(np.nanmedian(np.abs(ERR[n][m])) * 100) for n in KEY}309 h5[lev]["sd_median"] = float(np.nanmedian(het[m]))310H["H5"] = h5311res["hypotheses"] = H312313# ---------------------------------------------------------------- intervalles d'incertitude (couverture)314intervals = {}315for a in ["none", "hedo", "gbm"]:316 for (k, d, Tm) in [(5, 5, 12), (10, 10, 24)]:317 c = e2col("geo", k, d, Tm, a, "equal", "median"); e = est2[:, c].astype(float); s = sd2[:, c].astype(float)318 okv = VAL & ~np.isnan(e) & ~np.isnan(s) & (s > 0)319 z_needed = np.abs(lnP[okv] - e[okv]) / s[okv]; z80 = float(np.percentile(z_needed, 80))320 okt = TEST & ~np.isnan(e) & ~np.isnan(s) & (s > 0)321 cov_nominal = float(np.mean(np.abs(lnP[okt] - e[okt]) <= 1.2816 * s[okt]) * 100)322 cov_cal = float(np.mean(np.abs(lnP[okt] - e[okt]) <= z80 * s[okt]) * 100)323 width = float(np.median(2 * z80 * s[okt]) * 100)324 intervals[f"{a}|{k}|{d}|{Tm}"] = dict(z80_val=z80, couverture_nominale_80=cov_nominal, couverture_calibree=cov_cal, largeur_mediane_pct=width, n_test=int(okt.sum()))325res["intervalles"] = intervals326327# ---------------------------------------------------------------- règles conditionnelles (système de recommandation)328cond = {}; est_cond = np.full(len(T), np.nan)329for lev in order:330 mv_ = VAL & (T["dens"] == lev).to_numpy(); mt_ = TEST & (T["dens"] == lev).to_numpy()331 m1 = metrics(err1, mv_); m2 = metrics(err2, mv_)332 allm = pd.concat([pd.concat([E1, m1], axis=1).assign(src="E1"), pd.concat([E2, m2], axis=1).assign(src="E2")], ignore_index=True)333 b = allm[allm.couverture >= 95].sort_values("mdape").iloc[0]334 e = (err1 if b.src == "E1" else err2)[:, int(b.cfg)]335 est_cond[mt_] = e[mt_]336 cond[lev] = dict(b[["src", "sel", "d", "T", "k", "adj", "w", "recon", "mdape", "couverture"]].to_dict(), mdape_test=float(np.nanmedian(np.abs(e[mt_])) * 100), n_test=int(mt_.sum()))337res["regles_conditionnelles"] = cond338res["regles_conditionnelles_global_test"] = dict(mdape=float(np.nanmedian(np.abs(est_cond[TEST])) * 100), p10=float(np.nanmean(np.abs(est_cond[TEST]) <= 0.1) * 100), couverture=float(np.mean(~np.isnan(est_cond[TEST])) * 100))339pd.DataFrame(cond).T.to_csv(TAB / "regles_conditionnelles.csv")340341# ---------------------------------------------------------------- stabilité temporelle, indice, divers342tq = []343for q, g in T[TEST].groupby("trimestre"):344 m = np.zeros(len(T), bool); m[g.index] = True345 tq.append(dict(trimestre=q, n=int(m.sum()), **{n.split(" ")[0]: float(np.nanmedian(np.abs(ERR[n][m])) * 100) for n in KEY},346 **{"biais_" + n.split(" ")[0]: float(np.nanmedian(ERR[n][m]) * 100) for n in KEY}))347pd.DataFrame(tq).to_csv(TAB / "stabilite_trimestres.csv", index=False); res["stabilite"] = tq348fig, ax = plt.subplots(figsize=(7, 3.2))349idxf = pd.read_csv(OUT / "index_full.csv")350for j, (r, lab_) in enumerate([("QC", "Québec"), ("06", "Montréal"), ("03", "Capitale-Nationale"), ("07", "Outaouais"), ("05", "Estrie")]):351 g = idxf[idxf.region.astype(str).str.zfill(2) == r] if r != "QC" else idxf[idxf.region == "QC"]352 if len(g): ax.plot(pd.to_datetime(g.ym), np.exp(g.log_idx) * 100, color=PAL[j], lw=2 if r == "QC" else 1.4, label=lab_)353ax.set_ylabel("indice hédonique (2021-01 = 100)"); ax.legend(fontsize=8); ax.set_title("Indice de prix mensuel (plein échantillon, descriptif)")354savefig(fig, "f7_indice")355res["dataset"] = load_json("dataset_stats.json")356desc = {}357for per in ["train", "val", "test"]:358 d = S[S["period"] == per]359 desc[per] = dict(n=int(len(d)), prix_med=float(d.amount.median()), prix_moy=float(d.amount.mean()), aire_med=float(d.floor_area.median()),360 terrain_med=float(d.land_area.median()), age_med=float(d.age.median()), detache=float((d.lien == 0).mean() * 100),361 etages_moy=float(d.etages.mean()), ratio_role_med=float(d.ratio_role.median()))362for per, m in [("val", VAL), ("test", TEST)]:363 desc[per]["dens_med"] = float(np.median(T[m].n_2km_12m)); desc[per]["n_cand_med"] = float(np.median(T[m].n_cand))364res["descriptives"] = desc365res["n_targets"] = dict(val=int(VAL.sum()), test=int(TEST.sum()))366res["densite_distribution"] = T[TEST]["dens"].value_counts().sort_index().to_dict()367res["n_cand_median"] = dict(val=float(np.median(T[VAL].n_cand)), test=float(np.median(T[TEST].n_cand)))368save_json(res, "results.json")369log("analyse terminée")370