spb/wp2_uqo Public

1#!/usr/bin/env python32# Author: Simon-Pierre Boucher — contact@spboucher.ai3#4"""Step 4 — Robustness checks for the full model (D).56Faithful port of hedonic_maison.py step 12 (same seeds, same estimators):7  (a) VIF for all regressors8  (b) Breusch-Pagan heteroskedasticity test9  (c) quantile regressions (tau = 0.25, 0.50, 0.75)10  (d) winsorized re-estimation (price trimmed at 1%/99%)11  (e) bootstrap standard errors (B = 1,000, seed 42)12  (f) Lasso / Elastic Net variable selection (5-fold CV, seed 42)1314Inputs : data/processed/hedonic_maison_results.csv15Outputs: results/robustness_results.csv (long format: test, variable, value)16"""1718import sys19from pathlib import Path2021import numpy as np22import pandas as pd23import statsmodels.api as sm24from sklearn.linear_model import ElasticNetCV, LassoCV25from statsmodels.stats.diagnostic import het_breuschpagan26from statsmodels.stats.outliers_influence import variance_inflation_factor2728sys.path.insert(0, str(Path(__file__).resolve().parents[1]))2930from src import config31from src.models import fit_all_models, fit_ols, standardized_design32from src.references import SIM_COLS333435def main():36    df = pd.read_csv(config.ANALYSIS_CSV)37    y = df["log_price"]38    models, specs, sig_sims = fit_all_models(df)39    mD = models["D"]40    full_vars = specs["D"]4142    rows = []43    X_full = standardized_design(df, full_vars)4445    # (a) VIF ---------------------------------------------------------------46    print("VIF ...")47    for i, col in enumerate(X_full.columns):48        if col == "const":49            continue50        vif = variance_inflation_factor(X_full.values, i)51        if col.startswith("sim_"):52            rows.append({"test": "VIF", "variable": col, "statistic": "VIF",53                         "value": round(vif, 4), "p_value": np.nan,54                         "note": "VIF>10 indicates severe multicollinearity"})5556    # (b) Breusch-Pagan -----------------------------------------------------57    print("Breusch-Pagan ...")58    ols_plain = sm.OLS(y, X_full).fit()59    bp_stat, bp_pval, bp_f, bp_fp = het_breuschpagan(ols_plain.resid, X_full)60    rows.append({"test": "Breusch-Pagan", "variable": "Model D (global)",61                 "statistic": "LM", "value": round(bp_stat, 4),62                 "p_value": round(bp_pval, 6),63                 "note": "HC3 justified" if bp_pval < 0.05 else "Homoskedastic"})64    rows.append({"test": "Breusch-Pagan", "variable": "Model D (global)",65                 "statistic": "F", "value": round(bp_f, 4),66                 "p_value": round(bp_fp, 6), "note": ""})6768    # (c) Quantile regressions ----------------------------------------------69    print("Quantile regressions ...")70    for tau in (0.25, 0.50, 0.75):71        qr = sm.QuantReg(y, X_full).fit(q=tau, max_iter=5000)72        for col in SIM_COLS:73            rows.append({"test": f"QuantReg_tau{tau}", "variable": col,74                         "statistic": "coeff",75                         "value": round(qr.params.get(col, np.nan), 6),76                         "p_value": round(qr.pvalues.get(col, 1.0), 6),77                         "note": f"tau={tau}"})7879    # (d) Winsorized regression ----------------------------------------------80    print("Winsorized regression ...")81    p01, p99 = df["price"].quantile([0.01, 0.99])82    df_w = df[(df["price"] >= p01) & (df["price"] <= p99)]83    mD_w = fit_ols(df_w, full_vars)84    for col in SIM_COLS:85        c_full, c_wins = mD.params[col], mD_w.params[col]86        rows.append({"test": "Winsorized_1pct", "variable": col,87                     "statistic": "coeff", "value": round(c_wins, 6),88                     "p_value": round(mD_w.pvalues[col], 6),89                     "note": f"full={round(c_full, 6)}, delta={round(c_wins - c_full, 6)}"})90    rows.append({"test": "Winsorized_1pct", "variable": "Model D (global)",91                 "statistic": "R2", "value": round(mD_w.rsquared, 6),92                 "p_value": np.nan, "note": f"n={int(mD_w.nobs)}"})9394    # (e) Bootstrap SEs (seed and loop identical to the original) ------------95    print("Bootstrap (B=1000) ...")96    np.random.seed(config.SEED)97    B, n = 1000, len(df)98    boot_coefs = {col: [] for col in SIM_COLS}99    for _ in range(B):100        idx = np.random.choice(n, size=n, replace=True)101        mb = sm.OLS(y.iloc[idx], X_full.iloc[idx]).fit()102        for col in SIM_COLS:103            boot_coefs[col].append(mb.params[col])104    for col in SIM_COLS:105        bc = np.array(boot_coefs[col])106        hc3_se = mD.bse[col]107        ci_lo, ci_hi = np.percentile(bc, [2.5, 97.5])108        rows.append({"test": "Bootstrap_1000", "variable": col,109                     "statistic": "boot_SE", "value": round(bc.std(), 6),110                     "p_value": np.nan,111                     "note": (f"HC3_SE={round(hc3_se, 6)}, "112                              f"ratio={round(bc.std() / hc3_se, 4)}, "113                              f"CI=[{round(ci_lo, 5)},{round(ci_hi, 5)}]")})114115    # (f) Lasso / Elastic Net -------------------------------------------------116    print("Lasso / Elastic Net ...")117    X_lasso = X_full.drop(columns="const")118    lasso = LassoCV(cv=5, random_state=config.SEED, max_iter=10000).fit(X_lasso, y)119    lasso_coefs = pd.Series(lasso.coef_, index=X_lasso.columns)120    for col in SIM_COLS:121        c = lasso_coefs[col]122        rows.append({"test": "Lasso_CV", "variable": col, "statistic": "coeff",123                     "value": round(c, 6), "p_value": np.nan,124                     "note": "selected" if abs(c) > 0 else "dropped"})125126    enet = ElasticNetCV(l1_ratio=[0.1, 0.5, 0.7, 0.9, 0.95, 0.99],127                        cv=5, random_state=config.SEED, max_iter=10000).fit(X_lasso, y)128    enet_coefs = pd.Series(enet.coef_, index=X_lasso.columns)129    for col in SIM_COLS:130        c = enet_coefs[col]131        rows.append({"test": "ElasticNet_CV", "variable": col, "statistic": "coeff",132                     "value": round(c, 6), "p_value": np.nan,133                     "note": f"selected, l1={enet.l1_ratio_}" if abs(c) > 0 else "dropped"})134135    config.RESULTS_DIR.mkdir(parents=True, exist_ok=True)136    out = config.RESULTS_DIR / "robustness_results.csv"137    pd.DataFrame(rows).to_csv(out, index=False)138    print(f"-> {out} ({len(rows)} rows)")139140    n_lasso = int((lasso_coefs[SIM_COLS].abs() > 0).sum())141    print(f"Lasso alpha={lasso.alpha_:.6f}, semantic features kept: {n_lasso}/20")142    print(f"Winsorized R2={mD_w.rsquared:.4f} (full: {mD.rsquared:.4f})")143144145if __name__ == "__main__":146    main()147