spb/traite-evaluation-qc

1#!/usr/bin/env python32"""Chapitre 20 — réconciliation et incertitude : combinaison des indications du cas 1 (Gatineau, données synthétiques).34Reprend sans les recalculer « à la main » : la sélection du chapitre 16 (traite.cas1.selection), la grille hédonique du5chapitre 17 (traite.cas1.ajustements_hedoniques), le modèle de référence et la prédiction du chapitre 186(traite.cas1.hedonique / predire / intervalle_prediction) et le coût neuf du chapitre 21 (ch21_cout_neuf.cout_neuf_sujet).7Ajoute : bootstrap des comparables (candidats, ajustements, les deux), énumération exhaustive des sous-ensembles de cinq8comparables, indication par le coût sous hypothèses déclarées, pondération (égale, inverse de la variance, moindres carrés9généralisés avec corrélation, raisonnée), fourchette de la conclusion, et « vérité » de la loi de génération.10Imprime tous les chiffres cités dans le chapitre.11Usage (racine du dépôt) : PYTHONPATH=analysis/python python3 analysis/python/cases/ch20_reconciliation.py"""12from __future__ import annotations13import sys, pathlib, math, itertools, datetime14import numpy as np, pandas as pd15HERE = pathlib.Path(__file__).resolve()16sys.path.insert(0, str(HERE.parents[1]))17sys.path.insert(0, str(HERE.parent))18from traite import cas1, cout as ct, comparables as cp19import ch21_cout_neuf as c212021ROOT = HERE.parents[3]22GRAINE = 2026103123B = 2000                                   # réplications bootstrap24PLANCHER_BRUT = 0.01                       # plancher de l'inverse du brut (ch. 17, eq. 17.poids)25GROUPES = ["conditions_marche", "localisation", "superficie", "terrain", "age", "qualite", "garage", "sous_sol", "piscine", "type"]26COLS_GROUPES = {"conditions_marche": ["mois"], "localisation": ["dist_centre"] + [f"sect_{x}" for x in cas1.SECTEURS[1:]],27                "superficie": ["ln_S"], "terrain": ["ln_T"], "age": ["age", "age2_100"], "qualite": ["qualite"],28                "garage": ["garage"], "sous_sol": ["sous_sol"], "piscine": ["piscine"], "type": ["etages"]}29# Indication par le coût : reprise du chapitre 24 (ch24_depreciation.py : terrain extrait des ventes récentes de Hull, dépréciation30# extraite du marché à 28 ans + désuétude fonctionnelle). Les « hypothèses d'usage » ci-dessous ne servent que de contre-exemple.31ALPHA_TERRAIN_USAGE = 0.30  # allocation (part du terrain) — hypothèse d'usage, contre-exemple32VIE_ECONOMIQUE_USAGE = 60   # âge-vie — hypothèse d'usage, contre-exemple33AGE_EFFECTIF = 2834B_COUT = B                  # réplications du bootstrap de l'indication par le coût (mêmes tirages des 600 ventes que le modèle)35# Prédictions du sujet sous les quatre spécifications « en logarithme » du ch. 18 (tab. 18.variantes : semi-log, log-log,36# Box–Cox, indicatrices trimestrielles) : sert à l'incertitude de spécification du modèle hédonique.37PRED_SPEC_CH18 = [637953.0, 627260.0, 635171.0, 649980.0]38POIDS_RAISONNES = {"comparaison": 0.40, "hedonique": 0.45, "cout": 0.15}39ARRONDI_CONCLUSION = 5000404142# ---------------------------------------------------------------- 1. indications héritées des chapitres 16 à 1843def base() -> dict:44    t = cas1.tout()45    m = t["modele"]; sel = t["selection"]; s = t["sujet"]; v = t["ventes"]46    cand = sel["candidats"].copy()                      # 16 candidats, triés par distance d'attributs croissante47    comps = sel["comparables"]48    adj16 = cas1.ajustements_hedoniques(m, cand, s)49    adj16["distance"] = cand["distance"].values50    res = pd.Series(m["residus"], index=v.index)51    adj16["residu_log"] = res.loc[cand.index].values52    adj16["retenu"] = np.arange(len(cand)) < len(comps)53    adj5 = adj16[adj16["retenu"]].copy()54    # indication par comparaison, exactement comme au chapitre 17 (bibliothèque traite.comparables)55    ca = [cp.ComparableAjuste(r["id_vente"], r["prix_vente"], {g: r[g] for g in GROUPES}) for _, r in adj5.iterrows()]56    w = cp.poids_inverse_brut(ca, PLANCHER_BRUT)57    rec = cp.reconciliation(ca, w)58    cols = list(m["X"].columns)59    Xs = cas1.matrice(cas1.sujet_df(s))[cols].values[0]60    Xc16 = cas1.matrice(cand)[cols].values61    G = np.zeros((len(cols), len(GROUPES)))62    for j, g in enumerate(GROUPES):63        for c in COLS_GROUPES[g]:64            G[cols.index(c), j] = 1.065    return {**t, "cand": cand, "comps": comps, "adj16": adj16, "adj5": adj5, "poids_comp": np.array(w), "reconciliation": rec,66            "V_comp": rec["valeur_ponderee"], "sd_comp_pondere": rec["ecart_type_pondere"], "V_hed": t["prediction_hedonique"],67            "ip_hed": t["intervalle"], "cols": cols, "Xs": Xs, "Xc16": Xc16, "P16": cand["prix_vente"].values.astype(float), "G": G,68            "X": m["X"].values, "y": m["y"], "beta": m["beta"], "V_hc1": m["V_hc1"], "sigma": m["sigma"]}697071def grille_rapide(beta: np.ndarray, Xs: np.ndarray, Xc: np.ndarray, P: np.ndarray, G: np.ndarray, plancher: float = PLANCHER_BRUT) -> dict:72    """Cascade hédonique vectorisée (identique à traite.cas1.ajustements_hedoniques) : Δ_k par groupe, bases cumulées,73    ajustements A_k = B_{k-1}(e^{Δ_k} − 1), prix ajusté P e^{ΣΔ}, brut Σ|A_k|/P, poids inverse du brut, valeur pondérée."""74    D = ((Xs[None, :] - Xc) * beta[None, :]) @ G           # n_c × n_groupes75    cum = np.cumsum(D, axis=1)76    bases = P[:, None] * np.exp(cum)77    bases_prec = np.hstack([P[:, None], bases[:, :-1]])78    A = bases_prec * (np.exp(D) - 1.0)79    P_aj = bases[:, -1]80    brut = np.abs(A).sum(axis=1) / P81    w = 1.0 / np.maximum(brut, plancher); w = w / w.sum()82    V = float(w @ P_aj)83    return {"P_aj": P_aj, "brut": brut, "poids": w, "valeur": V, "sd_pondere": float(math.sqrt(w @ (P_aj - V) ** 2))}848586def valeur_ponderee(P_aj: np.ndarray, brut: np.ndarray, plancher: float = PLANCHER_BRUT) -> float:87    w = 1.0 / np.maximum(brut, plancher)88    return float((w @ P_aj) / w.sum())899091# ---------------------------------------------------------------- 2. bootstrap des comparables92def bootstrap(b: dict, B: int = B, graine: int = GRAINE, k: int | None = None) -> dict:93    """Trois rééchantillonnages : (A) les 16 candidats avec remise puis sélection des 5 plus proches (modèle fixe) ;94    (B) les 600 ventes avec remise, réestimation du modèle, grille des 5 comparables fixes et prédiction du sujet ;95    (C) les deux à la fois (même réplication b pour la corrélation) ; variante paramétrique β* ~ N(β̂, V_HC1)."""96    rng = np.random.default_rng(graine)97    n_c = len(b["P16"]); k = k or len(b["comps"])98    g0 = grille_rapide(b["beta"], b["Xs"], b["Xc16"], b["P16"], b["G"])99    P_aj16, brut16 = g0["P_aj"], g0["brut"]100    # (A) candidats101    idxA = np.sort(rng.integers(0, n_c, size=(B, n_c)), axis=1)[:, :k]     # candidats triés par distance : rang = indice102    VA = np.array([valeur_ponderee(P_aj16[i], brut16[i]) for i in idxA])103    # (B) modèle (bootstrap par paires)104    X, y = b["X"], b["y"]; n = len(y)105    idxB = rng.integers(0, n, size=(B, n))106    betas = np.empty((B, X.shape[1])); duan = np.empty(B)107    for r in range(B):108        Xb, yb = X[idxB[r]], y[idxB[r]]109        beta_b, *_ = np.linalg.lstsq(Xb, yb, rcond=None)110        betas[r] = beta_b; duan[r] = float(np.mean(np.exp(yb - Xb @ beta_b)))111    pred = np.exp(betas @ b["Xs"]) * duan112    retenus = np.arange(k)113    VB = np.empty(B); VC = np.empty(B)114    for r in range(B):115        g = grille_rapide(betas[r], b["Xs"], b["Xc16"], b["P16"], b["G"])116        VB[r] = valeur_ponderee(g["P_aj"][retenus], g["brut"][retenus])117        VC[r] = valeur_ponderee(g["P_aj"][idxA[r]], g["brut"][idxA[r]])118    # variante paramétrique (coefficients tirés de la normale asymptotique, comparables fixes)119    betasP = rng.multivariate_normal(b["beta"], b["V_hc1"], size=B)120    VP = np.array([valeur_ponderee(*(lambda g: (g["P_aj"][retenus], g["brut"][retenus]))(grille_rapide(bp, b["Xs"], b["Xc16"], b["P16"], b["G"]))) for bp in betasP])121    def resume(x):122        q = np.percentile(x, [2.5, 5, 25, 50, 75, 95, 97.5])123        return {"moyenne": float(x.mean()), "ecart_type": float(x.std(ddof=1)), "q025": q[0], "q05": q[1], "q25": q[2], "q50": q[3], "q75": q[4], "q95": q[5], "q975": q[6]}124    rho = float(np.corrcoef(VC, pred)[0, 1])125    rhoB = float(np.corrcoef(VB, pred)[0, 1])126    return {"P_aj16": P_aj16, "brut16": brut16, "idxA": idxA, "idxB": idxB, "VA": VA, "VB": VB, "VC": VC, "VP": VP, "pred": pred, "betas": betas,127            "resume": {"candidats": resume(VA), "ajustements": resume(VB), "conjoint": resume(VC), "parametrique": resume(VP), "prediction": resume(pred)},128            "rho_conjoint_pred": rho, "rho_ajustements_pred": rhoB,129            "part_A_sous_hed": float((VA < b["V_hed"]).mean()), "part_C_sous_hed": float((VC < b["V_hed"]).mean()),130            "part_C_au_dessus_comp": float((VC > b["V_comp"]).mean())}131132133def exhaustif(b: dict, P_aj16: np.ndarray, brut16: np.ndarray, k: int = 5) -> dict:134    """Tous les sous-ensembles de k candidats parmi les 16 (C(16,5) = 4368), modèle fixe."""135    n_c = len(P_aj16)136    vals = []; n_retenus = []137    for sub in itertools.combinations(range(n_c), k):138        i = np.array(sub)139        vals.append(valeur_ponderee(P_aj16[i], brut16[i])); n_retenus.append(int((i < k).sum()))140    vals = np.array(vals); n_retenus = np.array(n_retenus)141    q = np.percentile(vals, [2.5, 5, 50, 95, 97.5])142    par_n = pd.DataFrame({"n_retenus": n_retenus, "valeur": vals}).groupby("n_retenus")["valeur"].agg(["count", "mean", "min", "max"])143    return {"valeurs": vals, "n": len(vals), "moyenne": float(vals.mean()), "ecart_type": float(vals.std(ddof=1)), "min": float(vals.min()),144            "max": float(vals.max()), "q025": q[0], "q05": q[1], "q50": q[2], "q95": q[3], "q975": q[4], "par_n_retenus": par_n,145            "rang_selection": int((vals < b["V_comp"]).sum()), "part_sous_hed": float((vals < b["V_hed"]).mean())}146147148def elargissement(b: dict, P_aj16: np.ndarray, brut16: np.ndarray) -> pd.DataFrame:149    """Indication par comparaison lorsque l'on retient les k premiers candidats (k = 5 … 16) : valeur pondérée (inverse du brut),150    moyenne simple, résidu moyen (exp(e) − 1) des comparables retenus. Le modèle est fixe (arbitrage biais–variance du ch. 16)."""151    res = b["adj16"]["residu_log"].values152    rows = []153    for k in range(5, len(P_aj16) + 1):154        i = np.arange(k)155        rows.append({"k": k, "valeur_ponderee": valeur_ponderee(P_aj16[i], brut16[i]), "moyenne": float(P_aj16[i].mean()),156                     "residu_moyen_pct": float(np.mean(np.exp(res[i]) - 1)), "brut_moyen": float(brut16[i].mean()),157                     "distance_max": float(b["adj16"]["distance"].values[k - 1])})158    return pd.DataFrame(rows)159160161# ---------------------------------------------------------------- 3. indication par le coût (reprise du chapitre 24)162def _ch24() -> dict:163    import io, contextlib164    import ch24_depreciation as c24165    with contextlib.redirect_stdout(io.StringIO()):166        return c24.main()167168169def cout(b: dict) -> dict:170    """Indication par le coût du chapitre 24 : V = V_L + C_N − D_extraite(28 ans) − D_fonctionnelle. Ajoute les diagnostics171    d'extraction aux deux autres indications, la lecture « hypothèses d'usage » (allocation 30 %, âge-vie 60 ans) comme172    contre-exemple, et l'incertitude de spécification (variantes d'extraction du ch. 24)."""173    import ch24_depreciation as c24174    r24 = _ch24()175    CN = r24["cn"]["total"]; vs = r24["valeur"]; V_L = r24["vl_retenue"]["retenue"]; fonc = r24["salle_eau"]["desuetude"]176    D_ext = vs["d_extraite"]; D_tot = vs["B_extraction"]["depreciation_totale"]177    # variantes d'extraction du ch. 24 (fenêtre des ventes récentes, élasticité, qualité) → valeur sous chaque variante178    sx = r24["sens_extraction"].copy()179    sx["valeur"] = [ct.valeur_par_cout(float(vl), CN, float(t) * CN, fonc)["valeur"] for vl, t in zip(sx["vl_hull"], sx["taux_28"])]180    # contre-exemple : hypothèses d'usage (allocation sur l'indication par comparaison, âge-vie)181    V_L_usage = ALPHA_TERRAIN_USAGE * b["V_comp"]; D_usage = ct.depreciation_age_vie(CN, AGE_EFFECTIF, VIE_ECONOMIQUE_USAGE)182    usage = ct.valeur_par_cout(V_L_usage, CN, D_usage)183    return {"r24": r24, "CN": CN, "V_L": V_L, "V_L_mediane": r24["vl_retenue"]["mediane"], "n_terrain": r24["vl_retenue"]["n"],184            "part_terrain": r24["vl_retenue"]["part"], "D_ext": D_ext, "taux_ext": vs["taux_extrait_28"], "D_fonc": fonc, "D_tot": D_tot,185            "taux_tot": vs["B_extraction"]["taux_depreciation"], "contribution_batiment": vs["B_extraction"]["contribution_batiment"],186            "valeur": vs["valeur_retenue"], "arrondie": vs["valeur_arrondie"], "ventilee": vs["A_ventilee"]["valeur"],187            "taux_ventilee": vs["A_ventilee"]["taux_depreciation"], "tornade": r24["tornade"], "sens_extraction": sx,188            "sd_spec": float(sx["valeur"].std(ddof=1)), "cout_neuf_plus_terrain": CN + V_L,189            "extraction_hed": ct.extraction_depreciation_marche(b["V_hed"], V_L, CN, AGE_EFFECTIF),190            "extraction_comp": ct.extraction_depreciation_marche(b["V_comp"], V_L, CN, AGE_EFFECTIF),191            "usage": {"alpha": ALPHA_TERRAIN_USAGE, "vie": VIE_ECONOMIQUE_USAGE, "V_L": V_L_usage, "D": D_usage, "valeur": usage["valeur"],192                      "extraction_hed": ct.extraction_depreciation_marche(b["V_hed"], V_L_usage, CN, AGE_EFFECTIF)},193            "fonc_fn": fonc}194195196def bootstrap_cout(b: dict, bs: dict, co: dict, B: int = B_COUT) -> dict:197    """Bootstrap de l'indication par le coût du ch. 24 : mêmes rééchantillons des 600 ventes que le modèle (idxB) ; pour chaque198    réplication, prix ramenés à la date par la tendance du modèle rééstimé, terrain extrait des ventes récentes de Hull,199    dépréciation extraite à 28 ans (ajustement quadratique), désuétude fonctionnelle fixe. Les ventes d'un secteur sans vente200    récente dans la réplication (terrain non extractible) sont exclues de l'ajustement du profil."""201    import ch24_depreciation as c24202    v = b["ventes"]; CN = co["CN"]; fonc = co["D_fonc"]; i_mois = b["cols"].index("mois")203    VK = np.empty(B); VL = np.empty(B); T28 = np.empty(B)204    for r in range(B):205        vb = v.iloc[bs["idxB"][r]].reset_index(drop=True)206        mb = {"croissance_mensuelle": float(np.exp(bs["betas"][r][i_mois]) - 1.0)}207        w, vl = c24.extraction(vb, mb, CN)208        w = w.dropna(subset=["taux"])          # secteurs sans vente récente dans la réplication : pas de terrain, exclus du profil209        q = c24.ajuster_profil(w)["quadratique"]; t28 = float(c24.taux_quadratique(q, AGE_EFFECTIF))210        VL[r] = float(vl.loc["Hull", "vl_mediane"]); T28[r] = t28211        VK[r] = ct.valeur_par_cout(VL[r], CN, t28 * CN, fonc)["valeur"]212    qv = np.percentile(VK, [2.5, 5, 50, 95, 97.5])213    return {"VK": VK, "VL": VL, "T28": T28,214            "resume": {"moyenne": float(VK.mean()), "ecart_type": float(VK.std(ddof=1)), "q025": qv[0], "q05": qv[1], "q50": qv[2], "q95": qv[3], "q975": qv[4]},215            "sd_VL": float(VL.std(ddof=1)), "sd_T28": float(T28.std(ddof=1)),216            "rho_cout_pred": float(np.corrcoef(VK, bs["pred"][:B])[0, 1]), "rho_cout_comp": float(np.corrcoef(VK, bs["VC"][:B])[0, 1])}217218219# ---------------------------------------------------------------- 4. combinaison220def poids_optimaux_deux(s1: float, s2: float, rho: float) -> float:221    """Poids optimal de l'indication 1 parmi deux indications corrélées (Bates–Granger) : w1 = (σ2² − ρσ1σ2)/(σ1² + σ2² − 2ρσ1σ2)."""222    return (s2 ** 2 - rho * s1 * s2) / (s1 ** 2 + s2 ** 2 - 2 * rho * s1 * s2)223224225def poids_gls(Sigma: np.ndarray) -> np.ndarray:226    """w = Σ⁻¹1 / 1'Σ⁻¹1 (variance minimale sous Σw = 1, poids non contraints)."""227    inv = np.linalg.inv(Sigma); un = np.ones(Sigma.shape[0])228    w = inv @ un229    return w / w.sum()230231232def combinaison(b: dict, bs: dict, co: dict, bk: dict, graine: int = GRAINE) -> dict:233    V = np.array([b["V_comp"], b["V_hed"], co["valeur"]])234    sd_spec = float(np.std(PRED_SPEC_CH18, ddof=1))235    sd_pred_boot = bs["resume"]["prediction"]["ecart_type"]236    s = np.array([bs["resume"]["conjoint"]["ecart_type"], math.sqrt(sd_pred_boot ** 2 + sd_spec ** 2), math.sqrt(bk["resume"]["ecart_type"] ** 2 + co["sd_spec"] ** 2)])237    rho = bs["rho_conjoint_pred"]; rho_ck = bk["rho_cout_comp"]; rho_hk = bk["rho_cout_pred"]238    R = np.array([[1.0, rho, rho_ck], [rho, 1.0, rho_hk], [rho_ck, rho_hk, 1.0]])239    Sigma = R * np.outer(s, s)240    schemas = {241        "comparaison seule": np.array([1.0, 0.0, 0.0]),242        "hédonique seule": np.array([0.0, 1.0, 0.0]),243        "égale (un tiers chacune)": np.ones(3) / 3,244        "inverse de la variance (indépendance)": (1 / s ** 2) / (1 / s ** 2).sum(),245        "variance minimale avec corrélation (MCG)": poids_gls(Sigma),246        "raisonnée (retenue)": np.array([POIDS_RAISONNES["comparaison"], POIDS_RAISONNES["hedonique"], POIDS_RAISONNES["cout"]]),247    }248    tab = pd.DataFrame([{"schema": k, "w_comp": w[0], "w_hed": w[1], "w_cout": w[2], "valeur": float(w @ V), "ecart_type": float(math.sqrt(w @ Sigma @ w))}249                        for k, w in schemas.items()])250    # distribution bootstrap de la conclusion retenue (réplication conjointe : comparaison et modèle issus du même tirage ; coût tiré de la grille)251    n_b = len(bk["VK"]); cout_b = bk["VK"]252    w = schemas["raisonnée (retenue)"]253    Vb = w[0] * bs["VC"][:n_b] + w[1] * bs["pred"][:n_b] + w[2] * cout_b254    q = np.percentile(Vb, [2.5, 5, 50, 95, 97.5])255    V_raisonnee = float(w @ V)256    arrondi = ARRONDI_CONCLUSION257    conclusion = round(V_raisonnee / arrondi) * arrondi258    fourchette90 = (math.floor(q[1] / arrondi) * arrondi, math.ceil(q[3] / arrondi) * arrondi)259    # fourchette de valeur conclue : encadrée par les deux indications crédibles (comparaison, hédonique), arrondies vers l'extérieur260    fourchette_conclue = (math.floor(min(V[0], V[1]) / arrondi) * arrondi, math.ceil(max(V[0], V[1]) / arrondi) * arrondi)261    # courbes de sensibilité aux poids : w_hed de 0 à 1 − w_cout, pour trois poids du coût262    wh = np.linspace(0, 1, 101)263    courbes = {wc: np.array([ (1 - h - wc) * V[0] + h * V[1] + wc * V[2] if h <= 1 - wc + 1e-12 else np.nan for h in wh]) for wc in (0.0, 0.15, 0.30)}264    # moyenne simple sous les hypothèses d'usage du coût (contre-exemple)265    V_usage = np.array([b["V_comp"], b["V_hed"], co["usage"]["valeur"]])266    # variance de la combinaison de deux indications (comparaison, hédonique) selon ρ267    var_deux = {r: np.array([((1 - h) ** 2) * s[0] ** 2 + h ** 2 * s[1] ** 2 + 2 * h * (1 - h) * r * s[0] * s[1] for h in wh]) for r in (0.0, 0.3, 0.6, 0.9)}268    w_opt_deux = {r: poids_optimaux_deux(s[0], s[1], r) for r in (0.0, rho, 0.3, 0.6, 0.9)}269    return {"V": V, "s": s, "sd_spec_hed": sd_spec, "sd_pred_boot": sd_pred_boot, "rho": rho, "R": R, "Sigma": Sigma, "schemas": schemas, "tableau": tab,270            "moyenne_simple_usage": float(V_usage.mean()), "inverse_variance_usage": None,271            "Vb": Vb, "cout_b": cout_b, "q": {"q025": q[0], "q05": q[1], "q50": q[2], "q95": q[3], "q975": q[4]}, "V_raisonnee": V_raisonnee,272            "conclusion": conclusion, "fourchette90": fourchette90, "fourchette_conclue": fourchette_conclue,273            "largeur_relative_conclue": (fourchette_conclue[1] - fourchette_conclue[0]) / conclusion, "wh": wh, "courbes": courbes, "var_deux": var_deux, "w_opt_deux": w_opt_deux,274            "ecart_conclusion_comp": conclusion / b["V_comp"] - 1, "ecart_conclusion_hed": conclusion / b["V_hed"] - 1}275276277# ---------------------------------------------------------------- 5. ce que révèle la loi de génération278def verite(b: dict) -> dict:279    """Espérance conditionnelle « vraie » du prix du sujet selon la loi déclarée dans data/synthetic/README.md ;280    la constante (normalisation de la médiane) est estimée sur les 600 ventes (erreur-type σ/√n ≈ 0,35 %)."""281    from traite.synth import _mois_depuis282    v = b["ventes"]; s = b["sujet"]283    med = 475000.0; med_sect = {"Gatineau": 464900.0, "Aylmer": 530000.0, "Buckingham": 395000.0, "Masson-Angers": 395000.0}284    prime = {"Hull": 0.0, **{k: math.log(m / med) for k, m in med_sect.items()}}; prime["Masson-Angers"] -= 0.02285    g = 1.0906 ** 0.5 - 1286    ref = _mois_depuis([datetime.date(2024, 7, 1)])[0]287    def f(df, dates):288        mois = _mois_depuis([d.date() for d in pd.to_datetime(dates)])289        return (0.55 * np.log(df["superficie_habitable_m2"].values / 140) + 0.10 * np.log(df["superficie_terrain_m2"].values / 560)290                - 0.0055 * df["age"].values + 0.00003 * df["age"].values ** 2 + 0.06 * (df["qualite_1a5"].values - 3)291                + 0.035 * df["garage_places"].values + 0.03 * df["sous_sol_fini"].values + 0.02 * df["piscine"].values292                + np.array([prime[x] for x in df["secteur"].values]) - 0.006 * np.clip(df["dist_centre_km"].values - 3, 0, None)293                + 0.02 * (df["type"].values == "à étages") + math.log1p(g) * (mois - ref) / 12)294    fv = f(v, v["date_vente"]); e = np.log(v["prix_vente"].values) - fv295    c_hat = float(e.mean()); se_c = float(e.std(ddof=1) / math.sqrt(len(e)))296    sd = cas1.sujet_df(s); fs = float(f(sd, sd["date_vente"])[0])297    V_true = math.exp(fs + c_hat) * math.exp(0.085 ** 2 / 2)298    return {"c_hat": c_hat, "se_c": se_c, "sigma_e": float(e.std(ddof=1)), "V_true": V_true,299            "V_true_bas": V_true * math.exp(-1.96 * se_c), "V_true_haut": V_true * math.exp(1.96 * se_c),300            "ecart_comp": b["V_comp"] / V_true - 1, "ecart_hed": b["V_hed"] / V_true - 1}301302303# ---------------------------------------------------------------- 6. exercices (chiffres des solutions)304def exercices(b: dict, bs: dict, co: dict, cb: dict) -> dict:305    out = {}306    # ex. 1 : deux indications indépendantes307    V1, s1, V2, s2 = 582000.0, 30000.0, 636000.0, 12000.0308    w1 = (1 / s1 ** 2) / (1 / s1 ** 2 + 1 / s2 ** 2)309    out["ex1"] = {"w1": w1, "w2": 1 - w1, "V": w1 * V1 + (1 - w1) * V2, "sd": math.sqrt(1 / (1 / s1 ** 2 + 1 / s2 ** 2)),310                  "V_egal": 0.5 * (V1 + V2), "sd_egal": math.sqrt(0.25 * s1 ** 2 + 0.25 * s2 ** 2)}311    # ex. 2 : intervalle de confiance de la valeur (E[P|x]) contre intervalle de prédiction du prix (ch. 18)312    lnp, var_est, sig2 = 13.3627, 0.00010, 0.0830 ** 2313    duan = 1.00335314    out["ex2"] = {"ic_bas": math.exp(lnp - 1.96 * math.sqrt(var_est)) * duan, "ic_haut": math.exp(lnp + 1.96 * math.sqrt(var_est)) * duan,315                  "ip_bas": math.exp(lnp - 1.96 * math.sqrt(var_est + sig2)), "ip_haut": math.exp(lnp + 1.96 * math.sqrt(var_est + sig2)),316                  "demi_largeur_ic_pct": 1.96 * math.sqrt(var_est), "demi_largeur_ip_pct": 1.96 * math.sqrt(var_est + sig2),317                  "rapport_largeurs": math.sqrt(var_est + sig2) / math.sqrt(var_est)}318    # ex. 3 : poids optimaux avec corrélation319    out["ex3"] = {r: {"w1": poids_optimaux_deux(s1, s2, r),320                      "sd": math.sqrt(max(0.0, (lambda w: w ** 2 * s1 ** 2 + (1 - w) ** 2 * s2 ** 2 + 2 * w * (1 - w) * r * s1 * s2)(poids_optimaux_deux(s1, s2, r))))}321                  for r in (0.0, 0.3, 0.5, 0.6)}322    # ex. 4 : bootstrap à la main (percentiles de 20 réplications données)323    rep = np.array([561, 568, 571, 574, 576, 577, 579, 580, 582, 583, 585, 586, 588, 590, 592, 595, 598, 603, 611, 624]) * 1000.0324    out["ex4"] = {"moyenne": float(rep.mean()), "ecart_type": float(rep.std(ddof=1)), "q05": float(np.percentile(rep, 5)), "q95": float(np.percentile(rep, 95)),325                  "q10": float(np.percentile(rep, 10)), "q90": float(np.percentile(rep, 90))}326    # ex. 5 : coût sous hypothèses d'usage, puis extraction avec le terrain du ch. 24327    CN = co["CN"]328    out["ex5"] = {"a": ct.valeur_par_cout(0.35 * b["V_comp"], CN, ct.depreciation_age_vie(CN, 22, 60))["valeur"],329                  "b": ct.valeur_par_cout(0.35 * b["V_comp"], CN, ct.depreciation_age_vie(CN, 28, 75))["valeur"],330                  "extraction_hed": ct.extraction_depreciation_marche(b["V_hed"], co["V_L"], CN, 28),331                  "extraction_comp": ct.extraction_depreciation_marche(b["V_comp"], co["V_L"], CN, 28),332                  "V_L_pour_hed": b["V_hed"] - (CN - co["D_tot"])}333    # ex. 6 : autres pondérations et comparaison élargie (k = 16)334    V16 = float(r_el["valeur_ponderee"].iloc[-1]) if (r_el := _CACHE.get("el")) is not None else float(elargissement(b, bs["P_aj16"], bs["brut16"])["valeur_ponderee"].iloc[-1])335    w2 = np.array([0.30, 0.55, 0.15]); w1 = np.array([POIDS_RAISONNES["comparaison"], POIDS_RAISONNES["hedonique"], POIDS_RAISONNES["cout"]])336    V = np.array([b["V_comp"], b["V_hed"], co["valeur"]]); V16v = np.array([V16, b["V_hed"], co["valeur"]])337    a = ARRONDI_CONCLUSION338    bs8 = bootstrap(b, B=B, k=8)339    out["ex6"] = {"sd_conjoint_k8": bs8["resume"]["conjoint"]["ecart_type"], "V_k8": float(valeur_ponderee(bs["P_aj16"][:8], bs["brut16"][:8])),340                  "V_30_55_15": float(w2 @ V), "arrondi_30_55_15": round(float(w2 @ V) / a) * a, "V16": V16,341                  "V_k16_raisonnee": float(w1 @ V16v), "arrondi_k16": round(float(w1 @ V16v) / a) * a,342                  "V_55_30_15": float(np.array([0.55, 0.30, 0.15]) @ V), "arrondi_55_30_15": round(float(np.array([0.55, 0.30, 0.15]) @ V) / a) * a}343    # ex. 7 : incidence d'une hypothèse extraordinaire (eq. 10.incidence / 10.esperance)344    out["ex7"] = {"V_H": cb["conclusion"], "V_nonH": cb["conclusion"] - 14000.0, "p": 0.7}345    out["ex7"]["delta"] = out["ex7"]["V_H"] - out["ex7"]["V_nonH"]346    out["ex7"]["esperance"] = out["ex7"]["V_nonH"] + out["ex7"]["p"] * out["ex7"]["delta"]347    return out348349350# ---------------------------------------------------------------- tout351_CACHE: dict = {}352353def tout(verbose: bool = False) -> dict:354    if "r" in _CACHE:355        return _CACHE["r"]356    b = base()357    bs = bootstrap(b)358    ex = exhaustif(b, bs["P_aj16"], bs["brut16"])359    el = elargissement(b, bs["P_aj16"], bs["brut16"]); _CACHE["el"] = el360    co = cout(b)361    bk = bootstrap_cout(b, bs, co)362    cb = combinaison(b, bs, co, bk)363    vt = verite(b)364    exo = exercices(b, bs, co, cb)365    r = {"base": b, "bootstrap": bs, "exhaustif": ex, "elargissement": el, "cout": co, "bootstrap_cout": bk, "combinaison": cb, "verite": vt, "exercices": exo}366    _CACHE["r"] = r367    if verbose:368        imprimer(r)369    return r370371372def imprimer(r: dict) -> None:373    b, bs, ex, co, cb, vt = r["base"], r["bootstrap"], r["exhaustif"], r["cout"], r["combinaison"], r["verite"]374    f = lambda x: f"{x:,.0f}".replace(",", " ")375    print("=== Indications héritées ===")376    print(f"Comparaison (ch. 17) : pondérée {f(b['V_comp'])} $, écart-type pondéré {f(b['sd_comp_pondere'])} $, poids {np.round(b['poids_comp'], 3)}")377    ip = b["ip_hed"]378    print(f"Hédonique (ch. 18) : {f(b['V_hed'])} $ (Duan) ; IP 95 % [{f(ip['bas'])} ; {f(ip['haut'])}]")379    print("\n=== 16 candidats : prix ajusté, brut, résidu (%) ===")380    a = b["adj16"]381    print(a[["id_vente", "distance", "prix_vente", "prix_ajuste", "ajust_brut_pct", "residu_log", "retenu"]].assign(residu_pct=lambda d: np.exp(d["residu_log"]) - 1).round(4).to_string(index=False))382    print(f"résidu moyen (exp−1) des 5 retenus : {np.mean(np.exp(a.loc[a['retenu'], 'residu_log']) - 1):.4f} ; des 11 autres : {np.mean(np.exp(a.loc[~a['retenu'], 'residu_log']) - 1):.4f} ; des 16 : {np.mean(np.exp(a['residu_log']) - 1):.4f}")383    print("\n=== Bootstrap (B = %d) ===" % B)384    for k, v in bs["resume"].items():385        print(f"{k:14s} moy {f(v['moyenne'])} é.-t. {f(v['ecart_type'])} q2.5 {f(v['q025'])} q5 {f(v['q05'])} q50 {f(v['q50'])} q95 {f(v['q95'])} q97.5 {f(v['q975'])}")386    print(f"corrélation (conjoint, prédiction) = {bs['rho_conjoint_pred']:.3f} ; (ajustements, prédiction) = {bs['rho_ajustements_pred']:.3f}")387    print(f"part des rééchantillons de candidats sous la prédiction hédonique : {bs['part_A_sous_hed']:.3f} ; conjoint : {bs['part_C_sous_hed']:.3f}")388    print("\n=== Énumération exhaustive C(16,5) ===")389    print(f"n = {ex['n']} ; moyenne {f(ex['moyenne'])} ; é.-t. {f(ex['ecart_type'])} ; min {f(ex['min'])} ; max {f(ex['max'])} ; q5 {f(ex['q05'])} ; q50 {f(ex['q50'])} ; q95 {f(ex['q95'])}")390    print(f"rang de la sélection du ch. 16 : {ex['rang_selection']} sous-ensembles sur {ex['n']} donnent moins ; part sous la prédiction hédonique {ex['part_sous_hed']:.3f}")391    print(ex["par_n_retenus"].round(0).to_string())392    print("\n=== Élargissement k = 5 … 16 (modèle fixe) ===")393    print(r["elargissement"].round(4).to_string(index=False))394    bk = r["bootstrap_cout"]395    print("\n=== Coût (chapitre 24) ===")396    print(f"C_N {f(co['CN'])} ; V_L {f(co['V_L'])} (médiane {f(co['V_L_mediane'])} de {co['n_terrain']} ventes récentes de Hull, part {co['part_terrain']:.3f}) ; D extraite {f(co['D_ext'])} ({co['taux_ext']:.4f}) ; D fonctionnelle {f(co['D_fonc'])} ; D totale {f(co['D_tot'])} ({co['taux_tot']:.4f}) ; contribution bâtiment {f(co['contribution_batiment'])} ; valeur {f(co['valeur'])} → {f(co['arrondie'])} ; lecture ventilée {f(co['ventilee'])} ({co['taux_ventilee']:.4f})")397    print(f"C_N + V_L = {f(co['cout_neuf_plus_terrain'])} ; extraction à la prédiction hédonique : D = {f(co['extraction_hed']['depreciation'])} (taux {co['extraction_hed']['taux']:.4f}) ; à la comparaison : D = {f(co['extraction_comp']['depreciation'])} (taux {co['extraction_comp']['taux']:.4f})")398    u = co["usage"]399    print(f"hypothèses d'usage (allocation {u['alpha']:.2f} × comparaison = {f(u['V_L'])}, âge-vie {u['vie']} ans : D = {f(u['D'])}) : valeur {f(u['valeur'])} ; extraction au modèle avec ce terrain : D = {f(u['extraction_hed']['depreciation'])} ({u['extraction_hed']['taux']:.3f})")400    print("variantes d'extraction (ch. 24) :"); print(co["sens_extraction"][["variante", "vl_hull", "taux_28", "valeur"]].round(4).to_string(index=False)); print(f"écart-type de spécification {f(co['sd_spec'])}")401    print("tornade ch. 24 :"); [print(f"  {t['param']:50s} {f(t['bas'])} – {f(t['haut'])}") for t in co["tornade"]]402    rk = bk["resume"]403    print(f"bootstrap du coût (B = {len(bk['VK'])}) : moy {f(rk['moyenne'])} é.-t. {f(rk['ecart_type'])} q5 {f(rk['q05'])} q50 {f(rk['q50'])} q95 {f(rk['q95'])} ; é.-t. V_L {f(bk['sd_VL'])} ; é.-t. taux 28 ans {bk['sd_T28']:.4f} ; ρ(coût, prédiction) {bk['rho_cout_pred']:.3f} ; ρ(coût, comparaison) {bk['rho_cout_comp']:.3f}")404    print("\n=== Combinaison ===")405    print(f"indications : {[f(x) for x in cb['V']]} ; écarts-types : {[f(x) for x in cb['s']]} (hédonique : boot {f(cb['sd_pred_boot'])} + spécification {f(cb['sd_spec_hed'])} ; coût : boot {f(bk['resume']['ecart_type'])} + spécification {f(co['sd_spec'])})")406    print("matrice de corrélation :"); print(np.round(cb["R"], 3))407    print(f"moyenne simple avec le coût sous hypothèses d'usage : {f(cb['moyenne_simple_usage'])}")408    print(cb["tableau"].round(3).to_string(index=False))409    print(f"poids optimal de la comparaison à deux indications selon ρ : { {round(k,3): round(v,3) for k, v in cb['w_opt_deux'].items()} }")410    q = cb["q"]411    print(f"conclusion raisonnée {f(cb['V_raisonnee'])} → arrondie {f(cb['conclusion'])} ; bootstrap de la combinaison : q2.5 {f(q['q025'])} q5 {f(q['q05'])} q50 {f(q['q50'])} q95 {f(q['q95'])} q97.5 {f(q['q975'])} ; fourchette 90 % arrondie {f(cb['fourchette90'][0])}–{f(cb['fourchette90'][1])} ; fourchette conclue {f(cb['fourchette_conclue'][0])}–{f(cb['fourchette_conclue'][1])} ({cb['largeur_relative_conclue']:.3f})")412    print(f"écart conclusion / comparaison {cb['ecart_conclusion_comp']:.4f} ; / hédonique {cb['ecart_conclusion_hed']:.4f}")413    print("\n=== Loi de génération ===")414    print(f"c_hat {vt['c_hat']:.4f} (e.-t. {vt['se_c']:.4f}, σ_e {vt['sigma_e']:.4f}) ; V_vraie {f(vt['V_true'])} [{f(vt['V_true_bas'])} ; {f(vt['V_true_haut'])}] ; écart comparaison {vt['ecart_comp']:.4f} ; hédonique {vt['ecart_hed']:.4f} ; conclusion {cb['conclusion']/vt['V_true']-1:.4f}")415    print("\n=== Exercices ===")416    for k, v in r["exercices"].items():417        print(k, {kk: (round(vv, 4) if isinstance(vv, float) else vv) for kk, vv in v.items()} if isinstance(v, dict) else v)418419420if __name__ == "__main__":421    tout(verbose=True)422