UQO Working Paper No. 10 — The assessment gap in Quebec: vertical and horizontal inequity in municipal property assessment.
1#!/usr/bin/env python32# Author: Simon-Pierre Boucher — contact@spboucher.ai3"""Step 06 — Supplementary analyses added at the 2026-10-17 scientific revision.45Every number cited in the revised manuscript that is not produced by steps602–03 comes from this script. Outputs (``results/reproduced/rev_*.csv``) feed7``scripts/05_make_tables.py`` (tables ``results/tables/rev_*.tex``).89 rev_bounds.csv within-cell moments, Kochin–Parks reverse-regression10 bound, Clapp first stage, implied noise shares11 rev_spline.csv piecewise-linear (Sunderman et al. 1990) slopes by12 within-cell price quintile, cell FE, clustered SEs13 rev_taxshift_alt.csv excess burden by decile under alternative sorting14 variables and the revenue-neutral benchmark15 rev_cod.csv municipal COD under three time-control schemes16 rev_prb_blocks.csv municipality-year PRB by roll-lag band and by year17 rev_cities.csv Cheng FE / Clapp IV for each of the ten largest markets18 rev_trim.csv headline elasticities and tax shift after IAAO-style19 1.5 x IQR outlier trimming within cells20 rev_clusters.csv number of municipality clusters behind each subgroup21 rev_trim_grid.csv full battery under five outlier-trimming rules22 rev_het_trimmed.csv subgroup elasticities after IAAO 1.5 x IQR trimming23 rev_cities_trimmed.csv ten largest markets after IAAO 1.5 x IQR trimming24 rev_outlier_by_decile.csv share of IQR outliers by within-cell price decile25 rev_prb_blocks_trimmed.csv / rev_iaao_year_trimmed.csv PRB by block and year, trimmed sample2627Usage: python scripts/06_revision_analyses.py (≈ 5 min)28"""29import sys30import time31from pathlib import Path3233import numpy as np34import pandas as pd35import statsmodels.api as sm36from linearmodels.iv import IV2SLS37from linearmodels.iv.absorbing import AbsorbingLS3839sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src"))4041from wp10 import config, iaao, models, sample # noqa: E4024243OUT = config.REPRODUCED444546def _cells_ok(d: pd.DataFrame) -> pd.DataFrame:47 c = d.groupby("cell")["cell"].transform("size")48 return d[c >= config.CELL_MIN_OBS]495051def bounds(df: pd.DataFrame) -> pd.DataFrame:52 d = models._demean(df, ["ln_av", "ln_price"])53 x, y = d["ln_price_w"], d["ln_av_w"]54 vx, vy, cxy = x.var(), y.var(), np.cov(x, y)[0, 1]55 beta_f = cxy / vx # forward (Cheng) within slope56 beta_r = cxy / vy # reverse within slope57 rows = [("var_lnP_within", vx), ("var_lnAV_within", vy), ("cov_within", cxy),58 ("var_lnP_total", df["ln_price"].var()),59 ("beta_forward", beta_f), ("beta_reverse_inverse", 1 / beta_r),60 ("corr_within", cxy / np.sqrt(vx * vy)),61 ("iqr_lnP_within", x.quantile(.75) - x.quantile(.25)),62 # noise share lambda = var(u)/var(lnP_w) that would rationalise beta = 163 ("lambda_for_beta_1", 1 - beta_f),64 ("noise_sd_for_beta_1", np.sqrt((1 - beta_f) * vx))]65 # Clapp IV first stage (clustered by municipality)66 g_av = d.groupby("cell")["ln_av"].rank(pct=True)67 g_sp = d.groupby("cell")["ln_price"].rank(pct=True)68 z = np.zeros(len(d))69 z[(g_av <= 1 / 3) & (g_sp <= 1 / 3)] = -1.070 z[(g_av > 2 / 3) & (g_sp > 2 / 3)] = 1.071 d["z"] = z72 fs = sm.OLS(x.to_numpy(), sm.add_constant(z)).fit(73 cov_type="cluster", cov_kwds={"groups": d["muni"].to_numpy()})74 rows += [("first_stage_coef", fs.params[1]), ("first_stage_se", fs.bse[1]),75 ("first_stage_F_cluster", float(fs.tvalues[1] ** 2)),76 ("first_stage_r2", fs.rsquared),77 ("share_z_minus1", float((z == -1).mean())),78 ("share_z_plus1", float((z == 1).mean()))]79 iv = IV2SLS(y, None, d[["ln_price_w"]], d[["z"]]).fit(80 cov_type="clustered", clusters=d["muni"].astype("category"))81 b_iv = float(iv.params["ln_price_w"])82 rows += [("beta_iv", b_iv), ("se_iv", float(iv.std_errors["ln_price_w"])),83 ("lambda_implied_by_iv", 1 - beta_f / b_iv),84 ("noise_sd_implied_by_iv", np.sqrt(max(0, 1 - beta_f / b_iv) * vx))]85 # pooled (no FE) counterpart of the bound86 vxp, vyp = df["ln_price"].var(), df["ln_av"].var()87 cp = np.cov(df["ln_price"], df["ln_av"])[0, 1]88 rows += [("beta_forward_pooled", cp / vxp), ("beta_reverse_inverse_pooled", vyp / cp)]89 return pd.DataFrame(rows, columns=["quantity", "value"])909192def spline(df: pd.DataFrame, n_seg: int = 5) -> pd.DataFrame:93 """Continuous piecewise-linear within-cell regression; knots at the94 pooled quintiles of the within-cell demeaned log price."""95 d = models._demean(df, ["ln_price"])96 x = d["ln_price_w"]97 knots = [x.quantile(q) for q in np.linspace(0, 1, n_seg + 1)[1:-1]]98 X = pd.DataFrame({"x": x})99 for j, k in enumerate(knots):100 X[f"h{j+1}"] = np.clip(x - k, 0, None)101 mod = AbsorbingLS(d["ln_av"], X,102 absorb=pd.DataFrame({"cell": d["cell"].astype("category")}))103 res = mod.fit(cov_type="clustered", clusters=d["muni"].astype("category"))104 b, V = res.params.to_numpy(), res.cov.to_numpy()105 rows = []106 bounds_lo = [x.min()] + knots107 bounds_hi = knots + [x.max()]108 for j in range(n_seg):109 a = np.zeros(len(b)); a[: j + 1] = 1.0 # slope in segment j = b0 + h1 + … + hj110 slope = float(a @ b); se = float(np.sqrt(a @ V @ a))111 seg = x.between(bounds_lo[j], bounds_hi[j], inclusive="left" if j < n_seg - 1 else "both")112 rows.append({"segment": j + 1, "lo": bounds_lo[j], "hi": bounds_hi[j],113 "beta": slope, "se": se, "gamma": slope - 1, "n": int(seg.sum())})114 out = pd.DataFrame(rows)115 out.attrs["n"] = int(res.nobs)116 return out117118119def taxshift_alt(df: pd.DataFrame, n_bins: int = 10) -> pd.DataFrame:120 d = df.copy()121 med = d.groupby("cell")["ratio"].transform("median")122 agg = (d.groupby("cell")["role_valeur_immeuble"].transform("sum")123 / d.groupby("cell")["amount"].transform("sum"))124 d["rel_med"] = d["ratio"] / med - 1.0125 d["rel_agg"] = d["ratio"] / agg - 1.0126 d["proxy"] = 0.5 * d["amount"] + 0.5 * d["role_valeur_immeuble"] / med127 sorts = {"sale price": "amount", "IAAO blended value proxy": "proxy",128 "assessed value": "role_valeur_immeuble"}129 rows = []130 for label, col in sorts.items():131 dec = (d.groupby("cell")[col].rank(pct=True).mul(n_bins)132 .add(1 - 1e-9).astype(int).clip(1, n_bins))133 for bench, rel in (("cell median ratio", "rel_med"),134 ("revenue-neutral (aggregate ratio)", "rel_agg")):135 g = d.groupby(dec)[rel]136 t = pd.DataFrame({"decile": g.median().index, "median_rel": g.median().to_numpy(),137 "mean_rel": g.mean().to_numpy(), "n": g.size().to_numpy()})138 t["sort"] = label; t["benchmark"] = bench139 rows.append(t)140 return pd.concat(rows, ignore_index=True)141142143def cod_robust(df: pd.DataFrame) -> pd.DataFrame:144 """Municipal COD (median across blocks, munis with ≥100 sales) under145 (a) municipality × year blocks (baseline), (b) municipality × quarter146 blocks, (c) municipality × year blocks after a province-wide monthly147 time adjustment within roll vintage."""148 d = df.copy()149 d["quarter"] = d["sale_date"].dt.year.astype(str) + "Q" + d["sale_date"].dt.quarter.astype(str)150 d["month"] = d["sale_date"].dt.to_period("M").astype(str)151 m_month = d.groupby(["roll", "month"])["ratio"].transform("median")152 m_year = d.groupby(["roll", "sale_year"])["ratio"].transform("median")153 d["ratio_adj"] = d["ratio"] * m_year / m_month154 d["av_adj"] = d["ratio_adj"] * d["amount"]155156 def muni_table(data, by, av_col, min_n):157 data = data.assign(role_valeur_immeuble=data[av_col])158 blocks = iaao.group_metrics(data, by, min_n=min_n)159 blocks["muni"] = blocks["group"].str.split("_").str[0]160 agg = blocks.groupby("muni").agg(n=("n", "sum"), cod=("cod", "median"),161 prb=("prb", "median")).reset_index()162 return agg[agg["n"] >= config.MUNI_MIN_SALES]163164 schemes = {165 "Municipality x year blocks (baseline)": (d, ["muni", "sale_year"], "role_valeur_immeuble", 50),166 "Municipality x quarter blocks": (d, ["muni", "quarter"], "role_valeur_immeuble", 30),167 "Municipality x year, time-adjusted ratios": (d, ["muni", "sale_year"], "av_adj", 50),168 }169 rows = []170 for label, (data, by, col, mn) in schemes.items():171 t = muni_table(data, by, col, mn)172 rows.append({"scheme": label, "n_munis": len(t), "median_cod": t["cod"].median(),173 "p25_cod": t["cod"].quantile(.25), "p75_cod": t["cod"].quantile(.75),174 "share_cod_gt_15": (t["cod"] > 15).mean(),175 "share_cod_gt_20": (t["cod"] > 20).mean(),176 "median_prb": t["prb"].median(),177 "share_prb_neg": (t["prb"] < 0).mean(),178 "share_prb_lt_m05": (t["prb"] < -0.05).mean()})179 # province-wide COD on time-adjusted ratios180 r = d["ratio_adj"].to_numpy()181 rows.append({"scheme": "Province-wide, time-adjusted ratios", "n_munis": np.nan,182 "median_cod": iaao.cod(r), "p25_cod": np.nan, "p75_cod": np.nan,183 "share_cod_gt_15": np.nan, "share_cod_gt_20": np.nan,184 "median_prb": iaao.prb(d["av_adj"].to_numpy(), d["amount"].to_numpy())[0],185 "share_prb_neg": np.nan, "share_prb_lt_m05": np.nan})186 return pd.DataFrame(rows)187188189def prb_blocks(df: pd.DataFrame) -> pd.DataFrame:190 blocks = iaao.group_metrics(df, ["muni", "sale_year"], min_n=50)191 lag = df.groupby(["muni", "sale_year"])["lag_months"].median()192 lag.index = [f"{m}_{y}" for m, y in lag.index]193 blocks["lag"] = blocks["group"].map(lag)194 blocks["year"] = blocks["group"].str.split("_").str[1].astype(int)195 blocks["band"] = pd.cut(blocks["lag"], [0, 24, 48, 200],196 labels=["Roll lag < 24 m", "Roll lag 24–48 m", "Roll lag > 48 m"])197 rows = []198 for key, g in list(blocks.groupby("band", observed=True)) + list(blocks.groupby("year")):199 rows.append({"group": str(key), "n_blocks": len(g), "n_sales": int(g["n"].sum()),200 "median_prb": g["prb"].median(), "mean_prb": g["prb"].mean(),201 "share_prb_neg": (g["prb"] < 0).mean(),202 "median_ratio": g["median_ratio"].median(),203 "median_cod": g["cod"].median()})204 return pd.DataFrame(rows)205206207def cities(df: pd.DataFrame) -> pd.DataFrame:208 rows = []209 for city in config.BIG_CITIES:210 sub = _cells_ok(df[df["role_municipalite"] == city])211 if len(sub) < 2000:212 continue213 dep, exog = sub["ln_av"], sub[["ln_price"]]214 absorb = pd.DataFrame({"cell": sub["cell"].astype("category")})215 fe = AbsorbingLS(dep, exog, absorb=absorb).fit(cov_type="robust")216 d = models._demean(sub, ["ln_av", "ln_price"])217 g_av = d.groupby("cell")["ln_av"].rank(pct=True)218 g_sp = d.groupby("cell")["ln_price"].rank(pct=True)219 z = np.zeros(len(d)); z[(g_av <= 1/3) & (g_sp <= 1/3)] = -1; z[(g_av > 2/3) & (g_sp > 2/3)] = 1220 d["z"] = z221 iv = IV2SLS(d["ln_av_w"], None, d[["ln_price_w"]], d[["z"]]).fit(cov_type="robust")222 av, sp = sub["role_valeur_immeuble"].to_numpy(float), sub["amount"].to_numpy(float)223 share_condo = float((sub["prop_class"] == "condo").mean())224 rows.append({"city": city, "n": len(sub), "n_cells": sub["cell"].nunique(),225 "gamma_fe": float(fe.params["ln_price"]) - 1, "se_fe": float(fe.std_errors["ln_price"]),226 "gamma_iv": float(iv.params["ln_price_w"]) - 1, "se_iv": float(iv.std_errors["ln_price_w"]),227 "share_condo": share_condo})228 return pd.DataFrame(rows)229230231def trim_iqr(df: pd.DataFrame) -> tuple[pd.DataFrame, pd.DataFrame]:232 q1 = df.groupby("cell")["ratio"].transform(lambda s: s.quantile(.25))233 q3 = df.groupby("cell")["ratio"].transform(lambda s: s.quantile(.75))234 iqr = q3 - q1235 keep = df["ratio"].between(q1 - 1.5 * iqr, q3 + 1.5 * iqr)236 sub = _cells_ok(df[keep])237 fe, iv = models.cheng_fe(sub), models.clapp_iv(sub)238 ts = models.tax_shift(sub)239 av, sp = sub["role_valeur_immeuble"].to_numpy(float), sub["amount"].to_numpy(float)240 head = pd.DataFrame([{"variant": "IAAO 1.5 x IQR trim within cells", "n": len(sub),241 "share_dropped": 1 - keep.mean(),242 "gamma_fe": fe["gamma"], "se_fe": fe["se"],243 "gamma_iv": iv["gamma"], "se_iv": iv["se"],244 "cod": iaao.cod(av / sp), "prd": iaao.prd(av, sp),245 "prb": iaao.prb(av, sp)[0],246 "median_rel_d1": float(ts.loc[ts.decile == 1, "median_rel"].iloc[0]),247 "median_rel_d2": float(ts.loc[ts.decile == 2, "median_rel"].iloc[0]),248 "median_rel_d10": float(ts.loc[ts.decile == 10, "median_rel"].iloc[0])}])249 return head, ts250251252def clusters(df: pd.DataFrame) -> pd.DataFrame:253 muni_sales = df.groupby("muni")["muni"].transform("size")254 groups = {255 "Single-family": df["prop_class"] == "single_family",256 "Condominium": df["prop_class"] == "condo",257 "Plex (2–5 units)": df["prop_class"] == "plex",258 "Cottage": df["prop_class"] == "cottage",259 "Muni < 1k sales": muni_sales < 1_000,260 "Muni 1k–10k sales": muni_sales.between(1_000, 10_000),261 "Muni > 10k sales": muni_sales > 10_000,262 }263 rows = []264 for k, m in groups.items():265 sub = _cells_ok(df[m])266 rows.append({"group": k, "n": len(sub), "n_clusters": sub["muni"].nunique(),267 "n_cells": sub["cell"].nunique()})268 rows.append({"group": "Baseline", "n": len(df), "n_clusters": df["muni"].nunique(),269 "n_cells": df["cell"].nunique()})270 return pd.DataFrame(rows)271272273# ---------------------------------------------------------------- trimming grid274def _trim_mask(df: pd.DataFrame, rule: str) -> pd.Series:275 r = df["ratio"]276 g = df.groupby("cell")["ratio"]277 if rule == "baseline":278 return pd.Series(True, index=df.index)279 if rule.endswith("iqr"):280 k = float(rule.replace("iqr", ""))281 q1, q3 = g.transform(lambda s: s.quantile(.25)), g.transform(lambda s: s.quantile(.75))282 return r.between(q1 - k * (q3 - q1), q3 + k * (q3 - q1))283 if rule.endswith("pct"):284 a = float(rule.replace("pct", "")) / 100285 lo, hi = g.transform(lambda s: s.quantile(a)), g.transform(lambda s: s.quantile(1 - a))286 return r.between(lo, hi)287 raise ValueError(rule)288289290RULES = {"baseline": "Baseline (1/99 pct by roll vintage)",291 "3iqr": "Within-cell 3 x IQR (extreme outliers)",292 "1.5iqr": "Within-cell 1.5 x IQR (IAAO outlier rule)",293 "2.5pct": "Within-cell 2.5/97.5 pct",294 "5pct": "Within-cell 5/95 pct"}295296297def muni_shares(d: pd.DataFrame) -> dict:298 blocks = iaao.group_metrics(d, ["muni", "sale_year"], min_n=50)299 blocks["muni"] = blocks["group"].str.rsplit("_", n=1).str[0]300 agg = blocks.groupby("muni").agg(n=("n", "sum"), cod=("cod", "median"),301 prb=("prb", "median")).reset_index()302 agg = agg[agg["n"] >= config.MUNI_MIN_SALES]303 return {"n_munis": len(agg), "median_muni_cod": agg["cod"].median(),304 "share_muni_cod_gt_15": (agg["cod"] > 15).mean(),305 "share_muni_prb_neg": (agg["prb"] < 0).mean(),306 "share_muni_prb_lt_m05": (agg["prb"] < -0.05).mean(),307 "share_muni_prb_lt_m10": (agg["prb"] < -0.10).mean(),308 "median_muni_prb": agg["prb"].median()}309310311def trim_grid(df: pd.DataFrame) -> pd.DataFrame:312 rows = []313 for rule, label in RULES.items():314 keep = _trim_mask(df, rule)315 sub = _cells_ok(df[keep])316 fe, iv = models.cheng_fe(sub), models.clapp_iv(sub)317 ts = models.tax_shift(sub)318 av, sp = sub["role_valeur_immeuble"].to_numpy(float), sub["amount"].to_numpy(float)319 sp_ = spline(sub)320 row = {"rule": rule, "label": label, "n": len(sub), "share_dropped": 1 - keep.mean(),321 "median_ratio": float(np.median(av / sp)), "cod": iaao.cod(av / sp),322 "prd": iaao.prd(av, sp), "prb": iaao.prb(av, sp)[0], "prb_se": iaao.prb(av, sp)[1],323 "gamma_fe": fe["gamma"], "se_fe": fe["se"], "gamma_iv": iv["gamma"], "se_iv": iv["se"],324 "d1_median": float(ts.loc[ts.decile == 1, "median_rel"].iloc[0]),325 "d10_median": float(ts.loc[ts.decile == 10, "median_rel"].iloc[0]),326 "spline_seg1_beta": float(sp_.loc[0, "beta"]), "spline_seg5_beta": float(sp_.loc[4, "beta"])}327 row.update(muni_shares(sub))328 rows.append(row)329 print(f" {label:<44} drop={row['share_dropped']:.1%} COD={row['cod']:.1f} PRB={row['prb']:+.3f} "330 f"gFE={row['gamma_fe']:+.3f} gIV={row['gamma_iv']:+.3f} d1={row['d1_median']:+.1%} "331 f"muniPRB<-.05={row['share_muni_prb_lt_m05']:.0%} muniCOD>15={row['share_muni_cod_gt_15']:.0%}")332 return pd.DataFrame(rows)333334335def het_trimmed(df: pd.DataFrame) -> pd.DataFrame:336 sub = df[_trim_mask(df, "1.5iqr")]337 muni_sales = sub.groupby("muni")["muni"].transform("size")338 groups = {339 "Single-family": sub["prop_class"] == "single_family",340 "Condominium": sub["prop_class"] == "condo",341 "Plex (2–5 units)": sub["prop_class"] == "plex",342 "Age < 20 y": sub["age"] < 20, "Age > 60 y": sub["age"] > 60,343 "Land share < 0.2": sub["land_share"] < 0.2, "Land share > 0.4": sub["land_share"] > 0.4,344 "Muni < 1k sales": muni_sales < 1_000, "Muni > 10k sales": muni_sales > 10_000,345 "Sales 2021–2023": sub["sale_year"] <= 2023, "Sales 2024–2026": sub["sale_year"] >= 2024,346 }347 out = models.gamma_by_group(sub, groups)348 ivs = []349 for label, mask in groups.items():350 s2 = _cells_ok(sub[mask])351 ivs.append(models.clapp_iv(s2)["gamma"] if len(s2) >= 2000 else np.nan)352 out["gamma_iv"] = ivs353 return out354355356def cities_trimmed(df: pd.DataFrame) -> pd.DataFrame:357 sub = df[_trim_mask(df, "1.5iqr")]358 big = sub[sub["role_municipalite"].isin(config.BIG_CITIES)].copy()359 big["muni"] = big["role_municipalite"]360 blocks = iaao.group_metrics(big, ["muni", "sale_year"], min_n=200)361 blocks["muni"] = blocks["group"].str.rsplit("_", n=1).str[0]362 agg = blocks.groupby("muni").agg(n=("n", "sum"), median_ratio=("median_ratio", "median"),363 cod=("cod", "median"), prd=("prd", "median"),364 prb=("prb", "median")).reset_index()365 share_neg = blocks.assign(neg=blocks["prb"] < 0).groupby("muni")["neg"].mean().rename("share_years_prb_neg")366 agg = agg.merge(share_neg, on="muni")367 g = cities(sub).rename(columns={"city": "muni"})368 return agg.merge(g[["muni", "gamma_fe", "se_fe", "gamma_iv", "se_iv", "share_condo"]], on="muni")369370371def main() -> None:372 config.ensure_dirs()373 t0 = time.time()374 df = sample.load()375 print(f"sample {len(df):,}; share lag > 54 m: {(df['lag_months'] > 54).mean():.3%}; "376 f"muni-years with 2 rolls: {(df.groupby(['muni','sale_year'])['roll'].nunique() > 1).sum()}")377 pd.Series({"share_lag_gt_54": (df["lag_months"] > 54).mean(),378 "share_lag_gt_60": (df["lag_months"] > 60).mean(),379 "max_lag": df["lag_months"].max(),380 "muni_years": df.groupby(["muni", "sale_year"]).ngroups,381 "muni_years_two_rolls": int((df.groupby(["muni", "sale_year"])["roll"].nunique() > 1).sum()),382 "share_ratio_gt_1_5": (df["ratio"] > 1.5).mean(),383 "share_ratio_gt_2": (df["ratio"] > 2).mean()}).to_csv(OUT / "rev_misc.csv")384385 b = bounds(df); b.to_csv(OUT / "rev_bounds.csv", index=False); print(b.to_string())386 s = spline(df); s.to_csv(OUT / "rev_spline.csv", index=False); print(s.to_string())387 t = taxshift_alt(df); t.to_csv(OUT / "rev_taxshift_alt.csv", index=False)388 print(t[t.decile.isin([1, 2, 10])].to_string())389 c = cod_robust(df); c.to_csv(OUT / "rev_cod.csv", index=False); print(c.to_string())390 p = prb_blocks(df); p.to_csv(OUT / "rev_prb_blocks.csv", index=False); print(p.to_string())391 ci = cities(df); ci.to_csv(OUT / "rev_cities.csv", index=False); print(ci.to_string())392 head, ts = trim_iqr(df); head.to_csv(OUT / "rev_trim.csv", index=False)393 ts.to_csv(OUT / "rev_trim_taxshift.csv", index=False); print(head.T.to_string())394 cl = clusters(df); cl.to_csv(OUT / "rev_clusters.csv", index=False); print(cl.to_string())395 print("Trimming grid")396 tg = trim_grid(df); tg.to_csv(OUT / "rev_trim_grid.csv", index=False)397 ht = het_trimmed(df); ht.to_csv(OUT / "rev_het_trimmed.csv", index=False); print(ht.to_string())398 ct = cities_trimmed(df); ct.to_csv(OUT / "rev_cities_trimmed.csv", index=False); print(ct.to_string())399 # share of decile-1 sales that are ratio outliers400 d = df.copy()401 d["decile"] = d.groupby("cell")["amount"].rank(pct=True).mul(10).add(1 - 1e-9).astype(int).clip(1, 10)402 d["outlier"] = ~_trim_mask(df, "1.5iqr")403 (d.groupby("decile")["outlier"].mean().rename("share_outlier").reset_index()404 .to_csv(OUT / "rev_outlier_by_decile.csv", index=False))405 print(d.groupby("decile")["outlier"].mean().round(3).to_dict())406 pt = prb_blocks(df[_trim_mask(df, "1.5iqr")])407 pt.to_csv(OUT / "rev_prb_blocks_trimmed.csv", index=False); print(pt.to_string())408 # province-wide PRB by year on the trimmed sample409 tr = df[_trim_mask(df, "1.5iqr")]410 rows = [{"year": "All", "prb": iaao.prb(tr["role_valeur_immeuble"].to_numpy(float), tr["amount"].to_numpy(float))[0],411 "cod": iaao.cod((tr["role_valeur_immeuble"] / tr["amount"]).to_numpy()), "n": len(tr)}]412 for y, g in tr.groupby("sale_year"):413 rows.append({"year": str(y), "prb": iaao.prb(g["role_valeur_immeuble"].to_numpy(float), g["amount"].to_numpy(float))[0],414 "cod": iaao.cod((g["role_valeur_immeuble"] / g["amount"]).to_numpy()), "n": len(g)})415 pd.DataFrame(rows).to_csv(OUT / "rev_iaao_year_trimmed.csv", index=False); print(pd.DataFrame(rows).to_string())416 print(f"done in {time.time() - t0:.0f} s")417418419if __name__ == "__main__":420 main()421