SPB Git

spb/valoplex Public

ValoPlex — moteur d'évaluation spécialisé pour les plex au Québec, petit frère de Vrai-Prix.

TypeScript 90.3% Python 7.1% CSS 2.5%
16.6 KB · 364 lines python
Raw Blame History
1#!/usr/bin/env python32# Auteur : Simon-Pierre Boucher — contact@spboucher.ai3"""ValoPlex — modèle hédonique avancé spécialisé plex (2+ logements, CUBF 1000).45Avancées vs modèle généraliste Vrai-Prix :6  · CIBLE EN RATIO : le modèle prédit log(prix / valeur au rôle), pas le prix.7    Invariant d'échelle, il évalue aussi justement le duplex de 400 k$ que la8    tour au rôle de 700 M$ (un modèle en prix brut compresse les gros9    immeubles vers la masse des duplex — bogue observé et corrigé) ;10  · domaine restreint aux plex, features « par porte » (aire/porte,11    valeur au rôle/porte, bâtiment $/m²) et signaux de mixité ;12  · fourchette P10-P90 par régression quantile PUIS calibration conforme13    (split conformal) : la couverture 80 % est garantie empiriquement sur un14    ensemble de calibration jamais vu ;15  · validation IAAO segmentée (par nombre de portes, par année, temporel 2026).1617Usage : python3 hedonic_plex.py train | predict18"""19import json20import os21import sys2223import joblib24import lightgbm as lgb25import numpy as np26import pandas as pd2728QC = "/Users/simon-pierreboucher/Desktop/qc_house_eval"29BASE = "/Users/simon-pierreboucher/Desktop/valoplex"30MODELS = f"{BASE}/models"31ESTIM = f"{BASE}/data/estimations"32os.makedirs(MODELS, exist_ok=True)33os.makedirs(ESTIM, exist_ok=True)3435T0 = pd.Timestamp("2021-01-01")3637NUM = ["portes", "aire_etages_m2", "aire_par_porte", "superficie_terrain_m2",38       "front_terrain_m", "nb_etages", "nb_chambres_locatives",39       "nb_locaux_non_resid", "age", "lat", "lng", "n_adresses",40       "t_mois", "mois", "anciennete_role_mois",41       "log_val_terrain", "log_val_batiment", "log_val_immeuble",42       "val_par_porte", "bat_par_m2", "part_terrain", "croissance_role"]43CAT = ["lien_physique", "genre_construction", "code_mun_i", "region"]44FEATS = NUM + CAT45LGB = dict(objective="regression", metric="l1", num_leaves=192,46           learning_rate=0.045, n_estimators=4000, min_child_samples=30,47           colsample_bytree=0.85, subsample=0.85, subsample_freq=1,48           reg_lambda=1.5, n_jobs=16, verbose=-1)495051def c_tx(name):  # colonnes côté transactions enrichies52    m = {"rl0308a": "role_aire_etages_m2", "rl0302a": "role_superficie_terrain_m2",53         "rl0301a": "role_front_terrain_m", "rl0306a": "role_nb_etages",54         "rl0311a": "role_nb_logements", "rl0312a": "role_nb_locaux_non_resid",55         "rl0313a": "role_nb_chambres_locatives", "lat": "role_lat",56         "lng": "role_lng", "n_adresses": "role_n_adresses",57         "rl0309a": "role_lien_physique", "rl0310a": "role_genre_construction",58         "code_mun": "role_code_mun", "rl0404a": "role_valeur_immeuble",59         "rl0402a": "role_valeur_terrain", "rl0403a": "role_valeur_batiment",60         "rl0405a": "role_valeur_role_anterieur", "rl0307a": "role_annee_construction",61         "dat_cond_mrche": "role_date_cond_marche"}62    return m[name]636465def eng(df, col, when):66    """Features plex. `col` : fonction nom générique -> nom de colonne."""67    o = pd.DataFrame(index=df.index)68    num = lambda n: pd.to_numeric(df[col(n)], errors="coerce")69    portes = num("rl0311a")70    aire = num("rl0308a")71    val = num("rl0404a")72    terr = num("rl0402a")73    bat = num("rl0403a")74    prev = num("rl0405a")75    yb = num("rl0307a")7677    o["portes"] = portes78    o["aire_etages_m2"] = aire79    o["aire_par_porte"] = (aire / portes.replace(0, np.nan)).clip(10, 500)80    o["superficie_terrain_m2"] = num("rl0302a")81    o["front_terrain_m"] = num("rl0301a")82    o["nb_etages"] = num("rl0306a")83    o["nb_chambres_locatives"] = num("rl0313a")84    o["nb_locaux_non_resid"] = num("rl0312a")85    o["age"] = (when.dt.year - yb).where((yb > 1600) & (yb <= when.dt.year + 2))86    o["lat"] = num("lat")87    o["lng"] = num("lng")88    o["n_adresses"] = num("n_adresses")89    o["t_mois"] = (when - T0).dt.days / 30.4490    o["mois"] = when.dt.month.astype(float)91    # ancienneté du rôle : mois entre la date des conditions du marché du rôle92    # (dat_cond_mrche) et la date d'évaluation — LA variable d'un modèle en ratio :93    # plus le rôle est vieux, plus le marché s'en est éloigné94    dcm = pd.to_datetime(df[col("dat_cond_mrche")].astype(str).str[:10], errors="coerce")95    o["anciennete_role_mois"] = ((when - dcm).dt.days / 30.44).clip(-6, 72)96    o["log_val_terrain"] = np.log1p(terr.clip(lower=0))97    o["log_val_batiment"] = np.log1p(bat.clip(lower=0))98    o["log_val_immeuble"] = np.log1p(val.clip(lower=0))99    o["val_par_porte"] = (val / portes.replace(0, np.nan)).clip(5_000, 2_000_000)100    o["bat_par_m2"] = (bat / aire.replace(0, np.nan)).clip(50, 20_000)101    o["part_terrain"] = (terr / val.replace(0, np.nan)).clip(0, 1)102    o["croissance_role"] = (val / prev.replace(0, np.nan)).clip(0.2, 5)103    o["lien_physique"] = df[col("rl0309a")].astype(str).astype("category")104    o["genre_construction"] = df[col("rl0310a")].astype(str).astype("category")105    cm = pd.to_numeric(df[col("code_mun")].astype(str).str.strip(), errors="coerce")106    o["code_mun_i"] = cm.fillna(-1).astype(int).astype("category")107    o["region"] = (cm // 1000).fillna(-1).astype(int).astype("category")108    return o109110111def align(X, levels=None):112    if levels is None:113        return {k: X[k].cat.categories for k in CAT}114    for k in CAT:115        X[k] = pd.Categorical(X[k], categories=levels[k])116    return levels117118119def build_training():120    df = pd.read_parquet(f"{QC}/province_transactions_enrichi.parquet")121    lg_ = pd.to_numeric(df["role_nb_logements"], errors="coerce")122    cu = pd.to_numeric(df["role_cubf"], errors="coerce")123    val = pd.to_numeric(df["role_valeur_immeuble"], errors="coerce")124    ratio = df["amount"] / val.replace(0, np.nan)125    keep = (126        (cu == 1000) & (lg_ >= 2) & df["role_id_provinc"].notna() &127        (df["match_valeur_exacte"] | (df["match_dist_m"] <= 50)) &128        df["amount"].between(60_000, 30_000_000) & ratio.between(0.25, 4.0)129    )130    df = df[keep].copy()131    df["when"] = pd.to_datetime(df["date"])132    print(f"transactions plex retenues : {len(df)}", flush=True)133    val_ok = pd.to_numeric(df["role_valeur_immeuble"], errors="coerce")134    df = df[val_ok > 10_000].copy()  # cible en ratio : rôle > 0 requis135    X = eng(df, c_tx, df["when"])136    base = pd.to_numeric(df["role_valeur_immeuble"], errors="coerce").values137    y = np.log(df["amount"].values / base)  # log(prix / valeur au rôle)138    meta = df[["id", "date", "amount", "tx_year"]].reset_index(drop=True)139    meta["portes"] = pd.to_numeric(df["role_nb_logements"], errors="coerce").values140    meta["base"] = base141    meta["valeur_role"] = base142    return X.reset_index(drop=True), y, meta143144145def iaao(pred, price):146    r = pred / price147    med = np.median(r)148    cod = 100 * np.mean(np.abs(r - med)) / med149    prd = np.mean(r) / (np.sum(pred) / np.sum(price))150    proxy = 0.5 * (pred / med + price)151    z = np.log(proxy) / np.log(2)152    beta = np.linalg.lstsq(np.c_[np.ones(len(z)), z], r / med - 1, rcond=None)[0]153    return dict(ratio_median=round(float(med), 4), COD=round(float(cod), 2),154                PRD=round(float(prd), 4), PRB=round(float(beta[1]), 4))155156157def metrics(pred, price):158    if len(price) == 0:159        return {}160    ape = np.abs(pred - price) / price161    return dict(n=int(len(price)), MdAPE=round(float(np.median(ape)) * 100, 2),162                MAPE=round(float(np.mean(ape)) * 100, 2),163                dans_10pct=round(float((ape <= .10).mean()) * 100, 1),164                dans_20pct=round(float((ape <= .20).mean()) * 100, 1),165                R2_log=round(float(1 - np.var(np.log(pred) - np.log(price)) /166                                   np.var(np.log(price))), 4),167                **iaao(pred, price))168169170BANDS = [(2, 2), (3, 3), (4, 5), (6, 12), (13, 10**9)]171172173def band_of(doors):174    import numpy as _np175    b = _np.full(len(doors), len(BANDS) - 1)176    for i, (lo, hi) in enumerate(BANDS):177        b[(doors >= lo) & (doors <= hi)] = i178    return b179180181def shrink_weight(doors):182    """Poids du modèle selon le support des ventes : 1 dans le domaine,183    décroissant hors domaine (13+ portes = extrapolation)."""184    import numpy as _np185    w = _np.ones(len(doors))186    w[(doors >= 13) & (doors <= 24)] = 0.55187    w[(doors >= 25) & (doors <= 49)] = 0.35188    w[doors >= 50] = 0.20189    return w190191192def train():193    X, y, meta = build_training()194    rng = np.random.default_rng(20260809)195    n = len(X)196    u = rng.random(n)197    test = u < 0.15                    # évaluation finale198    calib = (u >= 0.15) & (u < 0.30)   # calibration conforme199    fit = u >= 0.30                    # entraînement200    report = {"n_total": int(n), "n_fit": int(fit.sum()),201              "n_calib": int(calib.sum()), "n_test": int(test.sum())}202203    print("=== modèle central (monotone) ===", flush=True)204    m = lgb.LGBMRegressor(**LGB)205    m.fit(X[FEATS][fit], y[fit], eval_set=[(X[FEATS][test], y[test])],206          callbacks=[lgb.early_stopping(200, verbose=False)])207    best = int(m.best_iteration_ or LGB["n_estimators"])208    report["best_iter"] = best209    sm = float(np.mean(np.exp(y[fit] - m.predict(X[FEATS][fit]))))210    report["smearing"] = round(sm, 4)211    base_test = meta.loc[test, "base"].values212    pred_test = base_test * np.exp(m.predict(X[FEATS][test])) * sm213    price_test = base_test * np.exp(y[test])214    report["holdout"] = metrics(pred_test, price_test)215216    # segmentation par portes217    pt = meta.loc[test, "portes"].values218    vr = meta.loc[test, "valeur_role"].values219    for lab, mask in [("2 portes", pt == 2), ("3 portes", pt == 3),220                      ("4-5 portes", (pt >= 4) & (pt <= 5)),221                      ("6-12 portes", (pt >= 6) & (pt <= 12)),222                      ("13+ portes", pt >= 13),223                      ("rôle > 2M$", vr > 2_000_000),224                      ("rôle > 5M$", vr > 5_000_000)]:225        report[f"seg_{lab}"] = metrics(pred_test[mask], price_test[mask])226    # temporel227    t26 = (meta["tx_year"] == 2026).values228    m26 = lgb.LGBMRegressor(**{**LGB, "n_estimators": best})229    m26.fit(X[FEATS][~t26], y[~t26])230    sm26 = float(np.mean(np.exp(y[~t26] - m26.predict(X[FEATS][~t26]))))231    b26 = meta.loc[t26, "base"].values232    report["test_2026"] = metrics(b26 * np.exp(m26.predict(X[FEATS][t26])) * sm26,233                                  b26 * np.exp(y[t26]))234235    print("=== quantiles + calibration conforme ===", flush=True)236    qmods = {}237    for a, nm in [(0.1, "p10"), (0.9, "p90")]:238        mq = lgb.LGBMRegressor(**{**LGB, "objective": "quantile", "alpha": a,239                                  "metric": "quantile", "n_estimators": best})240        mq.fit(X[FEATS][fit], y[fit])241        qmods[nm] = mq242    lo_c = qmods["p10"].predict(X[FEATS][calib])243    hi_c = qmods["p90"].predict(X[FEATS][calib])244    # score de non-conformité (intervalle symétrisé) et quantile 80 %245    s = np.maximum(lo_c - y[calib], y[calib] - hi_c)246    k = int(np.ceil(0.80 * (calib.sum() + 1))) - 1247    qhat = float(np.sort(s)[min(k, len(s) - 1)])248    report["conformal_qhat_log"] = round(qhat, 4)249    # couverture avant/après sur le test250    lo_t = qmods["p10"].predict(X[FEATS][test])251    hi_t = qmods["p90"].predict(X[FEATS][test])252    report["couverture_brute_pct"] = round(float(((y[test] >= lo_t) & (y[test] <= hi_t)).mean()) * 100, 1)253    report["couverture_conforme_pct"] = round(float(((y[test] >= lo_t - qhat) & (y[test] <= hi_t + qhat)).mean()) * 100, 1)254255    # --- calibration Mondrian (par bande de portes) + ancre de rétrécissement ---256    doors_cal = meta.loc[calib, "portes"].values257    bc = band_of(doors_cal)258    band_stats = []259    for i, (blo, bhi) in enumerate(BANDS):260        m_b = bc == i261        yb = y[calib][m_b]262        if m_b.sum() >= 30:263            sb = np.maximum(lo_c[m_b] - yb, yb - hi_c[m_b])264            kb = int(np.ceil(0.80 * (m_b.sum() + 1))) - 1265            qhat_b = float(np.sort(sb)[min(kb, len(sb) - 1)])266        else:267            qhat_b = qhat  # repli global si bande trop mince268        band_stats.append({269            "lo": blo, "hi": bhi, "n_calib": int(m_b.sum()),270            "qhat": qhat_b,271            "ratio_median": float(np.exp(np.median(yb))) if m_b.sum() >= 10 else 1.0,272            "ratio_q10": float(np.exp(np.quantile(yb, 0.10))) if m_b.sum() >= 10 else None,273            "ratio_q90": float(np.exp(np.quantile(yb, 0.90))) if m_b.sum() >= 10 else None,274        })275    report["band_stats"] = band_stats276    # couverture Mondrian sur le test277    bt = band_of(meta.loc[test, "portes"].values)278    q_t = np.array([band_stats[i]["qhat"] for i in bt])279    report["couverture_mondrian_pct"] = round(float(((y[test] >= lo_t - q_t) & (y[test] <= hi_t + q_t)).mean()) * 100, 1)280281    # modèle final : toutes les données282    print("=== réentraînement final ===", flush=True)283    mf = lgb.LGBMRegressor(**{**LGB, "n_estimators": int(best * 1.15)})284    mf.fit(X[FEATS], y)285    sm_f = float(np.mean(np.exp(y - mf.predict(X[FEATS]))))286    qf = {}287    for a, nm in [(0.1, "p10"), (0.9, "p90")]:288        mq = lgb.LGBMRegressor(**{**LGB, "objective": "quantile", "alpha": a,289                                  "metric": "quantile", "n_estimators": best})290        mq.fit(X[FEATS], y)291        qf[nm] = mq292293    imp = pd.Series(mf.feature_importances_, index=FEATS).sort_values(ascending=False)294    report["importance"] = {k: int(v) for k, v in imp.head(12).items()}295    joblib.dump({"model": mf, "smearing": sm_f, "p10": qf["p10"], "p90": qf["p90"],296                 "qhat": qhat, "band_stats": band_stats, "levels": align(X),297                 "feats": FEATS},298                f"{MODELS}/valoplex_lgb.joblib")299    with open(f"{MODELS}/rapport_validation.json", "w") as f:300        json.dump(report, f, indent=2, ensure_ascii=False, default=str)301    print(json.dumps({k: v for k, v in report.items() if k != "importance"},302                     indent=2, ensure_ascii=False))303    print("TRAIN_DONE")304305306def predict():307    b = joblib.load(f"{MODELS}/valoplex_lgb.joblib")308    mf, sm, qhat = b["model"], b["smearing"], b["qhat"]309    idx = pd.read_csv(f"{QC}/data/indexRole2026.csv")310    muns = dict(zip(idx.iloc[:, 0].astype(str).str.lstrip("0"), idx.iloc[:, 1]))311    ident = lambda n: n312    for Y in range(2021, 2027):313        df = pd.read_parquet(f"{QC}/data/parquet/role_{Y}.parquet")314        cu = pd.to_numeric(df["rl0105a"], errors="coerce")315        lg_ = pd.to_numeric(df["rl0311a"], errors="coerce")316        df = df[(cu == 1000) & (lg_ >= 2)].reset_index(drop=True)317        when = pd.Series(pd.Timestamp(f"{Y}-06-15"), index=df.index)318        X = eng(df, ident, when)319        align(X, b["levels"])320        base = pd.to_numeric(df["rl0404a"], errors="coerce").values321        ok = base > 10_000322        doors_arr = pd.to_numeric(df["rl0311a"], errors="coerce").fillna(2).values323        bands_arr = band_of(doors_arr)324        bs = b["band_stats"]325        pl_raw = mf.predict(X[b["feats"]])326        # rétrécissement vers l'ancre empirique de la bande hors domaine327        w = shrink_weight(doors_arr)328        anchor = np.array([np.log(max(bs[i]["ratio_median"], 0.05)) for i in bands_arr])329        pl = w * pl_raw + (1 - w) * anchor330        # intervalles : quantiles + qhat de bande dans le domaine ;331        # dispersion empirique de bande recentrée hors domaine (13+ portes)332        q_b = np.array([bs[i]["qhat"] for i in bands_arr])333        lo = b["p10"].predict(X[b["feats"]]) - q_b334        hi = b["p90"].predict(X[b["feats"]]) + q_b335        big = doors_arr >= 13336        if big.any():337            med_b = np.array([np.log(max(bs[i]["ratio_median"], 0.05)) for i in bands_arr])338            q10_b = np.array([np.log(bs[i]["ratio_q10"]) if bs[i]["ratio_q10"] else med_b[k] - 0.4339                              for k, i in enumerate(bands_arr)])340            q90_b = np.array([np.log(bs[i]["ratio_q90"]) if bs[i]["ratio_q90"] else med_b[k] + 0.4341                              for k, i in enumerate(bands_arr)])342            lo[big] = pl[big] + (q10_b[big] - med_b[big])343            hi[big] = pl[big] + (q90_b[big] - med_b[big])344        out = pd.DataFrame({345            "id_provinc": df["id_provinc"], "annee": Y,346            "code_mun": df["code_mun"], "arrond": df["arrond"],347            "adresse": df["adresse"], "lat": df["lat"], "lng": df["lng"],348            "portes": lg_[ (cu == 1000) & (lg_ >= 2) ].values,349            "valeur_role": base,350            "valeur_estimee": np.where(ok, np.round(base * np.exp(pl) * sm, -2), np.nan),351            "valeur_estimee_p10": np.where(ok, np.round(base * np.exp(lo) * sm, -2), np.nan),352            "valeur_estimee_p90": np.where(ok, np.round(base * np.exp(hi) * sm, -2), np.nan),353        })354        out["municipalite"] = out["code_mun"].astype(str).str.strip() \355            .str.lstrip("0").map(muns)356        out.to_parquet(f"{ESTIM}/plex_estimes_{Y}.parquet", index=False)357        print(f"[{Y}] {len(out)} plex | médiane {out['valeur_estimee'].median():,.0f} $",358              flush=True)359    print("PREDICT_DONE")360361362if __name__ == "__main__":363    {"train": train, "predict": predict}[sys.argv[1]]()364