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"], "cs_coef_variance": rc["coef_variance"]}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 rc = h.ventes_repetees_case_shiller(sans_flip["prix_1"], sans_flip["prix_2"], sans_flip["t_1"], sans_flip["t_2"], T)266 filt["sans les reventes rénovées (oracle)"] = {"n": len(sans_flip), "n_flips": 0, "bmn": precision(rb["indice"], vrai), "cs": precision(rc["indice"], vrai),267 "bmn_indice": rb["indice"], "cs_indice": rc["indice"], "cs_coef_variance": rc["coef_variance"]}268 out["filtre"] = filt269 # (c) taille de l'échantillon : sous-échantillons aléatoires de paires, erreur-type du dernier trimestre et EAM (moyenne de 200 tirages)270 rng = np.random.default_rng(graine); taille = {}271 for frac in (1.0, 0.5, 0.25):272 eam, fin = [], []273 for _ in range(200 if frac < 1 else 1):274 sub = p.sample(frac=frac, random_state=int(rng.integers(1e9))) if frac < 1 else p275 try:276 rb = h.ventes_repetees_bmn(sub["prix_1"], sub["prix_2"], sub["t_1"], sub["t_2"], T)277 except Exception:278 continue279 eam.append(precision(rb["indice"], vrai)["EAM_pct"]); fin.append(rb["indice"][-1])280 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")}281 out["taille"] = taille282 # (d) fenêtre d'estimation de l'indice hédonique : huit derniers trimestres seulement, rechaîné sur 2024T4283 sub = v[v["trim"] >= 7]284 X = cas1.matrice(sub).drop(columns=["mois"]); X = X[[c for c in X.columns if X[c].std() > 0 or c == "const"]]285 for k in range(8, T): X[f"trim_{k}"] = (sub["trim"] == k).astype(float)286 r = h.mco(X.values, np.log(sub["prix_vente"].astype(float)).values)287 tau = np.array([r["beta"][list(X.columns).index(f"trim_{k}")] for k in range(8, T)])288 idx_court = td["indice"][7] * np.exp(np.concatenate([[0.0], tau]))289 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),290 "sigma_court": float(np.sqrt(r["sigma2"]))}291 return out292293# ----------------------------------------------------------------------------- exercices294def exercices(tab: pd.DataFrame, b: dict, cs: dict, sp: dict, ref: dict) -> dict:295 e = {}296 # 35.1 : médiane et composition (deux trimestres fictifs)297 e["1"] = {"mediane_t0": 450_000, "mediane_t1": 480_000, "var_mediane_pct": (480_000 / 450_000 - 1) * 100,298 "composition_pct": (math.exp(0.539 * math.log(150 / 140)) - 1) * 100,299 "variation_pure_pct": (480_000 / 450_000 / math.exp(0.539 * math.log(150 / 140)) - 1) * 100}300 # 35.2 : indice à indicatrices301 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,302 "ic_t3": (100 * math.exp(0.102 - 1.96 * 0.021), 100 * math.exp(0.102 + 1.96 * 0.021))}303 # 35.3 : BMN à la main sur quatre paires et trois périodes304 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]))305 beta = np.linalg.solve(D.T @ D, D.T @ y)306 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}307 # 35.4 : Case–Shiller : variance et poids308 a0, a1 = 0.0145, 0.0012; gaps = np.array([1, 4, 8, 12])309 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)}310 # 35.5 : SPAR sur dix ventes fictives311 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])312 base_ratio = 1.045; r_ar = P.sum() / A.sum(); r_geo = math.exp(np.mean(np.log(P / A)))313 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}314 # 35.6 : ajustement de G0404 selon l'indice315 e["6"] = {"P": 506_500.0, "I_2025T2": {k: float(tab[k].iloc[9]) for k in ("hedonique", "bmn", "cs", "mediane")},316 "I_2026T3": {k: float(tab[k].iloc[14]) for k in ("hedonique", "bmn", "cs", "mediane")}}317 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"]}318 bm = ref["beta"][list(ref["X"].columns).index("mois")]; e["6"]["tendance"] = 506_500.0 * (math.exp(bm * 15.2) - 1)319 # 35.7 : révision et biais de sélection (chiffres du script : voir sensibilité) — valeurs conceptuelles320 # 35.8 : script (pas de chiffre fixe)321 return e322323# ----------------------------------------------------------------------------- principal324def main(verbose: bool = True) -> dict:325 v = charger(); ref = cas1.hedonique(v.drop(columns=["trim", "mois"]))326 vrai = verite()327 td = indice_td(v); ch = indice_chaine(v); imp = indice_imputation(v, td)328 naifs = indices_naifs(v, td)329 p = generer_reventes(v); rp = resume_reventes(p)330 b = bmn(p); cs = case_shiller(p); sp = spar(v)331 tab = tableau_indices(vrai, naifs, td, ch, b, cs, sp)332 prec = {k: precision(tab[k].to_numpy(), vrai) for k in ("mediane", "moy_geo", "hedonique", "chaine", "bmn", "cs", "spar", "spar_arith")}333 prec["hedonique"]["largeur_ic_finale_pct"] = float((td["haut"][-1] / td["bas"][-1] - 1) * 100)334 prec["bmn"]["largeur_ic_finale_pct"] = float((b["haut"][-1] / b["bas"][-1] - 1) * 100)335 prec["cs"]["largeur_ic_finale_pct"] = float((cs["haut"][-1] / cs["bas"][-1] - 1) * 100)336 prec["spar"]["largeur_ic_finale_pct"] = float((sp["haut"][-1] / sp["bas"][-1] - 1) * 100)337 # corrélation entre l'écart médiane − hédonique et la composition338 ecart_med = np.log(naifs["indice_mediane"].to_numpy() / td["indice"]); corr_comp = float(np.corrcoef(ecart_med, naifs["composition_log"])[0, 1])339 corr_surf = float(np.corrcoef(ecart_med, np.log(naifs["superficie_mediane"]))[0, 1])340 adj = ajustement_comparable(v, tab, ref)341 sens = sensibilite(v, p, vrai, td)342 ex = exercices(tab, b, cs, sp, ref)343 out = {"ventes": v, "ref": ref, "vrai": vrai, "td": td, "chaine": ch, "imputation": imp, "naifs": naifs, "paires": p, "resume_paires": rp,344 "bmn": b, "cs": cs, "spar": sp, "tableau": tab, "precision": prec, "corr_ecart_composition": corr_comp, "corr_ecart_superficie": corr_surf,345 "ajustement": adj, "sensibilite": sens, "exercices": ex, "beta_vrai": BETA_VRAI}346 if verbose:347 pd.set_option("display.width", 240); pd.set_option("display.max_columns", 40)348 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}")349 print("== Médiane, moyennes et composition :\n", naifs.round(3).to_string(index=False))350 print(f" corrélation(ln médiane/hédonique, composition) = {corr_comp:.3f} ; avec ln superficie médiane = {corr_surf:.3f}")351 print(f"== Hédonique à indicatrices : k = {td['k']}, σ = {td['sigma']:.4f}, R² aj. = {td['r2_ajuste']:.4f}")352 print(f"== Imputation géométrique : écart max à l'indice à indicatrices = {imp['ecart_max_vs_td']:.2e}")353 print("== Reventes synthétiques :", {k: (round(x, 4) if isinstance(x, float) else x) for k, x in rp.items()})354 print(f"== BMN : σ = {b['sigma']:.4f} ; écart max forme fermée vs lstsq = {np.max(np.abs(b['beta'] - b['beta_ferme'])):.1e}")355 print(f"== Case–Shiller : var(ε) = {cs['coef_variance'][0]:.5f} + {cs['coef_variance'][1]:.5f} × intervalle (trimestres) ; théorie {cs['variance_theorique']}")356 print(" poids relatifs :", {k: round(x, 3) for k, x in cs["poids_relatif_a_gap"].items()})357 print(f"== SPAR : sd ln(P/A) = {sp['sd_ln_ratio']:.4f} ; ratio médian global {sp['ratio_median_global']:.4f}")358 print("== Indices par trimestre :\n", tab.round(2).to_string(index=False))359 print("== Précision vs vérité :")360 for k, d in prec.items(): print(f" {k:11s}", {kk: round(x, 2) for kk, x in d.items()})361 print("== Ajustement de G0404 :", {k: (round(x, 1) if isinstance(x, float) else x) for k, x in adj.items()})362 print("== Sensibilité — révisions :", sens["revisions"])363 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} CS EAM={d['cs']['EAM_pct']:.2f} final={d['cs']['croissance_totale_pct']:.2f} var(ε) = {d['cs_coef_variance'][0]:.5f} + {d['cs_coef_variance'][1]:.5f}·gap")364 print(" taille :", sens["taille"])365 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()})366 print("== Exercices :")367 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()})368 return out369370if __name__ == "__main__":371 main()372