spb/traite-evaluation-qc

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>

added analysis/python/cases/ch37_ml.py +411 −0

		@@ -0,0 +1,411 @@
	1	+#!/usr/bin/env python3
	2	+"""Chapitre 37 — Apprentissage automatique et explicabilité : exemple intégral sur les 600 ventes du cas 1
	3	+(Gatineau, DONNÉES SYNTHÉTIQUES, graine 20261031).
	4	+
	5	+Le modèle hédonique de référence est `traite.cas1.hedonique` (semi-log, MCO, HC1) ; il n'est PAS redéfini ici. Ce script lui
	6	+oppose un arbre de régression, une forêt aléatoire (Breiman 2001), un gradient boosting (Friedman 2001 ; scikit-learn) et XGBoost
	7	+(Chen et Guestrin 2016), sous trois schémas de validation : plis aléatoires (les MÊMES dix plis que `ch18_hedonique.validation_croisee`,
	8	+graine 20261031), blocs spatiaux stratifiés par secteur (k-moyennes sur les coordonnées), et partition temporelle (ventes de 2026
	9	+prédites à partir de 2023-2025). Il ajoute la validation « un secteur exclu » (extrapolation hors support), un scénario où la loi de
	10	+génération est perturbée par une interaction et une rupture déclarées (le boosting y gagne), les valeurs SHAP (TreeSHAP) du boosting,
	11	+la dépendance partielle, l'importance par permutation, la prédiction du sujet par chaque modèle et les chiffres des exercices.
	12	+Usage (racine du dépôt) : PYTHONPATH=analysis/python python3 analysis/python/cases/ch37_ml.py"""
	13	+from __future__ import annotations
	14	+import sys, math, pathlib, warnings
	15	+import numpy as np, pandas as pd
	16	+sys.path.insert(0, str(pathlib.Path(__file__).resolve().parents[1]))
	17	+from traite import cas1
	18	+from traite import hedonique as h
	19	+
	20	+warnings.filterwarnings("ignore", category=FutureWarning)
	21	+ROOT = pathlib.Path(__file__).resolve().parents[3]
	22	+GRAINE = 20261031
	23	+N_ARBRES_FORET = 500
	24	+
	25	+# ----------------------------------------------------------------------------- paramètres de génération (README synthétique)
	26	+_g = (1 + 9.06 / 100) ** 0.5 - 1
	27	+VRAIS = {"ln_S": 0.55, "ln_T": 0.10, "age": -0.0055, "age2": 0.00003, "qualite": 0.06, "garage": 0.035, "sous_sol": 0.03, "piscine": 0.02,
	28	+ "etages": 0.02, "dist_centre": -0.006, "sigma": 0.085, "mois": math.log1p(_g) / 12,
	29	+ "sect": {"Hull": 0.0, "Gatineau": math.log(464_900 / 475_000), "Aylmer": math.log(530_000 / 475_000),
	30	+ "Buckingham": math.log(395_000 / 475_000), "Masson-Angers": math.log(395_000 / 475_000) - 0.02}}
	31	+# scénario non linéaire déclaré (section 37.6) : prime « patrimoine rénové » et « mur de la mi-vie »
	32	+SCENARIO = {"prime_qualite_age": 0.12, "qualite_min": 4, "age_min_prime": 40, "mur": -0.10, "age_mur": (25, 45)}
	33	+
	34	+# ----------------------------------------------------------------------------- régresseurs
	35	+VARIABLES_ML = ["superficie", "terrain", "age", "qualite", "garage", "sous_sol", "piscine", "etages", "dist_centre", "lat", "lon", "mois"] + \
	36	+ [f"sect_{s}" for s in cas1.SECTEURS]
	37	+ETIQUETTES = {"superficie": "Superficie habitable (m²)", "terrain": "Terrain (m²)", "age": "Âge (années)", "qualite": "Qualité (1-5)",
	38	+ "garage": "Garage (places)", "sous_sol": "Sous-sol fini", "piscine": "Piscine", "etages": "À étages", "dist_centre": "Distance au centre (km)",
	39	+ "lat": "Latitude", "lon": "Longitude", "mois": "Mois depuis 2024-07", "sect_Hull": "Secteur Hull", "sect_Gatineau": "Secteur Gatineau",
	40	+ "sect_Aylmer": "Secteur Aylmer", "sect_Buckingham": "Secteur Buckingham", "sect_Masson-Angers": "Secteur Masson-Angers"}
	41	+
	42	+def matrice_ml(df: pd.DataFrame, s: dict | None = None) -> pd.DataFrame:
	43	+ """Régresseurs bruts pour les méthodes à base d'arbres : aucune transformation (les arbres sont invariants aux transformations
	44	+ monotones), secteur en indicatrices complètes, coordonnées brutes, temps en mois depuis la date de référence du jeu."""
	45	+ X = pd.DataFrame(index=df.index)
	46	+ X["superficie"] = df["superficie_habitable_m2"].astype(float); X["terrain"] = df["superficie_terrain_m2"].astype(float)
	47	+ X["age"] = df["age"].astype(float); X["qualite"] = df["qualite_1a5"].astype(float); X["garage"] = df["garage_places"].astype(float)
	48	+ X["sous_sol"] = df["sous_sol_fini"].astype(float); X["piscine"] = df["piscine"].astype(float)
	49	+ X["etages"] = (df["type"] == "à étages").astype(float); X["dist_centre"] = df["dist_centre_km"].astype(float)
	50	+ X["lat"] = df["lat"].astype(float); X["lon"] = df["lon"].astype(float)
	51	+ d = df["date_vente"] if "date_vente" in df else pd.Series([pd.Timestamp(s["date_evaluation"])] * len(df), index=df.index)
	52	+ X["mois"] = ((pd.to_datetime(d) - cas1.REF_DATE).dt.days / 30.4375).astype(float)
	53	+ for sct in cas1.SECTEURS:
	54	+ X[f"sect_{sct}"] = (df["secteur"] == sct).astype(float)
	55	+ return X[VARIABLES_ML]
	56	+
	57	+# ----------------------------------------------------------------------------- modèles
	58	+class Hedonique:
	59	+ """Enveloppe du modèle hédonique de référence (traite.hedonique.mco sur traite.cas1.matrice) avec l'interface fit/predict."""
	60	+ def __init__(self, colonnes_exclues=()): self.excl = list(colonnes_exclues)
	61	+ def fit(self, X: pd.DataFrame, y):
	62	+ self.cols = [c for c in X.columns if c not in self.excl]
	63	+ self.r = h.mco(X[self.cols].values, np.asarray(y, float)); self.beta = self.r["beta"]; return self
	64	+ def predict(self, X: pd.DataFrame): return X[self.cols].values @ self.beta
	65	+
	66	+def _foret(**k):
	67	+ from sklearn.ensemble import RandomForestRegressor
	68	+ return RandomForestRegressor(**{**dict(n_estimators=N_ARBRES_FORET, min_samples_leaf=2, max_features=0.5, random_state=GRAINE, n_jobs=-1), **k})
	69	+
	70	+def _boosting(**k):
	71	+ from sklearn.ensemble import GradientBoostingRegressor
	72	+ return GradientBoostingRegressor(**{**dict(n_estimators=1200, learning_rate=0.03, max_depth=2, subsample=0.8, min_samples_leaf=5, random_state=GRAINE), **k})
	73	+
	74	+def _xgboost(**k):
	75	+ import xgboost as xgb
	76	+ return xgb.XGBRegressor(**{**dict(n_estimators=1200, learning_rate=0.03, max_depth=2, subsample=0.8, colsample_bytree=0.8, min_child_weight=5,
	77	+ reg_lambda=1.0, random_state=GRAINE, n_jobs=4, verbosity=0), **k})
	78	+
	79	+def _arbre(**k):
	80	+ from sklearn.tree import DecisionTreeRegressor
	81	+ return DecisionTreeRegressor(**{**dict(max_depth=3, min_samples_leaf=5, random_state=GRAINE), **k})
	82	+
	83	+HYPERPARAMETRES = {"Arbre (profondeur 3)": "profondeur 3, feuilles ≥ 5", "Forêt aléatoire": f"{N_ARBRES_FORET} arbres, feuilles ≥ 2, 50 % des variables par division, bootstrap",
	84	+ "Gradient boosting": "1 200 arbres de profondeur 2, ν = 0,03, sous-échantillon 80 %, feuilles ≥ 5",
	85	+ "XGBoost": "1 200 arbres de profondeur 2, ν = 0,03, sous-échantillon 80 %, 80 % des colonnes, poids min. 5, λ = 1"}
	86	+MODELES = {"Hédonique (référence)": ("hedo", Hedonique), "Arbre (profondeur 3)": ("ml", _arbre), "Forêt aléatoire": ("ml", _foret),
	87	+ "Gradient boosting": ("ml", _boosting), "XGBoost": ("ml", _xgboost)}
	88	+
	89	+# ----------------------------------------------------------------------------- plis
	90	+def plis_aleatoires(n: int, nfolds: int = 10, graine: int = GRAINE) -> list[np.ndarray]:
	91	+ """Les mêmes dix plis que ch18_hedonique.validation_croisee (même graine, même permutation, même découpage)."""
	92	+ rng = np.random.default_rng(graine); idx = rng.permutation(n); return list(np.array_split(idx, nfolds))
	93	+
	94	+def plis_spatiaux(v: pd.DataFrame, par_secteur: int = 2, graine: int = GRAINE) -> tuple[list[np.ndarray], np.ndarray]:
	95	+ """Blocs spatiaux stratifiés : dans chaque secteur, k-moyennes (k = 2) sur les coordonnées en km → 10 blocs contigus.
	96	+ Chaque secteur reste représenté à l'estimation (par l'autre moitié), mais le voisinage immédiat d'une vente testée ne l'est jamais."""
	97	+ from sklearn.cluster import KMeans
	98	+ km = h.coordonnees_km(v["lat"].to_numpy(), v["lon"].to_numpy()); bloc = np.zeros(len(v), int); b = 0
	99	+ for sct in cas1.SECTEURS:
	100	+ ii = np.where(v["secteur"].to_numpy() == sct)[0]
	101	+ lab = KMeans(par_secteur, n_init=10, random_state=graine).fit_predict(km[ii]); bloc[ii] = b + lab; b += par_secteur
	102	+ return [np.where(bloc == j)[0] for j in range(b)], bloc
	103	+
	104	+def plis_secteurs(v: pd.DataFrame) -> list[np.ndarray]:
	105	+ """Un secteur exclu à la fois (extrapolation hors support spatial)."""
	106	+ return [np.where(v["secteur"].to_numpy() == sct)[0] for sct in cas1.SECTEURS]
	107	+
	108	+def partition_temporelle(v: pd.DataFrame, date: str = "2026-01-01") -> tuple[np.ndarray, np.ndarray]:
	109	+ test = np.where(v["date_vente"] >= pd.Timestamp(date))[0]; return np.setdiff1d(np.arange(len(v)), test), test
	110	+
	111	+# ----------------------------------------------------------------------------- métriques
	112	+def metriques(pred_log: np.ndarray, y: np.ndarray, facteur: float = 1.0) -> dict:
	113	+ """RMSE du logarithme, R² hors échantillon, biais moyen (log), RMSE en dollars, MdAPE, statistiques de ratio IAAO
	114	+ (traite.hedonique.statistiques_ratio). `facteur` : correction de retransformation (Duan) appliquée aux prédictions en niveau."""
	115	+ m = ~np.isnan(pred_log); e = pred_log[m] - y[m]; P = np.exp(y[m]); V = np.exp(pred_log[m]) * facteur
	116	+ r = h.statistiques_ratio(V, P)
	117	+ return {"n": int(m.sum()), "rmse_log": float(np.sqrt(np.mean(e ** 2))), "biais_log": float(e.mean()),
	118	+ "r2_hors": float(1 - np.sum(e ** 2) / np.sum((y[m] - y[m].mean()) ** 2)), "rmse_dollars": float(np.sqrt(np.mean((V - P) ** 2))),
	119	+ "mdape": float(np.median(np.abs(V / P - 1)) * 100), "COD": r["COD"], "PRD": r["PRD"], "PRB": r["PRB"], "PRB_se": r["PRB_se"],
	120	+ "ratio_median": r["ratio_median"]}
	121	+
	122	+def valider(constructeur, genre: str, Xh: pd.DataFrame, Xm: pd.DataFrame, y: np.ndarray, plis: list[np.ndarray]) -> np.ndarray:
	123	+ """Prédictions hors échantillon (chaque pli prédit par un modèle estimé sans lui)."""
	124	+ X = Xh if genre == "hedo" else Xm; n = len(y); pred = np.full(n, np.nan); tous = np.arange(n)
	125	+ for f in plis:
	126	+ tr = np.setdiff1d(tous, f); m = constructeur(); m.fit(X.iloc[tr], y[tr]); pred[f] = np.asarray(m.predict(X.iloc[f]), float)
	127	+ return pred
	128	+
	129	+def facteur_duan(pred_log: np.ndarray, y: np.ndarray) -> float:
	130	+ """Facteur de retransformation estimé sur les résidus hors échantillon (honnête pour les modèles qui surajustent en échantillon)."""
	131	+ m = ~np.isnan(pred_log); return float(np.mean(np.exp(y[m] - pred_log[m])))
	132	+
	133	+def comparaison(v: pd.DataFrame, Xh: pd.DataFrame, Xm: pd.DataFrame, y: np.ndarray) -> dict:
	134	+ """Trois schémas de validation pour les cinq modèles. Retourne les prédictions hors échantillon et les métriques."""
	135	+ n = len(y); pa = plis_aleatoires(n); ps, bloc = plis_spatiaux(v); tr_t, te_t = partition_temporelle(v)
	136	+ out = {"plis": {"aleatoire": pa, "spatial": ps, "bloc_spatial": bloc, "temporel": (tr_t, te_t)}, "pred": {}, "metriques": {}, "duan": {}}
	137	+ for nom, (genre, ctor) in MODELES.items():
	138	+ p_a = valider(ctor, genre, Xh, Xm, y, pa); F = facteur_duan(p_a, y)
	139	+ p_s = valider(ctor, genre, Xh, Xm, y, ps)
	140	+ p_t = valider(ctor, genre, Xh, Xm, y, [te_t]) # un seul pli de test : les ventes de 2026
	141	+ out["pred"][nom] = {"aleatoire": p_a, "spatial": p_s, "temporel": p_t}; out["duan"][nom] = F
	142	+ out["metriques"][nom] = {"aleatoire": metriques(p_a, y, F), "spatial": metriques(p_s, y, F), "temporel": metriques(p_t, y, F)}
	143	+ # biais temporel par mois de 2026 pour le boosting et l'hédonique
	144	+ mois = Xm["mois"].to_numpy()
	145	+ out["biais_mensuel_2026"] = {nom: pd.Series(out["pred"][nom]["temporel"][te_t] - y[te_t]).groupby(np.floor(mois[te_t])).mean().to_dict()
	146	+ for nom in ("Hédonique (référence)", "Gradient boosting", "Forêt aléatoire")}
	147	+ return out
	148	+
	149	+def boosting_detendance(v: pd.DataFrame, Xh: pd.DataFrame, Xm: pd.DataFrame, y: np.ndarray) -> dict:
	150	+ """Remède à l'extrapolation temporelle : boosting estimé sur ln P − β̂_mois·mois (tendance de l'hédonique estimée sur l'entraînement),
	151	+ sans la variable mois, puis tendance rajoutée. Variante « feuilles ≥ 1 » pour montrer l'effet de la feuille extrême."""
	152	+ tr, te = partition_temporelle(v); n = len(y); cols = list(Xh.columns)
	153	+ hd = Hedonique().fit(Xh.iloc[tr], y[tr]); bm = hd.beta[cols.index("mois")]; mois = Xm["mois"].to_numpy()
	154	+ Xs = Xm.drop(columns=["mois"])
	155	+ g = _boosting().fit(Xs.iloc[tr], y[tr] - bm * mois[tr]); p = np.full(n, np.nan); p[te] = g.predict(Xs.iloc[te]) + bm * mois[te]
	156	+ g1 = _boosting(min_samples_leaf=1).fit(Xm.iloc[tr], y[tr]); p1 = np.full(n, np.nan); p1[te] = g1.predict(Xm.iloc[te])
	157	+ return {"beta_mois_entrainement": float(bm), "pred_detendance": p, "metriques_detendance": metriques(p, y), "pred_feuille1": p1,
	158	+ "metriques_feuille1": metriques(p1, y), "n_entrainement_apres_mois_17": int((mois[tr] > 17).sum())}
	159	+
	160	+def validation_secteurs(v: pd.DataFrame, Xh: pd.DataFrame, Xm: pd.DataFrame, y: np.ndarray) -> dict:
	161	+ """Un secteur exclu à la fois. L'hédonique n'a aucune information sur un secteur absent : son indicatrice est retirée et le secteur
	162	+ est prédit au niveau du secteur de référence (Hull) ; si Hull est exclu, Gatineau devient la référence. Ce choix est un JUGEMENT
	163	+ de l'analyste, explicite ; les arbres font implicitement le leur à partir des coordonnées et de la distance au centre."""
	164	+ n = len(y); plis = plis_secteurs(v); out = {"pred": {}, "metriques": {}, "biais_secteur": {}}
	165	+ for nom, (genre, ctor) in MODELES.items():
	166	+ pred = np.full(n, np.nan)
	167	+ for sct, f in zip(cas1.SECTEURS, plis):
	168	+ tr = np.setdiff1d(np.arange(n), f)
	169	+ if genre == "hedo":
	170	+ drop = f"sect_{sct}" if sct != "Hull" else "sect_Gatineau"; m = Hedonique([drop]).fit(Xh.iloc[tr], y[tr]); pred[f] = m.predict(Xh.iloc[f])
	171	+ else:
	172	+ m = ctor(); m.fit(Xm.iloc[tr], y[tr]); pred[f] = m.predict(Xm.iloc[f])
	173	+ out["pred"][nom] = pred; out["metriques"][nom] = metriques(pred, y)
	174	+ out["biais_secteur"][nom] = {sct: float(np.mean(pred[f] - y[f])) for sct, f in zip(cas1.SECTEURS, plis)}
	175	+ return out
	176	+
	177	+# ----------------------------------------------------------------------------- scénario non linéaire déclaré
	178	+def perturbation_scenario(v: pd.DataFrame) -> np.ndarray:
	179	+ """ln P† = ln P + 0,12·1[qualité ≥ 4 et âge ≥ 40] − 0,10·1[25 ≤ âge ≤ 45] (hypothèse pédagogique déclarée)."""
	180	+ q = v["qualite_1a5"].to_numpy(float); a = v["age"].to_numpy(float); S = SCENARIO
	181	+ return S["prime_qualite_age"] * ((q >= S["qualite_min"]) & (a >= S["age_min_prime"])) + S["mur"] * ((a >= S["age_mur"][0]) & (a <= S["age_mur"][1]))
	182	+
	183	+def scenario_nonlineaire(v: pd.DataFrame, Xh: pd.DataFrame, Xm: pd.DataFrame, y: np.ndarray) -> dict:
	184	+ """Mêmes plis aléatoires ; hédonique de référence (mal spécifié désormais), forêt, boosting, XGBoost, et hédonique respécifié
	185	+ avec les deux indicatrices (ce que l'analyste ferait après avoir lu la dépendance partielle de l'âge et les interactions SHAP)."""
	186	+ delta = perturbation_scenario(v); y2 = y + delta; pa = plis_aleatoires(len(y))
	187	+ q = v["qualite_1a5"].to_numpy(float); a = v["age"].to_numpy(float); S = SCENARIO
	188	+ Xh2 = Xh.copy(); Xh2["prime"] = ((q >= S["qualite_min"]) & (a >= S["age_min_prime"])).astype(float); Xh2["mur"] = ((a >= S["age_mur"][0]) & (a <= S["age_mur"][1])).astype(float)
	189	+ out = {"delta": delta, "y": y2, "n_prime": int(Xh2["prime"].sum()), "n_mur": int(Xh2["mur"].sum()), "pred": {}, "metriques": {}, "Xh2": Xh2}
	190	+ for nom, (genre, ctor) in MODELES.items():
	191	+ if nom == "Arbre (profondeur 3)": continue
	192	+ p = valider(ctor, genre, Xh, Xm, y2, pa); out["pred"][nom] = p; out["metriques"][nom] = metriques(p, y2, facteur_duan(p, y2))
	193	+ p = valider(Hedonique, "hedo", Xh2, Xm, y2, pa); out["pred"]["Hédonique respécifié"] = p; out["metriques"]["Hédonique respécifié"] = metriques(p, y2, facteur_duan(p, y2))
	194	+ r = Hedonique().fit(Xh2, y2); cols = list(Xh2.columns)
	195	+ out["coef_respecifie"] = {"prime": float(r.beta[cols.index("prime")]), "mur": float(r.beta[cols.index("mur")]), "se_prime": float(r.r["se_hc1"][cols.index("prime")]),
	196	+ "se_mur": float(r.r["se_hc1"][cols.index("mur")]), "sigma": float(np.sqrt(r.r["sigma2"])), "r2": float(r.r["r2"])}
	197	+ out["gbm_scenario"] = _boosting().fit(Xm, y2)
	198	+ out["hedo_scenario"] = Hedonique().fit(Xh, y2)
	199	+ return out
	200	+
	201	+# ----------------------------------------------------------------------------- explicabilité
	202	+def shap_boosting(gbm, Xm: pd.DataFrame, xs: pd.DataFrame) -> dict:
	203	+ """TreeSHAP (Lundberg et al. 2020) du gradient boosting : valeurs SHAP des 600 ventes, du sujet, valeur de base, importance
	204	+ moyenne |φ|, interactions SHAP moyennes. Vérifie l'axiome d'efficacité (somme des φ + base = prédiction)."""
	205	+ import shap
	206	+ ex = shap.TreeExplainer(gbm); sv = np.asarray(ex.shap_values(Xm)); base = float(np.ravel(ex.expected_value)[0])
	207	+ svs = np.asarray(ex.shap_values(xs))[0]
	208	+ si = np.asarray(ex.shap_interaction_values(Xm)); inter = np.abs(si).mean(0).copy(); np.fill_diagonal(inter, 0.0)
	209	+ inter_df = pd.DataFrame(inter, index=Xm.columns, columns=Xm.columns)
	210	+ paires = inter_df.where(np.triu(np.ones_like(inter, bool), 1)).stack().sort_values(ascending=False)
	211	+ return {"valeurs": sv, "base": base, "sujet": pd.Series(svs, index=Xm.columns), "importance": pd.Series(np.abs(sv).mean(0), index=Xm.columns).sort_values(ascending=False),
	212	+ "interactions": paires, "efficacite_max_err": float(np.abs(sv.sum(1) + base - gbm.predict(Xm)).max()),
	213	+ "efficacite_sujet_err": float(abs(svs.sum() + base - gbm.predict(xs)[0]))}
	214	+
	215	+def importance_permutation(Xh: pd.DataFrame, Xm: pd.DataFrame, y: np.ndarray, plis: list[np.ndarray], repetitions: int = 10, graine: int = GRAINE) -> pd.DataFrame:
	216	+ """Importance par permutation (Breiman 2001) mesurée HORS échantillon : pour chaque pli, la variable j est permutée dans les
	217	+ ventes testées et l'accroissement de l'erreur quadratique (log) est moyenné. Boosting et forêt."""
	218	+ rng = np.random.default_rng(graine); n = len(y); res = {}
	219	+ for nom, ctor in (("Gradient boosting", _boosting), ("Forêt aléatoire", _foret)):
	220	+ gain = np.zeros(Xm.shape[1]); base_err = 0.0
	221	+ for f in plis:
	222	+ tr = np.setdiff1d(np.arange(n), f); m = ctor(); m.fit(Xm.iloc[tr], y[tr]); Xt = Xm.iloc[f].to_numpy(); e0 = np.mean((m.predict(Xt) - y[f]) ** 2); base_err += e0 * len(f) / n
	223	+ for j in range(Xm.shape[1]):
	224	+ acc = 0.0
	225	+ for _ in range(repetitions):
	226	+ Xp = Xt.copy(); Xp[:, j] = rng.permutation(Xp[:, j]); acc += np.mean((m.predict(Xp) - y[f]) ** 2)
	227	+ gain[j] += (acc / repetitions - e0) * len(f) / n
	228	+ res[nom] = pd.Series(gain, index=Xm.columns); res[nom + " (EQM de base)"] = base_err
	229	+ return res
	230	+
	231	+def dependance_partielle(modele, X: pd.DataFrame, variable: str, grille: np.ndarray) -> np.ndarray:
	232	+ """PDP (Friedman 2001) : moyenne des prédictions lorsque la variable est fixée à chaque valeur de la grille, les autres variables
	233	+ gardant leurs valeurs observées. Retourne la courbe en logarithme du prix."""
	234	+ out = np.empty(len(grille))
	235	+ for i, g in enumerate(grille):
	236	+ Xg = X.copy(); Xg[variable] = g; out[i] = float(np.mean(modele.predict(Xg)))
	237	+ return out
	238	+
	239	+def pdp_vraie(v: pd.DataFrame, variable: str, grille: np.ndarray, scenario: bool = False) -> np.ndarray:
	240	+ """Dépendance partielle de la LOI DE GÉNÉRATION (connue) : utile pour juger les modèles sur données synthétiques."""
	241	+ a0 = v["age"].to_numpy(float); S0 = v["superficie_habitable_m2"].to_numpy(float); q = v["qualite_1a5"].to_numpy(float); out = np.empty(len(grille))
	242	+ for i, g in enumerate(grille):
	243	+ a = np.full_like(a0, g) if variable == "age" else a0; S = np.full_like(S0, g) if variable == "superficie" else S0
	244	+ f = VRAIS["ln_S"] * np.log(S / 140.0) + VRAIS["age"] * a + VRAIS["age2"] * a ** 2
	245	+ if scenario:
	246	+ f = f + SCENARIO["prime_qualite_age"] * ((q >= SCENARIO["qualite_min"]) & (a >= SCENARIO["age_min_prime"])) + SCENARIO["mur"] * ((a >= SCENARIO["age_mur"][0]) & (a <= SCENARIO["age_mur"][1]))
	247	+ out[i] = float(np.mean(f))
	248	+ return out
	249	+
	250	+def courbes_apprentissage(Xm: pd.DataFrame, y: np.ndarray, plis: list[np.ndarray], profondeurs=(1, 2, 4)) -> dict:
	251	+ """Erreur quadratique (log) d'entraînement et de validation croisée selon le nombre d'arbres du boosting (staged_predict), pour
	252	+ plusieurs profondeurs : surapprentissage et arrêt précoce."""
	253	+ n = len(y); out = {}
	254	+ for d in profondeurs:
	255	+ tr_err = np.zeros(1200); te_err = np.zeros(1200)
	256	+ for f in plis:
	257	+ tr = np.setdiff1d(np.arange(n), f); m = _boosting(max_depth=d).fit(Xm.iloc[tr], y[tr])
	258	+ for k, (pt, pv) in enumerate(zip(m.staged_predict(Xm.iloc[tr]), m.staged_predict(Xm.iloc[f]))):
	259	+ tr_err[k] += np.mean((pt - y[tr]) ** 2) / len(plis); te_err[k] += np.mean((pv - y[f]) ** 2) * len(f) / n
	260	+ out[d] = {"entrainement": np.sqrt(tr_err), "validation": np.sqrt(te_err), "optimum": int(np.argmin(te_err)) + 1, "rmse_min": float(np.sqrt(te_err.min()))}
	261	+ return out
	262	+
	263	+def sensibilite_hyperparametres(Xm: pd.DataFrame, y: np.ndarray, plis: list[np.ndarray]) -> pd.DataFrame:
	264	+ """Grille d'hyperparamètres du boosting (profondeur, taux d'apprentissage, nombre d'arbres, taille minimale des feuilles) et de la
	265	+ forêt (feuilles, variables par division) : RMSE log et COD en validation croisée aléatoire."""
	266	+ lignes = []
	267	+ grille_b = [dict(max_depth=1, n_estimators=1200), dict(max_depth=2, n_estimators=300), dict(max_depth=2, n_estimators=1200), dict(max_depth=2, n_estimators=3000),
	268	+ dict(max_depth=3, n_estimators=1200), dict(max_depth=4, n_estimators=1200), dict(max_depth=2, n_estimators=1200, learning_rate=0.1),
	269	+ dict(max_depth=2, n_estimators=1200, min_samples_leaf=1), dict(max_depth=2, n_estimators=1200, subsample=1.0)]
	270	+ for k in grille_b:
	271	+ p = valider(lambda: _boosting(**k), "ml", None, Xm, y, plis); m = metriques(p, y, facteur_duan(p, y))
	272	+ lignes.append({"modele": "Gradient boosting", "profondeur": k.get("max_depth", 2), "arbres": k.get("n_estimators", 1200), "nu": k.get("learning_rate", 0.03),
	273	+ "feuille_min": k.get("min_samples_leaf", 5), "sous_ech": k.get("subsample", 0.8), "variables": "toutes", **m})
	274	+ for k in [dict(min_samples_leaf=1, max_features=1.0), dict(min_samples_leaf=2, max_features=0.5), dict(min_samples_leaf=5, max_features=0.33), dict(min_samples_leaf=2, max_features=0.5, n_estimators=100)]:
	275	+ p = valider(lambda: _foret(**k), "ml", None, Xm, y, plis); m = metriques(p, y, facteur_duan(p, y))
	276	+ lignes.append({"modele": "Forêt aléatoire", "profondeur": "max", "arbres": k.get("n_estimators", N_ARBRES_FORET), "nu": np.nan, "feuille_min": k["min_samples_leaf"],
	277	+ "sous_ech": "bootstrap", "variables": f"{int(round(k['max_features'] * 100))} %", **m})
	278	+ return pd.DataFrame(lignes)
	279	+
	280	+# ----------------------------------------------------------------------------- sujet
	281	+def prediction_sujet(modeles_ajustes: dict, duan: dict, Xh_s: pd.DataFrame, Xm_s: pd.DataFrame, ref: dict, s: dict) -> pd.DataFrame:
	282	+ """Prédiction du sujet par chaque modèle estimé sur les 600 ventes : logarithme, niveau avec correction de retransformation
	283	+ (facteur hors échantillon ; pour l'hédonique, traite.cas1.predire donne la valeur de référence 637 953 $)."""
	284	+ rows = []
	285	+ for nom, m in modeles_ajustes.items():
	286	+ lp = float(m.predict(Xh_s if nom.startswith("Hédonique") else Xm_s)[0]); F = duan[nom]
	287	+ rows.append({"modele": nom, "ln_pred": lp, "naif": math.exp(lp), "facteur": F, "prediction": math.exp(lp) * F})
	288	+ df = pd.DataFrame(rows).set_index("modele")
	289	+ df.loc["Hédonique (référence)", "prediction_cas1"] = float(cas1.predire(ref, cas1.sujet_df(s))[0])
	290	+ return df
	291	+
	292	+# ----------------------------------------------------------------------------- exercices
	293	+def exercices(v: pd.DataFrame, comp: dict, secteurs: dict, scen: dict, sh: dict) -> dict:
	294	+ e = {}
	295	+ # 37.1 : division d'un nœud sur huit ventes (ln P, superficie) — réduction de la somme des carrés
	296	+ sub = v.sort_values("id_vente").iloc[:8]; lnP = np.log(sub["prix_vente"].to_numpy(float)); S = sub["superficie_habitable_m2"].to_numpy(float)
	297	+ ordre = np.argsort(S); S, lnP = S[ordre], lnP[ordre]; sse0 = float(np.sum((lnP - lnP.mean()) ** 2)); best = None; essais = []
	298	+ for k in range(1, 8):
	299	+ g, d = lnP[:k], lnP[k:]; sse = float(np.sum((g - g.mean()) ** 2) + np.sum((d - d.mean()) ** 2)); seuil = float((S[k - 1] + S[k]) / 2)
	300	+ essais.append({"seuil": seuil, "sse": sse, "gain": sse0 - sse})
	301	+ if best is None or sse < best["sse"]: best = {"seuil": seuil, "sse": sse, "gain": sse0 - sse, "k": k, "c_gauche": float(g.mean()), "c_droite": float(d.mean())}
	302	+ e["1"] = {"ids": list(sub.sort_values("superficie_habitable_m2")["id_vente"]), "S": S.tolist(), "lnP": lnP.tolist(), "prix": sub.sort_values("superficie_habitable_m2")["prix_vente"].tolist(),
	303	+ "sse0": sse0, "essais": essais, "best": best}
	304	+ # 37.2 : variance d'une moyenne d'arbres corrélés
	305	+ s2, rho = 0.12 ** 2, 0.30; e["2"] = {"sigma": 0.12, "rho": rho, **{f"B{B}": rho * s2 + (1 - rho) * s2 / B for B in (1, 10, 100, 500)}, "limite": rho * s2, "ecart_type_limite": math.sqrt(rho * s2),
	306	+ "B_pour_1pct": (1 - rho) / 0.01 / rho if rho else np.inf}
	307	+ # 37.3 : boosting à la main
	308	+ nu = 0.1; F = 12.80; yv = 13.00; traj = [F]
	309	+ for _ in range(3): F = F + nu * (yv - F); traj.append(F)
	310	+ e["3"] = {"nu": nu, "F": traj, "residus": [yv - t for t in traj], "m_1pct": math.log(0.01) / math.log(1 - nu), "F_nu1": yv}
	311	+ # 37.4 : valeurs de Shapley d'un jeu à trois joueurs (S, Q, L) ; v(S) en milliers de dollars au-dessus de la base
	312	+ vS = {(): 0.0, ("S",): 60.0, ("Q",): 30.0, ("L",): 20.0, ("S", "Q"): 100.0, ("S", "L"): 85.0, ("Q", "L"): 45.0, ("S", "Q", "L"): 130.0}
	313	+ import itertools
	314	+ joueurs = ["S", "Q", "L"]; phi = {}
	315	+ for j in joueurs:
	316	+ tot = 0.0
	317	+ for perm in itertools.permutations(joueurs):
	318	+ avant = tuple(sorted(perm[:perm.index(j)])); apres = tuple(sorted(avant + (j,)))
	319	+ tot += vS[apres] - vS[avant]
	320	+ phi[j] = tot / 6
	321	+ e["4"] = {"v": {"+".join(k) if k else "∅": val for k, val in vS.items()}, "phi": phi, "somme": sum(phi.values()), "total": vS[("S", "Q", "L")]}
	322	+ # 37.5 : dix prédictions d'un AVM fictif
	323	+ prix = np.array([312000, 358500, 401000, 445000, 489500, 532000, 598000, 655000, 742000, 895000], float)
	324	+ pred = np.array([335000, 349000, 388000, 468000, 471000, 561000, 575000, 702000, 718000, 839000], float)
	325	+ r = h.statistiques_ratio(pred, prix); el = np.log(pred / prix)
	326	+ e["5"] = {"prix": prix.tolist(), "pred": pred.tolist(), "ratios": (pred / prix).tolist(), "rmse_log": float(np.sqrt(np.mean(el ** 2))), "mdape": float(np.median(np.abs(pred / prix - 1)) * 100),
	327	+ "mape": float(np.mean(np.abs(pred / prix - 1)) * 100), **{k: float(val) for k, val in r.items()}}
	328	+ # 37.6 : validation par secteur vs aléatoire (chiffres du script)
	329	+ e["6"] = {nom: {"aleatoire": comp["metriques"][nom]["aleatoire"]["COD"], "spatial": comp["metriques"][nom]["spatial"]["COD"], "secteur": secteurs["metriques"][nom]["COD"]}
	330	+ for nom in ("Hédonique (référence)", "Forêt aléatoire", "Gradient boosting")}
	331	+ e["6"]["biais_secteur"] = secteurs["biais_secteur"]
	332	+ # 37.7 : scénario — respécification
	333	+ e["7"] = {nom: scen["metriques"][nom]["COD"] for nom in scen["metriques"]}; e["7"]["coef"] = scen["coef_respecifie"]; e["7"]["interactions_top"] = sh["scenario_interactions"].head(3).to_dict()
	334	+ return e
	335	+
	336	+# ----------------------------------------------------------------------------- orchestration
	337	+def main(verbose: bool = True, complet: bool = True) -> dict:
	338	+ d = cas1.charger(); v = d["ventes"].reset_index(drop=True); s = d["sujet"]
	339	+ ref = cas1.hedonique(v); Xh = cas1.matrice(v); Xm = matrice_ml(v); y = np.log(v["prix_vente"].astype(float).to_numpy())
	340	+ sd = cas1.sujet_df(s); Xh_s = cas1.matrice(sd)[Xh.columns]; Xm_s = matrice_ml(sd, s)
	341	+ comp = comparaison(v, Xh, Xm, y)
	342	+ det = boosting_detendance(v, Xh, Xm, y)
	343	+ sect = validation_secteurs(v, Xh, Xm, y)
	344	+ scen = scenario_nonlineaire(v, Xh, Xm, y)
	345	+ ajustes = {nom: (ctor().fit(Xh if genre == "hedo" else Xm, y)) for nom, (genre, ctor) in MODELES.items()}
	346	+ sujet = prediction_sujet(ajustes, comp["duan"], Xh_s, Xm_s, ref, s)
	347	+ sh = shap_boosting(ajustes["Gradient boosting"], Xm, Xm_s)
	348	+ sh_sc = shap_boosting(scen["gbm_scenario"], Xm, Xm_s); sh["scenario_interactions"] = sh_sc["interactions"]; sh["scenario_importance"] = sh_sc["importance"]
	349	+ pa = comp["plis"]["aleatoire"]
	350	+ perm = importance_permutation(Xh, Xm, y, pa) if complet else None
	351	+ grille_S = np.linspace(80, 260, 37); grille_A = np.linspace(1, 95, 48)
	352	+ pdp = {"superficie": {"grille": grille_S, "Gradient boosting": dependance_partielle(ajustes["Gradient boosting"], Xm, "superficie", grille_S),
	353	+ "Forêt aléatoire": dependance_partielle(ajustes["Forêt aléatoire"], Xm, "superficie", grille_S),
	354	+ "Hédonique (référence)": dependance_partielle(ajustes["Hédonique (référence)"], Xh, "ln_S", np.log(grille_S / 140.0)), "vraie": pdp_vraie(v, "superficie", grille_S)},
	355	+ "age": {"grille": grille_A, "Gradient boosting": dependance_partielle(ajustes["Gradient boosting"], Xm, "age", grille_A),
	356	+ "Forêt aléatoire": dependance_partielle(ajustes["Forêt aléatoire"], Xm, "age", grille_A), "vraie": pdp_vraie(v, "age", grille_A)}}
	357	+ # hédonique : l'âge entre par deux colonnes (age, age2_100) → PDP par substitution conjointe
	358	+ def pdp_hedo_age(modele, Xh_, grille):
	359	+ out = np.empty(len(grille))
	360	+ for i, g in enumerate(grille):
	361	+ Xg = Xh_.copy(); Xg["age"] = g; Xg["age2_100"] = g ** 2 / 100.0; out[i] = float(np.mean(modele.predict(Xg)))
	362	+ return out
	363	+ pdp["age"]["Hédonique (référence)"] = pdp_hedo_age(ajustes["Hédonique (référence)"], Xh, grille_A)
	364	+ pdp["age_scenario"] = {"grille": grille_A, "Gradient boosting": dependance_partielle(scen["gbm_scenario"], Xm, "age", grille_A),
	365	+ "Hédonique (référence)": pdp_hedo_age(scen["hedo_scenario"], Xh, grille_A), "vraie": pdp_vraie(v, "age", grille_A, scenario=True)}
	366	+ courbes = courbes_apprentissage(Xm, y, pa) if complet else None
	367	+ sens = sensibilite_hyperparametres(Xm, y, pa) if complet else None
	368	+ ex = exercices(v, comp, sect, scen, sh)
	369	+ out = {"ventes": v, "sujet": s, "modele_ref": ref, "Xh": Xh, "Xm": Xm, "y": y, "Xh_s": Xh_s, "Xm_s": Xm_s, "comparaison": comp, "detendance": det, "secteurs": sect,
	370	+ "scenario": scen, "ajustes": ajustes, "sujet_pred": sujet, "shap": sh, "permutation": perm, "pdp": pdp, "courbes": courbes, "sensibilite": sens, "exercices": ex}
	371	+ if verbose:
	372	+ pd.set_option("display.width", 220); pd.set_option("display.max_columns", 30)
	373	+ print("== Modèle hédonique de référence : R² %.3f, σ %.4f, n %d, k %d" % (ref["r2"], ref["sigma"], ref["n"], ref["k"]))
	374	+ for schema in ("aleatoire", "spatial", "temporel"):
	375	+ print(f"== Validation {schema} :")
	376	+ for nom, m in comp["metriques"].items():
	377	+ mm = m[schema]; print(f" {nom:24s} n={mm['n']:3d} RMSElog={mm['rmse_log']:.4f} R²hors={mm['r2_hors']:.3f} biais={mm['biais_log']:+.4f} RMSE$={mm['rmse_dollars']:,.0f} MdAPE={mm['mdape']:.2f} COD={mm['COD']:.2f} PRD={mm['PRD']:.3f} PRB={mm['PRB']:.3f} méd={mm['ratio_median']:.3f} Duan={comp['duan'][nom]:.4f}")
	378	+ print(" biais mensuel 2026 :", {k: {int(kk): round(vv, 3) for kk, vv in vals.items()} for k, vals in comp["biais_mensuel_2026"].items()})
	379	+ md = det["metriques_detendance"]; m1 = det["metriques_feuille1"]
	380	+ print(f"== Boosting détendancé (test 2026) : biais={md['biais_log']:+.4f} RMSElog={md['rmse_log']:.4f} COD={md['COD']:.2f} MdAPE={md['mdape']:.2f} ; β_mois entraînement {det['beta_mois_entrainement']:.5f} ; n(mois>17)={det['n_entrainement_apres_mois_17']}")
	381	+ print(f" Boosting feuilles ≥ 1 (test 2026) : biais={m1['biais_log']:+.4f} RMSElog={m1['rmse_log']:.4f} COD={m1['COD']:.2f}")
	382	+ print("== Un secteur exclu :")
	383	+ for nom, m in sect["metriques"].items():
	384	+ print(f" {nom:24s} RMSElog={m['rmse_log']:.4f} MdAPE={m['mdape']:.2f} COD={m['COD']:.2f} biais/secteur=", {k: round(vv, 3) for k, vv in sect["biais_secteur"][nom].items()})
	385	+ print(f"== Scénario non linéaire (n prime={scen['n_prime']}, n mur={scen['n_mur']}) :")
	386	+ for nom, m in scen["metriques"].items():
	387	+ print(f" {nom:24s} RMSElog={m['rmse_log']:.4f} MdAPE={m['mdape']:.2f} COD={m['COD']:.2f} PRB={m['PRB']:.3f}")
	388	+ print(" coefficients respécifiés :", {k: round(vv, 4) for k, vv in scen["coef_respecifie"].items()})
	389	+ print("== Sujet :"); print(sujet.round(4).to_string())
	390	+ print(f"== SHAP (boosting) : base {sh['base']:.4f} (exp {math.exp(sh['base']):,.0f} $) ; efficacité max err {sh['efficacite_max_err']:.2e}")
	391	+ print(" importance moyenne |φ| :", sh["importance"].round(4).to_dict())
	392	+ print(" sujet φ :", sh["sujet"].sort_values().round(4).to_dict())
	393	+ print(" interactions (base) :", sh["interactions"].head(4).round(4).to_dict()); print(" interactions (scénario) :", sh["scenario_interactions"].head(4).round(4).to_dict())
	394	+ if perm is not None:
	395	+ print("== Importance par permutation (ΔEQM log, hors échantillon) :")
	396	+ for nom in ("Gradient boosting", "Forêt aléatoire"): print(f" {nom} (EQM base {perm[nom + ' (EQM de base)']:.5f}) :", perm[nom].sort_values(ascending=False).round(5).head(8).to_dict())
	397	+ if courbes is not None:
	398	+ print("== Courbes d'apprentissage :", {d_: (c["optimum"], round(c["rmse_min"], 4), round(float(c["entrainement"][-1]), 4)) for d_, c in courbes.items()})
	399	+ if sens is not None:
	400	+ print("== Sensibilité aux hyperparamètres :"); print(sens[["modele", "profondeur", "arbres", "nu", "feuille_min", "sous_ech", "variables", "rmse_log", "mdape", "COD", "PRB"]].round(4).to_string(index=False))
	401	+ print("== Exercices :")
	402	+ print(" 37.1", {k: (round(vv, 4) if isinstance(vv, float) else vv) for k, vv in ex["1"]["best"].items()}, "SSE0", round(ex["1"]["sse0"], 4)); print(" ventes", ex["1"]["ids"], ex["1"]["S"], [round(x, 4) for x in ex["1"]["lnP"]])
	403	+ print(" 37.2", {k: (round(vv, 6) if isinstance(vv, float) else vv) for k, vv in ex["2"].items()})
	404	+ print(" 37.3", {k: (np.round(vv, 4).tolist() if isinstance(vv, list) else round(vv, 4)) for k, vv in ex["3"].items()})
	405	+ print(" 37.4", {k: (round(vv, 3) if isinstance(vv, float) else vv) for k, vv in ex["4"].items()})
	406	+ print(" 37.5", {k: (round(vv, 4) if isinstance(vv, float) else vv) for k, vv in ex["5"].items() if k not in ("prix", "pred", "ratios")}, [round(r_, 3) for r_ in ex["5"]["ratios"]])
	407	+ print(" 37.6", ex["6"]); print(" 37.7", ex["7"])
	408	+ return out
	409	+
	410	+if __name__ == "__main__":
	411	+ main()

