spb/wp10_uqo Public

1#!/usr/bin/env python32# Author: Simon-Pierre Boucher — contact@spboucher.ai3"""Step 03 — Vertical-inequity regressions, heterogeneity, horizontal4inequity and the implied tax shift.56Writes to ``results/reproduced/``:7  vertical.csv        Cheng OLS / Cheng FE / Clapp IV / Paglin–Fogarty8  quantile.csv        quantile-regression β(τ) on within-cell data9  heterogeneity.csv   Cheng-FE γ by property class, age, land share, year,10                      roll lag and municipality size11  binscatter.csv      within-cell mean ln ratio by price vigintile12  horizontal.csv      |deviation| regression (who gets noisy assessments)13  taxshift.csv        over/under-taxation by within-cell price decile14  robustness.csv      γ across sample and measurement variants1516Usage:  python scripts/03_estimate_regressions.py17"""18import sys19from pathlib import Path2021import numpy as np22import pandas as pd2324sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src"))2526from wp10 import config, models, sample  # noqa: E402272829def main() -> None:30    config.ensure_dirs()31    df = sample.load()32    out = config.REPRODUCED3334    # ------------------------------------------------------------ main table35    print("Vertical-inequity estimators")36    res = [models.cheng_pooled(df), models.cheng_fe(df), models.clapp_iv(df)]37    for r in res:38        print(f"  {r['estimator']:<28} beta={r['beta']:.4f} ({r['se']:.4f})  "39              f"gamma={r['gamma']:+.4f}  n={r['n']:,}")40    pf = models.paglin_fogarty(df)41    print(f"  {pf['estimator']:<28} a={pf['intercept']:,.0f} "42          f"({pf['intercept_se']:,.0f})  b={pf['slope']:.4f}")43    pd.DataFrame(res).to_csv(out / "vertical.csv", index=False)44    pd.Series(pf).to_csv(out / "paglin_fogarty.csv")4546    # ------------------------------------------------------------ quantiles47    qt = models.quantile_betas(df)48    qt.to_csv(out / "quantile.csv", index=False)49    print("Quantile betas:", {f"{r.tau:.2f}": round(r.beta, 3)50                              for r in qt.itertuples()})5152    # ------------------------------------------------------------ binscatter53    d = df.copy()54    d["lnr_w"] = d["ln_ratio"] - d.groupby("cell")["ln_ratio"].transform("mean")55    d["lnp_w"] = d["ln_price"] - d.groupby("cell")["ln_price"].transform("mean")56    d["bin"] = pd.qcut(d["lnp_w"], 20, labels=False)57    (d.groupby("bin")58       .agg(x=("lnp_w", "mean"), y=("lnr_w", "mean"),59            se=("lnr_w", lambda s: s.std() / np.sqrt(len(s))), n=("lnr_w", "size"))60       .reset_index()61       .to_csv(out / "binscatter.csv", index=False))6263    # ------------------------------------------------------------ heterogeneity64    muni_sales = df.groupby("muni")["muni"].transform("size")65    groups = {66        "Single-family": df["prop_class"] == "single_family",67        "Condominium": df["prop_class"] == "condo",68        "Plex (2–5 units)": df["prop_class"] == "plex",69        "Cottage": df["prop_class"] == "cottage",70        "Age < 20 y": df["age"] < 20,71        "Age 20–60 y": df["age"].between(20, 60),72        "Age > 60 y": df["age"] > 60,73        "Land share < 0.2": df["land_share"] < 0.2,74        "Land share 0.2–0.4": df["land_share"].between(0.2, 0.4),75        "Land share > 0.4": df["land_share"] > 0.4,76        "Roll lag < 24 m": df["lag_months"] < 24,77        "Roll lag 24–48 m": df["lag_months"].between(24, 48),78        "Roll lag > 48 m": df["lag_months"] > 48,79        "Muni < 1k sales": muni_sales < 1_000,80        "Muni 1k–10k sales": muni_sales.between(1_000, 10_000),81        "Muni > 10k sales": muni_sales > 10_000,82    }83    groups.update({f"Sales {y}": df["sale_year"] == y84                   for y in sorted(df["sale_year"].unique())})85    het = models.gamma_by_group(df, groups)86    het.to_csv(out / "heterogeneity.csv", index=False)87    print(f"Heterogeneity: {len(het)} subgroups estimated")8889    # ------------------------------------------------------------ horizontal90    tab, meta = models.horizontal_dispersion(df)91    tab.to_csv(out / "horizontal.csv")92    pd.Series(meta).to_csv(out / "horizontal_meta.csv")93    print("Horizontal-dispersion regression:", meta)9495    # ------------------------------------------------------------ tax shift96    ts = models.tax_shift(df)97    ts.to_csv(out / "taxshift.csv", index=False)98    print("Tax shift by decile (mean %):",99          {int(r.decile): f"{r.mean_rel:+.1%}" for r in ts.itertuples()})100101    # ------------------------------------------------------------ robustness102    variants = {103        "Baseline": df,104        "Condominiums only": df[df["prop_class"] == "condo"],105        "Excluding sales < $100k": df[df["amount"] >= 100_000],106        "Match score = 220 (max)": df[df["match_score"] >= 219.9],107        "Match distance <= 10 m": df[df["match_dist_m"] <= 10],108        "Single-family only": df[df["prop_class"] == "single_family"],109        "Ratio trim 5/95": df[df["ratio"].between(110            df["ratio"].quantile(.05), df["ratio"].quantile(.95))],111        "Sales 2021-2023": df[df["sale_year"] <= 2023],112        "Sales 2024-2026": df[df["sale_year"] >= 2024],113        "Munis >= 300 sales": df[df.groupby("muni")["muni"]114                                 .transform("size") >= 300],115        "Cells >= 50 sales": df[df.groupby("cell")["cell"]116                                .transform("size") >= 50],117    }118    rows = []119    for label, sub in variants.items():120        counts = sub.groupby("cell")["cell"].transform("size")121        sub = sub[counts >= config.CELL_MIN_OBS]122        est = models.cheng_fe(sub)123        iv = models.clapp_iv(sub)124        rows.append({"variant": label, "gamma_fe": est["gamma"],125                     "se_fe": est["se"], "gamma_iv": iv["gamma"],126                     "se_iv": iv["se"], "n": est["n"]})127        print(f"  {label:<26} gamma_FE={est['gamma']:+.4f} "128              f"gamma_IV={iv['gamma']:+.4f}  n={est['n']:,}")129    pd.DataFrame(rows).to_csv(out / "robustness.csv", index=False)130131132if __name__ == "__main__":133    main()134