#!/usr/bin/env python3 # Author: Simon-Pierre Boucher — contact@spboucher.ai """Step 06 — Supplementary analyses added at the 2026-10-17 scientific revision. Every number cited in the revised manuscript that is not produced by steps 02–03 comes from this script. Outputs (``results/reproduced/rev_*.csv``) feed ``scripts/05_make_tables.py`` (tables ``results/tables/rev_*.tex``). rev_bounds.csv within-cell moments, Kochin–Parks reverse-regression bound, Clapp first stage, implied noise shares rev_spline.csv piecewise-linear (Sunderman et al. 1990) slopes by within-cell price quintile, cell FE, clustered SEs rev_taxshift_alt.csv excess burden by decile under alternative sorting variables and the revenue-neutral benchmark rev_cod.csv municipal COD under three time-control schemes rev_prb_blocks.csv municipality-year PRB by roll-lag band and by year rev_cities.csv Cheng FE / Clapp IV for each of the ten largest markets rev_trim.csv headline elasticities and tax shift after IAAO-style 1.5 x IQR outlier trimming within cells rev_clusters.csv number of municipality clusters behind each subgroup rev_trim_grid.csv full battery under five outlier-trimming rules rev_het_trimmed.csv subgroup elasticities after IAAO 1.5 x IQR trimming rev_cities_trimmed.csv ten largest markets after IAAO 1.5 x IQR trimming rev_outlier_by_decile.csv share of IQR outliers by within-cell price decile rev_prb_blocks_trimmed.csv / rev_iaao_year_trimmed.csv PRB by block and year, trimmed sample Usage: python scripts/06_revision_analyses.py (≈ 5 min) """ import sys import time from pathlib import Path import numpy as np import pandas as pd import statsmodels.api as sm from linearmodels.iv import IV2SLS from linearmodels.iv.absorbing import AbsorbingLS sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src")) from wp10 import config, iaao, models, sample # noqa: E402 OUT = config.REPRODUCED def _cells_ok(d: pd.DataFrame) -> pd.DataFrame: c = d.groupby("cell")["cell"].transform("size") return d[c >= config.CELL_MIN_OBS] def bounds(df: pd.DataFrame) -> pd.DataFrame: d = models._demean(df, ["ln_av", "ln_price"]) x, y = d["ln_price_w"], d["ln_av_w"] vx, vy, cxy = x.var(), y.var(), np.cov(x, y)[0, 1] beta_f = cxy / vx # forward (Cheng) within slope beta_r = cxy / vy # reverse within slope rows = [("var_lnP_within", vx), ("var_lnAV_within", vy), ("cov_within", cxy), ("var_lnP_total", df["ln_price"].var()), ("beta_forward", beta_f), ("beta_reverse_inverse", 1 / beta_r), ("corr_within", cxy / np.sqrt(vx * vy)), ("iqr_lnP_within", x.quantile(.75) - x.quantile(.25)), # noise share lambda = var(u)/var(lnP_w) that would rationalise beta = 1 ("lambda_for_beta_1", 1 - beta_f), ("noise_sd_for_beta_1", np.sqrt((1 - beta_f) * vx))] # Clapp IV first stage (clustered by municipality) 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 fs = sm.OLS(x.to_numpy(), sm.add_constant(z)).fit( cov_type="cluster", cov_kwds={"groups": d["muni"].to_numpy()}) rows += [("first_stage_coef", fs.params[1]), ("first_stage_se", fs.bse[1]), ("first_stage_F_cluster", float(fs.tvalues[1] ** 2)), ("first_stage_r2", fs.rsquared), ("share_z_minus1", float((z == -1).mean())), ("share_z_plus1", float((z == 1).mean()))] iv = IV2SLS(y, None, d[["ln_price_w"]], d[["z"]]).fit( cov_type="clustered", clusters=d["muni"].astype("category")) b_iv = float(iv.params["ln_price_w"]) rows += [("beta_iv", b_iv), ("se_iv", float(iv.std_errors["ln_price_w"])), ("lambda_implied_by_iv", 1 - beta_f / b_iv), ("noise_sd_implied_by_iv", np.sqrt(max(0, 1 - beta_f / b_iv) * vx))] # pooled (no FE) counterpart of the bound vxp, vyp = df["ln_price"].var(), df["ln_av"].var() cp = np.cov(df["ln_price"], df["ln_av"])[0, 1] rows += [("beta_forward_pooled", cp / vxp), ("beta_reverse_inverse_pooled", vyp / cp)] return pd.DataFrame(rows, columns=["quantity", "value"]) def spline(df: pd.DataFrame, n_seg: int = 5) -> pd.DataFrame: """Continuous piecewise-linear within-cell regression; knots at the pooled quintiles of the within-cell demeaned log price.""" d = models._demean(df, ["ln_price"]) x = d["ln_price_w"] knots = [x.quantile(q) for q in np.linspace(0, 1, n_seg + 1)[1:-1]] X = pd.DataFrame({"x": x}) for j, k in enumerate(knots): X[f"h{j+1}"] = np.clip(x - k, 0, None) mod = AbsorbingLS(d["ln_av"], X, absorb=pd.DataFrame({"cell": d["cell"].astype("category")})) res = mod.fit(cov_type="clustered", clusters=d["muni"].astype("category")) b, V = res.params.to_numpy(), res.cov.to_numpy() rows = [] bounds_lo = [x.min()] + knots bounds_hi = knots + [x.max()] for j in range(n_seg): a = np.zeros(len(b)); a[: j + 1] = 1.0 # slope in segment j = b0 + h1 + … + hj slope = float(a @ b); se = float(np.sqrt(a @ V @ a)) seg = x.between(bounds_lo[j], bounds_hi[j], inclusive="left" if j < n_seg - 1 else "both") rows.append({"segment": j + 1, "lo": bounds_lo[j], "hi": bounds_hi[j], "beta": slope, "se": se, "gamma": slope - 1, "n": int(seg.sum())}) out = pd.DataFrame(rows) out.attrs["n"] = int(res.nobs) return out def taxshift_alt(df: pd.DataFrame, n_bins: int = 10) -> pd.DataFrame: d = df.copy() med = d.groupby("cell")["ratio"].transform("median") agg = (d.groupby("cell")["role_valeur_immeuble"].transform("sum") / d.groupby("cell")["amount"].transform("sum")) d["rel_med"] = d["ratio"] / med - 1.0 d["rel_agg"] = d["ratio"] / agg - 1.0 d["proxy"] = 0.5 * d["amount"] + 0.5 * d["role_valeur_immeuble"] / med sorts = {"sale price": "amount", "IAAO blended value proxy": "proxy", "assessed value": "role_valeur_immeuble"} rows = [] for label, col in sorts.items(): dec = (d.groupby("cell")[col].rank(pct=True).mul(n_bins) .add(1 - 1e-9).astype(int).clip(1, n_bins)) for bench, rel in (("cell median ratio", "rel_med"), ("revenue-neutral (aggregate ratio)", "rel_agg")): g = d.groupby(dec)[rel] t = pd.DataFrame({"decile": g.median().index, "median_rel": g.median().to_numpy(), "mean_rel": g.mean().to_numpy(), "n": g.size().to_numpy()}) t["sort"] = label; t["benchmark"] = bench rows.append(t) return pd.concat(rows, ignore_index=True) def cod_robust(df: pd.DataFrame) -> pd.DataFrame: """Municipal COD (median across blocks, munis with ≥100 sales) under (a) municipality × year blocks (baseline), (b) municipality × quarter blocks, (c) municipality × year blocks after a province-wide monthly time adjustment within roll vintage.""" d = df.copy() d["quarter"] = d["sale_date"].dt.year.astype(str) + "Q" + d["sale_date"].dt.quarter.astype(str) d["month"] = d["sale_date"].dt.to_period("M").astype(str) m_month = d.groupby(["roll", "month"])["ratio"].transform("median") m_year = d.groupby(["roll", "sale_year"])["ratio"].transform("median") d["ratio_adj"] = d["ratio"] * m_year / m_month d["av_adj"] = d["ratio_adj"] * d["amount"] def muni_table(data, by, av_col, min_n): data = data.assign(role_valeur_immeuble=data[av_col]) blocks = iaao.group_metrics(data, by, min_n=min_n) blocks["muni"] = blocks["group"].str.split("_").str[0] agg = blocks.groupby("muni").agg(n=("n", "sum"), cod=("cod", "median"), prb=("prb", "median")).reset_index() return agg[agg["n"] >= config.MUNI_MIN_SALES] schemes = { "Municipality x year blocks (baseline)": (d, ["muni", "sale_year"], "role_valeur_immeuble", 50), "Municipality x quarter blocks": (d, ["muni", "quarter"], "role_valeur_immeuble", 30), "Municipality x year, time-adjusted ratios": (d, ["muni", "sale_year"], "av_adj", 50), } rows = [] for label, (data, by, col, mn) in schemes.items(): t = muni_table(data, by, col, mn) rows.append({"scheme": label, "n_munis": len(t), "median_cod": t["cod"].median(), "p25_cod": t["cod"].quantile(.25), "p75_cod": t["cod"].quantile(.75), "share_cod_gt_15": (t["cod"] > 15).mean(), "share_cod_gt_20": (t["cod"] > 20).mean(), "median_prb": t["prb"].median(), "share_prb_neg": (t["prb"] < 0).mean(), "share_prb_lt_m05": (t["prb"] < -0.05).mean()}) # province-wide COD on time-adjusted ratios r = d["ratio_adj"].to_numpy() rows.append({"scheme": "Province-wide, time-adjusted ratios", "n_munis": np.nan, "median_cod": iaao.cod(r), "p25_cod": np.nan, "p75_cod": np.nan, "share_cod_gt_15": np.nan, "share_cod_gt_20": np.nan, "median_prb": iaao.prb(d["av_adj"].to_numpy(), d["amount"].to_numpy())[0], "share_prb_neg": np.nan, "share_prb_lt_m05": np.nan}) return pd.DataFrame(rows) def prb_blocks(df: pd.DataFrame) -> pd.DataFrame: blocks = iaao.group_metrics(df, ["muni", "sale_year"], min_n=50) lag = df.groupby(["muni", "sale_year"])["lag_months"].median() lag.index = [f"{m}_{y}" for m, y in lag.index] blocks["lag"] = blocks["group"].map(lag) blocks["year"] = blocks["group"].str.split("_").str[1].astype(int) blocks["band"] = pd.cut(blocks["lag"], [0, 24, 48, 200], labels=["Roll lag < 24 m", "Roll lag 24–48 m", "Roll lag > 48 m"]) rows = [] for key, g in list(blocks.groupby("band", observed=True)) + list(blocks.groupby("year")): rows.append({"group": str(key), "n_blocks": len(g), "n_sales": int(g["n"].sum()), "median_prb": g["prb"].median(), "mean_prb": g["prb"].mean(), "share_prb_neg": (g["prb"] < 0).mean(), "median_ratio": g["median_ratio"].median(), "median_cod": g["cod"].median()}) return pd.DataFrame(rows) def cities(df: pd.DataFrame) -> pd.DataFrame: rows = [] for city in config.BIG_CITIES: sub = _cells_ok(df[df["role_municipalite"] == city]) if len(sub) < 2000: continue dep, exog = sub["ln_av"], sub[["ln_price"]] absorb = pd.DataFrame({"cell": sub["cell"].astype("category")}) fe = AbsorbingLS(dep, exog, absorb=absorb).fit(cov_type="robust") d = models._demean(sub, ["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; z[(g_av > 2/3) & (g_sp > 2/3)] = 1 d["z"] = z iv = IV2SLS(d["ln_av_w"], None, d[["ln_price_w"]], d[["z"]]).fit(cov_type="robust") av, sp = sub["role_valeur_immeuble"].to_numpy(float), sub["amount"].to_numpy(float) share_condo = float((sub["prop_class"] == "condo").mean()) rows.append({"city": city, "n": len(sub), "n_cells": sub["cell"].nunique(), "gamma_fe": float(fe.params["ln_price"]) - 1, "se_fe": float(fe.std_errors["ln_price"]), "gamma_iv": float(iv.params["ln_price_w"]) - 1, "se_iv": float(iv.std_errors["ln_price_w"]), "share_condo": share_condo}) return pd.DataFrame(rows) def trim_iqr(df: pd.DataFrame) -> tuple[pd.DataFrame, pd.DataFrame]: q1 = df.groupby("cell")["ratio"].transform(lambda s: s.quantile(.25)) q3 = df.groupby("cell")["ratio"].transform(lambda s: s.quantile(.75)) iqr = q3 - q1 keep = df["ratio"].between(q1 - 1.5 * iqr, q3 + 1.5 * iqr) sub = _cells_ok(df[keep]) fe, iv = models.cheng_fe(sub), models.clapp_iv(sub) ts = models.tax_shift(sub) av, sp = sub["role_valeur_immeuble"].to_numpy(float), sub["amount"].to_numpy(float) head = pd.DataFrame([{"variant": "IAAO 1.5 x IQR trim within cells", "n": len(sub), "share_dropped": 1 - keep.mean(), "gamma_fe": fe["gamma"], "se_fe": fe["se"], "gamma_iv": iv["gamma"], "se_iv": iv["se"], "cod": iaao.cod(av / sp), "prd": iaao.prd(av, sp), "prb": iaao.prb(av, sp)[0], "median_rel_d1": float(ts.loc[ts.decile == 1, "median_rel"].iloc[0]), "median_rel_d2": float(ts.loc[ts.decile == 2, "median_rel"].iloc[0]), "median_rel_d10": float(ts.loc[ts.decile == 10, "median_rel"].iloc[0])}]) return head, ts def clusters(df: pd.DataFrame) -> pd.DataFrame: muni_sales = df.groupby("muni")["muni"].transform("size") groups = { "Single-family": df["prop_class"] == "single_family", "Condominium": df["prop_class"] == "condo", "Plex (2–5 units)": df["prop_class"] == "plex", "Cottage": df["prop_class"] == "cottage", "Muni < 1k sales": muni_sales < 1_000, "Muni 1k–10k sales": muni_sales.between(1_000, 10_000), "Muni > 10k sales": muni_sales > 10_000, } rows = [] for k, m in groups.items(): sub = _cells_ok(df[m]) rows.append({"group": k, "n": len(sub), "n_clusters": sub["muni"].nunique(), "n_cells": sub["cell"].nunique()}) rows.append({"group": "Baseline", "n": len(df), "n_clusters": df["muni"].nunique(), "n_cells": df["cell"].nunique()}) return pd.DataFrame(rows) # ---------------------------------------------------------------- trimming grid def _trim_mask(df: pd.DataFrame, rule: str) -> pd.Series: r = df["ratio"] g = df.groupby("cell")["ratio"] if rule == "baseline": return pd.Series(True, index=df.index) if rule.endswith("iqr"): k = float(rule.replace("iqr", "")) q1, q3 = g.transform(lambda s: s.quantile(.25)), g.transform(lambda s: s.quantile(.75)) return r.between(q1 - k * (q3 - q1), q3 + k * (q3 - q1)) if rule.endswith("pct"): a = float(rule.replace("pct", "")) / 100 lo, hi = g.transform(lambda s: s.quantile(a)), g.transform(lambda s: s.quantile(1 - a)) return r.between(lo, hi) raise ValueError(rule) RULES = {"baseline": "Baseline (1/99 pct by roll vintage)", "3iqr": "Within-cell 3 x IQR (extreme outliers)", "1.5iqr": "Within-cell 1.5 x IQR (IAAO outlier rule)", "2.5pct": "Within-cell 2.5/97.5 pct", "5pct": "Within-cell 5/95 pct"} def muni_shares(d: pd.DataFrame) -> dict: blocks = iaao.group_metrics(d, ["muni", "sale_year"], min_n=50) blocks["muni"] = blocks["group"].str.rsplit("_", n=1).str[0] agg = blocks.groupby("muni").agg(n=("n", "sum"), cod=("cod", "median"), prb=("prb", "median")).reset_index() agg = agg[agg["n"] >= config.MUNI_MIN_SALES] return {"n_munis": len(agg), "median_muni_cod": agg["cod"].median(), "share_muni_cod_gt_15": (agg["cod"] > 15).mean(), "share_muni_prb_neg": (agg["prb"] < 0).mean(), "share_muni_prb_lt_m05": (agg["prb"] < -0.05).mean(), "share_muni_prb_lt_m10": (agg["prb"] < -0.10).mean(), "median_muni_prb": agg["prb"].median()} def trim_grid(df: pd.DataFrame) -> pd.DataFrame: rows = [] for rule, label in RULES.items(): keep = _trim_mask(df, rule) sub = _cells_ok(df[keep]) fe, iv = models.cheng_fe(sub), models.clapp_iv(sub) ts = models.tax_shift(sub) av, sp = sub["role_valeur_immeuble"].to_numpy(float), sub["amount"].to_numpy(float) sp_ = spline(sub) row = {"rule": rule, "label": label, "n": len(sub), "share_dropped": 1 - keep.mean(), "median_ratio": float(np.median(av / sp)), "cod": iaao.cod(av / sp), "prd": iaao.prd(av, sp), "prb": iaao.prb(av, sp)[0], "prb_se": iaao.prb(av, sp)[1], "gamma_fe": fe["gamma"], "se_fe": fe["se"], "gamma_iv": iv["gamma"], "se_iv": iv["se"], "d1_median": float(ts.loc[ts.decile == 1, "median_rel"].iloc[0]), "d10_median": float(ts.loc[ts.decile == 10, "median_rel"].iloc[0]), "spline_seg1_beta": float(sp_.loc[0, "beta"]), "spline_seg5_beta": float(sp_.loc[4, "beta"])} row.update(muni_shares(sub)) rows.append(row) print(f" {label:<44} drop={row['share_dropped']:.1%} COD={row['cod']:.1f} PRB={row['prb']:+.3f} " f"gFE={row['gamma_fe']:+.3f} gIV={row['gamma_iv']:+.3f} d1={row['d1_median']:+.1%} " f"muniPRB<-.05={row['share_muni_prb_lt_m05']:.0%} muniCOD>15={row['share_muni_cod_gt_15']:.0%}") return pd.DataFrame(rows) def het_trimmed(df: pd.DataFrame) -> pd.DataFrame: sub = df[_trim_mask(df, "1.5iqr")] muni_sales = sub.groupby("muni")["muni"].transform("size") groups = { "Single-family": sub["prop_class"] == "single_family", "Condominium": sub["prop_class"] == "condo", "Plex (2–5 units)": sub["prop_class"] == "plex", "Age < 20 y": sub["age"] < 20, "Age > 60 y": sub["age"] > 60, "Land share < 0.2": sub["land_share"] < 0.2, "Land share > 0.4": sub["land_share"] > 0.4, "Muni < 1k sales": muni_sales < 1_000, "Muni > 10k sales": muni_sales > 10_000, "Sales 2021–2023": sub["sale_year"] <= 2023, "Sales 2024–2026": sub["sale_year"] >= 2024, } out = models.gamma_by_group(sub, groups) ivs = [] for label, mask in groups.items(): s2 = _cells_ok(sub[mask]) ivs.append(models.clapp_iv(s2)["gamma"] if len(s2) >= 2000 else np.nan) out["gamma_iv"] = ivs return out def cities_trimmed(df: pd.DataFrame) -> pd.DataFrame: sub = df[_trim_mask(df, "1.5iqr")] big = sub[sub["role_municipalite"].isin(config.BIG_CITIES)].copy() big["muni"] = big["role_municipalite"] blocks = iaao.group_metrics(big, ["muni", "sale_year"], min_n=200) blocks["muni"] = blocks["group"].str.rsplit("_", n=1).str[0] agg = blocks.groupby("muni").agg(n=("n", "sum"), median_ratio=("median_ratio", "median"), cod=("cod", "median"), prd=("prd", "median"), prb=("prb", "median")).reset_index() share_neg = blocks.assign(neg=blocks["prb"] < 0).groupby("muni")["neg"].mean().rename("share_years_prb_neg") agg = agg.merge(share_neg, on="muni") g = cities(sub).rename(columns={"city": "muni"}) return agg.merge(g[["muni", "gamma_fe", "se_fe", "gamma_iv", "se_iv", "share_condo"]], on="muni") def main() -> None: config.ensure_dirs() t0 = time.time() df = sample.load() print(f"sample {len(df):,}; share lag > 54 m: {(df['lag_months'] > 54).mean():.3%}; " f"muni-years with 2 rolls: {(df.groupby(['muni','sale_year'])['roll'].nunique() > 1).sum()}") pd.Series({"share_lag_gt_54": (df["lag_months"] > 54).mean(), "share_lag_gt_60": (df["lag_months"] > 60).mean(), "max_lag": df["lag_months"].max(), "muni_years": df.groupby(["muni", "sale_year"]).ngroups, "muni_years_two_rolls": int((df.groupby(["muni", "sale_year"])["roll"].nunique() > 1).sum()), "share_ratio_gt_1_5": (df["ratio"] > 1.5).mean(), "share_ratio_gt_2": (df["ratio"] > 2).mean()}).to_csv(OUT / "rev_misc.csv") b = bounds(df); b.to_csv(OUT / "rev_bounds.csv", index=False); print(b.to_string()) s = spline(df); s.to_csv(OUT / "rev_spline.csv", index=False); print(s.to_string()) t = taxshift_alt(df); t.to_csv(OUT / "rev_taxshift_alt.csv", index=False) print(t[t.decile.isin([1, 2, 10])].to_string()) c = cod_robust(df); c.to_csv(OUT / "rev_cod.csv", index=False); print(c.to_string()) p = prb_blocks(df); p.to_csv(OUT / "rev_prb_blocks.csv", index=False); print(p.to_string()) ci = cities(df); ci.to_csv(OUT / "rev_cities.csv", index=False); print(ci.to_string()) head, ts = trim_iqr(df); head.to_csv(OUT / "rev_trim.csv", index=False) ts.to_csv(OUT / "rev_trim_taxshift.csv", index=False); print(head.T.to_string()) cl = clusters(df); cl.to_csv(OUT / "rev_clusters.csv", index=False); print(cl.to_string()) print("Trimming grid") tg = trim_grid(df); tg.to_csv(OUT / "rev_trim_grid.csv", index=False) ht = het_trimmed(df); ht.to_csv(OUT / "rev_het_trimmed.csv", index=False); print(ht.to_string()) ct = cities_trimmed(df); ct.to_csv(OUT / "rev_cities_trimmed.csv", index=False); print(ct.to_string()) # share of decile-1 sales that are ratio outliers d = df.copy() d["decile"] = d.groupby("cell")["amount"].rank(pct=True).mul(10).add(1 - 1e-9).astype(int).clip(1, 10) d["outlier"] = ~_trim_mask(df, "1.5iqr") (d.groupby("decile")["outlier"].mean().rename("share_outlier").reset_index() .to_csv(OUT / "rev_outlier_by_decile.csv", index=False)) print(d.groupby("decile")["outlier"].mean().round(3).to_dict()) pt = prb_blocks(df[_trim_mask(df, "1.5iqr")]) pt.to_csv(OUT / "rev_prb_blocks_trimmed.csv", index=False); print(pt.to_string()) # province-wide PRB by year on the trimmed sample tr = df[_trim_mask(df, "1.5iqr")] rows = [{"year": "All", "prb": iaao.prb(tr["role_valeur_immeuble"].to_numpy(float), tr["amount"].to_numpy(float))[0], "cod": iaao.cod((tr["role_valeur_immeuble"] / tr["amount"]).to_numpy()), "n": len(tr)}] for y, g in tr.groupby("sale_year"): rows.append({"year": str(y), "prb": iaao.prb(g["role_valeur_immeuble"].to_numpy(float), g["amount"].to_numpy(float))[0], "cod": iaao.cod((g["role_valeur_immeuble"] / g["amount"]).to_numpy()), "n": len(g)}) pd.DataFrame(rows).to_csv(OUT / "rev_iaao_year_trimmed.csv", index=False); print(pd.DataFrame(rows).to_string()) print(f"done in {time.time() - t0:.0f} s") if __name__ == "__main__": main()