modified analysis/python/cases/ch39_bayes.py +60 −29

		@@ -35,18 +35,25 @@ LN_MEDIANE_HULL = math.log(475_000) # médiane APCIQ T3 2026 secteur H
35	35	# 1. Régression linéaire bayésienne conjuguée (normale–inverse-gamma)
36	36	# =============================================================================================================================
37	37
	38	+def prior_reference(colonnes) -> dict:
	39	+ """A priori de référence (impropre, Jeffreys) : p(β, σ²) ∝ 1/σ². Cas limite de la NIG avec V0⁻¹ = 0, a0 = −k/2, b0 = 0 : la
	40	+ postérieure est NIG(β̂, (X'X)⁻¹, (n−k)/2, SCR/2) et la prédictive coïncide avec l'intervalle de prédiction fréquentiste classique."""
	41	+ k = len(colonnes)
	42	+ return {"m0": np.zeros(k), "V0inv": np.zeros((k, k)), "a0": -k / 2, "b0": 0.0, "nom": "référence (impropre)"}
	43	+
38	44	def prior_faible(colonnes, v0: float = 1e4, a0: float = 1e-3, b0: float = 1e-3) -> dict:
39		− """A priori faiblement informatif : β | σ² ~ N(m0, σ² V0), σ² ~ IG(a0, b0). m0 = ln(médiane APCIQ de Hull) pour la constante,
40		− 0 ailleurs ; V0 = v0·I (écart-type a priori ≈ 100 σ ≈ 8 en logarithme : pratiquement plat)."""
	45	+ """A priori faiblement informatif PROPRE : β | σ² ~ N(m0, σ² V0), σ² ~ IG(a0, b0). m0 = ln(médiane APCIQ de Hull) pour la constante,
	46	+ 0 ailleurs ; V0 = v0·I (écart-type a priori ≈ 100 σ ≈ 8 en logarithme : pratiquement plat). Prix d'un a priori propre et vague :
	47	+ a_n = a0 + n/2 au lieu de (n−k)/2, donc E[σ²] ≈ SCR/(n−2) légèrement inférieur à s² = SCR/(n−k)."""
