Python 52.4%
TeX 47.6%
1#!/usr/bin/env python32# -*- coding: utf-8 -*-3"""Estimations complémentaires pour la révision JHE :4 - OLS A/B avec SE HC3 ET clusterisées par municipalité5 - ablation in-sample du modèle A (base / +texte / +texte+vision)6 - effets fixes de ville extrêmes7 - CV 5-fold OLS et RF pour les panels A (texte) et B (texte+vision)8 Résultats -> revision_extras.json9"""10import json, warnings11import numpy as np, pandas as pd12import statsmodels.formula.api as smf13from sklearn.linear_model import LinearRegression14from sklearn.ensemble import RandomForestRegressor, HistGradientBoostingRegressor15from sklearn.model_selection import KFold, cross_val_predict16from sklearn.metrics import r2_score17warnings.filterwarnings("ignore")1819# pipeline v5 jusqu'avant la CV (construit df, AMEN, rv, gv, vv, NUM_v4, NUM_v5, MISS, cats, y, kf)20exec(open("v5_ml.py").read().split("kf=KFold")[0])21amen_ = list(AMEN.keys())22out = {}2324# ---------- OLS A/B : HC3 vs cluster ----------25df["ln_price"] = np.log(df["price"])26df["ln_sqft"] = np.log(df["sqft_final"].fillna(df.groupby("ptype")["sqft_final"].transform("median")).fillna(df["sqft_final"].median()))27df["ln_lot"] = np.log(df["lot_sqft"].fillna(df["lot_sqft"].median())); df["has_lot"] = df["lot_sqft"].notna().astype(int)28df["has_age"] = df["age"].notna().astype(int)29for c in ["beds","baths","storeys","age","age2"]:30 df[c] = df[c].fillna(df.groupby("ptype")[c].transform("median")).fillna(df[c].median())31df["parking"] = df["parking"].fillna(0)3233form = ("ln_price ~ ln_sqft+ln_lot+has_lot+beds+baths+parking+storeys+age+age2+has_age+C(source)+C(ptype_fe)+"34 + "+".join(amen_) + "+n_images_c+virtual_tour+" + "+".join(rv) + "+" + "+".join(gv) + "+" + "+".join(vv) + "+C(city_fe)")35dfA = df[df["has_sqft"] == 1]36mA = smf.ols(form, data=dfA).fit(cov_type="HC3")37mB = smf.ols(form.replace("ln_sqft+",""), data=df).fit(cov_type="HC3")38mAc = smf.ols(form, data=dfA).fit(cov_type="cluster", cov_kwds={"groups": dfA["city_fe"]})39mBc = smf.ols(form.replace("ln_sqft+",""), data=df).fit(cov_type="cluster", cov_kwds={"groups": df["city_fe"]})4041keys = ["ln_sqft","ln_lot","beds","baths","parking","age","age2",42 "central_air","transit","waterfront",43 "g_cond","g_renovated_kitchen","g_hardwood","g_new_construction","g_income_property",44 "v_curb","v_reno","v_luxury","v_garage_vis",45 "v_mat_brique","v_mat_pierre","v_mat_vinyle","v_set_rural","v_land","v_cond"]46def grab(m, mc):47 o = {}48 for k in keys:49 if k in m.params.index:50 o[k] = {"b": float(m.params[k]), "se_hc3": float(m.bse[k]), "p_hc3": float(m.pvalues[k]),51 "se_cl": float(mc.bse[k]), "p_cl": float(mc.pvalues[k])}52 return o53out["modelA"] = {"r2": mA.rsquared, "adj": mA.rsquared_adj, "n": int(mA.nobs), "coef": grab(mA, mAc)}54out["modelB"] = {"r2": mB.rsquared, "adj": mB.rsquared_adj, "n": int(mB.nobs), "coef": grab(mB, mBc)}55print("A R2=%.4f N=%d | B R2=%.4f N=%d" % (mA.rsquared, mA.nobs, mB.rsquared, mB.nobs))5657# ---------- ablation in-sample modèle A ----------58base_form = ("ln_price ~ ln_sqft+ln_lot+has_lot+beds+baths+parking+storeys+age+age2+has_age+C(source)+C(ptype_fe)+"59 + "+".join(amen_) + "+n_images_c+virtual_tour+" + "+".join(rv) + "+C(city_fe)")60mA0 = smf.ols(base_form, data=dfA).fit()61mA1 = smf.ols(base_form + "+" + "+".join(gv), data=dfA).fit()62mA2 = smf.ols(base_form + "+" + "+".join(gv) + "+" + "+".join(vv), data=dfA).fit()63out["ablationA_insample"] = {"base": mA0.rsquared, "plus_text": mA1.rsquared, "plus_text_vision": mA2.rsquared}64print("ablation A in-sample:", out["ablationA_insample"])6566# ---------- effets fixes ville extrêmes ----------67fes = {k.split("T.")[1].rstrip("]"): float(v) for k, v in mA.params.items() if k.startswith("C(city_fe)")}68top = sorted(fes.items(), key=lambda kv: -kv[1])[:12]69bot = sorted(fes.items(), key=lambda kv: kv[1])[:8]70out["city_fe_top"] = [(k, (np.exp(v) - 1) * 100) for k, v in top]71out["city_fe_bot"] = [(k, (np.exp(v) - 1) * 100) for k, v in bot]72print("top villes:", out["city_fe_top"][:4])7374# ---------- CV : OLS et RF, panels texte / texte+vision ----------75kf = KFold(n_splits=5, shuffle=True, random_state=42)76Xcat = pd.get_dummies(df[cats], prefix=cats, dtype=float)77dfi = df.copy()78def stats(pred):79 e = np.abs(np.exp(pred - y) - 1)80 return {"r2": float(r2_score(y, pred)), "mdae": float(np.median(e) * 100),81 "w15": float(np.mean(e <= .15) * 100), "w25": float(np.mean(e <= .25) * 100)}82def cv_ols(numcols):83 X = pd.concat([dfi[numcols].reset_index(drop=True), Xcat.reset_index(drop=True)], axis=1)84 X = X.fillna(X.median())85 return stats(cross_val_predict(LinearRegression(), X, y, cv=kf, n_jobs=-1))86def cv_rf(numcols):87 X = pd.concat([df[numcols].reset_index(drop=True), Xcat.reset_index(drop=True)], axis=1)88 m = RandomForestRegressor(n_estimators=400, max_features=0.5, min_samples_leaf=3, n_jobs=-1, random_state=42)89 return stats(cross_val_predict(m, X, y, cv=kf, n_jobs=-1))90def cv_gbm(numcols):91 X = pd.concat([df[numcols].reset_index(drop=True), Xcat.reset_index(drop=True)], axis=1)92 m = HistGradientBoostingRegressor(max_iter=600, learning_rate=0.05, max_leaf_nodes=63, l2_regularization=1.0, random_state=42)93 return stats(cross_val_predict(m, X, y, cv=kf, n_jobs=-1))9495# matrice OLS imputée médiane par type (comme hedonic_ml)96for c in ["sqft_final","lot_sqft","beds","baths","storeys","age","age2"]:97 dfi[c] = dfi.groupby("ptype")[c].transform(lambda s: s.fillna(s.median()))98 dfi[c] = dfi[c].fillna(dfi[c].median())99dfi["parking"] = dfi["parking"].fillna(0)100101print(">>> CV OLS v4/v5")102out["cv"] = {"ols_v4": cv_ols(NUM_v4), "ols_v5": cv_ols(NUM_v5)}103print("OLS:", out["cv"]["ols_v4"], out["cv"]["ols_v5"])104print(">>> CV GBM v4/v5 (contrôle)")105out["cv"]["gbm_v4"] = cv_gbm(NUM_v4); out["cv"]["gbm_v5"] = cv_gbm(NUM_v5)106print("GBM:", out["cv"]["gbm_v4"], out["cv"]["gbm_v5"])107print(">>> CV RF v4/v5 (plus long)")108out["cv"]["rf_v4"] = cv_rf(NUM_v4); out["cv"]["rf_v5"] = cv_rf(NUM_v5)109print("RF:", out["cv"]["rf_v4"], out["cv"]["rf_v5"])110111json.dump(out, open("results/revision_extras.json", "w"), indent=1, default=float)112print("-> revision_extras.json")113