spb/wp10_uqo Public

1#!/usr/bin/env python32# Author: Simon-Pierre Boucher — contact@spboucher.ai3"""Step 04 — Publication figures (print-journal calibre).45Reads the analysis sample and the step-02/03 outputs, writes ten PNG figures6to ``figures/``. No titles are drawn inside single-panel figures — captions7in the manuscript carry the message; multi-panel figures use bold8"Panel A/B" headers. One accent hue per figure, doubled by linestyle or9marker so nothing is encoded by colour alone.1011Usage:  python scripts/04_make_figures.py12"""13import sys14from pathlib import Path1516import numpy as np17import pandas as pd1819sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src"))2021from wp10 import config, sample  # noqa: E40222from wp10.plotstyle import (BLUE, BLUES, GREY, INK, LIGHT, RED, TEXTWIDTH,  # noqa: E40223                            apply_style, panel_label, ygrid)2425import matplotlib.pyplot as plt  # noqa: E40226import matplotlib.dates as mdates  # noqa: E4022728OUT = config.FIGURES29RES = config.REPRODUCED303132def _save(fig, name):33    fig.savefig(OUT / name, bbox_inches="tight", pad_inches=0.02)34    plt.close(fig)353637# ---------------------------------------------------------------- fig 138def fig_ratio_dist(df):39    """Panel A: ratio histogram. Panel B: densities by roll lag."""40    fig, axes = plt.subplots(1, 2, figsize=(TEXTWIDTH, 2.9))4142    ax = axes[0]43    ax.hist(df["ratio"], bins=120, range=(0, 2.5), color=LIGHT,44            edgecolor=INK, linewidth=0.25)45    med = df["ratio"].median()46    ax.axvline(med, color=INK, lw=0.9)47    ax.axvline(1.0, color=GREY, lw=0.8, ls=(0, (4, 3)))48    ax.text(med - 0.05, ax.get_ylim()[1] * 0.97, f"median = {med:.2f}",49            ha="right", va="top", fontsize=8)50    ax.text(1.04, ax.get_ylim()[1] * 0.75, "AV = SP", fontsize=8, color=GREY)51    ax.set_xlabel("Assessment ratio $AV/SP$")52    ax.set_ylabel("Sales")53    ax.set_yticks([])54    ax.spines["left"].set_visible(False)55    panel_label(ax, "Panel A. All sales")5657    ax = axes[1]58    specs = [("Roll lag $<$ 24 m", df["lag_months"] < 24, BLUES[2], "-"),59             ("24–48 m", df["lag_months"].between(24, 48), BLUES[4], (0, (5, 2))),60             ("$>$ 48 m", df["lag_months"] > 48, BLUES[5], (0, (1, 1.2)))]61    for label, m, color, ls in specs:62        ax.hist(df.loc[m, "ratio"], bins=120, range=(0, 2.5), density=True,63                histtype="step", lw=1.2, color=color, ls=ls, label=label)64    ax.set_xlabel("Assessment ratio $AV/SP$")65    ax.set_ylabel("Density")66    ax.legend(loc="upper right", handlelength=2.4)67    panel_label(ax, "Panel B. By roll lag at sale")68    fig.subplots_adjust(wspace=0.25)69    _save(fig, "fig_ratio_dist.png")707172# ---------------------------------------------------------------- fig 273def fig_binscatter():74    """Within-cell binned scatter of ln ratio on ln price — the core fact."""75    b = pd.read_csv(RES / "binscatter.csv")76    fig, ax = plt.subplots(figsize=(0.72 * TEXTWIDTH, 3.3))77    ax.axhline(0, color=GREY, lw=0.7, ls=(0, (4, 3)))78    slope = np.polyfit(b["x"], b["y"], 1, w=b["n"])[0]79    icept = np.average(b["y"] - slope * b["x"], weights=b["n"])80    xs = np.linspace(b["x"].min(), b["x"].max(), 50)81    ax.plot(xs, slope * xs + icept, color=BLUE, lw=1.2, zorder=2)82    ax.errorbar(b["x"], b["y"], yerr=1.96 * b["se"], fmt="o", ms=4,83                mfc=INK, mec=INK, ecolor=INK, elinewidth=0.6, capsize=0,84                lw=0, zorder=3)85    ax.annotate(f"slope = ${slope:.3f}$", xy=(0.55, slope * 0.55 + icept),86                xytext=(0.32, 0.12), fontsize=8.5, color=BLUE,87                arrowprops=dict(arrowstyle="-", color=BLUE, lw=0.6,88                                shrinkA=2, shrinkB=2))89    ax.set_xlabel("ln sale price (demeaned within municipality × roll × year)")90    ax.set_ylabel("ln assessment ratio (demeaned)")91    _save(fig, "fig_binscatter.png")929394# ---------------------------------------------------------------- fig 395def fig_time(df):96    """Median ratio by sale month; sequential blues by roll vintage,97    each segment labelled directly (no legend)."""98    d = df.copy()99    d["month"] = d["sale_date"].dt.to_period("M").dt.to_timestamp()100    fig, ax = plt.subplots(figsize=(TEXTWIDTH, 2.9))101    rolls = sorted(d["roll"].unique())102    for i, roll in enumerate(rolls):103        gg = d[d["roll"] == roll].groupby("month")["ratio"]104        g = gg.median()[gg.size() >= 100]  # drop thin months (spurious spikes)105        g = g[g.index.notna()].sort_index()106        if len(g) < 3:107            continue108        color = BLUES[min(i, len(BLUES) - 1)]109        ax.plot(g.index, g.values, lw=1.3, color=color)110        ax.annotate(roll, xy=(g.index[-1], g.values[-1]),111                    xytext=(3, 0), textcoords="offset points",112                    fontsize=7.5, color=color, va="center")113    ax.axhline(1.0, color=GREY, lw=0.7, ls=(0, (4, 3)))114    ax.text(pd.Timestamp("2021-02-01"), 1.012, "parity", fontsize=7.5,115            color=GREY)116    ax.set_ylabel("Median assessment ratio")117    ax.xaxis.set_major_formatter(mdates.DateFormatter("%Y"))118    ygrid(ax)119    _save(fig, "fig_time.png")120121122# ---------------------------------------------------------------- fig 4123def fig_prb_muni():124    """Distribution of municipality-level PRB."""125    m = pd.read_csv(RES / "iaao_muni.csv")126    fig, ax = plt.subplots(figsize=(0.72 * TEXTWIDTH, 3.1))127    ax.hist(m["prb"], bins=40, color=LIGHT, edgecolor=INK, linewidth=0.3)128    lo, hi = config.IAAO_PRB_RANGE129    ax.axvspan(lo, hi, color="0.92", zorder=0)130    ax.axvline(0, color=GREY, lw=0.7, ls=(0, (4, 3)))131    med = m["prb"].median()132    ax.axvline(med, color=RED, lw=1.0)133    ymax = ax.get_ylim()[1]134    ax.text(med - 0.008, ymax * 0.97, f"median = {med:.2f}", ha="right",135            va="top", fontsize=8, color=RED)136    ax.text((lo + hi) / 2, ymax * 0.55, "IAAO\nband", ha="center",137            fontsize=7.5, color="0.35")138    share = (m["prb"] < 0).mean()139    ax.text(0.02, 0.97, f"{share:.0%} of municipalities\nhave PRB $<$ 0",140            transform=ax.transAxes, fontsize=8, va="top")141    ax.set_xlabel("Municipality-level PRB (median across sale years)")142    ax.set_ylabel("Municipalities")143    ygrid(ax)144    _save(fig, "fig_prb_muni.png")145146147# ---------------------------------------------------------------- fig 5148def fig_map():149    """Municipal PRB across the province (diverging hue, neutral midpoint)."""150    m = pd.read_csv(RES / "iaao_muni.csv")151    fig, ax = plt.subplots(figsize=(TEXTWIDTH, 4.4))152    v = m["prb"].clip(-0.30, 0.10)153    sc = ax.scatter(m["lng"], m["lat"], c=v, s=np.sqrt(m["n"]) * 0.9,154                    cmap="RdBu", vmin=-0.30, vmax=0.30,155                    edgecolor=INK, linewidth=0.25, alpha=0.9)156    cb = fig.colorbar(sc, ax=ax, shrink=0.75, pad=0.02)157    cb.set_label("PRB (negative = regressive)", fontsize=8)158    cb.ax.tick_params(labelsize=7.5)159    cb.outline.set_linewidth(0.5)160    offsets = {"Montréal": (6, -14), "Québec": (8, 4), "Gatineau": (-8, -14),161               "Sherbrooke": (8, -10), "Saguenay": (8, 4)}162    for _, r in m.nlargest(12, "n").iterrows():163        if r["name"] in offsets:164            ax.annotate(r["name"], (r["lng"], r["lat"]),165                        xytext=offsets[r["name"]], textcoords="offset points",166                        fontsize=7.5)167    ax.set_xlabel("Longitude")168    ax.set_ylabel("Latitude")169    ax.set_xlim(-80, -63)170    ax.set_ylim(44.9, 49.6)171    _save(fig, "fig_map.png")172173174# ---------------------------------------------------------------- fig 6175def fig_quantile():176    """β(τ) with a shaded 95% band; FE/IV benchmarks labelled directly."""177    q = pd.read_csv(RES / "quantile.csv")178    v = pd.read_csv(RES / "vertical.csv")179    fe = v.loc[v["estimator"].str.startswith("Cheng FE"), "beta"].iloc[0]180    iv = v.loc[v["estimator"].str.startswith("Clapp"), "beta"].iloc[0]181    fig, ax = plt.subplots(figsize=(0.72 * TEXTWIDTH, 3.3))182    ax.fill_between(q["tau"], q["beta"] - 1.96 * q["se"],183                    q["beta"] + 1.96 * q["se"], color=BLUE, alpha=0.18, lw=0)184    ax.plot(q["tau"], q["beta"], "o-", color=BLUE, ms=4, lw=1.3)185    x1, x0 = q["tau"].max(), q["tau"].min()186    for yv, lab, ls, xa, ha in [187            (1.0, r"$\beta = 1$ (proportional)", (0, (4, 3)), x1, "right"),188            (iv, f"Clapp IV ({iv:.2f})", (0, (1, 1.2)), x1, "right"),189            (fe, f"Cheng FE ({fe:.2f})", (0, (6, 2)), x0, "left")]:190        ax.axhline(yv, color=GREY, lw=0.8, ls=ls)191        ax.annotate(lab, xy=(xa, yv), xytext=(0, 3), textcoords="offset points",192                    ha=ha, fontsize=7.5, color="0.35")193    ax.set_xlabel(r"Quantile $\tau$ of the conditional $\ln AV$ distribution")194    ax.set_ylabel(r"$\beta(\tau)$")195    ax.set_xticks(q["tau"])196    _save(fig, "fig_quantile.png")197198199# ---------------------------------------------------------------- fig 7200def fig_heterogeneity():201    """Forest plot of Cheng-FE γ by subgroup, with panel groupings."""202    h = pd.read_csv(RES / "heterogeneity.csv").set_index("group")203    panels = [204        ("Property class", ["Single-family", "Condominium", "Plex (2–5 units)",205                            "Cottage"]),206        ("Building age", ["Age < 20 y", "Age 20–60 y", "Age > 60 y"]),207        ("Assessed land share", ["Land share < 0.2", "Land share 0.2–0.4",208                                 "Land share > 0.4"]),209        ("Roll lag at sale", ["Roll lag < 24 m", "Roll lag 24–48 m",210                              "Roll lag > 48 m"]),211        ("Municipality size", ["Muni < 1k sales", "Muni 1k–10k sales",212                               "Muni > 10k sales"]),213        ("Sale year", [f"Sales {y}" for y in range(2021, 2027)]),214    ]215    rows, ypos, headers, seps = [], [], [], []216    y = 0217    for name, keys in panels:218        headers.append((y, name))219        y -= 1220        for k in keys:221            if k in h.index:222                rows.append((y, k, h.loc[k]))223                ypos.append(y)224                y -= 1225        seps.append(y + 0.5)226        y -= 0.6227    fig, ax = plt.subplots(figsize=(0.85 * TEXTWIDTH, 5.6))228    for yy, k, r in rows:229        ax.errorbar(r["gamma"], yy, xerr=1.96 * r["se"], fmt="o", ms=3.8,230                    mfc=INK, mec=INK, ecolor=INK, elinewidth=0.8, capsize=1.5)231    ax.axvline(0, color=GREY, lw=0.7, ls=(0, (4, 3)))232    for yy, name in headers:233        ax.text(-0.72, yy, name, fontsize=8.5, fontweight="bold", va="center")234    labels = {yy: k.replace("<", "$<$").replace(">", "$>$")235              for yy, k, _ in rows}236    ax.set_yticks(list(labels.keys()))237    ax.set_yticklabels(labels.values(), fontsize=8)238    ax.set_ylim(y + 0.4, 1.0)239    ax.set_xlim(-0.72, 0.12)240    ax.set_xlabel(r"$\gamma$ = elasticity of the assessment ratio with respect"241                  " to price (negative = regressive)")242    ax.spines["left"].set_visible(False)243    ax.tick_params(axis="y", length=0)244    _save(fig, "fig_heterogeneity.png")245246247# ---------------------------------------------------------------- fig 8248def fig_taxshift():249    """Median excess tax burden by within-market price decile (polarity)."""250    t = pd.read_csv(RES / "taxshift.csv")251    fig, ax = plt.subplots(figsize=(0.85 * TEXTWIDTH, 3.1))252    vals = 100 * t["median_rel"]253    colors = [RED if v > 0 else BLUE for v in vals]254    ax.bar(t["decile"], vals, width=0.72, color=colors, edgecolor=INK,255           linewidth=0.4)256    ax.axhline(0, color=INK, lw=0.7)257    for d, v in zip(t["decile"], vals):258        va = "bottom" if v > 0 else "top"259        off = 1.2 if v > 0 else -1.2260        ax.text(d, v + off, f"{v:+.1f}", ha="center", va=va, fontsize=7.5)261    ax.text(2.5, 45, "over-taxed", fontsize=8, color=RED, ha="center")262    ax.text(8.5, 14, "under-taxed", fontsize=8, color=BLUE, ha="center")263    ax.set_xticks(t["decile"])264    ax.set_xlabel("Within-market sale-price decile")265    ax.set_ylabel("Median excess tax burden (%)")266    ax.set_ylim(min(vals) - 8, max(vals) + 9)267    _save(fig, "fig_taxshift.png")268269270# ---------------------------------------------------------------- fig 9271def fig_cod():272    """Panel A: COD distribution. Panel B: COD vs market size."""273    m = pd.read_csv(RES / "iaao_muni.csv")274    fig, axes = plt.subplots(1, 2, figsize=(TEXTWIDTH, 2.9))275276    ax = axes[0]277    ax.hist(m["cod"], bins=40, color=LIGHT, edgecolor=INK, linewidth=0.3)278    ax.axvline(config.IAAO_COD_MAX_SF, color=RED, lw=1.0)279    ymax = ax.get_ylim()[1]280    ax.text(config.IAAO_COD_MAX_SF + 0.6, ymax * 0.95,281            f"IAAO ceiling ({config.IAAO_COD_MAX_SF:.0f})", fontsize=7.5,282            color=RED, va="top")283    med = m["cod"].median()284    ax.axvline(med, color=INK, lw=0.9, ls=(0, (5, 2)))285    ax.text(med + 0.6, ymax * 0.72, f"median = {med:.1f}", fontsize=7.5)286    ax.set_xlabel("Municipality COD (median across sale years)")287    ax.set_ylabel("Municipalities")288    panel_label(ax, "Panel A. Distribution")289290    ax = axes[1]291    ax.scatter(m["n"], m["cod"], s=7, facecolor="none", edgecolor=INK,292               linewidth=0.5, alpha=0.65)293    ax.axhline(config.IAAO_COD_MAX_SF, color=RED, lw=0.9)294    ax.set_xscale("log")295    ax.set_xlabel("Sales in municipality, 2021–2026 (log scale)")296    ax.set_ylabel("COD")297    panel_label(ax, "Panel B. COD and market size")298    fig.subplots_adjust(wspace=0.25)299    _save(fig, "fig_cod.png")300301302# ---------------------------------------------------------------- fig 10303def fig_robustness():304    """γ across sample variants — FE (filled circles) vs IV (open squares)."""305    r = pd.read_csv(RES / "robustness.csv")306    fig, ax = plt.subplots(figsize=(0.85 * TEXTWIDTH, 3.6))307    y = np.arange(len(r))[::-1].astype(float)308    ax.errorbar(r["gamma_fe"], y + 0.16, xerr=1.96 * r["se_fe"], fmt="o",309                ms=4, mfc=BLUE, mec=BLUE, ecolor=BLUE, elinewidth=0.8,310                capsize=1.5, lw=0, label="Cheng FE")311    ax.errorbar(r["gamma_iv"], y - 0.16, xerr=1.96 * r["se_iv"], fmt="s",312                ms=4, mfc="white", mec=RED, ecolor=RED, elinewidth=0.8,313                capsize=1.5, lw=0, label="Clapp IV")314    ax.axvline(0, color=GREY, lw=0.7, ls=(0, (4, 3)))315    ax.set_yticks(y)316    labels = [str(v).replace("<=", "$\\leq$").replace(">=", "$\\geq$")317              .replace("<", "$<$").replace("$100k", "\\$100k")318              for v in r["variant"]]319    ax.set_yticklabels(labels, fontsize=8)320    ax.set_xlabel(r"$\gamma$ (negative = regressive)")321    ax.legend(loc="lower left", markerscale=1.1)322    ax.spines["left"].set_visible(False)323    ax.tick_params(axis="y", length=0)324    _save(fig, "fig_robustness.png")325326327def main() -> None:328    config.ensure_dirs()329    apply_style()330    df = sample.load()331    fig_ratio_dist(df)332    fig_binscatter()333    fig_time(df)334    fig_prb_muni()335    fig_map()336    fig_quantile()337    fig_heterogeneity()338    fig_taxshift()339    fig_cod()340    fig_robustness()341    made = sorted(p.name for p in OUT.glob("fig_*.png"))342    print(f"{len(made)} figures written:", ", ".join(made))343344345if __name__ == "__main__":346    main()347