41	48	k = len(colonnes); m0 = np.zeros(k)
42	49	if "const" in colonnes: m0[list(colonnes).index("const")] = LN_MEDIANE_HULL
43		− return {"m0": m0, "V0": np.eye(k) * v0, "a0": a0, "b0": b0, "nom": "faiblement informatif"}
	50	+ return {"m0": m0, "V0inv": np.eye(k) / v0, "a0": a0, "b0": b0, "nom": "faiblement informatif (propre)"}
44	51
45	52	def nig_posterieure(X: np.ndarray, y: np.ndarray, prior: dict) -> dict:
46		− """Postérieure NIG : V_n = (V0⁻¹ + X'X)⁻¹ ; m_n = V_n (V0⁻¹ m0 + X'y) ; a_n = a0 + n/2 ;
	53	+ """Postérieure NIG : V_n = (V0⁻¹ + X'X)⁻¹ ; m_n = V_n (V0⁻¹ m0 + X'y) ; a_n = a0 + n/2 ; (prior : m0, V0inv, a0, b0)
47	54	b_n = b0 + ½ (y'y + m0' V0⁻¹ m0 − m_n' V_n⁻¹ m_n). Marginale de β : Student-t à 2a_n ddl, échelle² = (b_n/a_n) diag(V_n)."""
48	55	X = np.asarray(X, float); y = np.asarray(y, float); n, k = X.shape
49		− V0inv = np.linalg.inv(prior["V0"]); m0 = prior["m0"]
	56	+ V0inv = prior["V0inv"]; m0 = prior["m0"]
50	57	Vn_inv = V0inv + X.T @ X; Vn = np.linalg.inv(Vn_inv)
51	58	mn = Vn @ (V0inv @ m0 + X.T @ y)
52	59	an = prior["a0"] + n / 2
		@@ -88,7 +95,7 @@ def prior_marche_voisin(post_voisin: dict, colonnes_voisin, colonnes_cible, kapp
88	95	m0[pos] = post_voisin["m_n"][idx]
89	96	V0[np.ix_(pos, pos)] = kappa * post_voisin["V_n"][np.ix_(idx, idx)]
90	97	ic = cc.index("const"); m0[ic] = LN_MEDIANE_HULL; V0[ic, ic] = v0_const
91		− return {"m0": m0, "V0": V0, "a0": post_voisin["a_n"], "b0": post_voisin["b_n"], "nom": f"marché voisin, κ = {kappa:g}"}
	98	+ return {"m0": m0, "V0inv": np.linalg.inv(V0), "a0": post_voisin["a_n"], "b0": post_voisin["b_n"], "nom": f"marché voisin, κ = {kappa:g}"}
92	99
93	100	def loo_exact_nig(X: np.ndarray, y: np.ndarray, prior: dict) -> dict:
94	101	"""LOO exact du modèle conjugué : Σ_i log p(y_i | y_{−i}) par n réestimations en forme fermée (prédictive Student-t)."""
		@@ -115,17 +122,21 @@ def matrice_hier(df: pd.DataFrame, moyennes: pd.Series | None = None):
115	122
116	123	def gibbs_hierarchique(y: np.ndarray, X: np.ndarray, z: np.ndarray, J: int = 5, n_iter: int = 3000, brulage: int = 1000,
117	124	chaines: int = 4, graine: int = GRAINE, prior_tau: str = "uniforme", tau_fixe: float | None = None,
118		− ig_a: float = 1e-3, ig_b: float = 1e-3) -> dict:
	125	+ ig_a: float = 1e-3, ig_b: float = 1e-3, bloc: bool = True) -> dict:
119	126	"""ln P_i = x_i'β + α_{s(i)} + ε_i, ε ~ N(0, σ²) ; α_s | α0, τ² ~ N(α0, τ²) ; p(β) ∝ 1, p(α0) ∝ 1, p(σ²) ∝ 1/σ² ;
120		− τ : p(τ) ∝ 1 (uniforme sur τ, BDA3 §5.7 ; « uniforme ») ou τ² ~ IG(ig_a, ig_b) (« ig ») ou τ fixé (« fixe », bayésien empirique).
	127	+ τ : p(τ) ∝ 1 (uniforme sur τ, BDA3 §5.7 ; « uniforme ») ou τ² ~ IG(ig_a, ig_b) (« ig ») ou τ fixé (bayésien empirique).
121	128	Conditionnelles complètes (toutes conjuguées) :
122	129	β | · ~ N((X'X)⁻¹X'(y − Zα), σ²(X'X)⁻¹)
123	130	α_s | · ~ N([n_s r̄_s/σ² + α0/τ²] / [n_s/σ² + 1/τ²], 1/[n_s/σ² + 1/τ²]), r̄_s = moyenne de (y − Xβ) dans le secteur s
124	131	α0 | · ~ N(moyenne des α_s, τ²/J)
125	132	σ² | · ~ IG(n/2, ½ Σ (y − Xβ − Zα)²)
126	133	τ² | · ~ IG((J−1)/2, ½ Σ (α_s − α0)²) [uniforme sur τ] ou IG(a + J/2, b + ½ Σ (α_s − α0)²) [IG]
	134	+ bloc=True : (β, α) sont tirés ENSEMBLE de leur conditionnelle jointe normale (Gibbs par blocs), ce qui supprime la lenteur de
	135	+ mélange due à la colinéarité entre la distance au centre et les niveaux des secteurs périphériques ; bloc=False : un à la fois.
127	136	Retourne les tirages après brûlage, de forme (chaînes, n_iter − brulage, ·)."""
128	137	y = np.asarray(y, float); X = np.asarray(X, float); n, k = X.shape
	138	+ Z = np.zeros((n, J)); Z[np.arange(n), z] = 1.0
	139	+ W = np.hstack([X, Z]); WtW = W.T @ W; Wty = W.T @ y
129	140	XtX_inv = np.linalg.inv(X.T @ X); L_beta = np.linalg.cholesky(XtX_inv)
130	141	ns = np.bincount(z, minlength=J).astype(float)
131	142	garder = n_iter - brulage
		@@ -133,21 +144,31 @@ def gibbs_hierarchique(y: np.ndarray, X: np.ndarray, z: np.ndarray, J: int = 5,
133	144	"sigma2": np.empty((chaines, garder)), "tau2": np.empty((chaines, garder))}
134	145	for c in range(chaines):
135	146	rng = np.random.default_rng([graine, c, 39])
136		− # valeurs initiales dispersées (BDA3 : points de départ surdispersés pour que R̂ ait un sens)
	147	+ # valeurs initiales surdispersées (BDA3 : points de départ dispersés pour que R̂ ait un sens)
137	148	beta = XtX_inv @ X.T @ y + rng.normal(0, 0.05, k)
138	149	alpha = np.full(J, float(np.mean(y))) + rng.normal(0, 0.3, J); alpha0 = float(np.mean(alpha))
139	150	sigma2 = float(np.exp(rng.normal(math.log(0.01), 0.5))); tau2 = float(np.exp(rng.normal(math.log(0.02), 1.0)))
140	151	if tau_fixe is not None: tau2 = tau_fixe ** 2
141	152	for it in range(n_iter):
142		− # β | α, σ²
143		− r = y - alpha[z]
144		− beta = XtX_inv @ X.T @ r + math.sqrt(sigma2) * (L_beta @ rng.standard_normal(k))
145		− # α_s | β, α0, σ², τ²
146		− e = y - X @ beta
147		− somme = np.bincount(z, weights=e, minlength=J)
148		− prec = ns / sigma2 + 1.0 / tau2
149		− moy = (somme / sigma2 + alpha0 / tau2) / prec
150		− alpha = moy + rng.standard_normal(J) / np.sqrt(prec)
	153	+ if bloc:
	154	+ # (β, α) | α0, σ², τ² : précision = W'W/σ² + diag(0_k, 1/τ² 1_J) ; moyenne = Préc⁻¹ (W'y/σ² + (0, α0/τ² 1_J))
	155	+ prec = WtW / sigma2; prec[k:, k:] += np.eye(J) / tau2
	156	+ b = Wty / sigma2; b[k:] += alpha0 / tau2
	157	+ Lp = np.linalg.cholesky(prec)
	158	+ moy = np.linalg.solve(Lp.T, np.linalg.solve(Lp, b))
	159	+ theta = moy + np.linalg.solve(Lp.T, rng.standard_normal(k + J))
	160	+ beta, alpha = theta[:k], theta[k:]
	161	+ e = y - X @ beta
	162	+ else:
	163	+ # β | α, σ²
	164	+ r = y - alpha[z]
	165	+ beta = XtX_inv @ X.T @ r + math.sqrt(sigma2) * (L_beta @ rng.standard_normal(k))
	166	+ # α_s | β, α0, σ², τ²
	167	+ e = y - X @ beta
	168	+ somme = np.bincount(z, weights=e, minlength=J)
	169	+ prec_a = ns / sigma2 + 1.0 / tau2
	170	+ moy_a = (somme / sigma2 + alpha0 / tau2) / prec_a
	171	+ alpha = moy_a + rng.standard_normal(J) / np.sqrt(prec_a)
151	172	# α0 | α, τ²
152	173	alpha0 = float(alpha.mean() + rng.standard_normal() * math.sqrt(tau2 / J))
153	174	# σ² | β, α
		@@ -164,7 +185,7 @@ def gibbs_hierarchique(y: np.ndarray, X: np.ndarray, z: np.ndarray, J: int = 5,
164	185	out["beta"][c, j] = beta; out["alpha"][c, j] = alpha; out["alpha0"][c, j] = alpha0
165	186	out["sigma2"][c, j] = sigma2; out["tau2"][c, j] = tau2
166	187	out["meta"] = {"n": n, "k": k, "J": J, "ns": ns, "n_iter": n_iter, "brulage": brulage, "chaines": chaines, "prior_tau": prior_tau,
167		− "tau_fixe": tau_fixe}
	188	+ "tau_fixe": tau_fixe, "bloc": bloc}
168	189	return out
169	190
170	191	def facteur_retrecissement(n_s: float, sigma2: float, tau2: float) -> float:
		@@ -284,12 +305,17 @@ def sous_echantillons(v: pd.DataFrame, s: dict) -> dict:
284	305	return {"tout": v, "hull": v[v["secteur"] == "Hull"], "voisin": v[v["secteur"] != "Hull"], "micro": micro}
285	306
286	307	def colonnes_actives(X: pd.DataFrame) -> list:
287		− return [c for c in X.columns if c == "const" or X[c].std() > 0]
	308	+ """Colonnes non constantes ; si le secteur de référence (Hull) est absent du sous-échantillon, les indicatrices restantes somment
	309	+ à 1 et l'une d'elles (la première) est retirée pour éviter la colinéarité avec la constante."""
	310	+ cols = [c for c in X.columns if c == "const" or X[c].std() > 0]
	311	+ sect = [c for c in cols if c.startswith("sect_")]
	312	+ if sect and np.allclose(X[sect].sum(axis=1), 1.0): cols.remove(sect[0])
	313	+ return cols
288	314
289	315	def ajuster_nig(df: pd.DataFrame, prior_fn=None, colonnes=None) -> dict:
290	316	X = cas1.matrice(df); cols = colonnes or colonnes_actives(X); X = X[cols]
291	317	y = np.log(df["prix_vente"].astype(float)).values
292		− prior = prior_fn(cols) if prior_fn else prior_faible(cols)
	318	+ prior = prior_fn(cols) if prior_fn else prior_reference(cols)
293	319	p = nig_posterieure(X.values, y, prior); p["colonnes"] = cols
294	320	return p
295	321
		@@ -304,7 +330,8 @@ def exemple_conjugue(v: pd.DataFrame, s: dict, ref: dict) -> dict:
304	330	"""Étapes 1-2 : postérieure NIG sur 600 ventes (a priori faible) ; petit échantillon (40 ventes à 3 km / 18 mois) avec a priori
305	331	faible vs informatif emprunté aux autres secteurs (κ = 1, 4, 16) ; prédictive du sujet."""
306	332	se = sous_echantillons(v, s)
307		− p600 = ajuster_nig(se["tout"]); cols = p600["colonnes"]
	333	+ p600 = ajuster_nig(se["tout"]); cols = p600["colonnes"] # a priori de référence
	334	+ p600_propre = ajuster_nig(se["tout"], prior_fn=prior_faible, colonnes=cols) # a priori propre vague (comparaison)
308	335	xs = x_sujet(s, cols)
309	336	tab = pd.DataFrame({"mco": ref["coef"]["coef"].reindex(cols), "se_hc1": ref["coef"]["se_hc1"].reindex(cols),
310	337	"post_moy": p600["m_n"], "post_sd": p600["sd"], "ic_bas": p600["ic95_bas"], "ic_haut": p600["ic95_haut"],
		@@ -314,7 +341,7 @@ def exemple_conjugue(v: pd.DataFrame, s: dict, ref: dict) -> dict:
314	341	# cohérence séquentielle : 2023-2024 puis 2025-2026
315	342	anc = v[v["date_vente"] < "2025-01-01"]; rec = v[v["date_vente"] >= "2025-01-01"]
316	343	p_anc = ajuster_nig(anc, colonnes=cols)
317		− prior_seq = {"m0": p_anc["m_n"], "V0": p_anc["V_n"], "a0": p_anc["a_n"], "b0": p_anc["b_n"], "nom": "séquentiel"}
	344	+ prior_seq = {"m0": p_anc["m_n"], "V0inv": np.linalg.inv(p_anc["V_n"]), "a0": p_anc["a_n"], "b0": p_anc["b_n"], "nom": "séquentiel"}
318	345	p_seq = ajuster_nig(rec, prior_fn=lambda c: prior_seq, colonnes=cols)
319	346	# petit échantillon : 40 ventes autour du sujet ; a priori emprunté aux autres secteurs
320	347	p_vois = ajuster_nig(se["voisin"])
		@@ -326,7 +353,7 @@ def exemple_conjugue(v: pd.DataFrame, s: dict, ref: dict) -> dict:
326	353	petits[nom] = {"post": p, "pred": pr, "largeur95_pct": (pr["haut95"] - pr["bas95"]) / pr["mediane"] * 100}
327	354	# Hull seul (186), a priori faible
328	355	p_hull = ajuster_nig(se["hull"]); pred_hull = predictive_t(p_hull, x_sujet(s, p_hull["colonnes"]))
329		− return {"p600": p600, "cols": cols, "xs": xs, "tab600": tab, "pred600": pred600, "p_anc": p_anc, "p_rec_seq": p_seq, "n_anc": len(anc),
	356	+ return {"p600": p600, "p600_propre": p600_propre, "pred600_propre": predictive_t(p600_propre, xs), "cols": cols, "xs": xs, "tab600": tab, "pred600": pred600, "p_anc": p_anc, "p_rec_seq": p_seq, "n_anc": len(anc),
330	357	"n_rec": len(rec), "p_vois": p_vois, "micro": micro, "cols_micro": cols_m, "petits": petits, "p_hull": p_hull,
331	358	"pred_hull": pred_hull, "n_micro": len(micro)}
332	359
		@@ -334,8 +361,10 @@ def exemple_hierarchique(v: pd.DataFrame, s: dict, ref: dict, n_iter: int = 3000
334	361	"""Étapes 3-4 : Gibbs sur les 600 ventes ; effets relatifs à Hull vs effets fixes (ch. 18) ; prédictive du sujet ;
335	362	expérience Masson-Angers réduit à 5, 10, 20 ventes ; sensibilité à l'a priori sur τ."""
336	363	Xc, moy, z = matrice_hier(v); y = np.log(v["prix_vente"].astype(float)).values; cols = list(Xc.columns)
337		− g = gibbs_hierarchique(y, Xc.values, z, n_iter=n_iter, brulage=brulage)
338		− res = resume_gibbs(g, cols); rel = effets_relatifs(g)
	364	+ g = gibbs_hierarchique(y, Xc.values, z, n_iter=n_iter, brulage=brulage) # Gibbs par blocs (résultats du chapitre)
	365	+ g_naif = gibbs_hierarchique(y, Xc.values, z, n_iter=n_iter, brulage=brulage, bloc=False) # un paramètre à la fois (diagnostics)
	366	+ res = resume_gibbs(g, cols); rel = effets_relatifs(g); res_naif = resume_gibbs(g_naif, cols)
	367	+ diag = pd.DataFrame({"rhat_bloc": res["rhat"], "ess_bloc": res["ess"], "rhat_naif": res_naif["rhat"], "ess_naif": res_naif["ess"]})
339	368	fixes = ref["coef"].loc[[f"sect_{x}" for x in SECTEURS[1:]], ["coef", "se_hc1"]]; fixes.index = SECTEURS[1:]
340	369	rel["fixe"] = fixes["coef"]; rel["fixe_se"] = fixes["se_hc1"]; rel["vrai"] = [VRAIS[f"sect_{x}"] for x in rel.index]
341	370	# prédictive du sujet (Hull) : ln P* = (x_s − x̄)'β + α_Hull + ε
		@@ -396,7 +425,7 @@ def exemple_hierarchique(v: pd.DataFrame, s: dict, ref: dict, n_iter: int = 3000
396	425	"ma_600": r600["moyenne"], "ma_600_sd": r600["sd"], "tau_med_5": float(np.median(np.sqrt(g5["tau2"]))),
397	426	"ma_5": r5["moyenne"], "ma_5_sd": r5["sd"], "ma_5_bas": r5["q025"], "ma_5_haut": r5["q975"]})
398	427	sensibilite = pd.DataFrame(sens_rows).set_index("prior")
399		− return {"g": g, "cols": cols, "moyennes": moy, "z": z, "y": y, "Xc": Xc, "resume": res, "relatifs": rel, "pred": pred, "B": B,
	428	+ return {"g": g, "g_naif": g_naif, "diag": diag, "cols": cols, "moyennes": moy, "z": z, "y": y, "Xc": Xc, "resume": res, "relatifs": rel, "pred": pred, "B": B,
400	429	"sigma2_med": s2m, "tau2_med": t2m, "experience": experience, "dist_ma": dist_ma, "sensibilite": sensibilite, "tau_eb": tau_eb,
401	430	"perm_ma": perm}
402	431
		@@ -405,7 +434,7 @@ def criteres(v: pd.DataFrame, hier: dict, ex: dict) -> pd.DataFrame:
405	434	y = hier["y"]
406	435	# pooling complet : matrice de référence sans indicatrices de secteur
407	436	X = cas1.matrice(v); cols_pool = [c for c in X.columns if not c.startswith("sect_")]
408		− p_pool = nig_posterieure(X[cols_pool].values, y, prior_faible(cols_pool)); p_pool["colonnes"] = cols_pool
	437	+ p_pool = nig_posterieure(X[cols_pool].values, y, prior_reference(cols_pool)); p_pool["colonnes"] = cols_pool
409	438	p_fixe = ex["p600"]
410	439	lignes = []
411	440	for nom, lp, loo in [("pooling complet (sans secteur)", log_vrais_pointwise_nig(p_pool), loo_exact_nig(p_pool["X"], y, p_pool["prior"])),
		@@ -487,6 +516,7 @@ def imprimer(r: dict) -> None:
487	516	print(t.round(4).to_string())
488	517	print(f"couverture des vrais paramètres (IC 95 %) : {int(t['couvre'].sum())}/{int(t['vrai'].notna().sum())}")
489	518	p = ex["p600"]; print(f"a_n = {p['a_n']:.1f} ; b_n = {p['b_n']:.4f} ; ν = {p['nu']:.0f} ; E[σ] ≈ {p['sigma_moy']:.4f} ; MCO σ = {r['modele']['sigma']:.4f}")
	519	+ pp = ex["p600_propre"]; print(f"a priori propre vague : a_n = {pp['a_n']:.3f} ; ν = {pp['nu']:.0f} ; E[σ] ≈ {pp['sigma_moy']:.4f} ; max|m_n − β̂| = {np.abs(pp['m_n'] - r['modele']['beta']).max():.2e} ; sujet 95 % [{_d(ex['pred600_propre']['bas95'])} ; {_d(ex['pred600_propre']['haut95'])}]")
490	520	pr = ex["pred600"]
491	521	print(f"prédictive du sujet : loc {pr['loc']:.4f}, échelle {pr['echelle']:.5f}, médiane {_d(pr['mediane'])}, part coefficients {pr['part_coefficients']*100:.2f} %")
492	522	print(f" crédibilité 90 % : [{_d(pr['bas90'])} ; {_d(pr['haut90'])}] 95 % : [{_d(pr['bas95'])} ; {_d(pr['haut95'])}] (ch. 18 : [{_d(IC18['bas95'])} ; {_d(IC18['haut95'])}])")
		@@ -497,8 +527,9 @@ def imprimer(r: dict) -> None:
497	527	print(f" {nom:9s}: ln_S {po['m_n'][i]:.3f} ({po['sd'][i]:.3f}) ; qualité {po['m_n'][iq]:.3f} ({po['sd'][iq]:.3f}) ; garage {po['m_n'][ig]:.3f} ({po['sd'][ig]:.3f}) ; "
498	528	f"E[σ] {po['sigma_moy']:.4f} ; sujet médiane {_d(pe['pred']['mediane'])} [{_d(pe['pred']['bas95'])} ; {_d(pe['pred']['haut95'])}] largeur {pe['largeur95_pct']:.1f} %")
499	529	ph = ex["pred_hull"]; print(f" Hull seul ({ex['p_hull']['n']}) : médiane {_d(ph['mediane'])} [{_d(ph['bas95'])} ; {_d(ph['haut95'])}]")
500		− print("\n=== modèle hiérarchique (Gibbs) : résumé ===")
	530	+ print("\n=== modèle hiérarchique (Gibbs par blocs) : résumé ===")
501	531	print(hi["resume"].round(4).to_string())
	532	+ print("\ndiagnostics : Gibbs par blocs vs un paramètre à la fois"); print(hi["diag"].round(3).to_string())
502	533	print("\neffets relatifs à Hull (hiérarchique vs effets fixes ch. 18 vs vrai) :"); print(hi["relatifs"].round(4).to_string())
503	534	print(f"σ² médian {hi['sigma2_med']:.5f} (σ {math.sqrt(hi['sigma2_med']):.4f}) ; τ² médian {hi['tau2_med']:.5f} (τ {math.sqrt(hi['tau2_med']):.4f})")
504	535	print("facteurs de rétrécissement B_s :", {k: round(v_, 4) for k, v_ in hi["B"].items()})

modified analysis/python/cases/ch40_montecarlo.py +20 −9

		@@ -196,7 +196,7 @@ def comparer_lhs(r: dict, n: int = 1000, k: int = 100) -> pd.DataFrame:
196	196	rows = []
197	197	for s in ["moyenne", "mediane", "q05", "q95", "ecart_type", "p_rcd1_120", "p_V_2M"]:
198	198	rows.append({"statistique": s, "moy_simple": simple[s].mean(), "moy_lhs": lhs[s].mean(), "se_simple": simple[s].std(ddof=1), "se_lhs": lhs[s].std(ddof=1),
199		− "ratio_variance": simple[s].var(ddof=1) / lhs[s].var(ddof=1)})
	199	+ "ratio_variance": (simple[s].var(ddof=1) / lhs[s].var(ddof=1)) if lhs[s].var(ddof=1) > 0 else float("inf")})
200	200	return pd.DataFrame(rows)
201	201
202	202	# ----------------------------------------------------------------------------- 40.5 mesures de risque, vérification analytique, stress tests
		@@ -272,17 +272,17 @@ NOMS_SOBOL = ["inoccupation_pertes", "croissance_loyers_marche", "rotation_annue
272	272	def _completer(x: dict) -> dict:
273	273	x = dict(x); x["taux_terminal"] = np.maximum(x["taux_actualisation"] - x["ecart_terminal"], R_T_PLANCHER); return x
274	274
275		−def sobol(r: dict, n: int = N_MC, graine: int = GRAINE) -> pd.DataFrame:
	275	+def sobol(r: dict, n: int = 2 * N_MC, graine: int = GRAINE) -> pd.DataFrame:
276	276	"""Indices de Sobol de premier ordre S_i et totaux S_Ti par échantillons appariés : deux matrices indépendantes A et B (n tirages
277	277	chacune), puis A_B^(i) = A où la seule colonne i vient de B. S_i ≈ (1/n) Σ f(B)·[f(A_B^(i)) − f(A)] / Var ; S_Ti ≈ (1/2n) Σ [f(A) − f(A_B^(i))]² / Var.
278	278	Les intrants de base (r et delta, non r et R_T) sont indépendants, condition de validité de la décomposition. Le facteur CAPEX est traité en groupe (3 facteurs)."""
279	279	nb = len(r["hyp"]["capex"])
280	280	A = tirer(n, graine, nb); B = tirer(n, graine + 1, nb)
281	281	fA = evaluer(r, A)["valeur"]; fB = evaluer(r, B)["valeur"]
282		− var = np.concatenate([fA, fB]).var(ddof=1); rows = []
	282	+ tous = np.concatenate([fA, fB]); var = tous.var(ddof=1); mu = tous.mean(); rows = []
283	283	for nom in NOMS_SOBOL:
284	284	AB = dict(A); AB[nom] = B[nom]; fAB = evaluer(r, _completer(AB))["valeur"]
285		− S1 = float(np.mean(fB * (fAB - fA)) / var); ST = float(0.5 * np.mean((fA - fAB) ** 2) / var)
	285	+ S1 = float(np.mean((fB - mu) * (fAB - fA)) / var); ST = float(0.5 * np.mean((fA - fAB) ** 2) / var) # centrage : même espérance, variance réduite
286	286	rows.append({"intrant": nom, "S1": S1, "ST": ST})
287	287	d = pd.DataFrame(rows); d.attrs["somme_S1"] = float(d["S1"].sum()); d.attrs["variance"] = float(var); d.attrs["n_evaluations"] = n * (2 + len(NOMS_SOBOL))
288	288	return d
		@@ -309,7 +309,7 @@ def regression_standardisee(mc: dict) -> dict:
309	309
310	310	def simuler_portefeuille(pf: pd.DataFrame, n: int = N_MC, graine: int = GRAINE, sigma_R: float = 0.005, rho_rmr: float = 0.60, rho_prov: float = 0.30,
311	311	croissance: pd.DataFrame | None = None, sigma_g_defaut: float = 0.03, choc_taux_renouvellement: float = 0.02, horizon_renouvellement: int = 2028,
312		− amortissement_mois: int = 300, seuil_rcd: float = 1.20) -> dict:
	312	+ amortissement_mois: int = 480, seuil_rcd: float = 1.20) -> dict:
313	313	"""Simulation à un an de la valeur et de la couverture de la dette d'un portefeuille d'immeubles à revenus avec corrélation régionale.
314	314
315	315	Modèle (hypothèses pédagogiques déclarées au ch. 40, §40.8) : pour l'immeuble i de la RMR k,
		@@ -339,6 +339,8 @@ def simuler_portefeuille(pf: pd.DataFrame, n: int = N_MC, graine: int = GRAINE,
339	339	sd = np.array([fin.Hypotheque(float(d), float(t) + (choc_taux_renouvellement if int(ech) <= horizon_renouvellement else 0.0), amortissement_mois).service_dette_annuel
340	340	for d, t, ech in zip(pf["dette"], pf["taux_dette"], pf["echeance_dette"])])
341	341	rcd = rne / sd
	342	+ sd0 = np.array([fin.Hypotheque(float(d), float(t), amortissement_mois).service_dette_annuel for d, t in zip(pf["dette"], pf["taux_dette"])])
	343	+ rcd_det = pf["rne_2026"].values / sd0
342	344	VP = V.sum(1); V0 = float(pf["valeur_2026"].sum()); VP_det = float((pf["rne_2026"].values * np.exp(g) / pf["taux_global"].values).sum())
343	345	q = lambda a, p: float(np.quantile(a, p))
344	346	stats_ = {"n": n, "N_immeubles": N, "valeur_2026": V0, "valeur_deterministe": VP_det, "mediane": q(VP, 0.5), "moyenne": float(VP.mean()), "ecart_type": float(VP.std(ddof=1)),
		@@ -346,7 +348,10 @@ def simuler_portefeuille(pf: pd.DataFrame, n: int = N_MC, graine: int = GRAINE,
346	348	"somme_ecarts_types": float(V.std(0, ddof=1).sum()), "nb_rcd_sous_seuil_moyen": float((rcd < seuil_rcd).sum(1).mean()), "nb_rcd_sous_1_moyen": float((rcd < 1.0).sum(1).mean()),
347	349	"P_au_moins_un_rcd_sous_1": float(((rcd < 1.0).any(1)).mean()), "P_au_moins_un_rcd_sous_seuil": float(((rcd < seuil_rcd).any(1)).mean()),
348	350	"rcd_portefeuille_mediane": q(rne.sum(1) / sd.sum(), 0.5), "P_rcd_portefeuille_sous_seuil": float((rne.sum(1) / sd.sum() < seuil_rcd).mean()),
349		− "nb_renouvellements_stress": int((pf["echeance_dette"] <= horizon_renouvellement).sum())}
	351	+ "nb_renouvellements_stress": int((pf["echeance_dette"] <= horizon_renouvellement).sum()),
	352	+ "nb_rcd_det_sous_seuil": int((rcd_det < seuil_rcd).sum()), "nb_rcd_det_sous_1": int((rcd_det < 1.0).sum()), "rcd_portefeuille_det": float(pf["rne_2026"].sum() / sd0.sum()),
	353	+ "nb_rcd_sous_seuil_q95": float(np.quantile((rcd < seuil_rcd).sum(1), 0.95)), "nb_rcd_sous_1_q95": float(np.quantile((rcd < 1.0).sum(1), 0.95)),
	354	+ "P_plus_de_10_sous_1": float(((rcd < 1.0).sum(1) > 10).mean()), "rpv_portefeuille": float(pf["dette"].sum() / V0)}
350	355	stats_["ratio_diversification"] = stats_["ecart_type"] / stats_["somme_ecarts_types"]
351	356	return {"valeurs": V, "valeur_portefeuille": VP, "rne": rne, "rcd": rcd, "service_dette": sd, "stats": stats_,
352	357	"par_rmr": pd.DataFrame({"rmr": rmrs.cat.categories, "ecart_type_valeur": [V[:, k == j].sum(1).std(ddof=1) for j in range(K)], "valeur_2026": [pf.loc[k == j, "valeur_2026"].sum() for j in range(K)]})}
		@@ -356,7 +361,7 @@ def comparer_structures_correlation(pf: pd.DataFrame, croissance: pd.DataFrame,
356	361	rows = []
357	362	for nom, rp, rr in [("Indépendance", 0.0, 0.0), ("Régionale (0,30 ; 0,60)", 0.30, 0.60), ("Parfaite", 1.0, 1.0)]:
358	363	s = simuler_portefeuille(pf, n, GRAINE, rho_rmr=rr, rho_prov=rp, croissance=croissance)["stats"]
359		− rows.append({"structure": nom, "rho_prov": rp, "rho_rmr": rr, **{k: s[k] for k in ["mediane", "ecart_type", "q05", "VaR05", "CVaR05", "ratio_diversification", "nb_rcd_sous_seuil_moyen", "nb_rcd_sous_1_moyen", "P_au_moins_un_rcd_sous_1", "P_rcd_portefeuille_sous_seuil"]}})
	364	+ rows.append({"structure": nom, "rho_prov": rp, "rho_rmr": rr, **{k: s[k] for k in ["mediane", "ecart_type", "q05", "VaR05", "CVaR05", "ratio_diversification", "nb_rcd_sous_seuil_moyen", "nb_rcd_sous_seuil_q95", "nb_rcd_sous_1_moyen", "nb_rcd_sous_1_q95", "P_plus_de_10_sous_1", "rcd_portefeuille_mediane"]}})
360	365	return pd.DataFrame(rows)
361	366
362	367	def variance_somme_correlee(sigmas: np.ndarray, rho: float) -> float:
		@@ -400,8 +405,13 @@ def tout(rapide: bool = False) -> dict:
400	405	D = {**r, "mc": mc, "tornade": torn, "lois": tableau_lois(x), "corr_rRT": correlation_r_RT(x), "erreurs": erreurs_theoriques(mc),
401	406	"risque": mesures_risque(mc, r["dcf"]["valeur"]), "rcd_analytique": rcd1_analytique(r), "stress": stress_deterministes(r, mc), "inverse": stress_inverse(r),
402	407	"sobol": sobol(r), "rangs": correlations_rang(mc), "reg": regression_standardisee(mc), "pert_rotation": pert(0.08, 0.15, 0.25)}
403		− k_rep = 20 if rapide else 100; k_se = 10 if rapide else 40
404		− D["repetitions"] = repetitions(r, N_MC, k_rep, graine0=1)
	408	+ k_rep = 20 if rapide else 100; k_se = 10 if rapide else 40; k_med = 40 if rapide else 300
	409	+ D["repetitions"] = repetitions(r, N_MC, k_med, graine0=1)
	410	+ rep = D["repetitions"]; D["graine_traite"] = {"mediane": mc["stats"]["mediane"], "moyenne_des_medianes": float(rep["mediane"].mean()), "se_mediane_empirique": float(rep["mediane"].std(ddof=1)),
	411	+ "z_mediane": float((mc["stats"]["mediane"] - rep["mediane"].mean()) / rep["mediane"].std(ddof=1)),
	412	+ "z_q95": float((mc["stats"]["q95"] - rep["q95"].mean()) / rep["q95"].std(ddof=1)), "z_q05": float((mc["stats"]["q05"] - rep["q05"].mean()) / rep["q05"].std(ddof=1)),
	413	+ "rang_mediane": int((rep["mediane"] < mc["stats"]["mediane"]).sum()), "k": len(rep),
	414	+ "moyenne_tiree_r": float(x["taux_actualisation"].mean()), "moyenne_tiree_gm": float(x["croissance_loyers_marche"].mean())}
405	415	D["convergence"] = convergence(r, N_max=20_000 if rapide else 100_000)
406	416	D["se_N"] = se_empirique_selon_N(r, k=k_se)
407	417	D["lhs"] = comparer_lhs(r, n=1000, k=k_rep)
		@@ -429,6 +439,7 @@ def imprimer(D: dict) -> None:
429	439	rep = D["repetitions"]
430	440	print(f"\nRépétitions ({len(rep)} graines, N = 10 000) : écart-type empirique — moyenne {rep['moyenne'].std(ddof=1):,.0f} | médiane {rep['mediane'].std(ddof=1):,.0f} | q05 {rep['q05'].std(ddof=1):,.0f} | q95 {rep['q95'].std(ddof=1):,.0f} | P(RCD1<1,20) {rep['p_rcd1_120'].std(ddof=1):.4f}")
431	441	print(f" moyenne des médianes {rep['mediane'].mean():,.0f} ; min {rep['mediane'].min():,.0f} ; max {rep['mediane'].max():,.0f} ; P(RCD1<1,20) moyenne {rep['p_rcd1_120'].mean():.4f}")
	442	+ print(" graine du traité (20261031) :", {k: (round(v, 4) if isinstance(v, float) else v) for k, v in D["graine_traite"].items()})
432	443	print("\nConvergence (une suite, graine 20261031) :\n", D["convergence"].to_string(index=False))
433	444	print(f"\nÉcart-type empirique selon N (pente log-log médiane = {D['se_N'].attrs['pente_log_mediane']:.3f}) :\n", D["se_N"].to_string(index=False))
434	445	print("\n=== 40.4 Hypercube latin vs simple (N = 1 000, répétitions) :\n", D["lhs"].to_string(index=False))

added analysis/python/cases/ch42_copro.py +326 −0

		@@ -0,0 +1,326 @@
	1	+#!/usr/bin/env python3
	2	+"""Chapitre 42 — Copropriétés et conversions : exemples chiffrés (DONNÉES SYNTHÉTIQUES, graine 20261031).
	3	+
	4	+Réutilise l'immeuble fictif de huit fractions du chapitre 7 (analysis/python/cases/ch07_servitudes_copro.py : quotes-parts,
	5	+répartition des charges, cotisation de seuil, déficit de capitalisation) et y ajoute quatre blocs, tous générés ici
	6	+(aucun fichier externe ; aucune statistique de marché n'est reproduite) :
	7	+ A. Évaluation de la fraction F4 par comparaison : deux ventes dans l'immeuble, trois ventes dans d'autres immeubles fictifs ;
	8	+ ajustements de date, d'étage, d'orientation, de stationnement, de superficie, de charges courantes (capitalisation) et de
	9	+ situation financière du syndicat (déficit du fonds de prévoyance, cotisation spéciale) ; réconciliation (traite.comparables).
	10	+ B. Effet d'une cotisation spéciale égale au déficit de capitalisation, répartie selon les quotes-parts (figure).
	11	+ C. Copropriété locative : revenu net d'exploitation et taux implicite d'une fraction louée (traite.finance).
	12	+ D. Conversion d'un immeuble locatif de huit logements en copropriété divise : valeur en bloc (capitalisation directe) contre
	13	+ somme des fractions nette des coûts, des délais et des locataires (rotation, préemption, vente occupée) ; sensibilités.
	14	+Toutes les fonctions financières viennent de analysis/python/traite/ (finance.py, comparables.py) ; rien n'est réécrit.
	15	+Usage (racine du dépôt) : PYTHONPATH=analysis/python python3 analysis/python/cases/ch42_copro.py
	16	+"""
	17	+from __future__ import annotations
	18	+import sys, math, pathlib
	19	+import numpy as np
	20	+HERE = pathlib.Path(__file__).resolve()
	21	+sys.path.insert(0, str(HERE.parents[1])); sys.path.insert(0, str(HERE.parent))
	22	+from traite import comparables as cp
	23	+from traite import finance as fin
	24	+import ch07_servitudes_copro as c7
	25	+
	26	+SEED = 20261031
	27	+DATE_EVAL = "2026-09-01"
	28	+
	29	+# ----------------------------------------------------------------------------- marché fictif : paramètres déclarés
	30	+# Prix implicites « du marché fictif » (hypothèses pédagogiques ; l'évaluateur les tient d'une analyse de ventes appariées / hédonique).
	31	+PRIX_M2_REF = 3_650.0 # $/m² pour une fraction de 98 m², 2e étage, orientation est/ouest, sans stationnement, immeuble bien capitalisé
	32	+ELASTICITE_S = 0.60 # élasticité prix-superficie (< 1) ; prix marginal au point de référence = 0,60 × 3 650 = 2 190 $/m²
	33	+S_REF = 98.0
	34	+COEF_ETAGE_MARCHE = {"RDC": -0.04, "2e": 0.00, "3e": 0.02} # proportion du prix (immeubles sans ascenseur)
	35	+COEF_ORIENT_MARCHE = {"nord": -0.03, "est": 0.00, "ouest": 0.00, "sud": 0.03}
	36	+STATIONNEMENT = 28_000.0 # $ ; case intérieure à usage restreint
	37	+G_MARCHE = 0.02 # croissance annuelle des prix des copropriétés sur le marché fictif (ajustement de conditions du marché)
	38	+R_CHARGES = 0.08 # taux de capitalisation des écarts de charges courantes (hypothèse pédagogique ; à extraire du marché)
	39	+KAPPA = 1.0 # coefficient de capitalisation du passif du syndicat (1 = dollar pour dollar, convention du ch. 7)
	40	+CHARGES_REF_M2 = 40_000.0 / 710.0 # 56,34 $/m²/an : charges courantes de l'immeuble fictif (hors fonds), référence du marché
	41	+SIGMA_BRUIT = 0.02 # bruit de transaction (log) dans la génération des prix des comparables
	42	+
	43	+# Situation financière des syndicats fictifs : déficit de capitalisation (ch. 7) ou cotisation spéciale votée non payée
	44	+QP = c7.quotes_parts()
	45	+INC = c7.incidence_valeur(QP) # déficit de l'immeuble fictif du ch. 7 : 150 908 $
	46	+DEFICIT_IMMEUBLE = float(INC["deficit"])
	47	+Q_F4 = float(QP["q"][3]); Q_F2 = float(QP["q"][1]); Q_F7 = float(QP["q"][6])
	48	+
	49	+SUJET = dict(nom="F4 (sujet)", immeuble="A (fictif, 8 fractions)", S=98.0, etage="2e", orient="sud", stat=True,
	50	+ charges_m2=CHARGES_REF_M2, passif=Q_F4 * DEFICIT_IMMEUBLE, mois=0.0)
	51	+
	52	+COMPARABLES = [ # (nom, immeuble, S, étage, orientation, stationnement, charges courantes $/m²/an, passif de la fraction ($), ancienneté (mois), date)
	53	+ dict(nom="C1 = F2", immeuble="A (même immeuble)", S=98.0, etage="RDC", orient="sud", stat=True, charges_m2=CHARGES_REF_M2, passif=Q_F2 * DEFICIT_IMMEUBLE, mois=6.0, date="2026-03-02"),
	54	+ dict(nom="C2 = F7", immeuble="A (même immeuble)", S=85.0, etage="2e", orient="est", stat=True, charges_m2=CHARGES_REF_M2, passif=Q_F7 * DEFICIT_IMMEUBLE, mois=10.0, date="2025-11-03"),
	55	+ dict(nom="C3", immeuble="B (24 fractions, ascenseur)", S=96.0, etage="2e", orient="sud", stat=True, charges_m2=78.0, passif=0.0, mois=3.0, date="2026-06-01"),
	56	+ dict(nom="C4", immeuble="C (12 fractions)", S=100.0, etage="3e", orient="ouest", stat=False, charges_m2=52.0, passif=0.09 * 60_000.0, mois=4.0, date="2026-05-04"),
	57	+ dict(nom="C5", immeuble="D (8 fractions)", S=92.0, etage="RDC", orient="sud", stat=True, charges_m2=58.0, passif=0.13 * 240_000.0, mois=8.0, date="2026-01-05"),
	58	+]
	59	+
	60	+
	61	+def valeur_modele(u: dict) -> float:
	62	+ """Valeur « vraie » d'une fraction dans le marché fictif à la date d'évaluation (hypothèse de génération).
	63	+ V = base(S) × (1 + c_étage + c_orient) + stationnement − κ × passif − (c_m² − c_ref) × S / R_c."""
	64	+ base = PRIX_M2_REF * S_REF * (u["S"] / S_REF) ** ELASTICITE_S
	65	+ v = base * (1.0 + COEF_ETAGE_MARCHE[u["etage"]] + COEF_ORIENT_MARCHE[u["orient"]])
	66	+ v += STATIONNEMENT if u["stat"] else 0.0
	67	+ v -= KAPPA * u["passif"]
	68	+ v -= (u["charges_m2"] - CHARGES_REF_M2) * u["S"] / R_CHARGES
	69	+ return v
	70	+
	71	+
	72	+def generer_prix(comparables=COMPARABLES, seed: int = SEED) -> list[dict]:
	73	+ """Prix observés = valeur du modèle déflatée à la date de vente (croissance g) × bruit, arrondis à 500 $."""
	74	+ rng = np.random.default_rng(seed)
	75	+ out = []
	76	+ for c in comparables:
	77	+ v = valeur_modele(c)
	78	+ p = v * (1.0 + G_MARCHE) ** (-c["mois"] / 12.0) * math.exp(rng.normal(0.0, SIGMA_BRUIT))
	79	+ d = dict(c); d["prix"] = float(np.round(p / 500.0) * 500.0); d["valeur_modele"] = v
	80	+ out.append(d)
	81	+ return out
	82	+
	83	+
	84	+def prix_marginal_m2(S: float = S_REF) -> float:
	85	+ """Prix implicite marginal d'un m² au voisinage de la superficie du sujet : dV/dS = ε × base(S) / S."""
	86	+ base = PRIX_M2_REF * S_REF * (S / S_REF) ** ELASTICITE_S
	87	+ return ELASTICITE_S * base / S
	88	+
	89	+
	90	+def ajuster(comp: dict, sujet: dict = SUJET, r_charges: float = R_CHARGES, kappa: float = KAPPA, g: float = G_MARCHE) -> cp.ComparableAjuste:
	91	+ """Grille : ajustement de conditions du marché (%, séquentiel) puis ajustements de propriété en dollars, signe du point de vue du sujet."""
	92	+ p_t = comp["prix"] * (1.0 + g) ** (comp["mois"] / 12.0) # prix ramené à la date d'évaluation
	93	+ A = {"marché": p_t - comp["prix"]}
	94	+ A["étage"] = p_t * (COEF_ETAGE_MARCHE[sujet["etage"]] - COEF_ETAGE_MARCHE[comp["etage"]])
	95	+ A["orientation"] = p_t * (COEF_ORIENT_MARCHE[sujet["orient"]] - COEF_ORIENT_MARCHE[comp["orient"]])
	96	+ A["stationnement"] = STATIONNEMENT * (int(sujet["stat"]) - int(comp["stat"]))
	97	+ A["superficie"] = prix_marginal_m2() * (sujet["S"] - comp["S"])
	98	+ A["charges"] = (comp["charges_m2"] - sujet["charges_m2"]) * sujet["S"] / r_charges
	99	+ A["syndicat"] = kappa * (comp["passif"] - sujet["passif"])
	100	+ return cp.ComparableAjuste(comp["nom"], comp["prix"], A)
	101	+
	102	+
	103	+def grille_fraction(r_charges: float = R_CHARGES, kappa: float = KAPPA, g: float = G_MARCHE, poids: str = "inverse_brut") -> dict:
	104	+ comps = generer_prix()
	105	+ ajustes = [ajuster(c, SUJET, r_charges, kappa, g) for c in comps]
	106	+ w = cp.poids_inverse_brut(ajustes) if poids == "inverse_brut" else None
	107	+ rec = cp.reconciliation(ajustes, w)
	108	+ return {"comparables": comps, "ajustes": ajustes, "reconciliation": rec, "poids": rec["poids"],
	109	+ "valeur_modele_sujet": valeur_modele(SUJET), "prix_marginal_m2": prix_marginal_m2()}
	110	+
	111	+
	112	+def sensibilite_fraction(r_charges_list=(0.05, 0.06, 0.08, 0.10, 0.12), kappas=(0.0, 0.5, 0.75, 1.0)) -> dict:
	113	+ grille = {k: [grille_fraction(r, k)["reconciliation"]["valeur_ponderee"] for r in r_charges_list] for k in kappas}
	114	+ return {"r_charges": list(r_charges_list), "kappas": list(kappas), "grille": grille}
	115	+
	116	+
	117	+# ----------------------------------------------------------------------------- B. cotisation spéciale
	118	+def cotisation_speciale(montant: float | None = None, qp: dict = QP) -> dict:
	119	+ """Répartition d'une cotisation spéciale (par défaut égale au déficit de capitalisation du ch. 7) selon les quotes-parts ;
	120	+ comparaison avec les charges annuelles normalisées (cotisation de seuil) de chaque fraction."""
	121	+ montant = DEFICIT_IMMEUBLE if montant is None else montant
	122	+ A_star = c7.cotisation_seuil()
	123	+ rep = c7.repartition_charges(A_star, qp)
	124	+ part = qp["q"] * montant
	125	+ return {"montant": montant, "par_fraction": part, "charges_annuelles": rep["par_fraction_total"],
	126	+ "ratio_annees": part / rep["par_fraction_total"], "A_star": A_star, "rep": rep}
	127	+
	128	+
	129	+# ----------------------------------------------------------------------------- C. copropriété locative (F4 louée)
	130	+LOYER_F4 = 1_700.0 # $/mois ; hypothèse pédagogique dans l'ordre de grandeur des loyers SCHL (oct. 2025) des copropriétés offertes en location, 2 ch.
	131	+INOCC_COPRO = 0.021 # SCHL, RMR de Montréal, copropriétés offertes en location, octobre 2025 (tableau 4.1.1)
	132	+TAXES_F4 = 3_300.0 # $/an ; ordre de grandeur du compte de taxes moyen d'une copropriété à Montréal en 2026 (La Presse, 2026-10-19, source secondaire)
	133	+ASSURANCE_COPRO = 600.0 # $/an ; hypothèse (BAC 2022 : 544 $ pour un copropriétaire, partiellement vérifié)
	134	+GESTION = 0.05 # part du revenu brut effectif (forfait usuel)
	135	+TAUX_PLEX = 0.0525 # taux global de référence des petits immeubles locatifs du marché fictif (dans la fourchette CBRE Québec Low Rise B T3 2026 : 5,00-6,00 %)
	136	+
	137	+
	138	+def copro_locative(valeur: float, loyer: float = LOYER_F4, inocc: float = INOCC_COPRO, taxes: float = TAXES_F4,
	139	+ assurance: float = ASSURANCE_COPRO, gestion: float = GESTION, frais_copro: float | None = None) -> dict:
	140	+ """RNE d'une fraction louée (convention de l'ouvrage : la cotisation au fonds de prévoyance normalisée tient lieu de réserve)."""
	141	+ if frais_copro is None:
	142	+ frais_copro = float(c7.repartition_charges(c7.cotisation_seuil(), QP)["par_fraction_total"][3]) # F4 : 10 218 $
	143	+ rbp = loyer * 12.0
	144	+ rbe = fin.revenu_brut_effectif(rbp, inocc)
	145	+ de = taxes + assurance + frais_copro + gestion * rbe
	146	+ noi = fin.rne(rbe, de)
	147	+ fonds = float(c7.repartition_charges(c7.cotisation_seuil(), QP)["par_fraction_fonds"][3]) # part « fonds de prévoyance » des frais (3 082 $)
	148	+ return {"rbp": rbp, "rbe": rbe, "depenses": de, "frais_copro": frais_copro, "noi": noi, "noi_sans_reserve": noi + fonds, "fonds": fonds,
	149	+ "taux_implicite_sans_reserve": fin.taux_global_extrait(noi + fonds, valeur),
	150	+ "taux_implicite": fin.taux_global_extrait(noi, valeur), "ratio_depenses": fin.ratio_depenses(de, rbe),
	151	+ "valeur_investissement_5_25": fin.capitalisation_directe(noi, TAUX_PLEX),
	152	+ "mrb": fin.multiplicateur_revenu_brut(valeur, rbp)}
	153	+
	154	+
	155	+# ----------------------------------------------------------------------------- D. conversion locatif → copropriété divise
	156	+# L'immeuble fictif de huit fractions, vu AVANT conversion comme un immeuble locatif de 8 logements (hypothèses déclarées).
	157	+LOYERS_LOG = {"F1": 1_250.0, "F2": 1_650.0, "F3": 1_300.0, "F4": 1_700.0, "F5": 1_450.0, "F6": 1_850.0, "F7": 1_500.0, "F8": 1_525.0} # $/mois, loyers contractuels
	158	+INOCC_LOCATIF = 0.03 # inoccupation + mauvaises créances (hypothèse ; SCHL RMR 2025 : 2,4 % Québec, 2,9 % Montréal, 3,8 % Gatineau)
	159	+RATIO_DEPENSES_BLOC = 0.42 # dépenses d'exploitation + réserve / RBE (hypothèse pédagogique pour un immeuble de 8 logements des années 1990)
	160	+LOCATAIRES = {"F1": True, "F2": False, "F3": True, "F4": False, "F5": True, "F6": True, "F7": True, "F8": False} # occupé par un locataire (bail en cours) ?
	161	+PREEMPTION = {"F6"} # locataire(s) qui acceptent l'offre de vente (art. 54.7) à la valeur libre
	162	+ROTATION = 0.12 # part des logements occupés libérés chaque année (départ volontaire ; aucune éviction : moratoire jusqu'au 6 juin 2027)
	163	+DECOTE_OCCUPE = 0.20 # décote d'une fraction vendue occupée (bail maintenu, art. 1937 C.c.Q.) à un investisseur — hypothèse
	164	+R_PROMOTEUR = 0.10 # taux d'actualisation exigé par le convertisseur
	165	+FRAIS_VENTE = 0.04 # courtage et frais de vente (proportion du prix) — hypothèse
	166	+COUTS_CONVERSION = { # dollars de 2026, hypothèses pédagogiques
	167	+ "Cadastre vertical et certificats de localisation (arpenteur-géomètre)": 28_000.0,
	168	+ "Déclaration de copropriété, publication et honoraires juridiques": 22_000.0,
	169	+ "Rapport d'expert et circulaire d'information (art. 54.6 Loi TAL)": 9_000.0,
	170	+ "Frais de demande (TAL, municipalité), budget prévisionnel, valeur relative": 6_000.0,
	171	+ "Individualisation des services et mise aux normes (8 logements)": 120_000.0,
	172	+ "Rafraîchissement des parties communes": 35_000.0,
	173	+ "Dotation initiale du fonds de prévoyance (solde d'ouverture requis, ch. 7)": 85_000.0,
	174	+}
	175	+HORIZON_CONV = 6 # années de projection des ventes
	176	+PORTAGE_PCT = 0.015 # coûts de portage annuels (taxes, assurance, intérêts nets des loyers résiduels) en proportion de la valeur des fractions invendues — hypothèse
	177	+
	178	+
	179	+def unites_fictives() -> list[dict]:
	180	+ """Les huit fractions de l'immeuble du ch. 7 comme unités du marché fictif (charges de référence, passif nul après dotation du fonds)."""
	181	+ out = []
	182	+ for i, (n, S, e, o, st) in enumerate(c7.FRACTIONS):
	183	+ out.append(dict(nom=n, S=S, etage=e, orient=o, stat=st, charges_m2=CHARGES_REF_M2, passif=0.0, loyer=LOYERS_LOG[n],
	184	+ locataire=LOCATAIRES[n], preemption=n in PREEMPTION, q=float(QP["q"][i])))
	185	+ return out
	186	+
	187	+
	188	+def valeur_bloc(taux: float = TAUX_PLEX, inocc: float = INOCC_LOCATIF, ratio_dep: float = RATIO_DEPENSES_BLOC) -> dict:
	189	+ u = unites_fictives()
	190	+ rbp = 12.0 * sum(x["loyer"] for x in u)
	191	+ rbe = fin.revenu_brut_effectif(rbp, inocc)
	192	+ de = ratio_dep * rbe
	193	+ noi = fin.rne(rbe, de)
	194	+ return {"rbp": rbp, "rbe": rbe, "depenses": de, "noi": noi, "taux": taux, "valeur": fin.capitalisation_directe(noi, taux),
	195	+ "mrb": fin.multiplicateur_revenu_brut(fin.capitalisation_directe(noi, taux), rbp), "n": len(u)}
	196	+
	197	+
	198	+def calendrier_ventes(rotation: float = ROTATION, decote_occ: float = DECOTE_OCCUPE, vendre_occupe: bool = False,
	199	+ horizon: int = HORIZON_CONV, unites=None) -> dict:
	200	+ """Produit brut attendu par année (fin d'année), avant frais de vente.
	201	+ An 1 : fractions libres et locataires qui exercent la préemption, à la valeur libre. Fractions occupées : si `vendre_occupe`,
	202	+ vendues à l'an 1 avec décote ; sinon vendues libres l'année où le locataire part (probabilité ρ par an), le reliquat à l'horizon
	203	+ étant vendu occupé avec décote. Valeurs libres = modèle du marché fictif (passif nul après dotation du fonds)."""
	204	+ u = unites or unites_fictives()
	205	+ produit = np.zeros(horizon + 1); invendu = np.zeros(horizon + 1)
	206	+ valeurs = {x["nom"]: valeur_modele(x) for x in u}
	207	+ libres_an1 = [x for x in u if (not x["locataire"]) or x["preemption"]]
	208	+ occupes = [x for x in u if x["locataire"] and not x["preemption"]]
	209	+ produit[1] += sum(valeurs[x["nom"]] for x in libres_an1)
	210	+ V_occ = sum(valeurs[x["nom"]] for x in occupes)
	211	+ if vendre_occupe:
	212	+ produit[1] += V_occ * (1.0 - decote_occ)
	213	+ else:
	214	+ reste = 1.0
	215	+ for t in range(1, horizon + 1):
	216	+ if t < horizon:
	217	+ produit[t] += V_occ * reste * rotation # fractions libérées pendant l'année t, vendues libres en fin d'année
	218	+ reste *= (1.0 - rotation)
	219	+ else:
	220	+ produit[t] += V_occ * reste * (1.0 - decote_occ) # reliquat vendu occupé à l'horizon
	221	+ reste = 0.0
	222	+ invendu[t] = V_occ * reste
	223	+ # valeur libre des fractions invendues au début de chaque année (assiette des coûts de portage)
	224	+ base_portage = np.zeros(horizon + 1)
	225	+ total = sum(valeurs.values())
	226	+ base_portage[1] = total
	227	+ if not vendre_occupe:
	228	+ reste = 1.0
	229	+ for t in range(2, horizon + 1):
	230	+ reste *= (1.0 - rotation)
	231	+ base_portage[t] = V_occ * reste
	232	+ return {"produit": produit, "valeurs": valeurs, "somme_libre": total, "V_occ": V_occ, "n_libres_an1": len(libres_an1),
	233	+ "n_occupes": len(occupes), "base_portage": base_portage}
	234	+
	235	+
	236	+def faisabilite(rotation: float = ROTATION, decote_occ: float = DECOTE_OCCUPE, r: float = R_PROMOTEUR, vendre_occupe: bool = False,
	237	+ frais_vente: float = FRAIS_VENTE, couts=COUTS_CONVERSION, portage: float = PORTAGE_PCT, horizon: int = HORIZON_CONV) -> dict:
	238	+ """Valeur de conversion V_conv = Σ_t [produit_t (1 − f) − portage_t] / (1+r)^t − coûts de conversion (t = 0) ; prime = V_conv − V_bloc."""
	239	+ cal = calendrier_ventes(rotation, decote_occ, vendre_occupe, horizon)
	240	+ flux = np.zeros(horizon + 1)
	241	+ flux[0] = -sum(couts.values())
	242	+ for t in range(1, horizon + 1):
	243	+ flux[t] = cal["produit"][t] * (1.0 - frais_vente) - portage * cal["base_portage"][t]
	244	+ v_conv = fin.van(list(flux), r)
	245	+ bloc = valeur_bloc()
	246	+ return {"flux": flux, "v_conv": v_conv, "v_bloc": bloc["valeur"], "prime": v_conv - bloc["valeur"],
	247	+ "somme_libre": cal["somme_libre"], "couts": sum(couts.values()), "frais_vente_total": float(sum(cal["produit"]) * frais_vente),
	248	+ "portage_total": float(sum(portage * cal["base_portage"][1:])), "produit": cal["produit"], "calendrier": cal, "bloc": bloc,
	249	+ "marge_sur_bloc": (v_conv - bloc["valeur"]) / bloc["valeur"],
	250	+ "somme_nette_non_actualisee": float(sum(cal["produit"]) * (1 - frais_vente) - sum(couts.values()) - sum(portage * cal["base_portage"][1:]))}
	251	+
	252	+
	253	+def sensibilite_conversion(rotations=(0.05, 0.08, 0.12, 0.20, 0.30), decotes=(0.10, 0.20, 0.30), r: float = R_PROMOTEUR) -> dict:
	254	+ grille = {d: [faisabilite(rho, d, r)["v_conv"] for rho in rotations] for d in decotes}
	255	+ occ = {d: faisabilite(ROTATION, d, r, vendre_occupe=True)["v_conv"] for d in decotes}
	256	+ return {"rotations": list(rotations), "decotes": list(decotes), "grille": grille, "vente_occupee_an1": occ, "v_bloc": valeur_bloc()["valeur"]}
	257	+
	258	+
	259	+def seuil_prix(r: float = R_PROMOTEUR) -> float:
	260	+ """Facteur multiplicatif uniforme des valeurs libres pour lequel V_conv = V_bloc (point mort de la conversion)."""
	261	+ global PRIX_M2_REF
	262	+ p0 = PRIX_M2_REF
	263	+ lo, hi = 0.3 * p0, 1.5 * p0
	264	+ for _ in range(100):
	265	+ mid = 0.5 * (lo + hi); PRIX_M2_REF = mid
	266	+ if faisabilite(r=r)["prime"] > 0: hi = mid
	267	+ else: lo = mid
	268	+ PRIX_M2_REF = p0
	269	+ return hi / p0
	270	+
	271	+
	272	+# ----------------------------------------------------------------------------- impression
	273	+def fmt(x): return f"{x:,.0f}".replace(",", " ")
	274	+
	275	+
	276	+def main() -> dict:
	277	+ out = {}
	278	+ g = grille_fraction(); out["grille"] = g
	279	+ print(f"== A. Fraction F4 (sujet) — marché fictif, graine {SEED}, date d'évaluation {DATE_EVAL}")
	280	+ print(f" passif du syndicat A (déficit {fmt(DEFICIT_IMMEUBLE)} $) : part de F4 = q4 {Q_F4:.4f} × déficit = {fmt(SUJET['passif'])} $ ; prix marginal {g['prix_marginal_m2']:,.0f} $/m²")
	281	+ print(f" valeur du modèle pour F4 (contrôle de génération) : {fmt(g['valeur_modele_sujet'])} $")
	282	+ for c, a, w in zip(g["comparables"], g["ajustes"], g["poids"]):
	283	+ print(f" {c['nom']:8s} {c['immeuble']:28s} S {c['S']:5.0f} {c['etage']:3s} {c['orient']:5s} stat {int(c['stat'])} charges {c['charges_m2']:5.1f} $/m² passif {fmt(c['passif']):>7} $ "
	284	+ f"{c['mois']:4.0f} mois prix {fmt(c['prix']):>8} $ (modèle {fmt(c['valeur_modele'])})")
	285	+ print(" ajustements : " + ", ".join(f"{k} {v:+,.0f}" for k, v in a.ajustements.items()).replace(",", " ").replace(" $", "$"))
	286	+ print(f" prix ajusté {fmt(a.prix_ajuste)} $ ; brut {a.brut_pct:.1%} ; net {a.net_pct:+.1%} ; poids {w:.3f}")
	287	+ rec = g["reconciliation"]
	288	+ print(f" réconciliation : pondérée {fmt(rec['valeur_ponderee'])} $ ; moyenne {fmt(rec['moyenne'])} ; médiane {fmt(rec['mediane'])} ; min {fmt(rec['min'])} ; max {fmt(rec['max'])} ; étendue {rec['etendue_pct']:.1%}")
	289	+ sf = sensibilite_fraction(); out["sens_fraction"] = sf
	290	+ print(" sensibilité (valeur pondérée) : lignes κ, colonnes R_c " + " ".join(f"{r:.0%}" for r in sf["r_charges"]))
	291	+ for k, row in sf["grille"].items():
	292	+ print(f" κ={k:.2f} " + " ".join(f"{v:9,.0f}" for v in row).replace(",", " "))
	293	+ cs = cotisation_speciale(); out["cotisation"] = cs
	294	+ print(f"\n== B. Cotisation spéciale de {fmt(cs['montant'])} $ (déficit du ch. 7) répartie selon q_i ; charges annuelles normalisées (A* = {fmt(cs['A_star'])} $)")
	295	+ for n, p, ch, r in zip(QP["noms"], cs["par_fraction"], cs["charges_annuelles"], cs["ratio_annees"]):
	296	+ print(f" {n} : cotisation {fmt(p):>7} $ ; charges {fmt(ch):>6} $/an ; = {r:.2f} années de charges")
	297	+ v_f4 = rec["valeur_ponderee"]
	298	+ cl = copro_locative(v_f4); out["locative"] = cl
	299	+ print(f"\n== C. F4 louée {fmt(LOYER_F4)} $/mois : RBP {fmt(cl['rbp'])} ; RBE {fmt(cl['rbe'])} ; dépenses {fmt(cl['depenses'])} (frais copro {fmt(cl['frais_copro'])}) ; RNE {fmt(cl['noi'])} $")
	300	+ print(f" RNE hors cotisation au fonds ({fmt(cl['fonds'])} $) = {fmt(cl['noi_sans_reserve'])} $ → taux {cl['taux_implicite_sans_reserve']:.2%}")
	301	+ print(f" taux implicite sur {fmt(v_f4)} $ = {cl['taux_implicite']:.2%} ; ratio de dépenses {cl['ratio_depenses']:.1%} ; MRB {cl['mrb']:.1f} ; valeur d'investissement à {TAUX_PLEX:.2%} = {fmt(cl['valeur_investissement_5_25'])} $")
	302	+ fb = faisabilite(); out["faisabilite"] = fb
	303	+ b = fb["bloc"]
	304	+ print(f"\n== D. Conversion — immeuble locatif de 8 logements : RBP {fmt(b['rbp'])} ; RBE {fmt(b['rbe'])} ; DE {fmt(b['depenses'])} ; RNE {fmt(b['noi'])} ; R {b['taux']:.2%} → V_bloc {fmt(b['valeur'])} $ ({fmt(b['valeur']/8)} $/log., MRB {b['mrb']:.1f})")
	305	+ cal = fb["calendrier"]
	306	+ print(" valeurs libres des fractions : " + ", ".join(f"{k} {fmt(v)}" for k, v in cal["valeurs"].items()) + f" ; somme {fmt(cal['somme_libre'])} $")
	307	+ print(f" libres ou préemption à l'an 1 : {cal['n_libres_an1']} fractions ; occupées : {cal['n_occupes']} (valeur libre {fmt(cal['V_occ'])} $) ; ρ = {ROTATION:.0%}, décote occupée {DECOTE_OCCUPE:.0%}, r = {R_PROMOTEUR:.0%}")
	308	+ print(" produit brut par année : " + ", ".join(f"an {t} {fmt(p)}" for t, p in enumerate(fb["produit"]) if t > 0))
	309	+ print(" flux : " + ", ".join(f"t{t} {f:,.0f}".replace(",", " ") for t, f in enumerate(fb["flux"])))
	310	+ print(f" coûts de conversion {fmt(fb['couts'])} ; frais de vente {fmt(fb['frais_vente_total'])} ; portage {fmt(fb['portage_total'])} ; somme nette non actualisée {fmt(fb['somme_nette_non_actualisee'])}")
	311	+ print(f" V_conv = {fmt(fb['v_conv'])} $ ; V_bloc = {fmt(fb['v_bloc'])} $ ; prime {fmt(fb['prime'])} $ ({fb['marge_sur_bloc']:+.1%})")
	312	+ fo = faisabilite(vendre_occupe=True); out["faisabilite_occ"] = fo
	313	+ print(f" variante : tout vendre à l'an 1 (occupées avec décote) → V_conv {fmt(fo['v_conv'])} $ ; prime {fmt(fo['prime'])} $")
	314	+ sc = sensibilite_conversion(); out["sens_conv"] = sc
	315	+ print(" sensibilité V_conv : lignes décote occupée, colonnes ρ " + " ".join(f"{r:.0%}" for r in sc["rotations"]) + " | vente occupée an 1")
	316	+ for d, row in sc["grille"].items():
	317	+ print(f" δ={d:.0%} " + " ".join(f"{v:10,.0f}" for v in row).replace(",", " ") + f" | {sc['vente_occupee_an1'][d]:10,.0f}".replace(",", " "))
	318	+ sp = seuil_prix(); out["seuil_prix"] = sp
	319	+ print(f" point mort : les valeurs libres peuvent baisser jusqu'à {sp:.1%} de leur niveau (−{1-sp:.1%}) avant que V_conv < V_bloc")
	320	+ for r in (0.08, 0.12):
	321	+ print(f" r = {r:.0%} : V_conv {fmt(faisabilite(r=r)['v_conv'])} $")
	322	+ return out
	323	+
	324	+
	325	+if __name__ == "__main__":
	326	+ main()

added analysis/python/cases/ch43_risques.py +386 −0

		@@ -0,0 +1,386 @@
	1	+#!/usr/bin/env python3
	2	+"""Chapitre 43 — Risques environnementaux et sinistres : DONNÉES SYNTHÉTIQUES et calculs de l'exemple intégral.
	3	+
	4	+Quatre blocs, tous reproductibles (graine 20261031) :
	5	+
	6	+1. Cas 4 (portefeuille, `traite.cas4`) : lecture du coefficient `inondable` et `sol` de `effets_regionaux`
	7	+ (prime de risque dans le log du taux global ; loi déclarée du générateur : +5 % du taux en zone inondable,
	8	+ +3 % si contamination), part de la valeur du portefeuille située en zone inondable.
	9	+2. Cas 4, doubles différences : on injecte sur les séries trimestrielles de RNE des trois immeubles
	10	+ `zone_inondable == 1` un choc d'inondation déclaré au printemps 2019 (chute transitoire puis perte
	11	+ permanente de 5 % du RNE), puis on l'estime par DiD 2 × 2 (naïf et par strate de RMR), par étude
	12	+ d'événement à effets fixes immeuble et RMR × trimestre (`traite.panel`) et par inférence de
	13	+ randomisation exacte (toutes les affectations de 3 immeubles traités parmi 48).
	14	+3. Municipalité riveraine fictive « Rivière-Fictive » : 625 lots dont 175 en zone inondable (classes très
	15	+ élevée / élevée / modérée / faible) le long d'une rivière et 75 lots riverains d'un lac non inondable, ventes
	16	+ 2016-2026 avec décote déclarée par classe, stigmate après une inondation le 1er mai 2019 et aménité « vue sur
	17	+ l'eau » (rivière et lac) qui brouille la mesure si on l'omet. Quatre méthodes de mesure :
	18	+ ventes appariées, hédonique à indicatrice (avec et sans contrôle de l'aménité), hédonique par classe avec
	19	+ interaction classe × après, doubles différences (cellules 2 × 2 et régression).
	20	+4. Modèles de valeur : terrain contaminé (coût de réhabilitation + délai + stigmate), immeuble sinistré
	21	+ (remise en état + perte d'usage + stigmate résiduel ; rôle de l'assurance), désuétude vs coût de correction.
	22	+
	23	+Écrit data/synthetic/ch43_riverain_ventes.csv et ch43_riverain_lots.csv s'ils n'existent pas.
	24	+Usage : python3 analysis/python/cases/ch43_risques.py
	25	+"""
	26	+from __future__ import annotations
	27	+import sys, pathlib, math, itertools
	28	+from datetime import date, timedelta
	29	+import numpy as np
	30	+import pandas as pd
	31	+
	32	+ROOT = pathlib.Path(__file__).resolve().parents[3]
	33	+sys.path.insert(0, str(ROOT / "analysis/python"))
	34	+from traite import cas4, comparables as cp, hedonique as h, panel # noqa: E402
	35	+
	36	+SEED = 20261031
	37	+D = ROOT / "data/synthetic"
	38	+
	39	+# =============================================================================== 1. cas 4 : prime de risque dans le taux
	40	+
	41	+def prime_taux_cas4() -> dict:
	42	+ """Coefficients `inondable` et `sol` de cas4.effets_regionaux : effet en % du taux, IC à 95 % (HC1),
	43	+ comparaison à la loi déclarée (+5 % / +3 %) et part de la valeur en zone inondable."""
	44	+ d = cas4.charger(); pf = d["portefeuille"]
	45	+ eff = cas4.effets_regionaux(pf); som = cas4.sommaire(pf)
	46	+ out = {"n": len(pf), "n_inondable": int(pf["zone_inondable"].sum()), "n_sol": int(pf["contamination_sol"].sum()),
	47	+ "part_valeur_inondable": som["total"]["part_zone_inondable"], "valeur_totale": som["total"]["valeur"]}
	48	+ for v, declare in [("inondable", 0.05), ("sol", 0.03)]:
	49	+ b, se = float(eff.loc[v, "coef"]), float(eff.loc[v, "se_hc1"])
	50	+ out[v] = {"coef": b, "se": se, "effet_pct": math.exp(b) - 1, "ic_bas_pct": math.exp(b - 1.96 * se) - 1,
	51	+ "ic_haut_pct": math.exp(b + 1.96 * se) - 1, "declare_pct": declare, "declare_log": math.log(1 + declare),
	52	+ "t_vs_declare": (b - math.log(1 + declare)) / se,
	53	+ # correction de Kennedy (variance non négligeable) : exp(b − V/2) − 1
	54	+ "effet_kennedy_pct": math.exp(b - 0.5 * se ** 2) - 1}
	55	+ # décote de valeur impliquée par la prime de taux, à RNE donné : V = RNE/R → ΔlnV = −Δln R
	56	+ out["decote_valeur_inondable_pct"] = math.exp(-out["inondable"]["coef"]) - 1
	57	+ out["decote_valeur_declaree_pct"] = 1 / 1.05 - 1
	58	+ out["immeubles_inondables"] = pf.loc[pf["zone_inondable"] == 1, ["id_immeuble", "rmr", "nb_logements", "valeur_2026", "rne_2026", "taux_global"]].reset_index(drop=True)
	59	+ out["immeubles_sol"] = pf.loc[pf["contamination_sol"] == 1, ["id_immeuble", "rmr", "nb_logements", "valeur_2026", "rne_2026", "taux_global"]].reset_index(drop=True)
	60	+ out["portefeuille"] = pf
	61	+ return out
	62	+
	63	+# =============================================================================== 2. cas 4 : DiD sur les séries de RNE
	64	+
	65	+CHOC_2019 = {"2019Q2": -0.25, "2019Q3": -0.12, "2019Q4": -0.06} # ln RNE, transitoire (perte d'usage, remise en état)
	66	+CHOC_PERMANENT = -0.05 # ln RNE à partir de 2020Q1 (loyers, inoccupation, assurance)
	67	+T_EVENEMENT = "2019Q2"
	68	+
	69	+def panel_choc(seed: int = SEED) -> pd.DataFrame:
	70	+ """Panel long (48 immeubles × 42 trimestres) de ln RNE avec le choc déclaré injecté sur les immeubles
	71	+ `zone_inondable == 1` ; colonnes : id, rmr, t (0..41), trimestre, traite, post, ln_rne, ln_rne_sans_choc, g."""
	72	+ d = cas4.charger(); pf, se = d["portefeuille"], d["series"].copy()
	73	+ se["trimestre"] = se["trimestre"].astype(str)
	74	+ se = se.merge(pf[["id_immeuble", "rmr", "zone_inondable"]], on="id_immeuble")
	75	+ trims = sorted(se["trimestre"].unique()); idx = {q: i for i, q in enumerate(trims)}
	76	+ se["t"] = se["trimestre"].map(idx); t_ev = idx[T_EVENEMENT]
	77	+ se["traite"] = se["zone_inondable"].astype(int); se["post"] = (se["t"] >= t_ev).astype(int)
	78	+ se["ln_rne_sans_choc"] = np.log(se["rne_trimestriel"].astype(float))
	79	+ choc = np.zeros(len(se))
	80	+ for q, c in CHOC_2019.items():
	81	+ choc[(se["trimestre"] == q).values & (se["traite"] == 1).values] = c
	82	+ choc[(se["t"] >= idx["2020Q1"]).values & (se["traite"] == 1).values] = CHOC_PERMANENT
	83	+ se["choc"] = choc; se["ln_rne"] = se["ln_rne_sans_choc"] + choc
	84	+ se["g"] = np.where(se["traite"] == 1, float(t_ev), np.nan)
	85	+ se["unite"] = se["id_immeuble"]
	86	+ return se.sort_values(["id_immeuble", "t"]).reset_index(drop=True)
	87	+
	88	+def did_cas4(pan: pd.DataFrame) -> dict:
	89	+ trims = sorted(pan["trimestre"].unique()); idx = {q: i for i, q in enumerate(trims)}
	90	+ avant = list(range(idx["2017Q1"], idx["2019Q1"] + 1)); apres = list(range(idx["2020Q1"], idx["2026Q2"] + 1))
	91	+ r = panel.did_2x2(pan, "ln_rne", "unite", "t", "traite", avant, apres)
	92	+ # version par strate : Δ_i centré sur la moyenne des témoins de la même RMR
	93	+ dl = r["deltas"].merge(pan[["unite", "rmr"]].drop_duplicates(), on="unite")
	94	+ m_rmr = dl[dl["traite"] == 0].groupby("rmr")["dy"].mean()
	95	+ dl["dy_strate"] = dl["dy"] - dl["rmr"].map(m_rmr)
	96	+ est_strate = float(dl.loc[dl["traite"] == 1, "dy_strate"].mean())
	97	+ # inférence de randomisation exacte : toutes les affectations de 3 traités parmi 48 (C(48,3) = 17 296)
	98	+ dy = dl["dy"].to_numpy(); ids = np.arange(len(dy)); n_tr = int(dl["traite"].sum())
	99	+ tot = dy.sum(); n = len(dy)
	100	+ combos = np.array(list(itertools.combinations(ids, n_tr)))
	101	+ s_tr = dy[combos].sum(axis=1)
	102	+ est_all = s_tr / n_tr - (tot - s_tr) / (n - n_tr)
	103	+ obs = r["estimateur"]
	104	+ p_ri = float(np.mean(np.abs(est_all) >= abs(obs) - 1e-12))
	105	+ # idem pour l'estimateur par strate
	106	+ dys = dl["dy_strate"].to_numpy() # témoins centrés par RMR ; pour les placebos on recentre à chaque tirage
	107	+ rmr_codes = dl["rmr"].astype("category").cat.codes.to_numpy()
	108	+ est_strate_all = np.empty(len(combos))
	109	+ for k, c in enumerate(combos):
	110	+ mask = np.zeros(n, bool); mask[c] = True
	111	+ # moyenne des témoins par RMR
	112	+ num = np.bincount(rmr_codes[~mask], weights=dy[~mask], minlength=rmr_codes.max() + 1)
	113	+ den = np.bincount(rmr_codes[~mask], minlength=rmr_codes.max() + 1)
	114	+ mr = np.where(den > 0, num / np.maximum(den, 1), 0.0)
	115	+ est_strate_all[k] = float(np.mean(dy[mask] - mr[rmr_codes[mask]]))
	116	+ p_ri_strate = float(np.mean(np.abs(est_strate_all) >= abs(est_strate) - 1e-12))
	117	+ # étude d'événement : effets fixes immeuble + RMR × trimestre, fenêtre −8..+12, référence −1
	118	+ es = panel.event_study(pan, "ln_rne", "unite", "t", "g", fenetre=(-8, 12), ref=-1, fe_supp="rmr")
	119	+ # TWFE statique avec RMR × trimestre
	120	+ pan2 = pan.copy(); pan2["traite_post"] = pan2["traite"] * pan2["post"]
	121	+ pan2["traite_permanent"] = pan2["traite"] * (pan2["t"] >= idx["2020Q1"]).astype(int)
	122	+ pan2["traite_transitoire"] = pan2["traite"] * pan2["trimestre"].isin(list(CHOC_2019)).astype(int)
	123	+ tw = panel.twfe(pan2, "ln_rne", ["traite_transitoire", "traite_permanent"], "unite", "t", fe_supp="rmr")
	124	+ # lien avec la valeur : Δln V = Δln RNE − Δln R
	125	+ pr = prime_taux_cas4()
	126	+ return {"cellules": r["cellules"], "estimateur_2x2": obs, "coef_reg": r["coef_reg"], "se_hc1": r["se_hc1"],
	127	+ "estimateur_strate": est_strate, "p_ri": p_ri, "p_ri_strate": p_ri_strate, "n_combinaisons": int(len(combos)),
	128	+ "rang_obs": int(np.sum(est_all <= obs)), "n_traites": r["n_traites"], "n_temoins": r["n_temoins"],
	129	+ "deltas": dl, "event_study": es, "twfe": tw["coef"], "placebo_dist": est_all, "placebo_dist_strate": est_strate_all,
	130	+ "verite_permanente": CHOC_PERMANENT, "verite_permanente_pct": math.exp(CHOC_PERMANENT) - 1,
	131	+ "verite_transitoire_moy": float(np.mean(list(CHOC_2019.values()))),
	132	+ "effet_valeur_estime_log": tw["coef"].loc["traite_permanent", "coef"] - pr["inondable"]["coef"],
	133	+ "effet_valeur_estime_pct": math.exp(tw["coef"].loc["traite_permanent", "coef"] - pr["inondable"]["coef"]) - 1,
	134	+ "effet_valeur_declare_log": CHOC_PERMANENT - math.log(1.05),
	135	+ "effet_valeur_declare_pct": math.exp(CHOC_PERMANENT - math.log(1.05)) - 1,
	136	+ "idx": idx, "trims": trims}
	137	+
	138	+# =============================================================================== 3. Rivière-Fictive : ventes riveraines
	139	+
	140	+NCOL, NLIG, PAS = 25, 25, 40.0 # 625 lots de 40 m × 40 m ; la rivière longe le bord y = 0, un lac (non inondable) le bord x = 24
	141	+COLS_LAC = {22, 23, 24} # lots riverains du lac : aménité « vue sur l'eau » sans exposition à l'inondation
	142	+CLASSES = ["très élevée", "élevée", "modérée", "faible", "hors zone"]
	143	+LIGNES_CLASSE = {0: "très élevée", 1: "élevée", 2: "élevée", 3: "modérée", 4: "modérée", 5: "faible", 6: "faible"}
	144	+DECOTE_AVANT = {"très élevée": -0.10, "élevée": -0.07, "modérée": -0.04, "faible": -0.015, "hors zone": 0.0} # ln P, hypothèses
	145	+STIGMATE = 0.5 # après l'inondation du 1er mai 2019, la décote est multipliée par (1 + STIGMATE)
	146	+AMENITE_RIVE = 0.04 # ln P, lots à moins de 120 m de la rivière OU riverains du lac (vue, accès) ; indépendant du risque
	147	+T_INOND = date(2019, 5, 1)
	148	+DATE_DEBUT, DATE_FIN, DATE_REF = date(2016, 1, 1), date(2026, 8, 31), date(2021, 1, 1)
	149	+PRIX_BASE, TAUX_MARCHE, SIGMA = 390_000.0, 0.035, 0.07
	150	+BETA = {"ln_hab": 0.55, "ln_terr": 0.10, "age": -0.005, "qualite": 0.06, "garage": 0.035}
	151	+
	152	+def generer_riverain(seed: int = SEED) -> tuple[pd.DataFrame, pd.DataFrame]:
	153	+ rng = np.random.default_rng(seed)
	154	+ lots = []
	155	+ for y in range(NLIG):
	156	+ for x in range(NCOL):
	157	+ cl = LIGNES_CLASSE.get(y, "hors zone")
	158	+ dist = (y + 0.5) * PAS + rng.uniform(-8, 8)
	159	+ lots.append({"id_lot": f"RF{y:02d}{x:02d}", "x": x, "y": y, "x_m": (x + 0.5) * PAS, "y_m": (y + 0.5) * PAS,
	160	+ "dist_riviere_m": round(dist, 1), "classe": cl, "expose": int(cl != "hors zone"), "rive": int(dist < 120 or x in COLS_LAC)})
	161	+ lots = pd.DataFrame(lots)
	162	+ # ventes : chaque lot se vend en moyenne 1,6 fois sur la période (Poisson), dates uniformes
	163	+ ventes = []
	164	+ n_jours = (DATE_FIN - DATE_DEBUT).days
	165	+ k = 0
	166	+ for _, L in lots.iterrows():
	167	+ nv = rng.poisson(1.6)
	168	+ hab = float(np.exp(rng.normal(np.log(140), 0.22))); terr = float(np.clip(rng.normal(560, 110), 300, 1200))
	169	+ annee = int(np.clip(rng.normal(1985, 18), 1930, 2016)); qual = int(np.clip(round(rng.normal(3, 0.8)), 1, 5)); gar = int(rng.random() < 0.55)
	170	+ for _ in range(nv):
	171	+ dv = DATE_DEBUT + timedelta(days=int(rng.integers(0, n_jours + 1)))
	172	+ age = dv.year - annee; post = int(dv >= T_INOND)
	173	+ decote = DECOTE_AVANT[L["classe"]] * (1 + STIGMATE * post)
	174	+ lnp = (math.log(PRIX_BASE) + BETA["ln_hab"] * math.log(hab / 140) + BETA["ln_terr"] * math.log(terr / 560) + BETA["age"] * age
	175	+ + BETA["qualite"] * (qual - 3) + BETA["garage"] * gar + AMENITE_RIVE * L["rive"] + decote
	176	+ + math.log(1 + TAUX_MARCHE) * ((dv - DATE_REF).days / 365.25) + rng.normal(0, SIGMA))
	177	+ k += 1
	178	+ ventes.append({"id_vente": f"V{k:04d}", "id_lot": L["id_lot"], "date_vente": dv.isoformat(), "prix_vente": round(math.exp(lnp) / 500) * 500,
	179	+ "sup_habitable_m2": round(hab, 1), "sup_terrain_m2": round(terr), "annee_construction": annee, "age": age, "qualite": qual,
	180	+ "garage": gar, "classe": L["classe"], "expose": L["expose"], "rive": L["rive"], "dist_riviere_m": L["dist_riviere_m"],
	181	+ "post": post, "decote_declaree": round(decote, 5)})
	182	+ ventes = pd.DataFrame(ventes).sort_values("date_vente").reset_index(drop=True)
	183	+ return lots, ventes
	184	+
	185	+def charger_ou_generer(ecrire: bool = True) -> tuple[pd.DataFrame, pd.DataFrame]:
	186	+ fl, fv = D / "ch43_riverain_lots.csv", D / "ch43_riverain_ventes.csv"
	187	+ lots, ventes = generer_riverain()
	188	+ if ecrire and not (fl.exists() and fv.exists()):
	189	+ lots.to_csv(fl, index=False); ventes.to_csv(fv, index=False)
	190	+ return lots, ventes
	191	+
	192	+def _X(v: pd.DataFrame, colonnes_sup: list[pd.Series]) -> np.ndarray:
	193	+ t = (pd.to_datetime(v["date_vente"]) - pd.Timestamp(DATE_REF)).dt.days / 365.25
	194	+ base = [np.ones(len(v)), np.log(v["sup_habitable_m2"] / 140), np.log(v["sup_terrain_m2"] / 560), v["age"].astype(float),
	195	+ (v["qualite"] - 3).astype(float), v["garage"].astype(float), t]
	196	+ return np.column_stack(base + [np.asarray(c, float) for c in colonnes_sup])
	197	+
	198	+def apparier(ventes: pd.DataFrame, fenetre_mois: float = 12.0) -> pd.DataFrame:
	199	+ """Ventes appariées après l'inondation : chaque vente exposée est appariée à la vente hors zone la plus proche
	200	+ (distance d'attributs : superficies, âge, qualité, garage, date ≤ 12 mois) ; témoin ajusté du temps et des écarts
	201	+ physiques par les coefficients hédoniques estimés sur les témoins (ch. 17)."""
	202	+ post = ventes[ventes["post"] == 1].copy(); post["d"] = pd.to_datetime(post["date_vente"])
	203	+ ex = post[post["expose"] == 1]; te = post[post["expose"] == 0].copy()
	204	+ # coefficients physiques sur les témoins seulement (hors zone)
	205	+ r = h.mco(_X(te, []), np.log(te["prix_vente"].values)); b = r["beta"]
	206	+ ech = {"sup_habitable_m2": 25, "sup_terrain_m2": 150, "age": 10, "qualite": 1, "garage": 1}
	207	+ rows = []; pris = set()
	208	+ for _, s in ex.iterrows():
	209	+ cand = te[(te["d"] - s["d"]).abs() <= pd.Timedelta(days=int(fenetre_mois * 30.4375))]
	210	+ cand = cand[~cand["id_vente"].isin(pris)]
	211	+ if cand.empty: continue
	212	+ dist = sum(((cand[k] - s[k]) / e) ** 2 for k, e in ech.items()) + ((cand["d"] - s["d"]).dt.days / 182.6) ** 2
	213	+ tem = cand.loc[dist.idxmin()]; pris.add(tem["id_vente"])
	214	+ mois = (s["d"] - tem["d"]).days / 30.4375
	215	+ aj_temps = cp.ajustement_temps_taux(tem["prix_vente"], TAUX_MARCHE, mois)
	216	+ fact = math.exp(b[1] * math.log(s["sup_habitable_m2"] / tem["sup_habitable_m2"]) + b[2] * math.log(s["sup_terrain_m2"] / tem["sup_terrain_m2"])
	217	+ + b[3] * (s["age"] - tem["age"]) + b[4] * (s["qualite"] - tem["qualite"]) + b[5] * (s["garage"] - tem["garage"]))
	218	+ p_aj = (tem["prix_vente"] + aj_temps) * fact
	219	+ rows.append({"id_expose": s["id_vente"], "id_temoin": tem["id_vente"], "classe": s["classe"], "rive": s["rive"], "prix_expose": s["prix_vente"],
	220	+ "prix_temoin": tem["prix_vente"], "ajust_temps": aj_temps, "facteur_physique": fact, "prix_temoin_ajuste": p_aj,
	221	+ "difference": s["prix_vente"] - p_aj, "difference_pct": s["prix_vente"] / p_aj - 1, "decote_declaree": s["decote_declaree"]})
	222	+ return pd.DataFrame(rows)
	223	+
	224	+def resume_paires(p: pd.DataFrame) -> dict:
	225	+ r = cp.ventes_appariees(p["prix_expose"].tolist(), p["prix_temoin_ajuste"].tolist())
	226	+ pct = p["difference_pct"]; n = len(p)
	227	+ out = {"n": n, "moyenne": r["moyenne"], "mediane": r["mediane"], "ecart_type": r["ecart_type"], "erreur_type": r["ecart_type"] / math.sqrt(n),
	228	+ "moyenne_pct": float(pct.mean()), "mediane_pct": float(pct.median()), "et_pct": float(pct.std(ddof=1) / math.sqrt(n)),
	229	+ "verite_pct": float(np.exp(p["decote_declaree"]).mean() - 1), "verite_log": float(p["decote_declaree"].mean())}
	230	+ out["ic95_pct"] = (out["moyenne_pct"] - 1.96 * out["et_pct"], out["moyenne_pct"] + 1.96 * out["et_pct"])
	231	+ out["par_classe"] = p.groupby("classe")["difference_pct"].agg(["count", "mean"]).rename(columns={"count": "n", "mean": "moyenne_pct"})
	232	+ out["par_rive"] = p.groupby("rive")["difference_pct"].agg(["count", "mean"]).rename(columns={"count": "n", "mean": "moyenne_pct"})
	233	+ return out
	234	+
	235	+def hedoniques(ventes: pd.DataFrame) -> dict:
	236	+ """(a) indicatrice `exposé` sans contrôle de l'aménité ; (b) avec `rive` ; (c) par classe avec interactions classe × après ;
	237	+ (d) DiD par régression exposé × après (avec rive)."""
	238	+ v = ventes; y = np.log(v["prix_vente"].values); out = {}
	239	+ ra = h.mco(_X(v, [v["expose"]]), y)
	240	+ out["a_sans_rive"] = {"coef": ra["beta"][7], "se": ra["se_hc1"][7], "r2": ra["r2"], "sigma": math.sqrt(ra["sigma2"])}
	241	+ rb = h.mco(_X(v, [v["expose"], v["rive"]]), y)
	242	+ out["b_avec_rive"] = {"coef": rb["beta"][7], "se": rb["se_hc1"][7], "coef_rive": rb["beta"][8], "se_rive": rb["se_hc1"][8], "r2": rb["r2"], "sigma": math.sqrt(rb["sigma2"])}
	243	+ cls = ["très élevée", "élevée", "modérée", "faible"]
	244	+ cols = [(v["classe"] == c).astype(float) for c in cls] + [(v["classe"] == c).astype(float) * v["post"] for c in cls] + [v["rive"]]
	245	+ rc = h.mco(_X(v, cols), y)
	246	+ noms = [f"{c}" for c in cls] + [f"{c} × après" for c in cls] + ["rive"]
	247	+ out["c_classes"] = pd.DataFrame({"coef": rc["beta"][7:], "se": rc["se_hc1"][7:]}, index=noms)
	248	+ out["c_classes"]["effet_pct"] = np.exp(out["c_classes"]["coef"]) - 1
	249	+ out["c_classes"]["verite"] = [DECOTE_AVANT[c] for c in cls] + [DECOTE_AVANT[c] * STIGMATE for c in cls] + [AMENITE_RIVE]
	250	+ out["c_classes"]["r2"] = rc["r2"]
	251	+ rd = h.mco(_X(v, [v["expose"], v["expose"] * v["post"], v["rive"]]), y) # `post` seul est colinéaire avec la tendance ? non : on l'ajoute
	252	+ rd = h.mco(_X(v, [v["expose"], v["post"], v["expose"] * v["post"], v["rive"]]), y)
	253	+ out["d_did_reg"] = {"coef_expose": rd["beta"][7], "se_expose": rd["se_hc1"][7], "coef_post": rd["beta"][8], "coef_did": rd["beta"][9], "se_did": rd["se_hc1"][9],
	254	+ "coef_rive": rd["beta"][10], "r2": rd["r2"], "sigma": math.sqrt(rd["sigma2"])}
	255	+ # cellules 2 × 2 sur ln P (sans contrôle)
	256	+ m = v.assign(lnp=y).groupby(["expose", "post"])["lnp"].mean()
	257	+ out["d_did_2x2"] = {"moyennes": m, "estimateur": float((m[(1, 1)] - m[(1, 0)]) - (m[(0, 1)] - m[(0, 0)])),
	258	+ "effectifs": v.groupby(["expose", "post"]).size().to_dict()}
	259	+ # vérités : décote moyenne (log) des ventes exposées avant / après, stigmate moyen
	260	+ ex = v[v["expose"] == 1]
	261	+ out["verite"] = {"decote_avant_log": float(ex.loc[ex["post"] == 0, "decote_declaree"].mean()), "decote_apres_log": float(ex.loc[ex["post"] == 1, "decote_declaree"].mean()),
	262	+ "decote_toutes_log": float(ex["decote_declaree"].mean())}
	263	+ out["verite"]["stigmate_log"] = out["verite"]["decote_apres_log"] - out["verite"]["decote_avant_log"]
	264	+ out["effectifs"] = {"ventes": len(v), "exposees": int(v["expose"].sum()), "rive": int(v["rive"].sum()), "apres": int(v["post"].sum())}
	265	+ return out
	266	+
	267	+def effets_par_methode(ventes: pd.DataFrame, paires: pd.DataFrame) -> pd.DataFrame:
	268	+ """Tableau des estimations de la décote (ln P) des ventes exposées après l'inondation, par méthode, avec IC à 95 %."""
	269	+ he = hedoniques(ventes); rp = resume_paires(paires)
	270	+ lignes = []
	271	+ lignes.append({"methode": "Ventes appariées (après 2019, sans contrôle de l'aménité)", "coef": math.log(1 + rp["moyenne_pct"]), "se": rp["et_pct"] / (1 + rp["moyenne_pct"]), "cible": "décote après"})
	272	+ lignes.append({"methode": "Hédonique, indicatrice exposé (sans « rive »)", "coef": he["a_sans_rive"]["coef"], "se": he["a_sans_rive"]["se"], "cible": "décote moyenne"})
	273	+ lignes.append({"methode": "Hédonique, indicatrice exposé + « rive »", "coef": he["b_avec_rive"]["coef"], "se": he["b_avec_rive"]["se"], "cible": "décote moyenne"})
	274	+ d = he["d_did_reg"]
	275	+ lignes.append({"methode": "DiD régression : exposé (avant)", "coef": d["coef_expose"], "se": d["se_expose"], "cible": "décote avant"})
	276	+ lignes.append({"methode": "DiD régression : exposé × après (stigmate)", "coef": d["coef_did"], "se": d["se_did"], "cible": "stigmate"})
	277	+ lignes.append({"methode": "DiD régression : décote après = somme", "coef": d["coef_expose"] + d["coef_did"], "se": float("nan"), "cible": "décote après"})
	278	+ lignes.append({"methode": "DiD 2 × 2 sur les moyennes de ln P (stigmate)", "coef": he["d_did_2x2"]["estimateur"], "se": float("nan"), "cible": "stigmate"})
	279	+ t = pd.DataFrame(lignes)
	280	+ ver = he["verite"]
	281	+ t["verite"] = t["cible"].map({"décote après": ver["decote_apres_log"], "décote moyenne": ver["decote_toutes_log"], "décote avant": ver["decote_avant_log"], "stigmate": ver["stigmate_log"]})
	282	+ t["effet_pct"] = np.exp(t["coef"]) - 1; t["ic_bas"] = t["coef"] - 1.96 * t["se"]; t["ic_haut"] = t["coef"] + 1.96 * t["se"]
	283	+ return t
	284	+
	285	+# =============================================================================== 4. modèles de valeur
	286	+
	287	+def valeur_terrain_contamine(v_propre: float, cout_rehab: float, delai_mois: float, r: float, stigmate: float, contingence: float = 0.0, perte_revenu_mensuelle: float = 0.0) -> dict:
	288	+ """V_c = V_p (1 − s)(1 + r)^(−T) − C(1 + k)(1 + r)^(−T/2) − Σ pertes de revenu actualisées pendant les travaux.
	289	+ T en années ; C payé à mi-parcours ; pertes mensuelles actualisées au taux périodique équivalent."""
	290	+ T = delai_mois / 12.0; rm = (1 + r) ** (1 / 12) - 1
	291	+ va_propre = v_propre * (1 - stigmate) / (1 + r) ** T
	292	+ va_cout = cout_rehab * (1 + contingence) / (1 + r) ** (T / 2)
	293	+ va_perte = sum(perte_revenu_mensuelle / (1 + rm) ** m for m in range(1, int(round(delai_mois)) + 1))
	294	+ v = va_propre - va_cout - va_perte
	295	+ return {"valeur": v, "va_propre_stigmatisee": va_propre, "va_cout": va_cout, "va_perte_revenu": va_perte, "cout_delai": v_propre * (1 - stigmate) - va_propre,
	296	+ "cout_stigmate": v_propre * stigmate, "ecart_total": v_propre - v, "ecart_pct": v / v_propre - 1}
	297	+
	298	+def valeur_apres_sinistre(v_avant: float, cout_remise: float, loyer_mensuel_perdu: float, mois_perte: float, r: float, stigmate: float, indemnite: float = 0.0) -> dict:
	299	+ """Valeur de l'immeuble sinistré non réparé : V_avant − coût de remise en état − VA de la perte d'usage − stigmate résiduel × V_avant ;
	300	+ position du propriétaire assuré : valeur + indemnité (l'indemnité couvre des coûts, jamais le stigmate)."""
	301	+ rm = (1 + r) ** (1 / 12) - 1
	302	+ va_perte = sum(loyer_mensuel_perdu / (1 + rm) ** m for m in range(1, int(round(mois_perte)) + 1))
	303	+ stig = stigmate * v_avant
	304	+ v = v_avant - cout_remise - va_perte - stig
	305	+ return {"valeur": v, "va_perte_usage": va_perte, "stigmate_dollars": stig, "perte_totale": v_avant - v, "perte_pct": v / v_avant - 1,
	306	+ "position_assure": v + indemnite, "perte_nette_assure": v_avant - (v + indemnite), "part_non_assuree": (v_avant - (v + indemnite)) / (v_avant - v)}
	307	+
	308	+def desuetude_vs_correction(cout_correction: float, perte_usage: float, perte_rne_annuelle: float, R: float, curable: bool) -> dict:
	309	+ """Pathologie curable : D = coût de correction + perte d'usage ; incurable : D = perte de RNE / R (capitalisation) — ch. 24."""
	310	+ d_cur = cout_correction + perte_usage; d_inc = perte_rne_annuelle / R
	311	+ return {"curable": d_cur, "incurable": d_inc, "retenu": d_cur if curable else d_inc, "regle": "curable si coût de correction ≤ gain de valeur"}
	312	+
	313	+def exemple_contamination() -> dict:
	314	+ """Hypothèses déclarées (fictives) : terrain d'un ancien garage avec immeuble de 8 logements à construire ou existant ;
	315	+ valeur propre 2 100 000 $ ; phase I 3 500 $, phase II 18 000 $, réhabilitation 185 000 $ (+ 25 % contingence) ; 15 mois ; 6,5 % ; stigmate 5 %."""
	316	+ hyp = {"v_propre": 2_100_000.0, "phase1": 3_500.0, "phase2": 18_000.0, "rehab": 185_000.0, "contingence": 0.25, "delai_mois": 15, "r": 0.065, "stigmate": 0.05, "perte_mensuelle": 0.0}
	317	+ c_total = hyp["phase1"] + hyp["phase2"] + hyp["rehab"]
	318	+ base = valeur_terrain_contamine(hyp["v_propre"], c_total, hyp["delai_mois"], hyp["r"], hyp["stigmate"], hyp["contingence"], hyp["perte_mensuelle"])
	319	+ sens = {}
	320	+ for s in [0.0, 0.05, 0.10, 0.15]:
	321	+ for T in [6, 15, 24]:
	322	+ sens[(s, T)] = valeur_terrain_contamine(hyp["v_propre"], c_total, T, hyp["r"], s, hyp["contingence"])["valeur"]
	323	+ sens_cout = {k: valeur_terrain_contamine(hyp["v_propre"], c_total * k, hyp["delai_mois"], hyp["r"], hyp["stigmate"], hyp["contingence"])["valeur"] for k in [0.7, 1.0, 1.5, 2.0]}
	324	+ return {"hyp": hyp, "cout_total": c_total, "base": base, "sensibilite": sens, "sensibilite_cout": sens_cout}
	325	+
	326	+def exemple_sinistre() -> dict:
	327	+ """Immeuble P41 du cas 4 (Montréal, zone inondable) : inondation fictive ; hypothèses déclarées."""
	328	+ pf = cas4.charger()["portefeuille"]; p = pf[pf["id_immeuble"] == "P41"].iloc[0]
	329	+ hyp = {"v_avant": float(p["valeur_2026"]), "rne": float(p["rne_2026"]), "nb_logements": int(p["nb_logements"]), "loyer": float(p["loyer_mensuel_moyen"]),
	330	+ "cout_remise": 240_000.0, "logements_touches": 4, "mois_perte": 5, "r": 0.065, "stigmate": 0.05, "franchise": 10_000.0, "part_couverte": 0.85}
	331	+ loyer_perdu = hyp["loyer"] * hyp["logements_touches"]
	332	+ indem = max(0.0, hyp["part_couverte"] * hyp["cout_remise"] - hyp["franchise"])
	333	+ base = valeur_apres_sinistre(hyp["v_avant"], hyp["cout_remise"], loyer_perdu, hyp["mois_perte"], hyp["r"], hyp["stigmate"], indem)
	334	+ sens = {s: valeur_apres_sinistre(hyp["v_avant"], hyp["cout_remise"], loyer_perdu, hyp["mois_perte"], hyp["r"], s, indem)["valeur"] for s in [0.0, 0.025, 0.05, 0.10, 0.15]}
	335	+ sans_assurance = valeur_apres_sinistre(hyp["v_avant"], hyp["cout_remise"], loyer_perdu, hyp["mois_perte"], hyp["r"], hyp["stigmate"], 0.0)
	336	+ return {"hyp": hyp, "loyer_perdu_mensuel": loyer_perdu, "indemnite": indem, "base": base, "sensibilite_stigmate": sens, "sans_assurance": sans_assurance, "immeuble": p}
	337	+
	338	+def exemple_pyrite() -> dict:
	339	+ """Pyrite (renvoi ch. 11 et 24) : excavation 25 000 $, perte d'usage 2 000 $ ; variante incurable : perte de RNE 1 800 $/an capitalisée à 5,5 %."""
	340	+ return {"curable": desuetude_vs_correction(25_000, 2_000, 0, 0.055, True), "incurable": desuetude_vs_correction(0, 0, 1_800, 0.055, False)}
	341	+
	342	+# =============================================================================== main
	343	+
	344	+def main(ecrire: bool = True, verbeux: bool = True) -> dict:
	345	+ pr = prime_taux_cas4(); pan = panel_choc(); dd = did_cas4(pan)
	346	+ lots, ventes = charger_ou_generer(ecrire); paires = apparier(ventes); rp = resume_paires(paires); he = hedoniques(ventes); epm = effets_par_methode(ventes, paires)
	347	+ ct = exemple_contamination(); sn = exemple_sinistre(); py = exemple_pyrite()
	348	+ res = {"prime": pr, "panel": pan, "did": dd, "lots": lots, "ventes": ventes, "paires": paires, "resume_paires": rp, "hedoniques": he, "effets_methodes": epm,
	349	+ "contamination": ct, "sinistre": sn, "pyrite": py}
	350	+ if not verbeux: return res
	351	+ P = print
	352	+ P("=== 1. Cas 4 — prime de risque dans le taux global (effets_regionaux) ===")
	353	+ P(f"n = {pr['n']} immeubles ; inondables = {pr['n_inondable']} ; contaminés = {pr['n_sol']} ; part de la valeur en zone inondable = {pr['part_valeur_inondable']*100:.2f} %")
	354	+ for v in ["inondable", "sol"]:
	355	+ e = pr[v]
	356	+ P(f"{v:10s}: coef = {e['coef']:.4f} (se HC1 {e['se']:.4f}) → +{e['effet_pct']*100:.2f} % du taux [IC95 {e['ic_bas_pct']*100:.2f} ; {e['ic_haut_pct']*100:.2f}] ; Kennedy {e['effet_kennedy_pct']*100:.2f} % ; déclaré +{e['declare_pct']*100:.0f} % (t = {e['t_vs_declare']:.2f})")
	357	+ P(f"décote de valeur impliquée (à RNE donné) : {pr['decote_valeur_inondable_pct']*100:.2f} % (déclarée {pr['decote_valeur_declaree_pct']*100:.2f} %)")
	358	+ P(pr["immeubles_inondables"].to_string(index=False))
	359	+ P("\n=== 2. Cas 4 — DiD sur ln RNE, choc 2019 injecté ===")
	360	+ P(f"vérité : transitoire {CHOC_2019} ; permanent {CHOC_PERMANENT} ({dd['verite_permanente_pct']*100:.2f} %)")
	361	+ P("cellules (ln RNE moyen) :"); P(dd["cellules"].to_string())
	362	+ P(f"DiD 2×2 naïf = {dd['estimateur_2x2']:.4f} ; régression Δ (se HC1) = {dd['coef_reg']:.4f} ({dd['se_hc1']:.4f}) ; par strate RMR = {dd['estimateur_strate']:.4f}")
	363	+ P(f"inférence de randomisation exacte : {dd['n_combinaisons']} affectations ; p (naïf) = {dd['p_ri']:.4f} ; p (strate) = {dd['p_ri_strate']:.4f}")
	364	+ P("TWFE immeuble + RMR×trimestre :"); P(dd["twfe"].to_string())
	365	+ P("étude d'événement (k, coef, se grappes) :"); P(dd["event_study"][["coef", "se_grappes", "ic_bas", "ic_haut"]].round(4).to_string())
	366	+ P(f"effet sur la valeur (Δln RNE − Δln R) : estimé {dd['effet_valeur_estime_pct']*100:.2f} % ; déclaré {dd['effet_valeur_declare_pct']*100:.2f} %")
	367	+ P("\n=== 3. Rivière-Fictive — ventes riveraines ===")
	368	+ P(f"lots {len(lots)} ; exposés {int(lots['expose'].sum())} ; par classe {lots['classe'].value_counts().to_dict()}")
	369	+ P(f"ventes {he['effectifs']} ; vérités (ln) : avant {he['verite']['decote_avant_log']:.4f}, après {he['verite']['decote_apres_log']:.4f}, stigmate {he['verite']['stigmate_log']:.4f}")
	370	+ P(f"paires : n = {rp['n']} ; moyenne {rp['moyenne']:.0f} $ ({rp['moyenne_pct']*100:.2f} %) ; médiane {rp['mediane_pct']*100:.2f} % ; IC95 [{rp['ic95_pct'][0]*100:.2f} ; {rp['ic95_pct'][1]*100:.2f}] ; vérité {rp['verite_pct']*100:.2f} %")
	371	+ P(rp["par_classe"].to_string()); P(rp["par_rive"].to_string())
	372	+ P(f"hédonique (a) exposé sans rive : {he['a_sans_rive']['coef']:.4f} ({he['a_sans_rive']['se']:.4f}) ; (b) avec rive : {he['b_avec_rive']['coef']:.4f} ({he['b_avec_rive']['se']:.4f}), rive {he['b_avec_rive']['coef_rive']:.4f} ({he['b_avec_rive']['se_rive']:.4f}) ; σ = {he['b_avec_rive']['sigma']:.4f}, R² = {he['b_avec_rive']['r2']:.3f}")
	373	+ P("(c) par classe :"); P(he["c_classes"].round(4).to_string())
	374	+ d = he["d_did_reg"]; P(f"(d) DiD régression : exposé {d['coef_expose']:.4f} ({d['se_expose']:.4f}) ; exposé×après {d['coef_did']:.4f} ({d['se_did']:.4f}) ; post {d['coef_post']:.4f} ; rive {d['coef_rive']:.4f}")
	375	+ P(f" DiD 2×2 : {he['d_did_2x2']['estimateur']:.4f} ; effectifs {he['d_did_2x2']['effectifs']}")
	376	+ P("effets par méthode :"); P(epm.round(4).to_string(index=False))
	377	+ P("\n=== 4. Modèles de valeur ===")
	378	+ b = ct["base"]; P(f"contamination : coût total {ct['cout_total']:.0f} $ ; V_c = {b['valeur']:.0f} $ ({b['ecart_pct']*100:.2f} %) ; VA propre stigmatisée {b['va_propre_stigmatisee']:.0f} ; VA coût {b['va_cout']:.0f} ; coût du délai {b['cout_delai']:.0f} ; stigmate {b['cout_stigmate']:.0f}")
	379	+ P("sensibilité (stigmate, mois) → valeur :"); P({k: round(v) for k, v in ct["sensibilite"].items()}); P("coût ×k :"); P({k: round(v) for k, v in ct["sensibilite_cout"].items()})
	380	+ s = sn["base"]; P(f"sinistre P41 : V_avant {sn['hyp']['v_avant']:.0f} ; loyer perdu {sn['loyer_perdu_mensuel']:.0f} $/mois × {sn['hyp']['mois_perte']} ; VA perte {s['va_perte_usage']:.0f} ; stigmate {s['stigmate_dollars']:.0f} ; V_après {s['valeur']:.0f} ({s['perte_pct']*100:.2f} %) ; indemnité {sn['indemnite']:.0f} ; position assuré {s['position_assure']:.0f} ; perte nette {s['perte_nette_assure']:.0f} ; part non assurée {s['part_non_assuree']*100:.1f} %")
	381	+ P("sensibilité stigmate :"); P({k: round(v) for k, v in sn["sensibilite_stigmate"].items()})
	382	+ P(f"pyrite : curable {py['curable']['retenu']:.0f} $ ; incurable {py['incurable']['retenu']:.0f} $")
	383	+ return res
	384	+
	385	+if __name__ == "__main__":
	386	+ main()

modified analysis/python/figures/ch20_figures.py +11 −10

		@@ -68,18 +68,19 @@ def fig_residus_elargissement(r):
68	68
69	69
70	70	def fig_fourchettes(r):
71		− b, bs, co, cb, vt = r["base"], r["bootstrap"], r["cout"], r["combinaison"], r["verite"]
	71	+ b, bs, co, cb, vt, bk = r["base"], r["bootstrap"], r["cout"], r["combinaison"], r["verite"], r["bootstrap_cout"]
72	72	ip = b["ip_hed"]
73	73	z90 = 1.6449
74	74	lignes = [
75	75	("Prix de vente probable\n(intervalle de prédiction à 95 %, ch. 18)", b["V_hed"], ip["bas"], ip["haut"], palette.cat[3], 0.5),
76		− ("Coût (hypothèses α, vie économique :\nétendue de la grille)", co["valeur"], co["min_grille"], co["max_grille"], palette.cat[2], 1.0),
	76	+ ("Coût, lecture ventilée (ch. 11), écartée\n(point)", co["ventilee"], np.nan, np.nan, palette.cat[2], 0.6),
	77	+ ("Coût, dépréciation extraite (ch. 24)\n(bootstrap, 90 %)", co["valeur"], bk["resume"]["q05"], bk["resume"]["q95"], palette.cat[2], 1.0),
77	78	("Comparaison, 5 comparables\n(bootstrap conjoint, 90 %)", b["V_comp"], bs["resume"]["conjoint"]["q05"], bs["resume"]["conjoint"]["q95"], palette.cat[1], 1.0),
78	79	("Comparaison élargie, 16 candidats\n(point)", r["elargissement"]["valeur_ponderee"].iloc[-1], np.nan, np.nan, palette.cat[1], 1.0),
79	80	("Hédonique (estimation + spécification,\n± 1,645 écart-type)", b["V_hed"], b["V_hed"] - z90 * cb["s"][1], b["V_hed"] + z90 * cb["s"][1], palette.cat[0], 1.0),
80		− ("Conclusion : valeur la plus probable,\nfourchette conclue et bootstrap 90 %", cb["conclusion"], cb["fourchette_conclue"][0], cb["fourchette_conclue"][1], palette.brand, 1.0),
	81	+ ("Conclusion : valeur la plus probable,\nfourchette conclue (bande : bootstrap 90 %)", cb["conclusion"], cb["fourchette_conclue"][0], cb["fourchette_conclue"][1], palette.brand, 1.0),
81	82	]
82		− fig, ax = plt.subplots(figsize=(LARGEUR, 3.6))
	83	+ fig, ax = plt.subplots(figsize=(LARGEUR, 3.9))
83	84	y = np.arange(len(lignes))[::-1]
84	85	for yi, (lab, pt, lo, hi, col, alpha) in zip(y, lignes):
85	86	if not np.isnan(lo):
		@@ -94,7 +95,7 @@ def fig_fourchettes(r):
94	95	ax.axvline(vt["V_true"], color=palette.ink2, lw=0.8, ls=":", zorder=0)
95	96	ax.text(vt["V_true"] + 3000, -0.55, "espérance de la loi\nde génération", fontsize=7, color=palette.ink2, ha="left", va="center")
96	97	ax.set_ylim(-0.9, len(lignes) - 0.3); ax.set_yticks(y); ax.set_yticklabels([l[0] for l in lignes], fontsize=7.2)
97		− ax.set_xlim(270000, 770000); ax.xaxis.set_major_formatter(FuncFormatter(fmt_k)); ax.set_xlabel("Dollars (synthétiques), 1er septembre 2026")
	98	+ ax.set_xlim(420000, 770000); ax.xaxis.set_major_formatter(FuncFormatter(fmt_k)); ax.set_xlabel("Dollars (synthétiques), 1er septembre 2026")
98	99	ax.grid(axis="y", visible=False)
99	100	fig.tight_layout()
100	101	savefig(fig, "ch20/ch20-fig-fourchettes-methodes.pdf")
		@@ -106,18 +107,18 @@ def fig_sensibilite_poids(r):
106	107	ax = axes[0]
107	108	for i, (wc, courbe) in enumerate(cb["courbes"].items()):
108	109	lab = "poids du coût = " + (f"{wc*100:.0f} %")
109		− ax.plot(cb["wh"], courbe, color=palette.cat[i], lw=1.6 if wc == 0.05 else 1.1, ls="-" if wc == 0.05 else "--", label=lab)
	110	+ ax.plot(cb["wh"], courbe, color=palette.cat[i], lw=1.6 if wc == 0.15 else 1.1, ls="-" if wc == 0.15 else "--", label=lab)
110	111	tab = cb["tableau"]
111	112	marqueurs = {"raisonnée (retenue)": ("o", palette.brand), "inverse de la variance (indépendance)": ("s", palette.cat[0]), "variance minimale avec corrélation (MCG)": ("^", palette.cat[0]), "égale (un tiers chacune)": ("D", palette.cat[2])}
112	113	for _, row in tab.iterrows():
113	114	if row["schema"] in marqueurs:
114	115	m, c = marqueurs[row["schema"]]
115	116	ax.plot(row["w_hed"], row["valeur"], m, color=c, ms=6, mec="white", mew=0.6, zorder=5)
116		− ax.annotate("raisonnée\n(0,40 / 0,55 / 0,05)", (0.55, tab.loc[tab["schema"] == "raisonnée (retenue)", "valeur"].iloc[0]), xytext=(0.62, 572000), fontsize=7, arrowprops=dict(arrowstyle="-", color=palette.ink2, lw=0.6))
117		− ax.annotate("inverse de la variance\net MCG", (0.70, tab.loc[tab["schema"] == "inverse de la variance (indépendance)", "valeur"].iloc[0]), xytext=(0.42, 628000), fontsize=7, arrowprops=dict(arrowstyle="-", color=palette.ink2, lw=0.6))
118		− ax.annotate("égale", (1 / 3, tab.loc[tab["schema"] == "égale (un tiers chacune)", "valeur"].iloc[0]), xytext=(0.12, 540000), fontsize=7, arrowprops=dict(arrowstyle="-", color=palette.ink2, lw=0.6))
	117	+ ax.annotate("raisonnée\n(0,40 / 0,45 / 0,15)", (0.45, tab.loc[tab["schema"] == "raisonnée (retenue)", "valeur"].iloc[0]), xytext=(0.62, 588000), fontsize=7, arrowprops=dict(arrowstyle="-", color=palette.ink2, lw=0.6))
	118	+ ax.annotate("inverse de la variance\net MCG", (0.616, tab.loc[tab["schema"] == "inverse de la variance (indépendance)", "valeur"].iloc[0]), xytext=(0.30, 630000), fontsize=7, arrowprops=dict(arrowstyle="-", color=palette.ink2, lw=0.6))
	119	+ ax.annotate("égale (un tiers au coût)", (1 / 3, tab.loc[tab["schema"] == "égale (un tiers chacune)", "valeur"].iloc[0]), xytext=(0.03, 622000), fontsize=7, arrowprops=dict(arrowstyle="-", color=palette.ink2, lw=0.6))
119	120	ax.set_xlabel("Poids de l'indication hédonique $w_h$ (le reste à la comparaison)"); ax.set_ylabel("Valeur combinée")
120		− ax.yaxis.set_major_formatter(FuncFormatter(fmt_k)); ax.set_xlim(0, 1); ax.legend(loc="lower right", fontsize=7, frameon=False)
	121	+ ax.yaxis.set_major_formatter(FuncFormatter(fmt_k)); ax.set_xlim(0, 1); ax.legend(loc="upper left", fontsize=7, frameon=False)
121	122	ax.set_title("(a) Conclusion selon les poids", loc="left", fontsize=8.5)
122	123	ax = axes[1]
123	124	s = cb["s"]

modified analysis/python/figures/ch36_figures.py +2 −2

		@@ -45,13 +45,13 @@ def fig_event_study(res):
45	45	ax.fill_between(es.index, es["ic_bas"], es["ic_haut"], color=palette.cat[0], alpha=0.18, lw=0)
46	46	ax.plot(es.index, es["coef"], "o-", color=palette.cat[0], ms=3.8, lw=1.4, label="Coefficient estimé (IC à 95 %, grappes)")
47	47	ax.step(es.index, es["vrai"], where="mid", color=palette.ink, lw=1.1, ls="--", label="Effet déclaré δ(k)")
48		− ax.set_title("(a) Étude d'événement TWFE (RMR × trimestre)", fontsize=8.5); ax.set_xlabel("Trimestres depuis la mise en service, k")
	48	+ ax.set_title("(a) Étude d'événement TWFE (RMR × trimestre)", fontsize=8); ax.set_xlabel("Trimestres depuis la mise en service, k")
49	49	ax.set_ylabel("Effet sur ln RNE (points de log)")
50	50	ax = axes[1]
51	51	ax.fill_between(cs.index, cs["ic_bas"], cs["ic_haut"], color=palette.cat[1], alpha=0.18, lw=0)
52	52	ax.plot(cs.index, cs["att"], "s-", color=palette.cat[1], ms=3.4, lw=1.4, label="ATT(k) de Callaway–Sant'Anna (IC bootstrap)")
53	53	kk = np.arange(-8, 16); ax.step(kk, c36.effet_vrai(kk), where="mid", color=palette.ink, lw=1.1, ls="--", label="Effet déclaré δ(k)")
54		− ax.set_title("(b) Callaway–Sant'Anna, témoins de la même RMR", fontsize=8.5); ax.set_xlabel("Trimestres depuis la mise en service, k")
	54	+ ax.set_title("(b) Callaway–Sant'Anna (témoins de la même RMR)", fontsize=8); ax.set_xlabel("Trimestres depuis la mise en service, k")
55	55	for ax in axes:
56	56	ax.axhline(0, color=palette.ink2, lw=0.8); ax.axvline(-0.5, color=palette.ink2, lw=0.8, ls=":"); ax.legend(fontsize=6.8, loc="upper left"); ax.grid(True, lw=0.4)
57	57	savefig(fig, "ch36/ch36-fig-event-study.pdf")

modified analysis/python/figures/ch38_figures.py +14 −12

		@@ -64,7 +64,7 @@ def fig_carte_cas1(r):
64	64	_km(ax, xy); ax.set_xlabel("Est, NAD83 / MTM zone 9 (km) — coordonnées synthétiques"); ax.set_ylabel("Nord (km)")
65	65	leg = [Patch(color=palette.seq[5], alpha=0.6, label="Zone d'exposition fictive, classe A (≤ 350 m)"), Patch(color=palette.seq[1], alpha=0.55, label="Classe B (350–900 m)"),
66	66	Line2D([], [], color=palette.cat[1], lw=1.4, label="Corridor de transport fictif et tampon de 300 m"),
67		− Line2D([], [], marker="^", color="none", markerfacecolor=palette.cat[2], markeredgecolor=palette.ink, markersize=7, label="Pôles fictifs (taille ∝ masse)"),
	67	+ Line2D([], [], marker="^", color="none", markerfacecolor=palette.cat[2], markeredgecolor=palette.ink, markersize=7, label="Pôles fictifs (taille selon la masse)"),
68	68	Line2D([], [], marker="o", color="none", markerfacecolor=palette.seq[9], markersize=5, label="Vente en classe A"),
69	69	Line2D([], [], marker="o", color="none", markerfacecolor=palette.seq[4], markersize=5, label="Vente en classe B"),
70	70	Line2D([], [], marker="o", color="none", markerfacecolor=palette.ink2, markersize=4, alpha=0.6, label="Autre vente"),
		@@ -80,25 +80,25 @@ def _carte_classes(ax, cells, z, bornes, cmap_colors, titre, fmt=lambda x: fmt_m
80	80	ax.fill(ext[:, 0], ext[:, 1], color=cmap(normc(val)), lw=0.3, ec=palette.surface)
81	81	ax.set_title(titre, fontsize=8.5)
82	82	handles = [Patch(color=cmap(i), label=f"{fmt(bornes[i])} – {fmt(bornes[i + 1])}") for i in range(len(bornes) - 1)]
83		− ax.legend(handles=handles, fontsize=6, loc="lower right", frameon=True, framealpha=0.9, edgecolor="none", facecolor=palette.surface, title_fontsize=6.5)
	83	+ ax.legend(handles=handles, fontsize=5.8, loc="upper left", frameon=True, framealpha=0.9, edgecolor="none", facecolor=palette.surface, handlelength=1.2, borderpad=0.3, labelspacing=0.25)
84	84
85	85
86	86	def fig_choroplethe(r):
87	87	g = r["grille"]; ok = g["cartographiees"]; occ = g["occupees"]; b = r["base"]; xy = b["xy_m"] / 1000
88	88	z = ok["prix_m2_med"].values; k = c38.K_CLASSES
89	89	seq = [palette.seq[i] for i in np.linspace(0, 9, k).round().astype(int)]
90		− fig, axes = plt.subplots(2, 2, figsize=(LARGEUR, 6.4), sharex=True, sharey=True)
91		− titres = {"intervalles égaux": "(a) Intervalles égaux", "quantiles": "(b) Quantiles", "Jenks": "(c) Ruptures naturelles (Fisher–Jenks)"}
	90	+ fig, axes = plt.subplots(2, 2, figsize=(LARGEUR, 4.0), sharex=True, sharey=True)
	91	+ titres = {"intervalles égaux": "(a) Intervalles égaux", "quantiles": "(b) Quantiles", "Jenks": "(c) Fisher–Jenks"}
92	92	for ax, (nom, cl) in zip(axes.flat[:3], g["classif"].items()):
93	93	for geom in occ.loc[occ["n"] < c38.MIN_VENTES_CELLULE].geometry:
94	94	ext = np.asarray(geom.exterior.coords) / 1000; ax.fill(ext[:, 0], ext[:, 1], facecolor="none", edgecolor=palette.ink2, hatch="////", lw=0.3)
95		− _carte_classes(ax, ok, z, cl["bornes"], seq, f"{titres[nom]} — GVF = {cl['gvf']:.3f}")
	95	+ _carte_classes(ax, ok, z, cl["bornes"], seq, f"{titres[nom]}, GVF = {cl['gvf']:.3f}".replace(".", ","))
96	96	ax.scatter(xy[:, 0], xy[:, 1], s=1.2, color=palette.ink, alpha=0.25, lw=0)
97	97	ax = axes.flat[3]; res = ok["residu_moy"].values; vmax = float(np.abs(res).max()); bornes = np.linspace(-vmax, vmax, 7)
98	98	div = [palette.div["neg"], "#7fa8e0", "#c5d4ea", "#ead0cf", "#eb9a99", palette.div["pos"]]
99	99	for geom in occ.loc[occ["n"] < c38.MIN_VENTES_CELLULE].geometry:
100	100	ext = np.asarray(geom.exterior.coords) / 1000; ax.fill(ext[:, 0], ext[:, 1], facecolor="none", edgecolor=palette.ink2, hatch="////", lw=0.3)
101		− _carte_classes(ax, ok, res, bornes, div, "(d) Résidu moyen du modèle de référence (divergent)", fmt=lambda x: f"{x:+.3f}".replace(".", ","))
	101	+ _carte_classes(ax, ok, res, bornes, div, "(d) Résidu moyen, modèle de référence", fmt=lambda x: f"{x:+.3f}".replace(".", ","))
102	102	ax.scatter(xy[:, 0], xy[:, 1], s=1.2, color=palette.ink, alpha=0.25, lw=0)
103	103	for ax in axes.flat:
104	104	_km(ax, xy, marge=1.0)
		@@ -106,8 +106,8 @@ def fig_choroplethe(r):
106	106	ax.set_xlabel("Est, MTM 9 (km)")
107	107	for ax in axes[:, 0]:
108	108	ax.set_ylabel("Nord (km)")
109		− fig.text(0.5, 0.005, f"Cellules hexagonales de {g['aire_cellule_km2']:.2f} km² ; hachures : cellules de moins de {c38.MIN_VENTES_CELLULE} ventes (non cartographiées) ; points : ventes", ha="center", fontsize=7, color=palette.ink2)
110		− fig.subplots_adjust(hspace=0.18, wspace=0.08)
	109	+ fig.text(0.5, 0.0, f"Cellules hexagonales de {g['aire_cellule_km2']:.2f} km² ; hachures : cellules de moins de {c38.MIN_VENTES_CELLULE} ventes (non cartographiées) ; points : ventes", ha="center", fontsize=7, color=palette.ink2)
	110	+ fig.subplots_adjust(hspace=0.28, wspace=0.06, bottom=0.12, top=0.95)
111	111	savefig(fig, "ch38/ch38-fig-choroplethe.pdf")
112	112
113	113
		@@ -122,7 +122,7 @@ def fig_faux_positifs(r):
122	122	for s in (-1.96, 1.96):
123	123	ax.axvline(s, color=palette.ink2, lw=0.7, ls=":")
124	124	ax.set_xlabel("Statistique $t$ de l'indicatrice d'un disque aléatoire de 1 km"); ax.set_ylabel("Densité"); ax.set_title("(a) Distribution des $t$ (346 disques)", fontsize=8.5)
125		− ax.legend(fontsize=6.8, loc="upper left")
	125	+ ax.legend(fontsize=6.8, loc="upper right")
126	126	ax = axes[1]
127	127	for i, spec in enumerate(["référence", "sans secteur", "sans secteur ni gradient"]):
128	128	w = tab[(tab["spec"] == spec) & (tab["erreur_type"] == "classique")].sort_values("rayon_km")
		@@ -161,13 +161,15 @@ def fig_cas4(r):
161	161	ax.scatter(x[m2], y[m2], marker="s", s=s[m2] + 70, facecolor="none", edgecolor=palette.cat[3], lw=1.3, zorder=4)
162	162	for rmr, w in g.groupby("rmr"):
163	163	cx, cy = w.geometry.x.mean() / 1000, w.geometry.y.mean() / 1000
164		− ax.annotate(f"{rmr} ({len(w)})", (cx, cy), xytext=(0, 16 if rmr != "Montréal" else -22), textcoords="offset points", fontsize=7, ha="center", color=palette.ink2)
	164	+ dx, dy, ha = (0, 14, "center") if rmr not in ("Montréal", "Sherbrooke") else ((0, -20, "center") if rmr == "Montréal" else (14, -4, "left"))
	165	+ ax.annotate(f"{rmr} ({len(w)})", (cx, cy), xytext=(dx, dy), textcoords="offset points", fontsize=7, ha=ha, color=palette.ink2)
165	166	ax.set_aspect("equal"); ax.grid(True, lw=0.4); ax.set_xlabel("Est, NAD83 / Québec Lambert (km)"); ax.set_ylabel("Nord (km)")
166		− ax.set_title("(a) 48 immeubles, taille ∝ valeur 2026", fontsize=8.5)
	167	+ ax.set_title("(a) 48 immeubles, taille selon la valeur 2026", fontsize=8.5)
167	168	leg = [Line2D([], [], marker="o", color="none", markerfacecolor=palette.cat[0], alpha=0.6, markersize=7, label="Immeuble"),
168	169	Line2D([], [], marker="o", color="none", markerfacecolor="none", markeredgecolor=palette.cat[5], markeredgewidth=1.6, markersize=9, label="Zone inondable (3)"),
169	170	Line2D([], [], marker="s", color="none", markerfacecolor="none", markeredgecolor=palette.cat[3], markeredgewidth=1.3, markersize=9, label="Contamination du sol (2)")]
170		− ax.legend(handles=leg, fontsize=6.8, loc="lower right")
	171	+ ax.legend(handles=leg, fontsize=6.8, loc="upper left")
	172	+ ax.set_ylim(y.min() - 55, y.max() + 60); ax.set_xlim(x.min() - 60, x.max() + 90)
171	173	ax = axes[1]; ordre = par.index[::-1]; yy = np.arange(len(ordre))
172	174	ax.barh(yy, par.loc[ordre, "part_inond_rmr"] * 100, color=palette.cat[5], alpha=0.85, height=0.6)
173	175	for yi, rmr in zip(yy, ordre):

added analysis/python/figures/ch39_figures.py +146 −0

Large diff (146 lines) is not loaded by default to keep the page fast.Load diff

added analysis/python/figures/ch40_figures.py +139 −0

Large diff (139 lines) is not loaded by default to keep the page fast.Load diff

added analysis/python/figures/ch42_bloc_piece.py +38 −0

Large diff (38 lines) is not loaded by default to keep the page fast.Load diff

added analysis/python/figures/ch42_cotisation_speciale.py +36 −0

Large diff (36 lines) is not loaded by default to keep the page fast.Load diff

added analysis/python/figures/ch42_grille.py +41 −0

Large diff (41 lines) is not loaded by default to keep the page fast.Load diff

added analysis/python/figures/ch42_medianes.py +42 −0

Large diff (42 lines) is not loaded by default to keep the page fast.Load diff

added analysis/python/figures/ch43_carte_cas4.py +38 −0

Large diff (38 lines) is not loaded by default to keep the page fast.Load diff

added analysis/python/figures/ch43_carte_riverain.py +52 −0

Large diff (52 lines) is not loaded by default to keep the page fast.Load diff

added analysis/python/figures/ch43_did_cas4.py +46 −0

Large diff (46 lines) is not loaded by default to keep the page fast.Load diff

added analysis/python/figures/ch43_effets_methodes.py +50 −0

Large diff (50 lines) is not loaded by default to keep the page fast.Load diff

added analysis/python/figures/ch43_placebos.py +32 −0

Large diff (32 lines) is not loaded by default to keep the page fast.Load diff

modified analysis/python/tests/test_ch20.py +45 −27

Large diff (105 lines) is not loaded by default to keep the page fast.Load diff

added analysis/python/tests/test_ch42.py +206 −0

Large diff (206 lines) is not loaded by default to keep the page fast.Load diff

modified assets/figures/ch20/ch20-fig-bootstrap-densite.pdf +0 −0

Binary file not shown.

modified assets/figures/ch20/ch20-fig-fourchettes-methodes.pdf +0 −0

Binary file not shown.

modified assets/figures/ch20/ch20-fig-residus-elargissement.pdf +0 −0

Binary file not shown.

modified assets/figures/ch20/ch20-fig-sensibilite-poids.pdf +0 −0

Binary file not shown.

modified assets/figures/ch38/ch38-fig-carte-cas1.pdf +0 −0

Binary file not shown.

modified assets/figures/ch38/ch38-fig-cas4.pdf +0 −0

Binary file not shown.

modified assets/figures/ch38/ch38-fig-choroplethe.pdf +0 −0

Binary file not shown.

modified assets/figures/ch38/ch38-fig-faux-positifs.pdf +0 −0

Binary file not shown.

modified assets/figures/ch38/ch38-fig-tampon.pdf +0 −0

Binary file not shown.

added assets/figures/ch39/ch39-fig-posterieures.pdf +0 −0

Binary file not shown.

added assets/figures/ch39/ch39-fig-predictive.pdf +0 −0

Binary file not shown.

added assets/figures/ch39/ch39-fig-prior-informatif.pdf +0 −0

Binary file not shown.

added assets/figures/ch39/ch39-fig-retrecissement.pdf +0 −0

Binary file not shown.

added assets/figures/ch39/ch39-fig-traces.pdf +0 −0

Binary file not shown.

added assets/figures/ch40/ch40-fig-convergence.pdf +0 −0

Binary file not shown.

added assets/figures/ch40/ch40-fig-distribution.pdf +0 −0

Binary file not shown.

added assets/figures/ch40/ch40-fig-lhs.pdf +0 −0

Binary file not shown.

added assets/figures/ch40/ch40-fig-nuage.pdf +0 −0

Binary file not shown.

added assets/figures/ch40/ch40-fig-portefeuille.pdf +0 −0

Binary file not shown.

added assets/figures/ch40/ch40-fig-rcd.pdf +0 −0

Binary file not shown.

added assets/figures/ch40/ch40-fig-sobol.pdf +0 −0

Binary file not shown.

added assets/figures/ch42/ch42-fig-bloc-piece.pdf +0 −0

Binary file not shown.

added assets/figures/ch42/ch42-fig-cotisation-speciale.pdf +0 −0

Binary file not shown.

added assets/figures/ch42/ch42-fig-grille-fraction.pdf +0 −0

Binary file not shown.

added assets/figures/ch42/ch42-fig-medianes-rmr.pdf +0 −0

Binary file not shown.

added assets/figures/ch43/ch43-fig-carte-cas4.pdf +0 −0

Binary file not shown.

added assets/figures/ch43/ch43-fig-carte-riverain.pdf +0 −0

Binary file not shown.

added assets/figures/ch43/ch43-fig-did-cas4.pdf +0 −0

Binary file not shown.

added assets/figures/ch43/ch43-fig-effets-methodes.pdf +0 −0

Binary file not shown.

added assets/figures/ch43/ch43-fig-placebos.pdf +0 −0

Binary file not shown.

added data/synthetic/ch43_riverain_lots.csv +626 −0

Large diff (626 lines) is not loaded by default to keep the page fast.Load diff

added data/synthetic/ch43_riverain_ventes.csv +765 −0

Large diff (765 lines) is not loaded by default to keep the page fast.Load diff

modified manuscript/chapters/ch34.tex +12 −12

Large diff (70 lines) is not loaded by default to keep the page fast.Load diff

modified manuscript/chapters/ch35.tex +30 −27

Large diff (142 lines) is not loaded by default to keep the page fast.Load diff

added references/lib-ch34.bib +16 −0

Large diff (16 lines) is not loaded by default to keep the page fast.Load diff

added references/lib-ch35.bib +23 −0

Large diff (23 lines) is not loaded by default to keep the page fast.Load diff

modified references/library.bib +18 −0

Large diff (21 lines) is not loaded by default to keep the page fast.Load diff

modified solutions/ch35.tex +38 −0

Large diff (39 lines) is not loaded by default to keep the page fast.Load diff

modified sources/registry.csv +6 −0

Large diff (10 lines) is not loaded by default to keep the page fast.Load diff

modified verification/calculations.csv +13 −0

Large diff (16 lines) is not loaded by default to keep the page fast.Load diff

added verification/chapters/ch34-claims.csv +24 −0

Large diff (24 lines) is not loaded by default to keep the page fast.Load diff

added verification/chapters/ch34-registry.csv +4 −0

Large diff (4 lines) is not loaded by default to keep the page fast.Load diff

added verification/chapters/ch34.md +35 −0

Large diff (35 lines) is not loaded by default to keep the page fast.Load diff

added verification/chapters/ch35-calculations.csv +14 −0

Large diff (14 lines) is not loaded by default to keep the page fast.Load diff

added verification/chapters/ch35-claims.csv +13 −0

Large diff (13 lines) is not loaded by default to keep the page fast.Load diff

added verification/chapters/ch35-figures.csv +6 −0

Large diff (6 lines) is not loaded by default to keep the page fast.Load diff

added verification/chapters/ch35-registry.csv +4 −0

Large diff (4 lines) is not loaded by default to keep the page fast.Load diff

modified verification/claims.csv +35 −0

Large diff (38 lines) is not loaded by default to keep the page fast.Load diff

modified verification/figures.csv +5 −0

Large diff (9 lines) is not loaded by default to keep the page fast.Load diff