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# Hypothèses déclarées de l'indication par le coût (à aligner avec les ch. 24 et 25 ; voir verification/chapters/ch20.md)30ALPHA_TERRAIN = 0.30        # part du terrain dans la valeur (méthode d'allocation) — hypothèse pédagogique, NON_VERIFIE31VIE_ECONOMIQUE = 60         # vie économique totale (ans), méthode âge-vie — hypothèse pédagogique, NON_VERIFIE32AGE_EFFECTIF = 28           # âge effectif = âge chronologique (entretien normal présumé, aucune inspection dans le dossier)33GRILLE_ALPHA = [0.20, 0.25, 0.30, 0.35, 0.40]34GRILLE_VIE = [50, 60, 75]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.55, "cout": 0.05}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) -> 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 = 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, "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 (hypothèses déclarées)162def cout(b: dict, alpha: float = ALPHA_TERRAIN, vie: float = VIE_ECONOMIQUE, age_eff: float = AGE_EFFECTIF) -> dict:163    devis, sujet = c21.charger()164    cn = c21.cout_neuf_sujet(devis)165    CN = cn["total"]166    V_L = alpha * b["V_comp"]                                   # allocation : part α de l'indication par comparaison167    D = ct.depreciation_age_vie(CN, age_eff, vie)168    val = ct.valeur_par_cout(V_L, CN, D)169    extraction_hed = ct.extraction_depreciation_marche(b["V_hed"], V_L, CN, age_eff)170    extraction_comp = ct.extraction_depreciation_marche(b["V_comp"], V_L, CN, age_eff)171    grille = pd.DataFrame([{"alpha": a, "vie": vv, "V_L": a * b["V_comp"], "D": ct.depreciation_age_vie(CN, age_eff, vv),172                            "valeur": ct.valeur_par_cout(a * b["V_comp"], CN, ct.depreciation_age_vie(CN, age_eff, vv))["valeur"]}173                           for a in GRILLE_ALPHA for vv in GRILLE_VIE])174    # âge effectif réduit (entretien supérieur) : variante175    age_var = {ae: ct.valeur_par_cout(V_L, CN, ct.depreciation_age_vie(CN, ae, vie))["valeur"] for ae in (20, 24, 28)}176    # terrain nécessaire pour que le coût rejoigne la comparaison ou le modèle (D âge-vie inchangée)177    return {"CN": CN, "V_L": V_L, "D": D, "taux_D": D / CN, "contribution_batiment": val["contribution_batiment"], "valeur": val["valeur"],178            "extraction_hed": extraction_hed, "extraction_comp": extraction_comp, "grille": grille,179            "sd_grille": float(grille["valeur"].std(ddof=0)), "min_grille": float(grille["valeur"].min()), "max_grille": float(grille["valeur"].max()),180            "age_effectif_variantes": age_var, "V_L_pour_comp": b["V_comp"] - (CN - D), "V_L_pour_hed": b["V_hed"] - (CN - D),181            "cout_neuf_plus_terrain": CN + V_L, "alpha": alpha, "vie": vie, "age_eff": age_eff}182183184# ---------------------------------------------------------------- 4. combinaison185def poids_optimaux_deux(s1: float, s2: float, rho: float) -> float:186    """Poids optimal de l'indication 1 parmi deux indications corrélées (Bates–Granger) : w1 = (σ2² − ρσ1σ2)/(σ1² + σ2² − 2ρσ1σ2)."""187    return (s2 ** 2 - rho * s1 * s2) / (s1 ** 2 + s2 ** 2 - 2 * rho * s1 * s2)188189190def poids_gls(Sigma: np.ndarray) -> np.ndarray:191    """w = Σ⁻¹1 / 1'Σ⁻¹1 (variance minimale sous Σw = 1, poids non contraints)."""192    inv = np.linalg.inv(Sigma); un = np.ones(Sigma.shape[0])193    w = inv @ un194    return w / w.sum()195196197def combinaison(b: dict, bs: dict, co: dict, graine: int = GRAINE) -> dict:198    V = np.array([b["V_comp"], b["V_hed"], co["valeur"]])199    sd_spec = float(np.std(PRED_SPEC_CH18, ddof=1))200    sd_pred_boot = bs["resume"]["prediction"]["ecart_type"]201    s = np.array([bs["resume"]["conjoint"]["ecart_type"], math.sqrt(sd_pred_boot ** 2 + sd_spec ** 2), co["sd_grille"]])202    rho = bs["rho_conjoint_pred"]203    Sigma = np.diag(s ** 2); Sigma[0, 1] = Sigma[1, 0] = rho * s[0] * s[1]204    schemas = {205        "comparaison seule": np.array([1.0, 0.0, 0.0]),206        "hédonique seule": np.array([0.0, 1.0, 0.0]),207        "égale (un tiers chacune)": np.ones(3) / 3,208        "inverse de la variance (indépendance)": (1 / s ** 2) / (1 / s ** 2).sum(),209        "variance minimale avec corrélation (MCG)": poids_gls(Sigma),210        "raisonnée (retenue)": np.array([POIDS_RAISONNES["comparaison"], POIDS_RAISONNES["hedonique"], POIDS_RAISONNES["cout"]]),211    }212    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))}213                        for k, w in schemas.items()])214    # distribution bootstrap de la conclusion retenue (réplication conjointe : comparaison et modèle issus du même tirage ; coût tiré de la grille)215    rng = np.random.default_rng(graine + 1)216    cout_b = rng.choice(co["grille"]["valeur"].values, size=len(bs["VC"]))217    w = schemas["raisonnée (retenue)"]218    Vb = w[0] * bs["VC"] + w[1] * bs["pred"] + w[2] * cout_b219    q = np.percentile(Vb, [2.5, 5, 50, 95, 97.5])220    V_raisonnee = float(w @ V)221    arrondi = ARRONDI_CONCLUSION222    conclusion = round(V_raisonnee / arrondi) * arrondi223    fourchette90 = (math.floor(q[1] / arrondi) * arrondi, math.ceil(q[3] / arrondi) * arrondi)224    # fourchette de valeur conclue : encadrée par les deux indications crédibles (comparaison, hédonique), arrondies vers l'extérieur225    fourchette_conclue = (math.floor(min(V[0], V[1]) / arrondi) * arrondi, math.ceil(max(V[0], V[1]) / arrondi) * arrondi)226    # courbes de sensibilité aux poids : w_hed de 0 à 1 − w_cout, pour trois poids du coût227    wh = np.linspace(0, 1, 101)228    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.05, 0.10)}229    # variance de la combinaison de deux indications (comparaison, hédonique) selon ρ230    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)}231    w_opt_deux = {r: poids_optimaux_deux(s[0], s[1], r) for r in (0.0, rho, 0.3, 0.6, 0.9)}232    return {"V": V, "s": s, "sd_spec_hed": sd_spec, "sd_pred_boot": sd_pred_boot, "rho": rho, "Sigma": Sigma, "schemas": schemas, "tableau": tab,233            "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,234            "conclusion": conclusion, "fourchette90": fourchette90, "fourchette_conclue": fourchette_conclue,235            "largeur_relative_conclue": (fourchette_conclue[1] - fourchette_conclue[0]) / conclusion, "wh": wh, "courbes": courbes, "var_deux": var_deux, "w_opt_deux": w_opt_deux,236            "ecart_conclusion_comp": conclusion / b["V_comp"] - 1, "ecart_conclusion_hed": conclusion / b["V_hed"] - 1}237238239# ---------------------------------------------------------------- 5. ce que révèle la loi de génération240def verite(b: dict) -> dict:241    """Espérance conditionnelle « vraie » du prix du sujet selon la loi déclarée dans data/synthetic/README.md ;242    la constante (normalisation de la médiane) est estimée sur les 600 ventes (erreur-type σ/√n ≈ 0,35 %)."""243    from traite.synth import _mois_depuis244    v = b["ventes"]; s = b["sujet"]245    med = 475000.0; med_sect = {"Gatineau": 464900.0, "Aylmer": 530000.0, "Buckingham": 395000.0, "Masson-Angers": 395000.0}246    prime = {"Hull": 0.0, **{k: math.log(m / med) for k, m in med_sect.items()}}; prime["Masson-Angers"] -= 0.02247    g = 1.0906 ** 0.5 - 1248    ref = _mois_depuis([datetime.date(2024, 7, 1)])[0]249    def f(df, dates):250        mois = _mois_depuis([d.date() for d in pd.to_datetime(dates)])251        return (0.55 * np.log(df["superficie_habitable_m2"].values / 140) + 0.10 * np.log(df["superficie_terrain_m2"].values / 560)252                - 0.0055 * df["age"].values + 0.00003 * df["age"].values ** 2 + 0.06 * (df["qualite_1a5"].values - 3)253                + 0.035 * df["garage_places"].values + 0.03 * df["sous_sol_fini"].values + 0.02 * df["piscine"].values254                + np.array([prime[x] for x in df["secteur"].values]) - 0.006 * np.clip(df["dist_centre_km"].values - 3, 0, None)255                + 0.02 * (df["type"].values == "à étages") + math.log1p(g) * (mois - ref) / 12)256    fv = f(v, v["date_vente"]); e = np.log(v["prix_vente"].values) - fv257    c_hat = float(e.mean()); se_c = float(e.std(ddof=1) / math.sqrt(len(e)))258    sd = cas1.sujet_df(s); fs = float(f(sd, sd["date_vente"])[0])259    V_true = math.exp(fs + c_hat) * math.exp(0.085 ** 2 / 2)260    return {"c_hat": c_hat, "se_c": se_c, "sigma_e": float(e.std(ddof=1)), "V_true": V_true,261            "V_true_bas": V_true * math.exp(-1.96 * se_c), "V_true_haut": V_true * math.exp(1.96 * se_c),262            "ecart_comp": b["V_comp"] / V_true - 1, "ecart_hed": b["V_hed"] / V_true - 1}263264265# ---------------------------------------------------------------- 6. exercices (chiffres des solutions)266def exercices(b: dict, bs: dict, co: dict, cb: dict) -> dict:267    out = {}268    # ex. 1 : deux indications indépendantes269    V1, s1, V2, s2 = 582000.0, 30000.0, 636000.0, 12000.0270    w1 = (1 / s1 ** 2) / (1 / s1 ** 2 + 1 / s2 ** 2)271    out["ex1"] = {"w1": w1, "w2": 1 - w1, "V": w1 * V1 + (1 - w1) * V2, "sd": math.sqrt(1 / (1 / s1 ** 2 + 1 / s2 ** 2)),272                  "V_egal": 0.5 * (V1 + V2), "sd_egal": math.sqrt(0.25 * s1 ** 2 + 0.25 * s2 ** 2)}273    # ex. 2 : intervalle de confiance de la valeur (E[P|x]) contre intervalle de prédiction du prix (ch. 18)274    lnp, var_est, sig2 = 13.3627, 0.00010, 0.0830 ** 2275    duan = 1.00335276    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,277                  "ip_bas": math.exp(lnp - 1.96 * math.sqrt(var_est + sig2)), "ip_haut": math.exp(lnp + 1.96 * math.sqrt(var_est + sig2)),278                  "demi_largeur_ic_pct": 1.96 * math.sqrt(var_est), "demi_largeur_ip_pct": 1.96 * math.sqrt(var_est + sig2),279                  "rapport_largeurs": math.sqrt(var_est + sig2) / math.sqrt(var_est)}280    # ex. 3 : poids optimaux avec corrélation281    out["ex3"] = {r: {"w1": poids_optimaux_deux(s1, s2, r),282                      "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))))}283                  for r in (0.0, 0.3, 0.5, 0.6)}284    # ex. 4 : bootstrap à la main (percentiles de 20 réplications données)285    rep = np.array([561, 568, 571, 574, 576, 577, 579, 580, 582, 583, 585, 586, 588, 590, 592, 595, 598, 603, 611, 624]) * 1000.0286    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)),287                  "q10": float(np.percentile(rep, 10)), "q90": float(np.percentile(rep, 90))}288    # ex. 5 : coût avec autres hypothèses289    CN = co["CN"]290    out["ex5"] = {"a": ct.valeur_par_cout(0.35 * b["V_comp"], CN, ct.depreciation_age_vie(CN, 22, 60))["valeur"],291                  "b": ct.valeur_par_cout(0.35 * b["V_comp"], CN, ct.depreciation_age_vie(CN, 28, 75))["valeur"],292                  "extraction": ct.extraction_depreciation_marche(b["V_hed"], 0.35 * b["V_comp"], CN, 28)}293    # ex. 6 : autres pondérations et comparaison élargie (k = 16)294    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])295    w2 = np.array([0.30, 0.65, 0.05]); w1 = np.array([POIDS_RAISONNES["comparaison"], POIDS_RAISONNES["hedonique"], POIDS_RAISONNES["cout"]])296    V = np.array([b["V_comp"], b["V_hed"], co["valeur"]]); V16v = np.array([V16, b["V_hed"], co["valeur"]])297    a = ARRONDI_CONCLUSION298    out["ex6"] = {"V_30_65_05": float(w2 @ V), "arrondi_30_65_05": round(float(w2 @ V) / a) * a, "V16": V16,299                  "V_k16_raisonnee": float(w1 @ V16v), "arrondi_k16": round(float(w1 @ V16v) / a) * a,300                  "V_60_35_05": float(np.array([0.60, 0.35, 0.05]) @ V)}301    # ex. 7 : incidence d'une hypothèse extraordinaire (eq. 10.incidence / 10.esperance)302    out["ex7"] = {"V_H": cb["conclusion"], "V_nonH": cb["conclusion"] - 14000.0, "p": 0.7}303    out["ex7"]["delta"] = out["ex7"]["V_H"] - out["ex7"]["V_nonH"]304    out["ex7"]["esperance"] = out["ex7"]["V_nonH"] + out["ex7"]["p"] * out["ex7"]["delta"]305    return out306307308# ---------------------------------------------------------------- tout309_CACHE: dict = {}310311def tout(verbose: bool = False) -> dict:312    if "r" in _CACHE:313        return _CACHE["r"]314    b = base()315    bs = bootstrap(b)316    ex = exhaustif(b, bs["P_aj16"], bs["brut16"])317    el = elargissement(b, bs["P_aj16"], bs["brut16"]); _CACHE["el"] = el318    co = cout(b)319    cb = combinaison(b, bs, co)320    vt = verite(b)321    exo = exercices(b, bs, co, cb)322    r = {"base": b, "bootstrap": bs, "exhaustif": ex, "elargissement": el, "cout": co, "combinaison": cb, "verite": vt, "exercices": exo}323    _CACHE["r"] = r324    if verbose:325        imprimer(r)326    return r327328329def imprimer(r: dict) -> None:330    b, bs, ex, co, cb, vt = r["base"], r["bootstrap"], r["exhaustif"], r["cout"], r["combinaison"], r["verite"]331    f = lambda x: f"{x:,.0f}".replace(",", " ")332    print("=== Indications héritées ===")333    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)}")334    ip = b["ip_hed"]335    print(f"Hédonique (ch. 18) : {f(b['V_hed'])} $ (Duan) ; IP 95 % [{f(ip['bas'])} ; {f(ip['haut'])}]")336    print("\n=== 16 candidats : prix ajusté, brut, résidu (%) ===")337    a = b["adj16"]338    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))339    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}")340    print("\n=== Bootstrap (B = %d) ===" % B)341    for k, v in bs["resume"].items():342        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'])}")343    print(f"corrélation (conjoint, prédiction) = {bs['rho_conjoint_pred']:.3f} ; (ajustements, prédiction) = {bs['rho_ajustements_pred']:.3f}")344    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}")345    print("\n=== Énumération exhaustive C(16,5) ===")346    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'])}")347    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}")348    print(ex["par_n_retenus"].round(0).to_string())349    print("\n=== Élargissement k = 5 … 16 (modèle fixe) ===")350    print(r["elargissement"].round(4).to_string(index=False))351    print("\n=== Coût (hypothèses : α = %.2f, vie = %d ans, âge effectif = %d) ===" % (co["alpha"], co["vie"], co["age_eff"]))352    print(f"C_N {f(co['CN'])} ; V_L {f(co['V_L'])} ; D âge-vie {f(co['D'])} ({co['taux_D']:.3f}) ; contribution bâtiment {f(co['contribution_batiment'])} ; valeur {f(co['valeur'])}")353    print(f"C_N + V_L = {f(co['cout_neuf_plus_terrain'])} ; extraction du marché à la prédiction hédonique : D = {f(co['extraction_hed']['depreciation'])} (taux {co['extraction_hed']['taux']:.3f}) ; à la comparaison : D = {f(co['extraction_comp']['depreciation'])} (taux {co['extraction_comp']['taux']:.3f})")354    print(f"terrain implicite pour rejoindre la comparaison : {f(co['V_L_pour_comp'])} ; le modèle : {f(co['V_L_pour_hed'])}")355    print(f"grille α × vie : min {f(co['min_grille'])}, max {f(co['max_grille'])}, é.-t. {f(co['sd_grille'])} ; âge effectif 20/24/28 : {[f(x) for x in co['age_effectif_variantes'].values()]}")356    print(co["grille"].pivot(index="alpha", columns="vie", values="valeur").round(0).to_string())357    print("\n=== Combinaison ===")358    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'])}) ; ρ = {cb['rho']:.3f}")359    print(cb["tableau"].round(3).to_string(index=False))360    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()} }")361    q = cb["q"]362    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})")363    print(f"écart conclusion / comparaison {cb['ecart_conclusion_comp']:.4f} ; / hédonique {cb['ecart_conclusion_hed']:.4f}")364    print("\n=== Loi de génération ===")365    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}")366    print("\n=== Exercices ===")367    for k, v in r["exercices"].items():368        print(k, {kk: (round(vv, 4) if isinstance(vv, float) else vv) for kk, vv in v.items()} if isinstance(v, dict) else v)369370371if __name__ == "__main__":372    tout(verbose=True)373