spb/wp10_uqo Public

1# Author: Simon-Pierre Boucher — contact@spboucher.ai2"""Econometric estimators for WP10.34Every vertical-inequity estimator ultimately measures how the assessment5ratio r = AV/SP moves with market value. The workhorse is the log-log6(Cheng 1974) regression78    ln AV_i = α_c(i) + β · ln SP_i + ε_i ,910where α_c(i) is a municipality × roll vintage × sale-year fixed effect that11absorbs the mechanical drift of ratios between triennial reference dates.12β < 1 ⇒ regressive assessment. We report γ ≡ β − 1, the elasticity of the13assessment *ratio* with respect to price (γ < 0 ⇒ regressive), estimated by1415  - pooled OLS (no fixed effects),16  - absorbing least squares with cell fixed effects,17  - Clapp's (1990) rank-based IV, which purges the attenuation/mean-reversion18    bias that pushes OLS toward spurious regressivity when sale prices carry19    idiosyncratic noise,20  - quantile regression across the conditional distribution.2122Inference is clustered by municipality throughout.23"""24from __future__ import annotations2526import numpy as np27import pandas as pd28import statsmodels.api as sm29from linearmodels.iv import IV2SLS30from linearmodels.iv.absorbing import AbsorbingLS3132from . import config333435# ---------------------------------------------------------------- helpers36def _demean(df: pd.DataFrame, cols: list[str], by: str = "cell") -> pd.DataFrame:37    """Within-transform ``cols`` inside groups of ``by`` (suffix ``_w``)."""38    out = df.copy()39    for c in cols:40        out[c + "_w"] = df[c] - df.groupby(by)[c].transform("mean")41    return out424344def _cluster_ols(y: np.ndarray, X: np.ndarray, clusters: np.ndarray,45                 dof_adjust: int = 0):46    """OLS with CRVE (cluster-robust) inference via statsmodels."""47    model = sm.OLS(y, X)48    return model.fit(cov_type="cluster", cov_kwds={"groups": clusters})495051# ---------------------------------------------------------------- Cheng / FE52def cheng_pooled(df: pd.DataFrame) -> dict:53    """Pooled Cheng regression: ln AV on ln SP, no fixed effects."""54    X = sm.add_constant(df["ln_price"].to_numpy())55    res = _cluster_ols(df["ln_av"].to_numpy(), X, df["muni"].to_numpy())56    return {"estimator": "Cheng OLS (pooled)", "beta": res.params[1],57            "se": res.bse[1], "gamma": res.params[1] - 1.0,58            "n": int(res.nobs), "r2": res.rsquared}596061def cheng_fe(df: pd.DataFrame, extra_controls: list[str] | None = None) -> dict:62    """Cheng regression absorbing municipality × roll × sale-year cells."""63    dep = df["ln_av"]64    exog_cols = ["ln_price"] + (extra_controls or [])65    exog = df[exog_cols]66    absorb = pd.DataFrame({"cell": df["cell"].astype("category")})67    mod = AbsorbingLS(dep, exog, absorb=absorb)68    res = mod.fit(cov_type="clustered", clusters=df["muni"].astype("category"))69    return {"estimator": "Cheng FE (cell)", "beta": float(res.params["ln_price"]),70            "se": float(res.std_errors["ln_price"]),71            "gamma": float(res.params["ln_price"]) - 1.0,72            "n": int(res.nobs), "r2": float(res.rsquared)}737475def paglin_fogarty(df: pd.DataFrame) -> dict:76    """Levels regression AV = a + b·SP (Paglin & Fogarty 1972).7778    A positive intercept with b below the overall assessment level indicates79    regressivity in levels.80    """81    X = sm.add_constant(df["amount"].to_numpy(float))82    res = _cluster_ols(df["role_valeur_immeuble"].to_numpy(float), X,83                       df["muni"].to_numpy())84    return {"estimator": "Paglin–Fogarty (levels)",85            "intercept": res.params[0], "intercept_se": res.bse[0],86            "slope": res.params[1], "slope_se": res.bse[1],87            "n": int(res.nobs), "r2": res.rsquared}888990# ---------------------------------------------------------------- Clapp IV91def clapp_iv(df: pd.DataFrame) -> dict:92    """Clapp (1990) rank-based IV on within-cell demeaned data.9394    The instrument z ∈ {−1, 0, +1} flags sales in the bottom/top third of95    *both* the ln AV and ln SP within-cell distributions. Because z carries96    only coarse rank information, it is (near-)orthogonal to the transitory97    component of either variable, undoing the attenuation that biases OLS98    toward regressivity (β̂ < β) when prices are noisy.99    """100    d = _demean(df, ["ln_av", "ln_price"])101    g_av = d.groupby("cell")["ln_av"].rank(pct=True)102    g_sp = d.groupby("cell")["ln_price"].rank(pct=True)103    z = np.zeros(len(d))104    z[(g_av <= 1 / 3) & (g_sp <= 1 / 3)] = -1.0105    z[(g_av > 2 / 3) & (g_sp > 2 / 3)] = 1.0106    d["z"] = z107    res = IV2SLS(d["ln_av_w"], None, d[["ln_price_w"]], d[["z"]]).fit(108        cov_type="clustered", clusters=d["muni"].astype("category"))109    beta = float(res.params["ln_price_w"])110    return {"estimator": "Clapp IV (rank instrument)", "beta": beta,111            "se": float(res.std_errors["ln_price_w"]), "gamma": beta - 1.0,112            "n": int(res.nobs), "r2": float(res.rsquared)}113114115# ---------------------------------------------------------------- quantiles116def quantile_betas(df: pd.DataFrame, taus: list[float] = config.QUANTILES,117                   max_n: int = 250_000, seed: int = config.SEED_BOOT) -> pd.DataFrame:118    """Quantile regressions of demeaned ln AV on demeaned ln SP.119120    Estimated on a seeded random subsample (IRLS on the full 600k sample is121    needlessly slow; the subsample SEs are already microscopic).122    """123    d = _demean(df, ["ln_av", "ln_price"])124    if len(d) > max_n:125        d = d.sample(max_n, random_state=seed)126    X = sm.add_constant(d["ln_price_w"].to_numpy())127    y = d["ln_av_w"].to_numpy()128    rows = []129    for tau in taus:130        r = sm.QuantReg(y, X).fit(q=tau, max_iter=2000)131        rows.append({"tau": tau, "beta": r.params[1], "se": r.bse[1],132                     "gamma": r.params[1] - 1.0, "n": len(d)})133    return pd.DataFrame(rows)134135136# ---------------------------------------------------------------- subgroups137def gamma_by_group(df: pd.DataFrame, groups: dict[str, pd.Series]) -> pd.DataFrame:138    """Cheng-FE γ estimated separately on each labelled subsample.139140    ``groups`` maps a label to a boolean mask over ``df``. Subsamples keep141    only cells that retain ≥ CELL_MIN_OBS sales after masking.142    """143    rows = []144    for label, mask in groups.items():145        sub = df[mask]146        counts = sub.groupby("cell")["cell"].transform("size")147        sub = sub[counts >= config.CELL_MIN_OBS]148        if len(sub) < 2_000:149            continue150        est = cheng_fe(sub)151        rows.append({"group": label, "gamma": est["gamma"], "se": est["se"],152                     "n": est["n"]})153    return pd.DataFrame(rows)154155156# ---------------------------------------------------------------- horizontal157def horizontal_dispersion(df: pd.DataFrame) -> tuple[pd.DataFrame, dict]:158    """Horizontal inequity: who gets the noisiest assessments?159160    Regresses the absolute log deviation of a sale's ratio from its cell161    median, |ln r_i − med_c ln r|, on property characteristics with cell162    fixed effects. Positive coefficients = less uniform assessment.163    Returns (coefficient table, fit metadata).164    """165    d = df.copy()166    d["abs_dev"] = (d["ln_ratio"]167                    - d.groupby("cell")["ln_ratio"].transform("median")).abs()168    d["age_dec"] = d["age"] / 10.0169    d["is_condo"] = (d["prop_class"] == "condo").astype(float)170    d["is_plex"] = (d["prop_class"] == "plex").astype(float)171    d["is_cottage"] = (d["prop_class"] == "cottage").astype(float)172    d = d.dropna(subset=["abs_dev", "age_dec", "land_share"])173    exog_cols = ["age_dec", "land_share", "is_condo", "is_plex", "is_cottage"]174    mod = AbsorbingLS(d["abs_dev"], d[exog_cols],175                      absorb=pd.DataFrame({"cell": d["cell"].astype("category")}))176    res = mod.fit(cov_type="clustered", clusters=d["muni"].astype("category"))177    tab = pd.DataFrame({"coef": res.params, "se": res.std_errors,178                        "tstat": res.tstats})179    return tab, {"n": int(res.nobs), "r2": float(res.rsquared),180                 "mean_dep": float(d["abs_dev"].mean())}181182183# ---------------------------------------------------------------- tax shift184def tax_shift(df: pd.DataFrame, n_bins: int = 10) -> pd.DataFrame:185    """Implied property-tax shift from differential assessment.186187    Within a taxing cell the levy is proportional to AV, so a property whose188    ratio exceeds the cell median by x% pays x% more tax than under uniform189    assessment. We compute rel_i = r_i / med_c(r) − 1 and average it by190    within-cell sale-price decile.191    """192    d = df.copy()193    d["rel"] = d["ratio"] / d.groupby("cell")["ratio"].transform("median") - 1.0194    d["decile"] = (d.groupby("cell")["amount"]195                   .rank(pct=True)196                   .mul(n_bins).add(1 - 1e-9).astype(int).clip(1, n_bins))197    out = (d.groupby("decile")198             .agg(n=("rel", "size"), mean_rel=("rel", "mean"),199                  median_rel=("rel", "median"),200                  se=("rel", lambda s: s.std() / np.sqrt(len(s))))201             .reset_index())202    return out203