#!/usr/bin/env python3 # Author: Simon-Pierre Boucher — contact@spboucher.ai # """Step 4 — Robustness checks for the full model (D). Faithful port of hedonic_maison.py step 12 (same seeds, same estimators): (a) VIF for all regressors (b) Breusch-Pagan heteroskedasticity test (c) quantile regressions (tau = 0.25, 0.50, 0.75) (d) winsorized re-estimation (price trimmed at 1%/99%) (e) bootstrap standard errors (B = 1,000, seed 42) (f) Lasso / Elastic Net variable selection (5-fold CV, seed 42) Inputs : data/processed/hedonic_maison_results.csv Outputs: results/robustness_results.csv (long format: test, variable, value) """ import sys from pathlib import Path import numpy as np import pandas as pd import statsmodels.api as sm from sklearn.linear_model import ElasticNetCV, LassoCV from statsmodels.stats.diagnostic import het_breuschpagan from statsmodels.stats.outliers_influence import variance_inflation_factor sys.path.insert(0, str(Path(__file__).resolve().parents[1])) from src import config from src.models import fit_all_models, fit_ols, standardized_design from src.references import SIM_COLS def main(): df = pd.read_csv(config.ANALYSIS_CSV) y = df["log_price"] models, specs, sig_sims = fit_all_models(df) mD = models["D"] full_vars = specs["D"] rows = [] X_full = standardized_design(df, full_vars) # (a) VIF --------------------------------------------------------------- print("VIF ...") for i, col in enumerate(X_full.columns): if col == "const": continue vif = variance_inflation_factor(X_full.values, i) if col.startswith("sim_"): rows.append({"test": "VIF", "variable": col, "statistic": "VIF", "value": round(vif, 4), "p_value": np.nan, "note": "VIF>10 indicates severe multicollinearity"}) # (b) Breusch-Pagan ----------------------------------------------------- print("Breusch-Pagan ...") ols_plain = sm.OLS(y, X_full).fit() bp_stat, bp_pval, bp_f, bp_fp = het_breuschpagan(ols_plain.resid, X_full) rows.append({"test": "Breusch-Pagan", "variable": "Model D (global)", "statistic": "LM", "value": round(bp_stat, 4), "p_value": round(bp_pval, 6), "note": "HC3 justified" if bp_pval < 0.05 else "Homoskedastic"}) rows.append({"test": "Breusch-Pagan", "variable": "Model D (global)", "statistic": "F", "value": round(bp_f, 4), "p_value": round(bp_fp, 6), "note": ""}) # (c) Quantile regressions ---------------------------------------------- print("Quantile regressions ...") for tau in (0.25, 0.50, 0.75): qr = sm.QuantReg(y, X_full).fit(q=tau, max_iter=5000) for col in SIM_COLS: rows.append({"test": f"QuantReg_tau{tau}", "variable": col, "statistic": "coeff", "value": round(qr.params.get(col, np.nan), 6), "p_value": round(qr.pvalues.get(col, 1.0), 6), "note": f"tau={tau}"}) # (d) Winsorized regression ---------------------------------------------- print("Winsorized regression ...") p01, p99 = df["price"].quantile([0.01, 0.99]) df_w = df[(df["price"] >= p01) & (df["price"] <= p99)] mD_w = fit_ols(df_w, full_vars) for col in SIM_COLS: c_full, c_wins = mD.params[col], mD_w.params[col] rows.append({"test": "Winsorized_1pct", "variable": col, "statistic": "coeff", "value": round(c_wins, 6), "p_value": round(mD_w.pvalues[col], 6), "note": f"full={round(c_full, 6)}, delta={round(c_wins - c_full, 6)}"}) rows.append({"test": "Winsorized_1pct", "variable": "Model D (global)", "statistic": "R2", "value": round(mD_w.rsquared, 6), "p_value": np.nan, "note": f"n={int(mD_w.nobs)}"}) # (e) Bootstrap SEs (seed and loop identical to the original) ------------ print("Bootstrap (B=1000) ...") np.random.seed(config.SEED) B, n = 1000, len(df) boot_coefs = {col: [] for col in SIM_COLS} for _ in range(B): idx = np.random.choice(n, size=n, replace=True) mb = sm.OLS(y.iloc[idx], X_full.iloc[idx]).fit() for col in SIM_COLS: boot_coefs[col].append(mb.params[col]) for col in SIM_COLS: bc = np.array(boot_coefs[col]) hc3_se = mD.bse[col] ci_lo, ci_hi = np.percentile(bc, [2.5, 97.5]) rows.append({"test": "Bootstrap_1000", "variable": col, "statistic": "boot_SE", "value": round(bc.std(), 6), "p_value": np.nan, "note": (f"HC3_SE={round(hc3_se, 6)}, " f"ratio={round(bc.std() / hc3_se, 4)}, " f"CI=[{round(ci_lo, 5)},{round(ci_hi, 5)}]")}) # (f) Lasso / Elastic Net ------------------------------------------------- print("Lasso / Elastic Net ...") X_lasso = X_full.drop(columns="const") lasso = LassoCV(cv=5, random_state=config.SEED, max_iter=10000).fit(X_lasso, y) lasso_coefs = pd.Series(lasso.coef_, index=X_lasso.columns) for col in SIM_COLS: c = lasso_coefs[col] rows.append({"test": "Lasso_CV", "variable": col, "statistic": "coeff", "value": round(c, 6), "p_value": np.nan, "note": "selected" if abs(c) > 0 else "dropped"}) enet = ElasticNetCV(l1_ratio=[0.1, 0.5, 0.7, 0.9, 0.95, 0.99], cv=5, random_state=config.SEED, max_iter=10000).fit(X_lasso, y) enet_coefs = pd.Series(enet.coef_, index=X_lasso.columns) for col in SIM_COLS: c = enet_coefs[col] rows.append({"test": "ElasticNet_CV", "variable": col, "statistic": "coeff", "value": round(c, 6), "p_value": np.nan, "note": f"selected, l1={enet.l1_ratio_}" if abs(c) > 0 else "dropped"}) config.RESULTS_DIR.mkdir(parents=True, exist_ok=True) out = config.RESULTS_DIR / "robustness_results.csv" pd.DataFrame(rows).to_csv(out, index=False) print(f"-> {out} ({len(rows)} rows)") n_lasso = int((lasso_coefs[SIM_COLS].abs() > 0).sum()) print(f"Lasso alpha={lasso.alpha_:.6f}, semantic features kept: {n_lasso}/20") print(f"Winsorized R2={mD_w.rsquared:.4f} (full: {mD.rsquared:.4f})") if __name__ == "__main__": main()