spb/qwhpi Public

1#!/usr/bin/env python32# =============================================================================3# QWHPI — Quebec Weekly Housing Price Index4# Author  : Simon-Pierre Boucher5# Contact : contact@spboucher.ai6# File    : engine/scripts/08_figures.py7# Purpose : Publication figures — aggregate, regions, cities, median-vs-8#           hedonic, repeat sales, downsampling, volumes, YoY map, gap.9# =============================================================================10"""Research figures (Execution Order step 10)."""1112from __future__ import annotations1314import sys15from pathlib import Path1617sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src"))1819import numpy as np20import polars as pl2122from qwhpi import plotting as P23from qwhpi.config import FIGURES_DIR, PROCESSED_DIR, TABLES_DIR, ensure_dirs2425import matplotlib.dates as mdates  # noqa: E40226import matplotlib.pyplot as plt  # noqa: E402272829def dts(weeks: pl.Series) -> np.ndarray:30    return weeks.str.to_date("%Y-%m-%d").to_numpy()313233def cell(c: pl.DataFrame, gid: str, t: str) -> pl.DataFrame:34    return c.filter((pl.col("geography_id") == gid)35                    & (pl.col("property_type") == t)).sort("week")363738def main() -> None:39    ensure_dirs()40    c = pl.read_parquet(PROCESSED_DIR / "qwhpi_weekly.parquet")4142    # 1 — Quebec aggregate with CI band43    q = cell(c, "quebec", "all")44    fig, ax = plt.subplots(figsize=(9, 4.5))45    x = dts(q["week"])46    ax.fill_between(x, q["lower_95"], q["upper_95"], color=P.BAND, label="95% CI")47    ax.plot(x, q["index_smoothed"], color=P.SERIES[0], label="QWHPI-QC (smoothed)")48    ax.plot(x, q["index"], color=P.SERIES[0], alpha=0.35, lw=1, label="raw weekly")49    ax.set_title("Quebec weekly housing price index — all types (2021 = 100)")50    ax.legend(loc="upper left")51    ax.xaxis.set_major_formatter(mdates.DateFormatter("%Y"))52    P.save(fig, FIGURES_DIR / "quebec_aggregate.png")5354    # 2 — 17 regions, small multiples55    regs = sorted(c.filter(pl.col("geography_level") == "region")["geography_id"].unique())56    fig, axes = plt.subplots(3, 6, figsize=(15, 7), sharex=True, sharey=True)57    for ax, rid in zip(axes.flat, regs):58        r = cell(c, rid, "all")59        ax.plot(dts(r["week"]), r["index_smoothed"], color=P.SERIES[0], lw=1.4)60        ax.set_title(r["geography_name"][0], fontsize=8)61        ax.xaxis.set_major_formatter(mdates.DateFormatter("%y"))62    for ax in axes.flat[len(regs):]:63        ax.axis("off")64    fig.suptitle("Regional indexes — all types (2021 = 100)")65    P.save(fig, FIGURES_DIR / "region_comparison.png")6667    # 3 — Big-4 cities68    fig, ax = plt.subplots(figsize=(9, 4.5))69    for i, (gid, label) in enumerate([("montreal", "Montréal"), ("quebec-city", "Québec"),70                                       ("laval", "Laval"), ("gatineau", "Gatineau")]):71        s = cell(c, gid, "all")72        ax.plot(dts(s["week"]), s["index_smoothed"], color=P.SERIES[i], label=label)73    ax.set_title("Major cities — all types (2021 = 100)")74    ax.legend(loc="upper left")75    P.save(fig, FIGURES_DIR / "big4_cities.png")7677    # 4 — Condo across cities78    fig, ax = plt.subplots(figsize=(9, 4.5))79    for i, (gid, label) in enumerate([("montreal", "Montréal"), ("quebec-city", "Québec"),80                                       ("gatineau", "Gatineau")]):81        s = cell(c, gid, "condo")82        ax.plot(dts(s["week"]), s["index_smoothed"], color=P.SERIES[i], label=label)83    ax.set_title("Condo indexes by city (2021 = 100)")84    ax.legend(loc="upper left")85    P.save(fig, FIGURES_DIR / "condo_by_city.png")8687    # 5 — Raw median vs hedonic (province, unifamilial)88    bi = pl.read_parquet("data/interim/baseline_indexes.parquet")89    b = bi.filter((pl.col("geography_id") == "quebec")90                  & (pl.col("property_type") == "unifamilial")).sort("week")91    med = b["raw_median"].to_numpy()92    base_med = np.nanmean(med[:52])93    fig, ax = plt.subplots(figsize=(9, 4.5))94    ax.plot(dts(b["week"]), med / base_med * 100, color=P.SERIES[1], lw=1.2,95            label="Raw weekly median (indexed)")96    h = cell(c, "quebec", "unifamilial")97    ax.plot(dts(h["week"]), h["index_smoothed"], color=P.SERIES[0],98            label="QWHPI (quality-adjusted)")99    ax.set_title("Composition noise: raw median vs hedonic index — Quebec unifamilial")100    ax.legend(loc="upper left")101    P.save(fig, FIGURES_DIR / "raw_median_vs_hedonic.png")102103    # 6 — Weekly volumes104    vol = pl.read_parquet(PROCESSED_DIR / "weekly_liquidity.parquet")105    v = vol.filter((pl.col("geography_id") == "quebec")106                   & (pl.col("property_type") == "all")).sort("week")107    fig, ax = plt.subplots(figsize=(9, 3.2))108    ax.bar(dts(v["week"]), v["transactions"], width=6, color=P.SERIES[0], alpha=0.6)109    ax.set_title("Weekly transaction volume — Quebec (research sample)")110    P.save(fig, FIGURES_DIR / "volumes.png")111112    # 7 — Downsampling stability113    ds = pl.read_csv(TABLES_DIR / "downsampling_results.csv")114    fig, ax = plt.subplots(figsize=(7, 4))115    ax.plot(ds["target_tx_per_week"], ds["rmse_log_pct"], "o-", color=P.SERIES[0],116            label="RMSE vs full-sample path (%)")117    ax.plot(ds["target_tx_per_week"], (1 - ds["ci95_coverage"]) * 10, "s--",118            color=P.SERIES[1], label="(1 − CI coverage) × 10")119    ax.set_xlabel("target transactions per week")120    ax.set_title("Downsampling experiment — Montréal condo")121    ax.legend()122    ax.invert_xaxis()123    P.save(fig, FIGURES_DIR / "downsampling_stability.png")124125    # 8 — YoY map (choropleth)126    import geopandas as gpd127    latest = c.filter((pl.col("geography_level") == "region")128                      & (pl.col("property_type") == "all")129                      & (~pl.col("is_partial_week")))130    lw = latest.filter(pl.col("week") == latest["week"].max())131    g = gpd.read_file("web/public/regions.geojson")132    g = g.merge(lw.select("geography_id", "yoy_pct").to_pandas(), on="geography_id")133    fig, ax = plt.subplots(figsize=(8, 7))134    g.plot(column="yoy_pct", cmap="Blues", legend=True, ax=ax,135           edgecolor="#fcfcfb", linewidth=0.8,136           legend_kwds={"label": "YoY %", "shrink": 0.6})137    ax.set_axis_off()138    ax.set_title(f"Year-over-year appreciation by region — week of {lw['week'][0]}")139    P.save(fig, FIGURES_DIR / "yoy_map.png")140141    # 9 — Assessment gap evolution142    gap = pl.read_parquet(PROCESSED_DIR / "assessment_gap.parquet")143    fig, ax = plt.subplots(figsize=(9, 4))144    for i, t in enumerate(["unifamilial", "condo", "plex"]):145        s = gap.filter((pl.col("geography_id") == "quebec")146                       & (pl.col("property_type") == t)).sort("week")147        # 13-week rolling median for readability148        vals = s["median_ratio"].rolling_median(window_size=13)149        ax.plot(dts(s["week"]), vals, color=P.SERIES[i], label=t)150    ax.axhline(1.0, color=P.TEXT2, lw=1, ls="--")151    ax.set_title("Assessment gap — median sale price / municipal assessment (13-wk median)")152    ax.legend(loc="upper left")153    P.save(fig, FIGURES_DIR / "assessment_gap.png")154155    # 10 — Reliability heat table (cells × grade)156    summary = (157        c.group_by(["geography_id", "property_type"])158        .agg(pl.col("reliability_grade").mode().first().alias("grade"))159    )160    counts = summary.group_by(["property_type", "grade"]).len().sort(["property_type", "grade"])161    fig, ax = plt.subplots(figsize=(7, 3.4))162    grades = ["A", "B", "C", "D", "E"]163    types = ["all", "unifamilial", "condo", "plex"]164    mat = np.zeros((len(types), len(grades)))165    for r in counts.iter_rows(named=True):166        mat[types.index(r["property_type"]), grades.index(r["grade"])] = r["len"]167    im = ax.imshow(mat, cmap="Blues", aspect="auto")168    ax.set_xticks(range(len(grades)), grades)169    ax.set_yticks(range(len(types)), types)170    for i in range(len(types)):171        for j in range(len(grades)):172            if mat[i, j]:173                ax.text(j, i, int(mat[i, j]), ha="center", va="center",174                        color=P.TEXT if mat[i, j] < mat.max() * 0.6 else "white")175    ax.set_title("Published series by reliability grade")176    ax.grid(False)177    P.save(fig, FIGURES_DIR / "reliability_matrix.png")178179    print("All figures written to outputs/figures/")180181182if __name__ == "__main__":183    main()184