spb/traite-evaluation-qc
Public
1#!/usr/bin/env python32"""Chapitre 35 — Indices immobiliers et ventes répétées : exemple intégral sur le cas 1 (Gatineau, DONNÉES SYNTHÉTIQUES, graine 20261031).34Construit et compare, sur les 600 ventes du cas 1 (data/synthetic/cas1_gatineau_ventes.csv, 15 trimestres 2023T1-2026T3) :5 * l'indice de la médiane et de la moyenne géométrique des prix (biais de composition, décomposition exacte) ;6 * l'indice hédonique à indicatrices temporelles (traite.hedonique.indice_dummies_temps, modèle de référence traite.cas1.matrice7 avec les indicatrices trimestrielles du ch. 18), l'indice chaîné par périodes adjacentes et l'indice d'imputation géométrique8 (égal à l'indice à indicatrices quand les coefficients sont communs) ;9 * un jeu de REVENTES SYNTHÉTIQUES dérivé des 600 ventes (chaque vente observée est la première vente ; la revente est simulée10 avec la loi de génération du jeu : tendance 4,43 %/an, bruit de transaction σ = 0,085 à chaque vente, marche aléatoire11 idiosyncrasique δ² = 0,0004 par mois, 15 % de reventes rapides « rénovées » avec prime) ;12 * les estimateurs de ventes répétées de Bailey-Muth-Nourse (traite.hedonique.ventes_repetees_bmn) et de Case-Shiller pondéré en13 trois étapes (traite.hedonique.ventes_repetees_case_shiller), avec erreurs-types ;14 * l'indice SPAR (prix / valeur au rôle, Bourassa, Hoesli et Sun 2006) construit avec valeur_role_2025 ;15 * la comparaison à la vérité de génération, l'ajustement de conditions du marché du comparable G0404 selon l'indice,16 la sensibilité (millésimes et révision, filtre des reventes rapides, taille de l'échantillon, fenêtre d'estimation) et les17 chiffres des solutions d'exercices.18Usage (racine du dépôt) : PYTHONPATH=analysis/python python3 analysis/python/cases/ch35_indices.py"""19from __future__ import annotations20import sys, math, pathlib21import numpy as np, pandas as pd22sys.path.insert(0, str(pathlib.Path(__file__).resolve().parents[1]))23sys.path.insert(0, str(pathlib.Path(__file__).resolve().parent))24from traite import cas125from traite import hedonique as h26import ch18_hedonique as c182728ROOT = pathlib.Path(__file__).resolve().parents[3]29GRAINE = 2026103130T = 15 # trimestres 2023T1 … 2026T331LABELS = [f"{2023 + k // 4}T{k % 4 + 1}" for k in range(T)]32BETA_VRAI = c18.VRAIS["mois"] # ln(1,0443)/12 : tendance mensuelle de génération33SIGMA = c18.VRAIS["sigma"] # bruit de transaction de génération34DELTA2 = 0.0004 # variance mensuelle de la marche aléatoire idiosyncrasique des reventes (hypothèse)35P_REVENTE = 0.60 # probabilité qu'une vente observée soit suivie d'une revente dans la fenêtre36PART_FLIP = 0.15 # part des reventes rapides « rénovées »37PRIME_FLIP = (0.12, 0.04) # prime moyenne et écart-type (log) des reventes rapides38FIN_FENETRE = pd.Timestamp("2026-08-31")3940# ----------------------------------------------------------------------------- données et vérité41def charger() -> pd.DataFrame:42 v = cas1.charger()["ventes"].copy()43 v["trim"] = c18.trimestre(v["date_vente"])44 v["mois"] = ((v["date_vente"] - cas1.REF_DATE).dt.days / 30.4375).astype(float)45 return v4647def mois_mi_trimestre() -> np.ndarray:48 mid = pd.to_datetime([f"{2023 + k // 4}-{3 * (k % 4) + 2:02d}-15" for k in range(T)])49 return ((mid - cas1.REF_DATE).days / 30.4375).to_numpy()5051def verite(beta: float = BETA_VRAI) -> np.ndarray:52 """Indice vrai (loi de génération) au milieu de chaque trimestre, base 2023T1 = 100."""53 m = mois_mi_trimestre()54 return 100 * np.exp(beta * (m - m[0]))5556# ----------------------------------------------------------------------------- médiane, moyenne, composition57def indices_naifs(v: pd.DataFrame, td: dict) -> pd.DataFrame:58 """Indice de la médiane, de la moyenne arithmétique, de la moyenne géométrique et du prix médian au m² par trimestre ;59 décomposition exacte de la moyenne géométrique : mean(ln P)_t = α + x̄_t'β + τ_t (les résidus MCO somment à zéro par trimestre)."""60 g = v.groupby("trim")61 out = pd.DataFrame({"trimestre": LABELS, "n": g.size().reindex(range(T)).to_numpy(),62 "mediane": g["prix_vente"].median().to_numpy(), "moyenne": g["prix_vente"].mean().to_numpy(),63 "moy_geo": np.exp(g["prix_vente"].apply(lambda s: np.log(s).mean())).to_numpy(),64 "mediane_m2": g.apply(lambda d: (d["prix_vente"] / d["superficie_habitable_m2"]).median()).to_numpy(),65 "superficie_mediane": g["superficie_habitable_m2"].median().to_numpy(),66 "part_aylmer": g["secteur"].apply(lambda s: (s == "Aylmer").mean()).to_numpy(),67 "part_hull": g["secteur"].apply(lambda s: (s == "Hull").mean()).to_numpy(),68 "age_median": g["age"].median().to_numpy()})69 for c in ("mediane", "moyenne", "moy_geo", "mediane_m2"):70 out[f"indice_{c}"] = 100 * out[c] / out[c].iloc[0]71 # composition : x̄_t'β̂ (hors indicatrices temporelles) relativement à 2023T172 X = td["X"]; beta = pd.Series(td["mco"]["beta"], index=X.columns)73 cols = [c for c in X.columns if not c.startswith("trim_")]74 xb = (X[cols].values @ beta[cols].values)75 comp = pd.Series(xb, index=v.index).groupby(v["trim"]).mean().to_numpy()76 out["composition_log"] = comp - comp[0]77 out["indice_composition"] = 100 * np.exp(out["composition_log"])78 out["indice_moy_geo_corrigee"] = 100 * np.exp(np.log(out["moy_geo"] / out["moy_geo"].iloc[0]) - out["composition_log"])79 return out8081# ----------------------------------------------------------------------------- hédonique : indicatrices, chaînage, imputation82def indice_td(v: pd.DataFrame) -> dict:83 """Indice hédonique à indicatrices trimestrielles (2023T1 = référence) sur le modèle de référence du cas 1."""84 X = c18.matrice_trimestres(v); y = np.log(v["prix_vente"].astype(float)).values85 r = h.mco(X.values, y)86 cols = [c for c in X.columns if c.startswith("trim_")]87 ix = [list(X.columns).index(c) for c in cols]88 tau = r["beta"][ix]; se = r["se_hc1"][ix]89 idx = h.indice_dummies_temps(tau)90 return {"X": X, "y": y, "mco": r, "tau": np.concatenate([[0.0], tau]), "se": np.concatenate([[0.0], se]), "indice": idx,91 "bas": np.concatenate([[100.0], 100 * np.exp(tau - 1.96 * se)]), "haut": np.concatenate([[100.0], 100 * np.exp(tau + 1.96 * se)]),92 "sigma": float(np.sqrt(r["sigma2"])), "r2_ajuste": r["r2_ajuste"], "k": r["k"]}9394def indice_chaine(v: pd.DataFrame) -> dict:95 """Indice hédonique chaîné : pour chaque paire de trimestres adjacents (t−1, t), MCO sur les ventes des deux trimestres avec une96 indicatrice de t (coefficients des caractéristiques propres à la paire) ; I_t = I_{t−1} exp(δ_t)."""97 delta = np.zeros(T); se = np.zeros(T)98 for t in range(1, T):99 sub = v[v["trim"].isin([t - 1, t])]100 X = cas1.matrice(sub).drop(columns=["mois"])101 X = X[[c for c in X.columns if X[c].std() > 0 or c == "const"]]102 X["d_t"] = (sub["trim"] == t).astype(float)103 r = h.mco(X.values, np.log(sub["prix_vente"].astype(float)).values)104 j = list(X.columns).index("d_t"); delta[t] = r["beta"][j]; se[t] = r["se_hc1"][j]105 ln_idx = np.cumsum(delta)106 return {"delta": delta, "se_delta": se, "indice": 100 * np.exp(ln_idx), "se_cumul": np.sqrt(np.cumsum(se ** 2))}107108def indice_imputation(v: pd.DataFrame, td: dict) -> dict:109 """Imputation géométrique (Laspeyres : panier de la période de base ; Paasche : panier courant ; Fisher = moyenne géométrique)110 avec le modèle à indicatrices (coefficients communs) : chaque vente de la période 0 (resp. t) est « revendue » à la période t (resp. 0)111 au prix prédit. Avec des coefficients communs, les trois indices coïncident avec exp(τ_t) ; la fonction le vérifie numériquement."""112 X = td["X"]; beta = td["mco"]["beta"]; cols = list(X.columns)113 tcols = [c for c in cols if c.startswith("trim_")]114 def predire_a(Xsub: pd.DataFrame, t: int) -> np.ndarray:115 Z = Xsub.copy()116 for c in tcols: Z[c] = 0.0117 if t > 0: Z[f"trim_{t}"] = 1.0118 return Z.values @ beta119 lasp = np.zeros(T); paas = np.zeros(T)120 X0 = X[v["trim"] == 0]121 for t in range(T):122 Xt = X[v["trim"] == t]123 lasp[t] = np.mean(predire_a(X0, t) - predire_a(X0, 0))124 paas[t] = np.mean(predire_a(Xt, t) - predire_a(Xt, 0))125 fisher = 0.5 * (lasp + paas)126 return {"laspeyres": 100 * np.exp(lasp), "paasche": 100 * np.exp(paas), "fisher": 100 * np.exp(fisher),127 "ecart_max_vs_td": float(np.max(np.abs(100 * np.exp(fisher) - td["indice"])))}128129# ----------------------------------------------------------------------------- reventes synthétiques130def generer_reventes(v: pd.DataFrame, graine: int = GRAINE, p_revente: float = P_REVENTE, part_flip: float = PART_FLIP,131 delta2: float = DELTA2, beta: float = BETA_VRAI, sigma: float = SIGMA) -> pd.DataFrame:132 """Paires de reventes synthétiques. Pour chaque vente observée (première vente, prix P1 à la date t1), avec probabilité p_revente :133 intervalle Δ (mois) tiré d'une loi déclarée (reventes « rapides » : 3 à 9 mois, prime log N(0,08 ; 0,03) ; autres : 6 + Exp(16), ≥ 6) ;134 ln P2 = ln P1 + β Δ + prime·flip + w + ε2 − ε1', avec ε1', ε2 ~ N(0, σ²) (bruit de transaction à chaque vente) et w ~ N(0, δ² Δ).135 Seules les reventes datées au plus tard le 31 août 2026 sont conservées ; prix arrondis à 500 $."""136 rng = np.random.default_rng(graine)137 rows = []138 for _, s in v.iterrows():139 if rng.random() > p_revente:140 continue141 flip = rng.random() < part_flip142 if flip:143 gap = int(rng.integers(3, 10)); prime = rng.normal(*PRIME_FLIP)144 else:145 gap = int(max(6, round(6 + rng.exponential(16.0)))); prime = 0.0146 d2 = s["date_vente"] + pd.DateOffset(months=gap)147 if d2 > FIN_FENETRE:148 continue149 w = rng.normal(0.0, math.sqrt(delta2 * gap)); e1 = rng.normal(0.0, sigma); e2 = rng.normal(0.0, sigma)150 lnp2 = math.log(s["prix_vente"]) + beta * gap + prime + w + e2 - e1151 p2 = 500 * round(math.exp(lnp2) / 500)152 rows.append({"id_vente": s["id_vente"], "secteur": s["secteur"], "date_1": s["date_vente"], "date_2": d2, "prix_1": float(s["prix_vente"]),153 "prix_2": float(p2), "t_1": int(s["trim"]), "t_2": int(c18.trimestre(pd.Series([d2]))[0]), "intervalle_mois": gap,154 "flip": bool(flip), "prime": prime})155 df = pd.DataFrame(rows)156 df["ln_ratio"] = np.log(df["prix_2"] / df["prix_1"])157 df["id_paire"] = [f"R{i + 1:04d}" for i in range(len(df))]158 return df159160def resume_reventes(p: pd.DataFrame) -> dict:161 return {"n_paires": len(p), "n_flips": int(p["flip"].sum()), "intervalle_median": float(p["intervalle_mois"].median()),162 "intervalle_moyen": float(p["intervalle_mois"].mean()), "intervalle_min": int(p["intervalle_mois"].min()), "intervalle_max": int(p["intervalle_mois"].max()),163 "part_moins_12": float((p["intervalle_mois"] < 12).mean()), "variation_mediane": float(np.exp(p["ln_ratio"].median()) - 1),164 "variation_moyenne_log": float(p["ln_ratio"].mean()), "sd_ln_ratio": float(p["ln_ratio"].std(ddof=1)),165 "paires_par_annee_1": p["date_1"].dt.year.value_counts().sort_index().to_dict(),166 "paires_par_annee_2": p["date_2"].dt.year.value_counts().sort_index().to_dict(),167 "meme_trimestre": int((p["t_1"] == p["t_2"]).sum())}168169# ----------------------------------------------------------------------------- ventes répétées : BMN et Case–Shiller avec erreurs-types170def _se_mco_sans_constante(D: np.ndarray, y: np.ndarray, beta: np.ndarray, w: np.ndarray | None = None) -> np.ndarray:171 if w is None: w = np.ones(len(y))172 Dw = D * w[:, None]; yw = y * w173 res = yw - Dw @ beta; n, k = D.shape174 s2 = res @ res / (n - k)175 return np.sqrt(np.diag(s2 * np.linalg.pinv(Dw.T @ Dw)))176177def bmn(p: pd.DataFrame) -> dict:178 r = h.ventes_repetees_bmn(p["prix_1"], p["prix_2"], p["t_1"], p["t_2"], T)179 se = _se_mco_sans_constante(r["D"], r["y"], r["beta"])180 r["se"] = np.concatenate([[0.0], se]); r["bas"] = 100 * np.exp(np.concatenate([[0.0], r["beta"]]) - 1.96 * r["se"])181 r["haut"] = 100 * np.exp(np.concatenate([[0.0], r["beta"]]) + 1.96 * r["se"]); r["sigma"] = float(np.sqrt(r["residus"] @ r["residus"] / (len(r["y"]) - T + 1)))182 # vérification de la forme fermée : β = (D'D)^{-1} D'y183 r["beta_ferme"] = np.linalg.solve(r["D"].T @ r["D"], r["D"].T @ r["y"])184 return r185186def case_shiller(p: pd.DataFrame) -> dict:187 r = h.ventes_repetees_case_shiller(p["prix_1"], p["prix_2"], p["t_1"], p["t_2"], T)188 e1 = h.ventes_repetees_bmn(p["prix_1"], p["prix_2"], p["t_1"], p["t_2"], T)189 gap_t = (p["t_2"] - p["t_1"]).to_numpy(float)190 se = _se_mco_sans_constante(e1["D"], e1["y"], r["beta"], r["poids"])191 r["se"] = np.concatenate([[0.0], se]); r["bas"] = 100 * np.exp(np.concatenate([[0.0], r["beta"]]) - 1.96 * r["se"])192 r["haut"] = 100 * np.exp(np.concatenate([[0.0], r["beta"]]) + 1.96 * r["se"])193 r["gap_trimestres"] = gap_t; r["residus_bmn2"] = e1["residus"] ** 2194 a = r["coef_variance"] # var(ε) = a0 + a1·(intervalle en trimestres)195 r["variance_a_gap"] = {g: float(a[0] + a[1] * g) for g in (1, 2, 4, 8, 12)}196 r["poids_relatif_a_gap"] = {g: float(math.sqrt((a[0] + a[1] * 1) / (a[0] + a[1] * g))) for g in (1, 2, 4, 8, 12)}197 # pente théorique : var = 2σ² + δ²·(3 mois par trimestre)·gap, en unités de trimestre198 r["variance_theorique"] = (2 * SIGMA ** 2, DELTA2 * 3)199 return r200201# ----------------------------------------------------------------------------- SPAR202def spar(v: pd.DataFrame) -> dict:203 """Indice SPAR (sale price appraisal ratio) : arithmétique [(ΣP_t/ΣA_t)/(ΣP_0/ΣA_0)] et géométrique [exp(mean ln(P/A)_t − mean ln(P/A)_0)],204 A = valeur au rôle 2025 (date de référence 2024-07-01, fictive)."""205 g = v.groupby("trim")206 ratio_arith = (g["prix_vente"].sum() / g["valeur_role_2025"].sum()).to_numpy()207 lr = np.log(v["prix_vente"] / v["valeur_role_2025"])208 ratio_geo = lr.groupby(v["trim"]).mean().to_numpy(); sd_geo = lr.groupby(v["trim"]).std(ddof=1).to_numpy(); n = g.size().to_numpy()209 se_geo = np.sqrt(sd_geo ** 2 / n + sd_geo[0] ** 2 / n[0]); se_geo[0] = 0.0210 return {"ratio_arith": ratio_arith, "ratio_geo": np.exp(ratio_geo), "indice_arith": 100 * ratio_arith / ratio_arith[0],211 "indice_geo": 100 * np.exp(ratio_geo - ratio_geo[0]), "se_log": se_geo,212 "bas": 100 * np.exp(ratio_geo - ratio_geo[0] - 1.96 * se_geo), "haut": 100 * np.exp(ratio_geo - ratio_geo[0] + 1.96 * se_geo),213 "sd_ln_ratio": float(lr.std(ddof=1)), "ratio_median_global": float(np.median(v["prix_vente"] / v["valeur_role_2025"]))}214215# ----------------------------------------------------------------------------- comparaison216def precision(idx: np.ndarray, vrai: np.ndarray) -> dict:217 e = 100 * (idx / vrai - 1)218 return {"EAM_pct": float(np.mean(np.abs(e[1:]))), "RMSE_pct": float(np.sqrt(np.mean(e[1:] ** 2))), "max_abs_pct": float(np.max(np.abs(e))),219 "croissance_totale_pct": float(idx[-1] - 100), "ecart_final_pct": float(e[-1]),220 "variations_trim_sd_pct": float(np.std(np.diff(np.log(idx))) * 100)}221222def tableau_indices(vrai, naifs, td, ch, b, cs, sp) -> pd.DataFrame:223 return pd.DataFrame({"trimestre": LABELS, "n_ventes": naifs["n"], "vrai": vrai, "mediane": naifs["indice_mediane"].to_numpy(),224 "moy_geo": naifs["indice_moy_geo"].to_numpy(), "hedonique": td["indice"], "hed_bas": td["bas"], "hed_haut": td["haut"],225 "chaine": ch["indice"], "bmn": b["indice"], "bmn_bas": b["bas"], "bmn_haut": b["haut"], "cs": cs["indice"],226 "cs_bas": cs["bas"], "cs_haut": cs["haut"], "spar": sp["indice_geo"], "spar_arith": sp["indice_arith"]})227228def ajustement_comparable(v: pd.DataFrame, tab: pd.DataFrame, ref: dict, id_vente: str = "G0404", date_eval: str = "2026-09-01") -> dict:229 """Ajustement de conditions du marché d'un comparable du ch. 16 selon l'indice retenu : A = P × (I_eval / I_vente − 1)."""230 c = v[v["id_vente"] == id_vente].iloc[0]; P = float(c["prix_vente"]); t1 = int(c["trim"]); t2 = T - 1231 dm = (pd.Timestamp(date_eval) - c["date_vente"]).days / 30.4375232 bm = ref["beta"][list(ref["X"].columns).index("mois")]233 out = {"prix": P, "trimestre_vente": LABELS[t1], "mois_ecoules": float(dm), "tendance": P * (math.exp(bm * dm) - 1)}234 for k in ("mediane", "hedonique", "chaine", "bmn", "cs", "spar"):235 out[k] = P * (tab[k].iloc[t2] / tab[k].iloc[t1] - 1)236 out["vrai"] = P * (math.exp(BETA_VRAI * dm) - 1)237 return out238239# ----------------------------------------------------------------------------- sensibilité240def sensibilite(v: pd.DataFrame, p: pd.DataFrame, vrai: np.ndarray, td: dict, graine: int = GRAINE) -> dict:241 out = {}242 # (a) millésimes : révision des indices de ventes répétées quand de nouvelles paires s'ajoutent243 mill = {}244 for fin, lab in ((9, "2025T2"), (11, "2025T4"), (T - 1, "2026T3")):245 sub = p[p["t_2"] <= fin]246 rb = h.ventes_repetees_bmn(sub["prix_1"], sub["prix_2"], sub["t_1"], sub["t_2"], fin + 1)247 rc = h.ventes_repetees_case_shiller(sub["prix_1"], sub["prix_2"], sub["t_1"], sub["t_2"], fin + 1)248 mill[lab] = {"n": len(sub), "bmn": rb["indice"], "cs": rc["indice"]}249 out["millesimes"] = mill250 rev = {}251 for t, lab in ((9, "2025T2"), (7, "2024T4"), (5, "2024T2")):252 rev[lab] = {"bmn_v2025T2": float(mill["2025T2"]["bmn"][t]), "bmn_v2025T4": float(mill["2025T4"]["bmn"][t]), "bmn_v2026T3": float(mill["2026T3"]["bmn"][t]),253 "cs_v2025T2": float(mill["2025T2"]["cs"][t]), "cs_v2026T3": float(mill["2026T3"]["cs"][t])}254 out["revisions"] = rev255 # (b) filtre des reventes rapides256 filt = {}257 for seuil, lab in ((0, "toutes"), (12, "intervalle ≥ 12 mois"), (18, "intervalle ≥ 18 mois")):258 sub = p[p["intervalle_mois"] >= seuil] if seuil else p259 rb = h.ventes_repetees_bmn(sub["prix_1"], sub["prix_2"], sub["t_1"], sub["t_2"], T)260 rc = h.ventes_repetees_case_shiller(sub["prix_1"], sub["prix_2"], sub["t_1"], sub["t_2"], T)261 filt[lab] = {"n": len(sub), "n_flips": int(sub["flip"].sum()), "bmn": precision(rb["indice"], vrai), "cs": precision(rc["indice"], vrai),262 "bmn_indice": rb["indice"], "cs_indice": rc["indice"]}263 sans_flip = p[~p["flip"]]264 rb = h.ventes_repetees_bmn(sans_flip["prix_1"], sans_flip["prix_2"], sans_flip["t_1"], sans_flip["t_2"], T)265 filt["sans les reventes rénovées (oracle)"] = {"n": len(sans_flip), "n_flips": 0, "bmn": precision(rb["indice"], vrai), "bmn_indice": rb["indice"]}266 out["filtre"] = filt267 # (c) taille de l'échantillon : sous-échantillons aléatoires de paires, erreur-type du dernier trimestre et EAM (moyenne de 200 tirages)268 rng = np.random.default_rng(graine); taille = {}269 for frac in (1.0, 0.5, 0.25):270 eam, fin = [], []271 for _ in range(200 if frac < 1 else 1):272 sub = p.sample(frac=frac, random_state=int(rng.integers(1e9))) if frac < 1 else p273 try:274 rb = h.ventes_repetees_bmn(sub["prix_1"], sub["prix_2"], sub["t_1"], sub["t_2"], T)275 except Exception:276 continue277 eam.append(precision(rb["indice"], vrai)["EAM_pct"]); fin.append(rb["indice"][-1])278 taille[frac] = {"n": int(round(frac * len(p))), "EAM_moyen_pct": float(np.mean(eam)), "sd_indice_final": float(np.std(fin)) if len(fin) > 1 else float("nan")}279 out["taille"] = taille280 # (d) fenêtre d'estimation de l'indice hédonique : huit derniers trimestres seulement, rechaîné sur 2024T4281 sub = v[v["trim"] >= 7]282 X = cas1.matrice(sub).drop(columns=["mois"]); X = X[[c for c in X.columns if X[c].std() > 0 or c == "const"]]283 for k in range(8, T): X[f"trim_{k}"] = (sub["trim"] == k).astype(float)284 r = h.mco(X.values, np.log(sub["prix_vente"].astype(float)).values)285 tau = np.array([r["beta"][list(X.columns).index(f"trim_{k}")] for k in range(8, T)])286 idx_court = td["indice"][7] * np.exp(np.concatenate([[0.0], tau]))287 out["fenetre"] = {"n": len(sub), "indice_court": idx_court, "indice_long": td["indice"][7:], "ecart_max_pct": float(np.max(np.abs(idx_court / td["indice"][7:] - 1)) * 100),288 "sigma_court": float(np.sqrt(r["sigma2"]))}289 return out290291# ----------------------------------------------------------------------------- exercices292def exercices(tab: pd.DataFrame, b: dict, cs: dict, sp: dict, ref: dict) -> dict:293 e = {}294 # 35.1 : médiane et composition (deux trimestres fictifs)295 e["1"] = {"mediane_t0": 450_000, "mediane_t1": 480_000, "var_mediane_pct": (480_000 / 450_000 - 1) * 100,296 "composition_pct": (math.exp(0.539 * math.log(150 / 140)) - 1) * 100,297 "variation_pure_pct": (480_000 / 450_000 / math.exp(0.539 * math.log(150 / 140)) - 1) * 100}298 # 35.2 : indice à indicatrices299 tau = np.array([0.0, 0.031, 0.058, 0.102]); e["2"] = {"indice": 100 * np.exp(tau), "variation_t3_t2_pct": (math.exp(0.102 - 0.058) - 1) * 100,300 "ic_t3": (100 * math.exp(0.102 - 1.96 * 0.021), 100 * math.exp(0.102 + 1.96 * 0.021))}301 # 35.3 : BMN à la main sur quatre paires et trois périodes302 D = np.array([[1.0, 0.0], [-1.0, 1.0], [0.0, 1.0], [1.0, 0.0]]); y = np.log(np.array([462_000 / 430_000, 512_000 / 470_000, 540_000 / 495_000, 455_000 / 436_000]))303 beta = np.linalg.solve(D.T @ D, D.T @ y)304 e["3"] = {"y": y, "DtD": D.T @ D, "Dty": D.T @ y, "beta": beta, "indice": 100 * np.exp(np.concatenate([[0.0], beta])), "residus": y - D @ beta}305 # 35.4 : Case–Shiller : variance et poids306 a0, a1 = 0.0145, 0.0012; gaps = np.array([1, 4, 8, 12])307 var = a0 + a1 * gaps; e["4"] = {"var": var, "poids": 1 / np.sqrt(var), "poids_rel": np.sqrt(var[0] / var), "sd_pct": 100 * np.sqrt(var)}308 # 35.5 : SPAR sur dix ventes fictives309 P = np.array([512_000, 468_500, 601_000, 445_000, 530_000, 487_500, 575_000, 455_000, 498_000, 522_000.0]); A = np.array([470_000, 452_000, 545_000, 425_000, 480_000, 468_000, 510_000, 440_000, 455_000, 471_000.0])310 base_ratio = 1.045; r_ar = P.sum() / A.sum(); r_geo = math.exp(np.mean(np.log(P / A)))311 e["5"] = {"ratio_arith": r_ar, "ratio_geo": r_geo, "spar_arith": 100 * r_ar / base_ratio, "spar_geo": 100 * r_geo / base_ratio, "ratios": P / A}312 # 35.6 : ajustement de G0404 selon l'indice313 e["6"] = {"P": 506_500.0, "I_2025T2": {k: float(tab[k].iloc[9]) for k in ("hedonique", "bmn", "cs", "mediane")},314 "I_2026T3": {k: float(tab[k].iloc[14]) for k in ("hedonique", "bmn", "cs", "mediane")}}315 e["6"]["ajustements"] = {k: 506_500.0 * (e["6"]["I_2026T3"][k] / e["6"]["I_2025T2"][k] - 1) for k in e["6"]["I_2025T2"]}316 bm = ref["beta"][list(ref["X"].columns).index("mois")]; e["6"]["tendance"] = 506_500.0 * (math.exp(bm * 15.2) - 1)317 # 35.7 : révision et biais de sélection (chiffres du script : voir sensibilité) — valeurs conceptuelles318 # 35.8 : script (pas de chiffre fixe)319 return e320321# ----------------------------------------------------------------------------- principal322def main(verbose: bool = True) -> dict:323 v = charger(); ref = cas1.hedonique(v.drop(columns=["trim", "mois"]))324 vrai = verite()325 td = indice_td(v); ch = indice_chaine(v); imp = indice_imputation(v, td)326 naifs = indices_naifs(v, td)327 p = generer_reventes(v); rp = resume_reventes(p)328 b = bmn(p); cs = case_shiller(p); sp = spar(v)329 tab = tableau_indices(vrai, naifs, td, ch, b, cs, sp)330 prec = {k: precision(tab[k].to_numpy(), vrai) for k in ("mediane", "moy_geo", "hedonique", "chaine", "bmn", "cs", "spar", "spar_arith")}331 prec["hedonique"]["largeur_ic_finale_pct"] = float((td["haut"][-1] / td["bas"][-1] - 1) * 100)332 prec["bmn"]["largeur_ic_finale_pct"] = float((b["haut"][-1] / b["bas"][-1] - 1) * 100)333 prec["cs"]["largeur_ic_finale_pct"] = float((cs["haut"][-1] / cs["bas"][-1] - 1) * 100)334 prec["spar"]["largeur_ic_finale_pct"] = float((sp["haut"][-1] / sp["bas"][-1] - 1) * 100)335 # corrélation entre l'écart médiane − hédonique et la composition336 ecart_med = np.log(naifs["indice_mediane"].to_numpy() / td["indice"]); corr_comp = float(np.corrcoef(ecart_med, naifs["composition_log"])[0, 1])337 corr_surf = float(np.corrcoef(ecart_med, np.log(naifs["superficie_mediane"]))[0, 1])338 adj = ajustement_comparable(v, tab, ref)339 sens = sensibilite(v, p, vrai, td)340 ex = exercices(tab, b, cs, sp, ref)341 out = {"ventes": v, "ref": ref, "vrai": vrai, "td": td, "chaine": ch, "imputation": imp, "naifs": naifs, "paires": p, "resume_paires": rp,342 "bmn": b, "cs": cs, "spar": sp, "tableau": tab, "precision": prec, "corr_ecart_composition": corr_comp, "corr_ecart_superficie": corr_surf,343 "ajustement": adj, "sensibilite": sens, "exercices": ex, "beta_vrai": BETA_VRAI}344 if verbose:345 pd.set_option("display.width", 240); pd.set_option("display.max_columns", 40)346 print(f"== Vérité de génération : β = {BETA_VRAI:.6f}/mois ({(math.exp(12 * BETA_VRAI) - 1) * 100:.2f} %/an) ; indice final {vrai[-1]:.2f}")347 print("== Médiane, moyennes et composition :\n", naifs.round(3).to_string(index=False))348 print(f" corrélation(ln médiane/hédonique, composition) = {corr_comp:.3f} ; avec ln superficie médiane = {corr_surf:.3f}")349 print(f"== Hédonique à indicatrices : k = {td['k']}, σ = {td['sigma']:.4f}, R² aj. = {td['r2_ajuste']:.4f}")350 print(f"== Imputation géométrique : écart max à l'indice à indicatrices = {imp['ecart_max_vs_td']:.2e}")351 print("== Reventes synthétiques :", {k: (round(x, 4) if isinstance(x, float) else x) for k, x in rp.items()})352 print(f"== BMN : σ = {b['sigma']:.4f} ; écart max forme fermée vs lstsq = {np.max(np.abs(b['beta'] - b['beta_ferme'])):.1e}")353 print(f"== Case–Shiller : var(ε) = {cs['coef_variance'][0]:.5f} + {cs['coef_variance'][1]:.5f} × intervalle (trimestres) ; théorie {cs['variance_theorique']}")354 print(" poids relatifs :", {k: round(x, 3) for k, x in cs["poids_relatif_a_gap"].items()})355 print(f"== SPAR : sd ln(P/A) = {sp['sd_ln_ratio']:.4f} ; ratio médian global {sp['ratio_median_global']:.4f}")356 print("== Indices par trimestre :\n", tab.round(2).to_string(index=False))357 print("== Précision vs vérité :")358 for k, d in prec.items(): print(f" {k:11s}", {kk: round(x, 2) for kk, x in d.items()})359 print("== Ajustement de G0404 :", {k: (round(x, 1) if isinstance(x, float) else x) for k, x in adj.items()})360 print("== Sensibilité — révisions :", sens["revisions"])361 for lab, d in sens["filtre"].items(): print(f" filtre {lab:36s} n={d['n']:3d} flips={d['n_flips']:2d} BMN EAM={d['bmn']['EAM_pct']:.2f} final={d['bmn']['croissance_totale_pct']:.2f}", f"CS EAM={d['cs']['EAM_pct']:.2f}" if "cs" in d else "")362 print(" taille :", sens["taille"])363 print(" fenêtre hédonique :", {k: (np.round(x, 2) if isinstance(x, np.ndarray) else round(x, 3)) for k, x in sens["fenetre"].items()})364 print("== Exercices :")365 for k, e in ex.items(): print(" ", k, {kk: (np.round(vv, 4) if isinstance(vv, (np.ndarray, float)) else vv) for kk, vv in e.items()})366 return out367368if __name__ == "__main__":369 main()370