spb/wp3_uqo Public

1#!/usr/bin/env python32# Author: Simon-Pierre Boucher — contact@spboucher.ai3#4"""Quantile-regression stability, inter-quantile tests, imputation and5winsorization sensitivity (refactor of the original v3_qr_imputation_analysis.py).67Produces results/v3_qr_imputation_results.pkl with:8  1. QR coefficient stability at tau=0.50 over 10 subsamples of 150,0009  2. Inter-quantile Wald tests (tau=0.10 vs tau=0.90)10  3. Imputation sensitivity: OLS re-estimated after dropping observations with11     originally missing year_built / lot_size_sqft (missingness read from DuckDB)12  4. Winsorization sensitivity: lot size capped at the 99.5th percentile,13     bike score capped at 1001415Usage:  python scripts/02_qr_imputation.py [--output results/v3_qr_imputation_results.pkl]16"""1718import argparse19import pickle20import sys21import time22import warnings23from pathlib import Path2425import duckdb26import numpy as np27import pandas as pd28import statsmodels.api as sm29from sklearn.linear_model import LinearRegression3031sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src"))32from wp3 import config, data3334warnings.filterwarnings("ignore")3536OLS_FEATURES = ["ln_sqft", "bedrooms", "bathrooms", "age", "age_sq",37                "ln_lot", "has_pool", "has_garage", "luxury_score",38                "tag_foreclosure", "on_waterfront"]394041def qr_stability(X_const, y_clean):42    """Median-regression coefficient stability across 10 random subsamples."""43    coef_matrix = {v: [] for v in config.KEY_VARS}44    for i, seed in enumerate(config.QR_STABILITY_SEEDS):45        rng = np.random.RandomState(seed)46        idx = rng.choice(len(X_const), size=config.QR_SUBSAMPLE_SIZE, replace=False)47        print(f"  QR subsample {i + 1}/{len(config.QR_STABILITY_SEEDS)} (seed={seed})...")48        res = sm.QuantReg(y_clean.iloc[idx], X_const.iloc[idx]).fit(q=0.50, max_iter=1000)49        for v in config.KEY_VARS:50            coef_matrix[v].append(res.params[v])5152    stability = {}53    for v, vals_list in coef_matrix.items():54        vals = np.array(vals_list)55        m, s = vals.mean(), vals.std()56        stability[v] = {57            "mean": m, "std": s,58            "cv": abs(s / m) if abs(m) > 1e-10 else np.inf,59            "min": vals.min(), "max": vals.max(),60            "sign_stability_pct": float(np.mean(np.sign(vals) == np.sign(m)) * 100),61        }62    return stability636465def inter_quantile_test(X_const, y_clean):66    """Wald z-tests for coefficient equality between tau=0.10 and tau=0.90."""67    rng = np.random.RandomState(config.RANDOM_STATE)68    idx = rng.choice(len(X_const), size=config.QR_SUBSAMPLE_SIZE, replace=False)69    X_sub, y_sub = X_const.iloc[idx], y_clean.iloc[idx]7071    print("  Fitting QR at tau=0.10 and tau=0.90...")72    res10 = sm.QuantReg(y_sub, X_sub).fit(q=0.10, max_iter=1000)73    res90 = sm.QuantReg(y_sub, X_sub).fit(q=0.90, max_iter=1000)7475    out = {}76    for v in config.KEY_VARS:77        diff = res10.params[v] - res90.params[v]78        z = diff / np.sqrt(res10.bse[v] ** 2 + res90.bse[v] ** 2)79        out[v] = {80            "beta_010": res10.params[v], "beta_090": res90.params[v],81            "difference": diff, "z_stat": z,82            "significant_5pct": bool(abs(z) > 1.96),83        }84    return out858687def load_missingness_flags():88    """Read original missingness of year_built / lot_size_sqft from the raw DB."""89    con = duckdb.connect(str(config.DUCKDB_PATH), read_only=True)90    flags = con.execute("""91        SELECT zpid,92               CASE WHEN year_built IS NULL THEN 1 ELSE 0 END AS yb_missing,93               CASE WHEN lot_size_sqft IS NULL THEN 1 ELSE 0 END AS lot_missing94        FROM properties95    """).fetchdf()96    con.close()97    flags["zpid"] = flags["zpid"].astype(np.int64)98    return flags99100101def build_unstandardized_matrix(df_merged, X_const):102    """Rebuild the OLS design matrix on the raw (unstandardized) scale."""103    xcols = [c for c in X_const.columns if c != "const"]104    missing = [c for c in xcols if c not in df_merged.columns]105    if missing:106        for cat_col in config.CATEGORICAL_COLS:107            if cat_col not in df_merged.columns:108                continue109            dummies = pd.get_dummies(df_merged[cat_col], prefix=cat_col)110            for dc in dummies.columns:111                if dc in missing:112                    df_merged[dc] = dummies[dc].astype(float)113    return df_merged[xcols].values.astype(np.float64), xcols114115116def imputation_sensitivity(df_merged, X_unstd, xcols):117    """OLS on subsets that drop originally-missing year_built / lot_size rows."""118    y = df_merged["ln_price"].values119120    def run(mask, label):121        lr = LinearRegression()122        lr.fit(X_unstd[mask], y[mask])123        coefs = dict(zip(xcols, lr.coef_))124        return {"label": label, "N": int(mask.sum()),125                "R2": lr.score(X_unstd[mask], y[mask]),126                "key_coefficients": {v: coefs[v] for v in OLS_FEATURES}}127128    no_yb = df_merged["yb_missing"].values == 0129    no_lot = df_merged["lot_missing"].values == 0130    scenarios = [131        (np.ones(len(df_merged), dtype=bool), "Full sample (baseline)"),132        (no_yb, "Drop missing year_built"),133        (no_lot, "Drop missing lot_size_sqft"),134        (no_yb & no_lot, "Drop both missing"),135    ]136    out = []137    for mask, label in scenarios:138        r = run(mask, label)139        out.append(r)140        print(f"  {label:<30} N={r['N']:>9,}  R2={r['R2']:.4f}")141    return out142143144def winsorization_sensitivity(df_merged, X_unstd, xcols):145    """Re-estimate OLS with lot size winsorized at p99.5 and bike score capped at 100."""146    df_w = df_merged.copy()147    y = df_merged["ln_price"].values148149    lot_p995 = df_w["lot_size_sqft"].quantile(0.995)150    n_lot = int((df_w["lot_size_sqft"] > lot_p995).sum())151    df_w["lot_size_sqft"] = df_w["lot_size_sqft"].clip(upper=lot_p995)152    df_w["ln_lot"] = np.where(df_w["lot_size_sqft"] > 0, np.log(df_w["lot_size_sqft"]), 0.0)153154    n_bike = int((df_w["bike_score"] > 100).sum())155    df_w["bike_score"] = df_w["bike_score"].clip(upper=100)156157    def fit(X):158        lr = LinearRegression()159        lr.fit(X, y)160        coefs = dict(zip(xcols, lr.coef_))161        return {"R2": lr.score(X, y),162                "key_coefficients": {v: coefs[v] for v in OLS_FEATURES}}163164    X_w = df_w[xcols].values.astype(np.float64)165    return {166        "baseline": fit(X_unstd),167        "winsorized": fit(X_w),168        "lot_p995": lot_p995,169        "n_lot_winsorized": n_lot,170        "n_bike_capped": n_bike,171    }172173174def main():175    parser = argparse.ArgumentParser(description=__doc__)176    parser.add_argument("--output", type=Path,177                        default=config.RESULTS_DIR / "v3_qr_imputation_results.pkl")178    args = parser.parse_args()179180    t0 = time.time()181    print("Loading model data...")182    md = data.load_model_data()183    X_const, y_clean, df_clean = md["X_const"], md["y_clean"], md["df_clean"]184185    results = {}186187    print("\n1. QR stability (tau=0.50, 10 subsamples)...")188    results["qr_stability"] = qr_stability(X_const, y_clean)189190    print("\n2. Inter-quantile difference test (tau=0.10 vs 0.90)...")191    results["iqr_test"] = inter_quantile_test(X_const, y_clean)192193    print("\n3. Imputation sensitivity...")194    flags = load_missingness_flags()195    df_merged = df_clean.merge(flags, on="zpid", how="left")196    X_unstd, xcols = build_unstandardized_matrix(df_merged, X_const)197    results["imputation_sensitivity"] = imputation_sensitivity(df_merged, X_unstd, xcols)198199    print("\n4. Winsorization sensitivity...")200    results["winsorization_sensitivity"] = winsorization_sensitivity(201        df_merged, X_unstd, xcols)202203    args.output.parent.mkdir(parents=True, exist_ok=True)204    with open(args.output, "wb") as f:205        pickle.dump(results, f)206    print(f"\nSaved to {args.output}  ({time.time() - t0:.0f}s total)")207208209if __name__ == "__main__":210    main()211