spb/traite-evaluation-qc

1#!/usr/bin/env python32"""Chapitre 1 — Théorie de la valeur et prix immobiliers : tous les chiffres cités dans le texte (DONNÉES SYNTHÉTIQUES, cas 1, graine 20261031).34Sections du script :5  1. modèle hédonique de référence du cas 1 (traite.cas1.hedonique) : dispersion des prix à caractéristiques données (sigma),6     prix implicites, contributions, rendements décroissants, anticipation (tendance), conformité (secteurs) ;7  2. valeur = espérance conditionnelle : prédiction pour le sujet, médiane/moyenne/mode log-normaux, facteur de Duan,8     intervalle de prédiction, probabilité d'un écart prix/valeur supérieur à un seuil ;9  3. prix des comparables du chapitre 16 (prix observés vs valeur estimée) ;10  4. modèle de recherche du vendeur (simulation étiquetée) : prix de réserve, délai de vente, dispersion des prix de biens identiques ;11  5. plus haute et meilleure utilisation : valeur résiduelle du terrain selon trois usages (hypothèses déclarées) ;12  6. sensibilité à la dispersion sigma.13Exécuter depuis la racine : python3 analysis/python/cases/ch01_valeur_prix.py14"""15from __future__ import annotations16import sys, pathlib, math17import numpy as np, pandas as pd1819ROOT = pathlib.Path(__file__).resolve().parents[3]20sys.path.insert(0, str(ROOT / "analysis/python"))21from traite import cas122from scipy.stats import norm2324GRAINE = 202610312526# ---------------------------------------------------------------------------- 1. modèle hédonique et principes2728LIBELLES = {29    "ln_S": "Superficie habitable (élasticité, réf. 140 m²)", "ln_T": "Superficie du terrain (élasticité, réf. 560 m²)",30    "age": "Âge (par année)", "age2_100": "Âge au carré (/100)", "qualite": "Qualité (par échelon, réf. 3)",31    "garage": "Garage (par place)", "sous_sol": "Sous-sol fini", "piscine": "Piscine", "etages": "Bâtiment à étages",32    "dist_centre": "Distance au centre au-delà de 3 km (par km)", "sect_Gatineau": "Secteur Gatineau (réf. Hull)",33    "sect_Aylmer": "Secteur Aylmer (réf. Hull)", "sect_Buckingham": "Secteur Buckingham (réf. Hull)",34    "sect_Masson-Angers": "Secteur Masson-Angers (réf. Hull)", "mois": "Tendance du marché (par mois)",35}3637def modele() -> dict:38    d = cas1.charger(); m = cas1.hedonique(d["ventes"])39    return {**d, "modele": m}4041def prediction_sujet(r: dict) -> dict:42    """Valeur du sujet comme espérance conditionnelle : médiane (exp(x'b)), moyenne (Duan), intervalle de prédiction à 95 %."""43    m, s = r["modele"], r["sujet"]44    sd = cas1.sujet_df(s)45    X = cas1.matrice(sd)[m["X"].columns]46    lnp = float((X.values @ m["beta"])[0])47    duan = float(np.mean(np.exp(m["residus"])))48    sigma = m["sigma"]49    ip = cas1.intervalle_prediction(m, sd).iloc[0]50    return {"ln_pred": lnp, "mediane": math.exp(lnp), "moyenne_duan": math.exp(lnp) * duan, "facteur_duan": duan,51            "moyenne_lognormale": math.exp(lnp + sigma ** 2 / 2), "mode_lognormal": math.exp(lnp - sigma ** 2),52            "ip_bas": float(ip["bas"]), "ip_haut": float(ip["haut"]), "sigma": sigma,53            "se_pred_log": float(math.sqrt(sigma ** 2 + (X.values @ m["V_hc1"] @ X.values.T)[0, 0])),54            "x_sujet": X.iloc[0]}5556def contributions(r: dict, pred: dict) -> pd.DataFrame:57    """Prix implicites et contributions marginales des attributs du sujet : C_k = P̂ [1 − exp(−β_k x_k)] (retrait de l'attribut)."""58    m = r["modele"]; beta = pd.Series(m["beta"], index=m["X"].columns); x = pred["x_sujet"]; P = pred["mediane"]59    rows = []60    for k in beta.index:61        if k == "const": continue62        bx = float(beta[k] * x[k])63        rows.append({"variable": k, "libelle": LIBELLES.get(k, k), "x_sujet": float(x[k]), "beta": float(beta[k]),64                     "se_hc1": float(m["coef"].loc[k, "se_hc1"]), "effet_pct_unitaire": math.exp(beta[k]) - 1,65                     "beta_x": bx, "effet_pct_sujet": math.exp(bx) - 1, "contribution_dollars": P * (1 - math.exp(-bx))})66    return pd.DataFrame(rows).set_index("variable")6768def rendements_decroissants(r: dict, pred: dict) -> dict:69    """Prix moyen et prix marginal du mètre carré habitable pour le sujet : dP/dS = β_S P / S."""70    m, s = r["modele"], r["sujet"]; bS = float(m["coef"].loc["ln_S", "coef"]); P = pred["mediane"]; S = s["superficie_habitable_m2"]71    out = {"elasticite_S": bS, "prix_moyen_m2": P / S, "prix_marginal_m2": bS * P / S}72    # valeur ajoutée par 10 m² supplémentaires, puis par 10 m² de plus (décroissance)73    out["gain_10m2"] = P * ((1 + 10 / S) ** bS - 1)74    out["gain_10m2_suivants"] = P * ((1 + 20 / S) ** bS - (1 + 10 / S) ** bS)75    return out7677def anticipation(r: dict, pred: dict) -> dict:78    """Tendance temporelle : croissance mensuelle, annuelle composée ; prix du sujet à la date de référence du rôle vs date d'évaluation."""79    m, s = r["modele"], r["sujet"]; g_m = m["croissance_mensuelle"]80    mois_ref = float(pred["x_sujet"]["mois"])  # mois entre 2024-07-01 et 2026-09-0181    return {"croissance_mensuelle": g_m, "croissance_annuelle": (1 + g_m) ** 12 - 1, "mois_depuis_reference": mois_ref,82            "facteur_cumule": (1 + g_m) ** mois_ref, "valeur_a_la_reference": pred["mediane"] / (1 + g_m) ** mois_ref,83            "valeur_role": s["valeur_role_2025"]}8485def conformite(r: dict, pred: dict) -> pd.DataFrame:86    """Valeur du sujet (médiane conditionnelle) s'il était situé dans chacun des cinq secteurs, toutes choses égales d'ailleurs."""87    m = r["modele"]; beta = pd.Series(m["beta"], index=m["X"].columns); P = pred["mediane"]88    rows = [{"secteur": "Hull", "effet_pct": 0.0, "valeur": P}]89    for sct in cas1.SECTEURS[1:]:90        b = float(beta[f"sect_{sct}"]); rows.append({"secteur": sct, "effet_pct": math.exp(b) - 1, "valeur": P * math.exp(b)})91    return pd.DataFrame(rows).set_index("secteur")9293# ---------------------------------------------------------------------------- 2. dispersion et espérance conditionnelle9495def dispersion(r: dict) -> dict:96    res = r["modele"]["residus"]; sigma = r["modele"]["sigma"]97    s = pd.Series(res)98    out = {"n": len(res), "sigma": sigma, "ecart_type_residus": float(s.std(ddof=1)), "asymetrie": float(s.skew()), "aplatissement_exces": float(s.kurt()),99           "q025": float(np.quantile(res, 0.025)), "q975": float(np.quantile(res, 0.975)),100           "part_dans_pm10pct": float(np.mean(np.abs(res) <= math.log(1.10))),101           "part_dans_pm5pct": float(np.mean(np.abs(res) <= math.log(1.05)))}102    out["r2"] = r["modele"]["r2"]; out["r2_ajuste"] = r["modele"]["r2_ajuste"]103    return out104105def proba_ecart(sigma: float, d: float) -> float:106    """P(|P/V − 1| > d) sous ln(P/V) ~ N(0, σ²) : 2[1 − Φ(ln(1+d)/σ)] (approximation symétrique en log)."""107    return 2 * (1 - norm.cdf(math.log(1 + d) / sigma))108109def tableau_sensibilite(sigmas=(0.04, 0.06, 0.083, 0.12, 0.16)) -> pd.DataFrame:110    rows = []111    for sg in sigmas:112        rows.append({"sigma": sg, "moyenne_sur_mediane": math.exp(sg ** 2 / 2) - 1, "mode_sur_mediane": math.exp(-sg ** 2) - 1,113                     "ip95_bas_pct": math.exp(-1.96 * sg) - 1, "ip95_haut_pct": math.exp(1.96 * sg) - 1,114                     "proba_ecart_10pct": proba_ecart(sg, 0.10), "proba_ecart_5pct": proba_ecart(sg, 0.05),115                     "n_pour_se_2pct": math.ceil((1.96 * sg / math.log(1.02)) ** 2)})116    return pd.DataFrame(rows).set_index("sigma")117118# ---------------------------------------------------------------------------- 3. comparables du chapitre 16 : prix observés vs valeur estimée119120def comparables_ch16(r: dict) -> pd.DataFrame:121    sel = cas1.selection(); comps = sel["comparables"].copy()122    m = r["modele"]123    idx = comps.index124    comps["valeur_estimee"] = np.exp(m["X"].loc[idx].values @ m["beta"])125    comps["residu_log"] = m["residus"][[list(r["ventes"].index).index(i) for i in idx]]126    comps["ecart_pct"] = comps["prix_vente"] / comps["valeur_estimee"] - 1127    return comps[["id_vente", "date_vente", "superficie_habitable_m2", "qualite_1a5", "garage_places", "prix_vente", "valeur_estimee", "ecart_pct"]]128129# ---------------------------------------------------------------------------- 4. recherche du vendeur (simulation étiquetée)130131def prix_reserve(centre: float, sigma_o: float, c_mois: float, taux_annuel: float, p_offre: float = 1.0, n_quad: int = 4001) -> float:132    """Point fixe W = −c + δ [ p E max(O, W) + (1 − p) W ], offres O ~ log-normale(ln centre, σ_o²) ; prix de réserve r = W.133    centre = médiane des enchères des acheteurs (≠ valeur marchande, qui est l'espérance du prix ACCEPTÉ) ;134    c_mois = coût de détention mensuel (proportion de centre) ; δ = 1/(1 + i/12)."""135    delta = 1 / (1 + taux_annuel / 12); c = c_mois * centre136    z = np.linspace(-6, 6, n_quad); w = norm.pdf(z); w /= w.sum()137    O = centre * np.exp(sigma_o * z)138    W = centre * 0.9139    for _ in range(10_000):140        EmaxO = float(np.sum(w * np.maximum(O, W)))141        Wn = -c + delta * (p_offre * EmaxO + (1 - p_offre) * W)142        if abs(Wn - W) < 1e-6: W = Wn; break143        W = Wn144    return float(W)145146def simuler_ventes(centre: float, sigma_o: float, c_mois: float, taux_annuel: float, p_offre: float = 1.0, n: int = 20_000, graine: int = GRAINE) -> dict:147    """n vendeurs de biens identiques de valeur centre suivent la règle du prix de réserve ; prix accepté et délai de vente."""148    rng = np.random.default_rng(graine)149    r = prix_reserve(centre, sigma_o, c_mois, taux_annuel, p_offre)150    prix = np.empty(n); delai = np.empty(n)151    for i in range(n):152        t = 0153        while True:154            t += 1155            if rng.random() <= p_offre:156                o = centre * math.exp(sigma_o * rng.standard_normal())157                if o >= r: prix[i] = o; delai[i] = t; break158    return {"prix_reserve": r, "prix": prix, "delai": delai, "prix_moyen": float(prix.mean()), "prix_median": float(np.median(prix)),159            "ecart_type_prix": float(prix.std(ddof=1)), "cv_prix": float(prix.std(ddof=1) / prix.mean()),160            "delai_moyen": float(delai.mean()), "delai_median": float(np.median(delai)),161            "delai_theorique": 1 / (p_offre * (1 - norm.cdf(math.log(r / centre) / sigma_o)))}162163PARAMS_RECHERCHE = {"sigma_o": 0.06, "c_mois": 0.004, "taux_annuel": 0.05, "p_offre": 0.6}164165def tableau_recherche(centre: float) -> pd.DataFrame:166    rows = []167    for c in (0.002, 0.004, 0.008, 0.015):168        s = simuler_ventes(centre, PARAMS_RECHERCHE["sigma_o"], c, PARAMS_RECHERCHE["taux_annuel"], PARAMS_RECHERCHE["p_offre"], n=10_000)169        rows.append({"cout_detention_mensuel": c, "prix_reserve_sur_V": s["prix_reserve"] / centre, "delai_moyen_mois": s["delai_moyen"],170                     "delai_theorique_mois": s["delai_theorique"], "prix_moyen_sur_V": s["prix_moyen"] / centre, "cv_prix": s["cv_prix"]})171    return pd.DataFrame(rows).set_index("cout_detention_mensuel")172173# ---------------------------------------------------------------------------- 5. PHMU : valeur résiduelle du terrain (hypothèses déclarées)174175USAGES_PHMU = pd.DataFrame({176    "usage": ["Unifamiliale (usage actuel)", "Duplex", "Triplex (si le zonage le permet)"],177    "valeur_produit_fini": [640_000, 880_000, 1_150_000],178    "cout_total_construction": [470_000, 720_000, 1_010_000],179}).set_index("usage")180181def phmu() -> pd.DataFrame:182    """Valeur résiduelle du terrain = valeur du produit fini − coût total (construction, frais, profit du promoteur inclus).183    Chiffres HYPOTHÉTIQUES, pédagogiques, déclarés dans le texte ; non calibrés."""184    t = USAGES_PHMU.copy(); t["valeur_residuelle_terrain"] = t["valeur_produit_fini"] - t["cout_total_construction"]185    return t186187# ---------------------------------------------------------------------------- 6. ancrage : modèle descriptif188189def estimation_ancree(valeur: float, ancre: float, lam: float) -> float:190    """Estimation = λ·ancre + (1 − λ)·valeur (modèle descriptif de l'ancrage-ajustement)."""191    return lam * ancre + (1 - lam) * valeur192193# ---------------------------------------------------------------------------- tout194195def tout() -> dict:196    r = modele(); pred = prediction_sujet(r)197    out = {**r, "pred": pred, "contrib": contributions(r, pred), "rendements": rendements_decroissants(r, pred),198           "anticipation": anticipation(r, pred), "conformite": conformite(r, pred), "dispersion": dispersion(r),199           "sensibilite": tableau_sensibilite(), "comparables": comparables_ch16(r), "phmu": phmu()}200    V = pred["mediane"]201    out["recherche"] = simuler_ventes(V, **PARAMS_RECHERCHE)202    out["recherche_tableau"] = tableau_recherche(V)203    return out204205def main():206    pd.set_option("display.width", 220); pd.set_option("display.max_columns", 30)207    t = tout(); m = t["modele"]; p = t["pred"]; d = t["dispersion"]208    print(f"== Modèle hédonique cas 1 : n = {m['n']}, k = {m['k']}, R² = {m['r2']:.4f}, R² ajusté = {m['r2_ajuste']:.4f}, sigma = {m['sigma']:.4f}")209    print(f"== Résidus : écart-type {d['ecart_type_residus']:.4f}, asymétrie {d['asymetrie']:.3f}, aplatissement (excès) {d['aplatissement_exces']:.3f}, "210          f"q2,5 % {d['q025']:.4f}, q97,5 % {d['q975']:.4f} ; part des ventes à ±10 % de la valeur estimée {d['part_dans_pm10pct']:.3f}, à ±5 % {d['part_dans_pm5pct']:.3f}")211    print(f"== Sujet : ln V = {p['ln_pred']:.4f} ; médiane {p['mediane']:,.0f} $ ; moyenne (Duan ×{p['facteur_duan']:.5f}) {p['moyenne_duan']:,.0f} $ ; "212          f"moyenne log-normale {p['moyenne_lognormale']:,.0f} $ ; mode {p['mode_lognormal']:,.0f} $ ; IP 95 % [{p['ip_bas']:,.0f} ; {p['ip_haut']:,.0f}] ; se(log) {p['se_pred_log']:.4f}")213    print(f"== P(|P/V − 1| > 10 %) = {proba_ecart(m['sigma'], 0.10):.3f} ; > 5 % = {proba_ecart(m['sigma'], 0.05):.3f} ; > 20 % = {proba_ecart(m['sigma'], 0.20):.4f}")214    print("== Prix implicites et contributions pour le sujet :\n", t["contrib"].round(5).to_string())215    rd = t["rendements"]216    print(f"== Rendements décroissants : élasticité {rd['elasticite_S']:.4f} ; prix moyen {rd['prix_moyen_m2']:,.0f} $/m² ; prix marginal {rd['prix_marginal_m2']:,.0f} $/m² ; "217          f"+10 m² = {rd['gain_10m2']:,.0f} $ ; 10 m² suivants = {rd['gain_10m2_suivants']:,.0f} $")218    a = t["anticipation"]219    print(f"== Anticipation : {a['croissance_mensuelle']*100:.4f} %/mois = {a['croissance_annuelle']*100:.2f} %/an ; {a['mois_depuis_reference']:.1f} mois depuis 2024-07-01, "220          f"facteur {a['facteur_cumule']:.4f} ; valeur du sujet à la référence {a['valeur_a_la_reference']:,.0f} $ vs rôle fictif {a['valeur_role']:,.0f} $")221    print("== Conformité (secteurs) :\n", t["conformite"].round(4).to_string())222    c = t["comparables"].copy(); c["date_vente"] = c["date_vente"].dt.strftime("%Y-%m-%d")223    print("== Comparables du chapitre 16 : prix observés vs valeur estimée :\n", c.round(3).to_string(index=False))224    s = t["recherche"]225    print(f"== Recherche (σ_o {PARAMS_RECHERCHE['sigma_o']}, c {PARAMS_RECHERCHE['c_mois']}/mois, i {PARAMS_RECHERCHE['taux_annuel']}, p {PARAMS_RECHERCHE['p_offre']}) : "226          f"prix de réserve {s['prix_reserve']:,.0f} $ ({s['prix_reserve']/p['mediane']:.4f} V) ; prix moyen {s['prix_moyen']:,.0f} $ ; CV {s['cv_prix']:.4f} ; "227          f"délai moyen {s['delai_moyen']:.2f} mois (théorique {s['delai_theorique']:.2f}), médian {s['delai_median']:.0f}")228    print("== Recherche : sensibilité au coût de détention :\n", t["recherche_tableau"].round(4).to_string())229    print("== PHMU (hypothèses déclarées) :\n", t["phmu"].to_string())230    print("== Sensibilité à sigma :\n", t["sensibilite"].round(4).to_string())231    return t232233if __name__ == "__main__":234    main()235