# Author: Simon-Pierre Boucher — contact@spboucher.ai """Econometric estimators for WP10. Every vertical-inequity estimator ultimately measures how the assessment ratio r = AV/SP moves with market value. The workhorse is the log-log (Cheng 1974) regression ln AV_i = α_c(i) + β · ln SP_i + ε_i , where α_c(i) is a municipality × roll vintage × sale-year fixed effect that absorbs the mechanical drift of ratios between triennial reference dates. β < 1 ⇒ regressive assessment. We report γ ≡ β − 1, the elasticity of the assessment *ratio* with respect to price (γ < 0 ⇒ regressive), estimated by - pooled OLS (no fixed effects), - absorbing least squares with cell fixed effects, - Clapp's (1990) rank-based IV, which purges the attenuation/mean-reversion bias that pushes OLS toward spurious regressivity when sale prices carry idiosyncratic noise, - quantile regression across the conditional distribution. Inference is clustered by municipality throughout. """ from __future__ import annotations import numpy as np import pandas as pd import statsmodels.api as sm from linearmodels.iv import IV2SLS from linearmodels.iv.absorbing import AbsorbingLS from . import config # ---------------------------------------------------------------- helpers def _demean(df: pd.DataFrame, cols: list[str], by: str = "cell") -> pd.DataFrame: """Within-transform ``cols`` inside groups of ``by`` (suffix ``_w``).""" out = df.copy() for c in cols: out[c + "_w"] = df[c] - df.groupby(by)[c].transform("mean") return out def _cluster_ols(y: np.ndarray, X: np.ndarray, clusters: np.ndarray, dof_adjust: int = 0): """OLS with CRVE (cluster-robust) inference via statsmodels.""" model = sm.OLS(y, X) return model.fit(cov_type="cluster", cov_kwds={"groups": clusters}) # ---------------------------------------------------------------- Cheng / FE def cheng_pooled(df: pd.DataFrame) -> dict: """Pooled Cheng regression: ln AV on ln SP, no fixed effects.""" X = sm.add_constant(df["ln_price"].to_numpy()) res = _cluster_ols(df["ln_av"].to_numpy(), X, df["muni"].to_numpy()) return {"estimator": "Cheng OLS (pooled)", "beta": res.params[1], "se": res.bse[1], "gamma": res.params[1] - 1.0, "n": int(res.nobs), "r2": res.rsquared} def cheng_fe(df: pd.DataFrame, extra_controls: list[str] | None = None) -> dict: """Cheng regression absorbing municipality × roll × sale-year cells.""" dep = df["ln_av"] exog_cols = ["ln_price"] + (extra_controls or []) exog = df[exog_cols] absorb = pd.DataFrame({"cell": df["cell"].astype("category")}) mod = AbsorbingLS(dep, exog, absorb=absorb) res = mod.fit(cov_type="clustered", clusters=df["muni"].astype("category")) return {"estimator": "Cheng FE (cell)", "beta": float(res.params["ln_price"]), "se": float(res.std_errors["ln_price"]), "gamma": float(res.params["ln_price"]) - 1.0, "n": int(res.nobs), "r2": float(res.rsquared)} def paglin_fogarty(df: pd.DataFrame) -> dict: """Levels regression AV = a + b·SP (Paglin & Fogarty 1972). A positive intercept with b below the overall assessment level indicates regressivity in levels. """ X = sm.add_constant(df["amount"].to_numpy(float)) res = _cluster_ols(df["role_valeur_immeuble"].to_numpy(float), X, df["muni"].to_numpy()) return {"estimator": "Paglin–Fogarty (levels)", "intercept": res.params[0], "intercept_se": res.bse[0], "slope": res.params[1], "slope_se": res.bse[1], "n": int(res.nobs), "r2": res.rsquared} # ---------------------------------------------------------------- Clapp IV def clapp_iv(df: pd.DataFrame) -> dict: """Clapp (1990) rank-based IV on within-cell demeaned data. The instrument z ∈ {−1, 0, +1} flags sales in the bottom/top third of *both* the ln AV and ln SP within-cell distributions. Because z carries only coarse rank information, it is (near-)orthogonal to the transitory component of either variable, undoing the attenuation that biases OLS toward regressivity (β̂ < β) when prices are noisy. """ d = _demean(df, ["ln_av", "ln_price"]) g_av = d.groupby("cell")["ln_av"].rank(pct=True) g_sp = d.groupby("cell")["ln_price"].rank(pct=True) z = np.zeros(len(d)) z[(g_av <= 1 / 3) & (g_sp <= 1 / 3)] = -1.0 z[(g_av > 2 / 3) & (g_sp > 2 / 3)] = 1.0 d["z"] = z res = IV2SLS(d["ln_av_w"], None, d[["ln_price_w"]], d[["z"]]).fit( cov_type="clustered", clusters=d["muni"].astype("category")) beta = float(res.params["ln_price_w"]) return {"estimator": "Clapp IV (rank instrument)", "beta": beta, "se": float(res.std_errors["ln_price_w"]), "gamma": beta - 1.0, "n": int(res.nobs), "r2": float(res.rsquared)} # ---------------------------------------------------------------- quantiles def quantile_betas(df: pd.DataFrame, taus: list[float] = config.QUANTILES, max_n: int = 250_000, seed: int = config.SEED_BOOT) -> pd.DataFrame: """Quantile regressions of demeaned ln AV on demeaned ln SP. Estimated on a seeded random subsample (IRLS on the full 600k sample is needlessly slow; the subsample SEs are already microscopic). """ d = _demean(df, ["ln_av", "ln_price"]) if len(d) > max_n: d = d.sample(max_n, random_state=seed) X = sm.add_constant(d["ln_price_w"].to_numpy()) y = d["ln_av_w"].to_numpy() rows = [] for tau in taus: r = sm.QuantReg(y, X).fit(q=tau, max_iter=2000) rows.append({"tau": tau, "beta": r.params[1], "se": r.bse[1], "gamma": r.params[1] - 1.0, "n": len(d)}) return pd.DataFrame(rows) # ---------------------------------------------------------------- subgroups def gamma_by_group(df: pd.DataFrame, groups: dict[str, pd.Series]) -> pd.DataFrame: """Cheng-FE γ estimated separately on each labelled subsample. ``groups`` maps a label to a boolean mask over ``df``. Subsamples keep only cells that retain ≥ CELL_MIN_OBS sales after masking. """ rows = [] for label, mask in groups.items(): sub = df[mask] counts = sub.groupby("cell")["cell"].transform("size") sub = sub[counts >= config.CELL_MIN_OBS] if len(sub) < 2_000: continue est = cheng_fe(sub) rows.append({"group": label, "gamma": est["gamma"], "se": est["se"], "n": est["n"]}) return pd.DataFrame(rows) # ---------------------------------------------------------------- horizontal def horizontal_dispersion(df: pd.DataFrame) -> tuple[pd.DataFrame, dict]: """Horizontal inequity: who gets the noisiest assessments? Regresses the absolute log deviation of a sale's ratio from its cell median, |ln r_i − med_c ln r|, on property characteristics with cell fixed effects. Positive coefficients = less uniform assessment. Returns (coefficient table, fit metadata). """ d = df.copy() d["abs_dev"] = (d["ln_ratio"] - d.groupby("cell")["ln_ratio"].transform("median")).abs() d["age_dec"] = d["age"] / 10.0 d["is_condo"] = (d["prop_class"] == "condo").astype(float) d["is_plex"] = (d["prop_class"] == "plex").astype(float) d["is_cottage"] = (d["prop_class"] == "cottage").astype(float) d = d.dropna(subset=["abs_dev", "age_dec", "land_share"]) exog_cols = ["age_dec", "land_share", "is_condo", "is_plex", "is_cottage"] mod = AbsorbingLS(d["abs_dev"], d[exog_cols], absorb=pd.DataFrame({"cell": d["cell"].astype("category")})) res = mod.fit(cov_type="clustered", clusters=d["muni"].astype("category")) tab = pd.DataFrame({"coef": res.params, "se": res.std_errors, "tstat": res.tstats}) return tab, {"n": int(res.nobs), "r2": float(res.rsquared), "mean_dep": float(d["abs_dev"].mean())} # ---------------------------------------------------------------- tax shift def tax_shift(df: pd.DataFrame, n_bins: int = 10) -> pd.DataFrame: """Implied property-tax shift from differential assessment. Within a taxing cell the levy is proportional to AV, so a property whose ratio exceeds the cell median by x% pays x% more tax than under uniform assessment. We compute rel_i = r_i / med_c(r) − 1 and average it by within-cell sale-price decile. """ d = df.copy() d["rel"] = d["ratio"] / d.groupby("cell")["ratio"].transform("median") - 1.0 d["decile"] = (d.groupby("cell")["amount"] .rank(pct=True) .mul(n_bins).add(1 - 1e-9).astype(int).clip(1, n_bins)) out = (d.groupby("decile") .agg(n=("rel", "size"), mean_rel=("rel", "mean"), median_rel=("rel", "median"), se=("rel", lambda s: s.std() / np.sqrt(len(s)))) .reset_index()) return out