spb/traite-evaluation-qc
Public
1#!/usr/bin/env python32"""Chapitre 18 — Modèles hédoniques : exemple intégral sur les 600 ventes du cas 1 (Gatineau, DONNÉES SYNTHÉTIQUES, graine 20261031).34Le modèle de référence est `traite.cas1.hedonique` (semi-log, HC1) ; il n'est PAS redéfini ici (les ch. 17, 19, 20, 35, 37, 39, 45 s'y réfèrent).5Ce script ajoute : la comparaison des coefficients estimés aux paramètres de génération (data/synthetic/README.md), les variantes de6spécification (linéaire, log-log, Box–Cox, sans secteurs, indicatrices trimestrielles), les diagnostics (Breusch–Pagan, Jarque–Bera, VIF,7distances de Cook, I de Moran des résidus), les prix implicites au sujet (en $), la prédiction avec correction de Duan et son intervalle,8les statistiques de ratio (COD, PRD, PRB) du modèle et du rôle fictif, et l'analyse de sensibilité (spécification, sous-échantillons).9Usage (racine du dépôt) : PYTHONPATH=analysis/python python3 analysis/python/cases/ch18_hedonique.py"""10from __future__ import annotations11import sys, math, pathlib12import numpy as np, pandas as pd13sys.path.insert(0, str(pathlib.Path(__file__).resolve().parents[1]))14from traite import cas115from traite import hedonique as h16from scipy import stats, optimize1718ROOT = pathlib.Path(__file__).resolve().parents[3]19GRAINE = 202610312021# ----------------------------------------------------------------------------- paramètres de génération (README synthétique)22_g = (1 + 9.06 / 100) ** 0.5 - 1 # croissance annuelle composée du jeu (variation du rôle 2027-2029 sur deux ans)23VRAIS = {"ln_S": 0.55, "ln_T": 0.10, "age": -0.0055, "age2_100": 0.0030, "qualite": 0.06, "garage": 0.035, "sous_sol": 0.03,24 "piscine": 0.02, "etages": 0.02, "dist_centre": -0.006,25 "sect_Gatineau": math.log(464_900 / 475_000), "sect_Aylmer": math.log(530_000 / 475_000),26 "sect_Buckingham": math.log(395_000 / 475_000), "sect_Masson-Angers": math.log(395_000 / 475_000) - 0.02,27 "mois": math.log1p(_g) / 12, "sigma": 0.085}28ETIQUETTES = {"const": "Constante", "ln_S": "ln(superficie habitable / 140 m²)", "ln_T": "ln(terrain / 560 m²)", "age": "Âge (années)",29 "age2_100": "Âge² / 100", "qualite": "Qualité (écart à 3)", "garage": "Garage (places)", "sous_sol": "Sous-sol fini",30 "piscine": "Piscine", "etages": "À étages", "dist_centre": "Distance au centre au-delà de 3 km (km)",31 "sect_Gatineau": "Secteur Gatineau", "sect_Aylmer": "Secteur Aylmer", "sect_Buckingham": "Secteur Buckingham",32 "sect_Masson-Angers": "Secteur Masson-Angers", "mois": "Mois depuis le 1er juillet 2024"}3334# ----------------------------------------------------------------------------- utilitaires35def _duan(res: np.ndarray) -> float:36 return float(np.mean(np.exp(res)))3738def ratios(valeurs: np.ndarray, prix: np.ndarray) -> dict:39 return h.statistiques_ratio(valeurs, prix)4041def _rmse(a, b) -> float:42 a = np.asarray(a, float); b = np.asarray(b, float)43 return float(np.sqrt(np.mean((a - b) ** 2)))4445def trimestre(d: pd.Series) -> np.ndarray:46 d = pd.to_datetime(d)47 return ((d.dt.year - 2023) * 4 + (d.dt.month - 1) // 3).to_numpy()4849# ----------------------------------------------------------------------------- variantes de spécification50def matrice_lineaire(v: pd.DataFrame) -> pd.DataFrame:51 """Forme linéaire : mêmes variables, superficie et terrain en niveaux (m²)."""52 X = cas1.matrice(v).copy()53 X["S"] = v["superficie_habitable_m2"].astype(float); X["T"] = v["superficie_terrain_m2"].astype(float)54 X = X.drop(columns=["ln_S", "ln_T"])55 return X[["const", "S", "T"] + [c for c in X.columns if c not in ("const", "S", "T")]]5657def matrice_loglog(v: pd.DataFrame) -> pd.DataFrame:58 """Forme log-log : l'âge entre en logarithme (ln(1 + âge)) à la place du polynôme."""59 X = cas1.matrice(v).copy()60 X["ln_age"] = np.log1p(v["age"].astype(float))61 X = X.drop(columns=["age", "age2_100"])62 return X[["const", "ln_S", "ln_T", "ln_age"] + [c for c in X.columns if c not in ("const", "ln_S", "ln_T", "ln_age")]]6364def matrice_sans_secteurs(v: pd.DataFrame) -> pd.DataFrame:65 X = cas1.matrice(v)66 return X[[c for c in X.columns if not c.startswith("sect_")]]6768def matrice_trimestres(v: pd.DataFrame) -> pd.DataFrame:69 """Indicatrices trimestrielles (2023T1 = référence) à la place de la tendance linéaire `mois`."""70 X = cas1.matrice(v).drop(columns=["mois"])71 q = trimestre(v["date_vente"])72 for k in range(1, int(q.max()) + 1):73 X[f"trim_{k}"] = (q == k).astype(float)74 return X7576def boxcox_profil(P: np.ndarray, X: np.ndarray) -> dict:77 """Transformation de Box–Cox (1964) de la variable dépendante : y(λ) = (P^λ − 1)/λ, λ = 0 → ln P.78 Log-vraisemblance concentrée : ℓ(λ) = −(n/2) ln(SCR_λ / n) + (λ − 1) Σ ln P. IC à 95 % : ℓ(λ) ≥ ℓ_max − χ²₁(0,95)/2."""79 P = np.asarray(P, float); n = len(P); slog = np.log(P).sum()80 def y_l(l): return np.log(P) if abs(l) < 1e-8 else (P ** l - 1) / l81 def ll(l):82 y = y_l(l); b = np.linalg.lstsq(X, y, rcond=None)[0]; r = y - X @ b83 return -(n / 2) * np.log(r @ r / n) + (l - 1) * slog84 opt = optimize.minimize_scalar(lambda l: -ll(l), bounds=(-1.0, 1.5), method="bounded")85 lam = float(opt.x); llmax = ll(lam); seuil = llmax - stats.chi2.ppf(0.95, 1) / 286 bas = optimize.brentq(lambda l: ll(l) - seuil, -1.0, lam); haut = optimize.brentq(lambda l: ll(l) - seuil, lam, 1.5)87 grille = np.linspace(-0.6, 1.2, 73)88 return {"lambda": lam, "ic95": (float(bas), float(haut)), "ll_max": float(llmax), "ll_0": float(ll(0.0)), "ll_1": float(ll(1.0)),89 "grille": grille, "ll_grille": np.array([ll(l) for l in grille]), "y": y_l(lam)}9091def _predire_lineaire(beta, X): return X @ beta9293def ajuster_variantes(v: pd.DataFrame, s: dict, ref: dict) -> dict:94 """Estime les variantes et renvoie, pour chacune : R², σ, RMSE en $, COD, prédiction du sujet."""95 P = v["prix_vente"].astype(float).to_numpy(); lnP = np.log(P)96 sd = cas1.sujet_df(s)97 out = {}98 # --- référence (semi-log)99 out["semi-log (référence)"] = {"r2": ref["r2"], "r2_ajuste": ref["r2_ajuste"], "sigma": ref["sigma"], "k": ref["k"],100 "ajuste": np.exp(ref["X"].values @ ref["beta"]) * _duan(ref["residus"]), "sujet": float(cas1.predire(ref, sd)[0])}101 # --- linéaire102 XL = matrice_lineaire(v); rL = h.mco(XL.values, P)103 XLs = matrice_lineaire(sd)[XL.columns].values104 out["linéaire"] = {"r2": rL["r2"], "r2_ajuste": rL["r2_ajuste"], "sigma": float(np.sqrt(rL["sigma2"])), "k": rL["k"],105 "ajuste": XL.values @ rL["beta"], "sujet": float((XLs @ rL["beta"])[0]), "mco": rL, "X": XL}106 # --- log-log107 XG = matrice_loglog(v); rG = h.mco(XG.values, lnP)108 XGs = matrice_loglog(sd)[XG.columns].values109 out["log-log"] = {"r2": rG["r2"], "r2_ajuste": rG["r2_ajuste"], "sigma": float(np.sqrt(rG["sigma2"])), "k": rG["k"],110 "ajuste": np.exp(XG.values @ rG["beta"]) * _duan(rG["residus"]), "sujet": float(np.exp(XGs @ rG["beta"])[0] * _duan(rG["residus"])),111 "mco": rG, "X": XG}112 # --- Box–Cox113 X0 = ref["X"].values; bc = boxcox_profil(P, X0)114 rB = h.mco(X0, bc["y"]); lam = bc["lambda"]115 def inv(yl): return np.exp(yl) if abs(lam) < 1e-8 else np.power(np.maximum(lam * yl + 1, 1e-9), 1 / lam)116 Xs0 = cas1.matrice(sd)[ref["X"].columns].values117 # retransformation par lissage de Duan généralisé : moyenne des inv(ŷ + e_i)118 aj_bc = np.array([np.mean(inv(yh + rB["residus"])) for yh in (X0 @ rB["beta"])])119 out["Box–Cox"] = {"r2": rB["r2"], "r2_ajuste": rB["r2_ajuste"], "sigma": float(np.sqrt(rB["sigma2"])), "k": rB["k"], "lambda": lam,120 "ic95": bc["ic95"], "ajuste": aj_bc, "sujet": float(np.mean(inv((Xs0 @ rB["beta"])[0] + rB["residus"]))), "mco": rB, "boxcox": bc}121 # --- sans secteurs122 XS = matrice_sans_secteurs(v); rS = h.mco(XS.values, lnP)123 XSs = matrice_sans_secteurs(sd)[XS.columns].values124 out["semi-log sans secteurs"] = {"r2": rS["r2"], "r2_ajuste": rS["r2_ajuste"], "sigma": float(np.sqrt(rS["sigma2"])), "k": rS["k"],125 "ajuste": np.exp(XS.values @ rS["beta"]) * _duan(rS["residus"]), "sujet": float(np.exp(XSs @ rS["beta"])[0] * _duan(rS["residus"])),126 "mco": rS, "X": XS, "coef": pd.Series(rS["beta"], index=XS.columns), "se": pd.Series(rS["se_hc1"], index=XS.columns)}127 # --- indicatrices trimestrielles128 XT = matrice_trimestres(v); rT = h.mco(XT.values, lnP)129 # sujet : date d'évaluation 2026-09-01 → trimestre 2026T3 (= dernier trimestre observé)130 sdT = cas1.matrice(sd).drop(columns=["mois"])131 for c in XT.columns:132 if c.startswith("trim_"): sdT[c] = 0.0133 sdT[f"trim_{int(trimestre(sd['date_vente']).max())}"] = 1.0134 XTs = sdT[XT.columns].values135 out["semi-log, indicatrices trimestrielles"] = {"r2": rT["r2"], "r2_ajuste": rT["r2_ajuste"], "sigma": float(np.sqrt(rT["sigma2"])), "k": rT["k"],136 "ajuste": np.exp(XT.values @ rT["beta"]) * _duan(rT["residus"]), "sujet": float(np.exp(XTs @ rT["beta"])[0] * _duan(rT["residus"])),137 "mco": rT, "X": XT}138 for nom, d in out.items():139 d["rmse"] = _rmse(d["ajuste"], P); d["ratios"] = ratios(d["ajuste"], P)140 return out141142# ----------------------------------------------------------------------------- diagnostics143def breusch_pagan(X: np.ndarray, res: np.ndarray) -> dict:144 """Test de Breusch–Pagan (1979), version studentisée de Koenker : LM = n·R² de la régression de e² sur X ; χ²(k − 1)."""145 e2 = res ** 2146 r = h.mco(X, e2, robuste=False)147 n, k = X.shape; LM = n * r["r2"]148 return {"LM": float(LM), "ddl": k - 1, "p": float(1 - stats.chi2.cdf(LM, k - 1))}149150def vif(X: pd.DataFrame) -> pd.Series:151 cols = [c for c in X.columns if c != "const"]152 out = {}153 for c in cols:154 autres = X[[k for k in X.columns if k != c]].values155 r = h.mco(autres, X[c].values, robuste=False)156 out[c] = 1.0 / max(1.0 - r["r2"], 1e-12)157 return pd.Series(out)158159def cook(X: np.ndarray, res: np.ndarray, sigma2: float) -> dict:160 n, k = X.shape161 H = X @ np.linalg.pinv(X.T @ X) @ X.T162 hii = np.diag(H)163 D = res ** 2 * hii / (k * sigma2 * (1 - hii) ** 2)164 return {"D": D, "h": hii, "seuil": 4.0 / n, "n_sup": int((D > 4.0 / n).sum())}165166def moran_residus(v: pd.DataFrame, res: np.ndarray, k: int = 8) -> dict:167 coords = np.column_stack([(v["lat"] - 45.428) * 111.0, (v["lon"] + 75.714) * 78.0])168 W = h.matrice_poids_knn(coords, k)169 return h.moran_i(res, W)170171def diagnostics(v: pd.DataFrame, ref: dict) -> dict:172 X = ref["X"]; res = ref["residus"]173 jb = stats.jarque_bera(res)174 return {"breusch_pagan": breusch_pagan(X.values, res), "jarque_bera": {"stat": float(jb.statistic), "p": float(jb.pvalue)},175 "skew": float(stats.skew(res)), "kurtosis": float(stats.kurtosis(res, fisher=False)),176 "vif": vif(X), "cook": cook(X.values, res, ref["sigma"] ** 2), "moran": moran_residus(v, res),177 "res_min": float(res.min()), "res_max": float(res.max()), "pct_abs_sup_15": float(np.mean(np.abs(res) > 0.15) * 100)}178179# ----------------------------------------------------------------------------- prix implicites au sujet180def prix_implicites_sujet(ref: dict, s: dict) -> pd.DataFrame:181 """Prix implicites en dollars au sujet, dérivés du semi-log autour de la prédiction corrigée P̂_s (Duan).182 Continues : ∂P/∂x = β P̂ (× dérivée de la transformation) ; indicatrices : P̂ [exp(β) − 1] (Halvorsen–Palmquist) et Kennedy."""183 beta = pd.Series(ref["beta"], index=ref["X"].columns); se = pd.Series(np.sqrt(np.diag(ref["V_hc1"])), index=ref["X"].columns)184 sd = cas1.sujet_df(s); Ps = float(cas1.predire(ref, sd)[0]); xs = cas1.matrice(sd).iloc[0]185 S = float(s["superficie_habitable_m2"]); T = float(s["superficie_terrain_m2"]); age = float(s["age"])186 rows = []187 def cont(nom, unite, b, scale, se_b):188 rows.append({"caracteristique": nom, "unite": unite, "type": "continue", "beta": b, "prix_implicite": b * Ps * scale, "se": abs(scale) * Ps * se_b})189 cont("Superficie habitable", "$/m²", beta["ln_S"], 1.0 / S, se["ln_S"])190 cont("Superficie du terrain", "$/m²", beta["ln_T"], 1.0 / T, se["ln_T"])191 d_age = beta["age"] + 2 * beta["age2_100"] * age / 100.0192 se_age = math.sqrt(ref["V_hc1"][list(beta.index).index("age"), list(beta.index).index("age")]193 + (2 * age / 100) ** 2 * ref["V_hc1"][list(beta.index).index("age2_100"), list(beta.index).index("age2_100")]194 + 2 * (2 * age / 100) * ref["V_hc1"][list(beta.index).index("age"), list(beta.index).index("age2_100")])195 rows.append({"caracteristique": "Âge (une année de plus, à 28 ans)", "unite": "$/an", "type": "continue", "beta": d_age, "prix_implicite": d_age * Ps, "se": Ps * se_age})196 cont("Conditions du marché (un mois)", "$/mois", beta["mois"], 1.0, se["mois"])197 rows.append({"caracteristique": "Distance au centre (sujet à 2,1 km : dans le plateau)", "unite": "$/km", "type": "continue", "beta": beta["dist_centre"], "prix_implicite": 0.0, "se": 0.0})198 def disc(nom, col, direction=+1):199 b = beta[col] * direction; v = se[col] ** 2200 hp = Ps * (math.exp(b) - 1); ken = Ps * (math.exp(b - v / 2) - 1)201 rows.append({"caracteristique": nom, "unite": "$", "type": "indicatrice", "beta": beta[col], "prix_implicite": hp, "kennedy": ken, "se": Ps * math.exp(b) * se[col]})202 rows.append({"caracteristique": "Qualité (un échelon de plus)", "unite": "$", "type": "indicatrice", "beta": beta["qualite"],203 "prix_implicite": Ps * (math.exp(beta["qualite"]) - 1), "kennedy": Ps * (math.exp(beta["qualite"] - se["qualite"] ** 2 / 2) - 1), "se": Ps * math.exp(beta["qualite"]) * se["qualite"]})204 disc("Garage (une place de plus)", "garage")205 disc("Sous-sol fini (contribution : fini vs non fini)", "sous_sol")206 disc("Piscine (ajout)", "piscine")207 disc("Cottage à étages vs plain-pied", "etages")208 disc("Secteur Aylmer vs Hull", "sect_Aylmer")209 disc("Secteur Gatineau vs Hull", "sect_Gatineau")210 df = pd.DataFrame(rows); df["P_sujet"] = Ps211 return df212213# ----------------------------------------------------------------------------- prédiction et intervalle214def prediction_sujet(ref: dict, s: dict) -> dict:215 sd = cas1.sujet_df(s); X = cas1.matrice(sd)[ref["X"].columns].values216 lnp = float((X @ ref["beta"])[0]); xVx = float((X @ ref["V_hc1"] @ X.T)[0, 0])217 duan = _duan(ref["residus"]); logn = math.exp(ref["sigma"] ** 2 / 2)218 ip = cas1.intervalle_prediction(ref, sd).iloc[0]219 ip80 = cas1.intervalle_prediction(ref, sd, 0.80).iloc[0]220 # prédiction à la date de référence du rôle fictif (2024-07-01) pour la comparaison avec la valeur au rôle221 sd_ref = cas1.sujet_df(s, "2024-07-01"); p_ref = float(cas1.predire(ref, sd_ref)[0])222 return {"ln_pred": lnp, "naif": math.exp(lnp), "duan": math.exp(lnp) * duan, "facteur_duan": duan, "facteur_lognormal": logn,223 "lognormal": math.exp(lnp) * logn, "xVx": xVx, "sigma2": ref["sigma"] ** 2, "sd_pred": math.sqrt(ref["sigma"] ** 2 + xVx),224 "ic95": (float(ip["bas"]), float(ip["haut"])), "ic80": (float(ip80["bas"]), float(ip80["haut"])),225 "pred_2024_07": p_ref, "ratio_role": float(s["valeur_role_2025"]) / p_ref}226227# ----------------------------------------------------------------------------- validation croisée et ratios228def validation_croisee(v: pd.DataFrame, nfolds: int = 10, graine: int = GRAINE) -> dict:229 rng = np.random.default_rng(graine)230 X = cas1.matrice(v).values; lnP = np.log(v["prix_vente"].astype(float).to_numpy()); P = np.exp(lnP)231 idx = rng.permutation(len(P)); folds = np.array_split(idx, nfolds)232 pred = np.empty_like(P)233 for f in folds:234 tr = np.setdiff1d(idx, f); r = h.mco(X[tr], lnP[tr])235 pred[f] = np.exp(X[f] @ r["beta"]) * _duan(r["residus"])236 return {"pred": pred, "rmse": _rmse(pred, P), "ratios": ratios(pred, P), "erreur_abs_pct_mediane": float(np.median(np.abs(pred / P - 1)) * 100)}237238def ratios_role_fictif(v: pd.DataFrame, ref: dict) -> dict:239 """Étude de ratios du rôle fictif : valeur au rôle / prix ramené à la date de référence (2024-07-01) par la tendance du modèle."""240 bm = ref["beta"][list(ref["X"].columns).index("mois")]241 prix_ref = v["prix_vente"].astype(float).to_numpy() * np.exp(-bm * ref["X"]["mois"].to_numpy())242 return {"ratios": ratios(v["valeur_role_2025"].to_numpy(float), prix_ref), "prix_ref": prix_ref,243 "r": v["valeur_role_2025"].to_numpy(float) / prix_ref}244245def indice_trimestriel(variantes: dict, ref: dict, v: pd.DataFrame) -> pd.DataFrame:246 rT = variantes["semi-log, indicatrices trimestrielles"]["mco"]; XT = variantes["semi-log, indicatrices trimestrielles"]["X"]247 cols = [c for c in XT.columns if c.startswith("trim_")]248 bq = np.array([rT["beta"][list(XT.columns).index(c)] for c in cols]); seq = np.array([rT["se_hc1"][list(XT.columns).index(c)] for c in cols])249 idx = h.indice_dummies_temps(bq)250 q = np.arange(len(cols) + 1)251 # tendance : mois au milieu de chaque trimestre depuis 2024-07-01, normalisée à 2023T1252 mid = pd.to_datetime([f"{2023 + k // 4}-{3 * (k % 4) + 2:02d}-15" for k in q])253 mois_mid = (mid - cas1.REF_DATE).days / 30.4375254 bm = ref["beta"][list(ref["X"].columns).index("mois")]255 tend = 100 * np.exp(bm * (mois_mid - mois_mid[0]))256 lab = [f"{2023 + k // 4}T{k % 4 + 1}" for k in q]257 return pd.DataFrame({"trimestre": lab, "indice_dummies": idx, "bas": np.concatenate([[100], 100 * np.exp(bq - 1.96 * seq)]),258 "haut": np.concatenate([[100], 100 * np.exp(bq + 1.96 * seq)]), "indice_tendance": tend,259 "n_ventes": [int((trimestre(v["date_vente"]) == k).sum()) for k in q]})260261# ----------------------------------------------------------------------------- sensibilité : sous-échantillons262def sensibilite_echantillon(v: pd.DataFrame, s: dict, ref: dict, diag: dict) -> pd.DataFrame:263 sd = cas1.sujet_df(s); Xs_full = cas1.matrice(sd)264 def fit(sub: pd.DataFrame, nom: str):265 X = cas1.matrice(sub); y = np.log(sub["prix_vente"].astype(float)).values266 garder = [c for c in X.columns if X[c].std() > 0 or c == "const"]267 X = X[garder]; r = h.mco(X.values, y)268 p = float(np.exp((Xs_full[garder].values @ r["beta"])[0]) * _duan(r["residus"]))269 bS = r["beta"][garder.index("ln_S")]270 return {"echantillon": nom, "n": len(sub), "r2": r["r2"], "sigma": float(np.sqrt(r["sigma2"])), "beta_ln_S": bS, "sujet": p}271 rows = [fit(v, "600 ventes (référence)")]272 rows.append(fit(v[v["secteur"] == "Hull"], "Secteur Hull seulement"))273 mois_av = (pd.Timestamp(s["date_evaluation"]) - v["date_vente"]).dt.days / 30.4375274 rows.append(fit(v[mois_av <= 18], "Ventes des 18 derniers mois"))275 rows.append(fit(v[v["type"] == "à étages"], "Cottages à étages seulement"))276 rows.append(fit(v[diag["cook"]["D"] <= diag["cook"]["seuil"]], "Sans les observations influentes (Cook > 4/n)"))277 rows.append(fit(v[(v["superficie_habitable_m2"].between(110, 230))], "Superficie de 110 à 230 m²"))278 return pd.DataFrame(rows)279280# ----------------------------------------------------------------------------- exercices (chiffres des solutions)281def exercices(ref: dict, pred: dict, vari: dict, idx: pd.DataFrame) -> dict:282 ex = {}283 # 18.1 : interprétation des coefficients284 ex["1"] = {"elasticite_S": 0.539, "garage_hp": math.exp(0.0324) - 1, "buck_hp": math.exp(-0.250) - 1, "buck_kennedy": math.exp(-0.250 - 0.072 ** 2 / 2) - 1,285 "mois_annuel": (1 + 0.0036) ** 12 - 1}286 # 18.2 : prix implicites à partir de coefficients arrondis287 P, S = 638_000.0, 165.0288 ex["2"] = {"pS_marginal": 0.539 * P / S, "pS_10m2_exact": P * ((175 / 165) ** 0.539 - 1), "pS_10m2_lineaire": 0.539 * P / S * 10,289 "age28": (-0.0051 + 2 * 0.0026 * 28 / 100) * P, "age60": (-0.0051 + 2 * 0.0026 * 60 / 100) * P, "age_min": 0.0051 / (2 * 0.0026 / 100),290 "piscine_hp": P * (math.exp(0.0263) - 1)}291 # 18.3 : Box–Cox292 ex["3"] = {"bc_1": (500_000 ** 1 - 1) / 1, "bc_05": (500_000 ** 0.5 - 1) / 0.5, "bc_0": math.log(500_000), "lambda_hat": vari["Box–Cox"]["lambda"], "ic": vari["Box–Cox"]["ic95"]}293 # 18.4 : Duan et intervalle avec chiffres arrondis294 lnp, sig, xVx, duan = 13.3627, 0.0830, round(pred["xVx"], 5), 1.0034295 sd = math.sqrt(sig ** 2 + xVx); z = 1.959964296 ex["4"] = {"naif": math.exp(lnp), "duan": math.exp(lnp) * duan, "lognormal": math.exp(lnp) * math.exp(sig ** 2 / 2), "sd": sd,297 "bas": math.exp(lnp - z * sd), "haut": math.exp(lnp + z * sd), "xVx": xVx, "part_xVx": xVx / (sig ** 2 + xVx),298 "demi_largeur_pct": (math.exp(z * sd) - 1) * 100}299 # 18.5 : étude de ratios (dix ventes fictives)300 prix = np.array([312_000, 358_500, 401_000, 445_000, 489_500, 532_000, 598_000, 655_000, 742_000, 895_000], float)301 val = np.array([298_000, 329_000, 371_000, 405_000, 452_000, 478_000, 539_000, 576_000, 645_000, 752_000], float)302 r5 = ratios(val, prix); ex["5"] = {"ratios": val / prix, **r5}303 # 18.6 : variables omises (coefficients sans secteurs)304 cS = vari["semi-log sans secteurs"]["coef"]; cR = pd.Series(ref["beta"], index=ref["X"].columns)305 ex["6"] = {"dist_ref": cR["dist_centre"], "dist_sans": cS["dist_centre"], "lnT_ref": cR["ln_T"], "lnT_sans": cS["ln_T"],306 "sigma_ref": ref["sigma"], "sigma_sans": vari["semi-log sans secteurs"]["sigma"], "sujet_ref": vari["semi-log (référence)"]["sujet"],307 "sujet_sans": vari["semi-log sans secteurs"]["sujet"], "r2_sans": vari["semi-log sans secteurs"]["r2"]}308 # variante demandée à l'exercice : omission de la qualité (tirée indépendamment des autres attributs dans la loi de génération)309 XQ = ref["X"].drop(columns=["qualite"]); rQ = h.mco(XQ.values, ref["y"])310 d = cas1.charger(); sdQ = cas1.matrice(cas1.sujet_df(d["sujet"]))[XQ.columns].values311 ex["6"].update({"sigma_sans_qualite": float(np.sqrt(rQ["sigma2"])), "r2_sans_qualite": rQ["r2"], "lnS_sans_qualite": rQ["beta"][list(XQ.columns).index("ln_S")],312 "sujet_sans_qualite": float(np.exp(sdQ @ rQ["beta"])[0] * _duan(rQ["residus"]))})313 # 18.7 : Rosen — prix implicite marginal et enchère314 ex["7"] = {"p_S": 0.539 * 638_000 / 165, "p_S_120": 0.539 * 638_000 * (120 / 165) ** 0.539 / 120, "p_S_220": 0.539 * 638_000 * (220 / 165) ** 0.539 / 220}315 # 18.8 : indice trimestriel vs tendance316 bm = ref["beta"][list(ref["X"].columns).index("mois")]317 ecart = (idx["indice_dummies"] / idx["indice_tendance"] - 1) * 100318 ex["8"] = {"g_mensuel": math.exp(bm) - 1, "g_annuel": math.exp(12 * bm) - 1, "ecart_abs_moyen_pct": float(ecart.abs().mean()), "ecart_max_pct": float(ecart.abs().max()),319 "G0404_tendance": 506_500 * (math.exp(bm * 15.2) - 1), "indice_fin": float(idx["indice_dummies"].iloc[-1]), "tendance_fin": float(idx["indice_tendance"].iloc[-1]),320 "indice_2025T2": float(idx.set_index("trimestre").loc["2025T2", "indice_dummies"]),321 "G0404_indices": 506_500 * (float(idx["indice_dummies"].iloc[-1]) / float(idx.set_index("trimestre").loc["2025T2", "indice_dummies"]) - 1)}322 return ex323324# ----------------------------------------------------------------------------- principal325def main(verbose: bool = True) -> dict:326 d = cas1.charger(); v = d["ventes"]; s = d["sujet"]327 ref = cas1.hedonique(v)328 coef = ref["coef"].copy(); coef["vrai"] = [VRAIS.get(c, np.nan) for c in coef.index]329 coef["ecart_t"] = (coef["coef"] - coef["vrai"]) / coef["se_hc1"]330 coef["ic_bas"] = coef["coef"] - 1.96 * coef["se_hc1"]; coef["ic_haut"] = coef["coef"] + 1.96 * coef["se_hc1"]331 coef["couvre"] = (coef["ic_bas"] <= coef["vrai"]) & (coef["vrai"] <= coef["ic_haut"])332 coef["etiquette"] = [ETIQUETTES[c] for c in coef.index]333 diag = diagnostics(v, ref)334 vari = ajuster_variantes(v, s, ref)335 pi = prix_implicites_sujet(ref, s)336 pred = prediction_sujet(ref, s)337 cv = validation_croisee(v)338 role = ratios_role_fictif(v, ref)339 idx = indice_trimestriel(vari, ref, v)340 sens = sensibilite_echantillon(v, s, ref, diag)341 ex = exercices(ref, pred, vari, idx)342 # comparables du ch. 16 : résidus du modèle343 sel = cas1.selection(); comps = sel["comparables"]344 res_comps = pd.Series(ref["residus"], index=v.index).loc[comps.index]345 P = v["prix_vente"].astype(float).to_numpy()346 out = {"ventes": v, "sujet": s, "modele": ref, "coef": coef, "diag": diag, "variantes": vari, "prix_implicites": pi, "prediction": pred,347 "cv": cv, "role": role, "indice": idx, "sensibilite": sens, "exercices": ex, "comparables": comps, "residus_comparables": res_comps,348 "ratios_modele": ratios(vari["semi-log (référence)"]["ajuste"], P)}349 if verbose:350 pd.set_option("display.width", 220); pd.set_option("display.max_columns", 30)351 print("== Modèle de référence (traite.cas1.hedonique) : n =", ref["n"], "k =", ref["k"], f"R² = {ref['r2']:.4f} R²aj = {ref['r2_ajuste']:.4f} σ = {ref['sigma']:.4f}")352 print(coef[["coef", "se_hc1", "t", "effet_pct", "vrai", "ecart_t", "couvre"]].round(4).to_string())353 print(f" couverture des IC 95 % : {int(coef['couvre'].sum())}/{int(coef['vrai'].notna().sum())} paramètres de génération couverts")354 print(f"== Croissance mensuelle {ref['croissance_mensuelle']*100:.3f} % ; annuelle {((1+ref['croissance_mensuelle'])**12-1)*100:.2f} % ; vraie {_g*100:.2f} %")355 bp = diag["breusch_pagan"]; print(f"== Breusch–Pagan (Koenker) LM = {bp['LM']:.2f}, ddl = {bp['ddl']}, p = {bp['p']:.3f}")356 print(f"== Jarque–Bera = {diag['jarque_bera']['stat']:.2f} (p = {diag['jarque_bera']['p']:.3f}) ; asymétrie {diag['skew']:.3f} ; aplatissement {diag['kurtosis']:.3f}")357 print(f"== Résidus : min {diag['res_min']:.3f}, max {diag['res_max']:.3f}, |e| > 0,15 : {diag['pct_abs_sup_15']:.1f} %")358 print("== VIF :\n", diag["vif"].round(2).to_string())359 ck = diag["cook"]; top = np.argsort(-ck["D"])[:3]360 print(f"== Cook : seuil 4/n = {ck['seuil']:.4f} ; {ck['n_sup']} observations au-dessus ; max {ck['D'].max():.4f} ; 3 plus grandes :", list(v['id_vente'].iloc[top]), np.round(ck['D'][top], 4))361 m = diag["moran"]; print(f"== Moran (résidus, 8 voisins) I = {m['I']:.4f}, z = {m['z']:.2f}, p = {m['p']:.3f}")362 print("== Variantes :")363 for nom, dd in vari.items():364 extra = f" λ = {dd['lambda']:.3f} IC [{dd['ic95'][0]:.3f}; {dd['ic95'][1]:.3f}]" if "lambda" in dd else ""365 print(f" {nom:42s} k={dd['k']:2d} R²={dd['r2']:.4f} R²aj={dd['r2_ajuste']:.4f} σ={dd['sigma']:.4f} RMSE={dd['rmse']:,.0f} $ COD={dd['ratios']['COD']:.2f} PRD={dd['ratios']['PRD']:.3f} PRB={dd['ratios']['PRB']:.3f} sujet={dd['sujet']:,.0f} ${extra}")366 cS = vari["semi-log sans secteurs"]367 print(f" sans secteurs : dist_centre {cS['coef']['dist_centre']:.4f} (réf {ref['beta'][list(ref['X'].columns).index('dist_centre')]:.4f}) ; ln_T {cS['coef']['ln_T']:.4f} (réf {ref['beta'][list(ref['X'].columns).index('ln_T')]:.4f})")368 print("== Prix implicites au sujet (P̂ = %.0f $) :" % pi["P_sujet"].iloc[0])369 print(pi[["caracteristique", "unite", "beta", "prix_implicite", "kennedy", "se"]].round(4).to_string(index=False))370 print("== Prédiction du sujet :", {k: (round(x, 5) if isinstance(x, float) and abs(x) < 10 else (tuple(round(y) for y in x) if isinstance(x, tuple) else round(x))) for k, x in pred.items()})371 print(f"== Ratios du modèle (in-sample) : COD {out['ratios_modele']['COD']:.2f}, PRD {out['ratios_modele']['PRD']:.4f}, PRB {out['ratios_modele']['PRB']:.4f}, médiane {out['ratios_modele']['ratio_median']:.4f}")372 print(f"== Validation croisée 10 plis : RMSE {cv['rmse']:,.0f} $, COD {cv['ratios']['COD']:.2f}, PRD {cv['ratios']['PRD']:.4f}, PRB {cv['ratios']['PRB']:.4f}, erreur absolue médiane {cv['erreur_abs_pct_mediane']:.2f} %")373 rr = role["ratios"]; print(f"== Rôle fictif vs prix ramenés au 2024-07-01 : médiane {rr['ratio_median']:.4f}, COD {rr['COD']:.2f}, PRD {rr['PRD']:.4f}, PRB {rr['PRB']:.4f}")374 print("== Indice trimestriel :\n", idx.round(2).to_string(index=False))375 print("== Sensibilité à l'échantillon :\n", sens.round(4).to_string(index=False))376 print("== Résidus des 5 comparables du ch. 16 :", dict(zip(comps["id_vente"], res_comps.round(4))), "moyenne", round(res_comps.mean(), 4))377 print("== Exercices :")378 for k, e in ex.items(): print(" ", k, {kk: (np.round(vv, 4) if not isinstance(vv, tuple) else tuple(round(x, 4) for x in vv)) for kk, vv in e.items()})379 return out380381if __name__ == "__main__":382 main()383