spb/wp10_uqo Public

1#!/usr/bin/env python32# Author: Simon-Pierre Boucher — contact@spboucher.ai3"""Step 02 — Descriptive statistics and IAAO ratio-study diagnostics.45Writes to ``results/reproduced/``:6  summary_stats.csv   sample descriptives used in Table 17  iaao_overall.csv    province-wide median ratio / COD / PRD / PRB with8                      bootstrap CIs, overall and by sale year9  iaao_cities.csv     the ten largest markets10  iaao_muni.csv       every municipality with ≥ 100 sales (maps + histograms)1112Usage:  python scripts/02_iaao_stats.py13"""14import sys15from pathlib import Path1617import numpy as np18import pandas as pd1920sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src"))2122from wp10 import config, iaao, sample  # noqa: E402232425def main() -> None:26    config.ensure_dirs()27    df = sample.load()28    out = config.REPRODUCED2930    # ------------------------------------------------------------ Table 131    desc_vars = {32        "amount": "Sale price ($)",33        "role_valeur_immeuble": "Assessed value ($)",34        "ratio": "Assessment ratio AV/SP",35        "lag_months": "Roll lag (months)",36        "land_share": "Assessed land share",37        "age": "Building age (years)",38        "role_superficie_terrain_m2": "Lot area (m2)",39        "role_aire_etages_m2": "Floor area (m2)",40    }41    rows = []42    for var, label in desc_vars.items():43        s = df[var].dropna()44        rows.append({"variable": label, "n": len(s), "mean": s.mean(),45                     "sd": s.std(), "p10": s.quantile(.10), "p50": s.median(),46                     "p90": s.quantile(.90)})47    pd.DataFrame(rows).to_csv(out / "summary_stats.csv", index=False)4849    counts = {"n_sales": len(df), "n_munis": df["muni"].nunique(),50              "n_cells": df["cell"].nunique()}51    for k, v in df.groupby("prop_class").size().items():52        counts[f"n_{k}"] = int(v)53    pd.Series(counts).to_csv(out / "sample_counts.csv")5455    # ------------------------------------------------------------ overall + by year56    blocks = [("All sales 2021–2026", df)]57    blocks += [(str(y), g) for y, g in df.groupby("sale_year")]58    rows = []59    for label, g in blocks:60        av = g["role_valeur_immeuble"].to_numpy(float)61        sp = g["amount"].to_numpy(float)62        r = av / sp63        b, se = iaao.prb(av, sp)64        row = {"group": label, "n": len(g),65               "median_ratio": float(np.median(r)), "cod": iaao.cod(r),66               "prd": iaao.prd(av, sp), "prb": b, "prb_se": se}67        cis = iaao.bootstrap_ci(av, sp, n_boot=200)68        for stat, (lo, hi) in cis.items():69            row[f"{stat}_lo"], row[f"{stat}_hi"] = lo, hi70        rows.append(row)71        print(f"  {label:<22} n={row['n']:>8,}  med={row['median_ratio']:.3f}  "72              f"COD={row['cod']:.1f}  PRD={row['prd']:.3f}  PRB={row['prb']:+.4f}")73    pd.DataFrame(rows).to_csv(out / "iaao_overall.csv", index=False)7475    # ---------------------------------------------------------------------76    # Municipality-level statistics. Within a municipality × sale-year block77    # a single roll is in force, so the roll lag is (nearly) constant and the78    # COD/PRB are not inflated by market-time drift. Annual blocks are then79    # aggregated to one row per municipality (median across years, total n).80    def annual_then_aggregate(data: pd.DataFrame, min_n: int) -> pd.DataFrame:81        blocks = iaao.group_metrics(data, ["muni", "sale_year"], min_n=min_n)82        blocks["muni"] = blocks["group"].str.rsplit("_", n=1).str[0]83        agg = (blocks.groupby("muni")84               .agg(n=("n", "sum"), n_years=("n", "size"),85                    median_ratio=("median_ratio", "median"),86                    cod=("cod", "median"), prd=("prd", "median"),87                    prb=("prb", "median"))88               .reset_index())89        share_neg = (blocks.assign(neg=blocks["prb"] < 0)90                     .groupby("muni")["neg"].mean().rename("share_years_prb_neg"))91        return agg.merge(share_neg, on="muni")9293    # ------------------------------------------------------------ ten largest markets94    big = df[df["role_municipalite"].isin(config.BIG_CITIES)].copy()95    big["muni"] = big["role_municipalite"]  # aggregate by display name96    tab = annual_then_aggregate(big, min_n=200)97    tab.to_csv(out / "iaao_cities.csv", index=False)9899    # ------------------------------------------------------------ every muni ≥ 100 sales100    tab = annual_then_aggregate(df, min_n=50)101    tab = tab[tab["n"] >= config.MUNI_MIN_SALES]102    coords = df.groupby("muni")[["lat", "lng"]].median()103    names = df.groupby("muni")["role_municipalite"].first()104    tab = tab.merge(coords, left_on="muni", right_index=True)105    tab = tab.merge(names.rename("name"), left_on="muni", right_index=True)106    tab.to_csv(out / "iaao_muni.csv", index=False)107    print(f"\nMunicipality-level metrics: {len(tab)} municipalities "108          f"(median within-year COD {tab['cod'].median():.1f}, "109          f"share PRB<0: {(tab['prb'] < 0).mean():.1%})")110111112if __name__ == "__main__":113    main()114