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/11_monthly_canonical.py7# Purpose : v2.1 canonical dataset — qhpi_monthly.parquet with growth8#           horizons, reliability grades, dollar values, vintages, coverage.9# =============================================================================10"""Canonical monthly dataset build (v2.1).1112Reliability grades are anchored on the monthly downsampling experiment13(outputs/tables/downsampling_monthly.csv): smoothed-path RMSE ≈ 1.0–1.2% at14150+ tx/month, ≈1.5% at 75, ≈1.9% at 40, ≥2.6% at 10.15  A: median ≥ 150 tx/month and se_log < 1.5%16  B: median ≥ 75 and se_log < 2.5%17  C: median ≥ 40 (direct estimation floor)18  D: median ≥ 15 (hierarchical, heavy shrinkage)19  E: < 15 (model-implied, transparency only)20"""2122from __future__ import annotations2324import datetime as dt25import sys26import unicodedata27from pathlib import Path2829sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src"))3031import numpy as np32import polars as pl3334from qwhpi.clean import CLEAN_PARQUET, research_sample35from qwhpi.config import (36    CITY_TARGETS, INTERIM_DIR, PROCESSED_DIR, PROPERTY_TYPES, TABLES_DIR,37    ensure_dirs,38)3940MODEL_VERSION = "2.1.0"41CANONICAL = PROCESSED_DIR / "qhpi_monthly.parquet"42FIRST_RELEASE = PROCESSED_DIR / "first_release_monthly.parquet"43VINTAGE_DIR = PROCESSED_DIR / "vintages_monthly"44KEY = ["geography_id", "property_type", "period"]4546GROWTH = {"monthly_pct": 1, "three_month_pct": 3, "six_month_pct": 6, "yoy_pct": 12}474849def slug(name: str) -> str:50    s = unicodedata.normalize("NFKD", name).encode("ascii", "ignore").decode()51    return s.lower().replace(" ", "-").replace("'", "")525354def city_slug(name: str) -> str:55    return "quebec-city" if name == "Québec" else slug(name)565758def main() -> None:59    ensure_dirs()60    data_vintage = dt.date.today().isoformat()61    m = pl.read_parquet(INTERIM_DIR / "monthly_indexes.parquet")62    cl = research_sample(pl.read_parquet(CLEAN_PARQUET)).with_columns(63        pl.col("date").dt.strftime("%Y-%m").alias("period"))6465    # ---------------- transactions per published cell ---------------- #66    cell_frames = [67        cl.select(pl.lit("quebec").alias("geography_id"),68                  pl.col("propertyType").alias("property_type"),69                  "period", "log_amount", "amount"),70        cl.filter(pl.col("region_code").is_not_null()).select(71            pl.concat_str([pl.lit("region-"), pl.col("region_code")]).alias("geography_id"),72            pl.col("propertyType").alias("property_type"),73            "period", "log_amount", "amount"),74    ]75    for city in CITY_TARGETS:76        sub = cl.filter(pl.col("municipality") == city)77        if sub.height:78            cell_frames.append(sub.select(79                pl.lit(city_slug(city)).alias("geography_id"),80                pl.col("propertyType").alias("property_type"),81                "period", "log_amount", "amount"))82    tx = pl.concat(cell_frames)83    tx_all = pl.concat([tx, tx.with_columns(pl.lit("all").alias("property_type"))])8485    # ---------------- representative dollar value ---------------- #86    basket = (87        tx_all.join(m.select(*KEY, "index_smoothed"), on=KEY, how="inner")88        .group_by(["geography_id", "property_type"])89        .agg((pl.col("log_amount") - (pl.col("index_smoothed") / 100.0).log())90             .mean().exp().round(0).alias("basket_value"))91    )92    canon = m.join(basket, on=["geography_id", "property_type"], how="left") \93        .with_columns((pl.col("basket_value") * pl.col("index_smoothed") / 100.0)94                      .round(0).alias("representative_value"))9596    # ---------------- growth horizons ---------------- #97    keys = ["geography_id", "property_type"]98    canon = canon.sort([*keys, "period"]).with_columns([99        ((pl.col("index_smoothed") / pl.col("index_smoothed").shift(h).over(keys) - 1) * 100)100        .round(3).alias(name)101        for name, h in GROWTH.items()102    ])103104    # ---------------- reliability grades ---------------- #105    med = canon.group_by(keys).agg(106        pl.col("transactions").fill_nan(None).median().alias("cell_median_tx"))107    canon = canon.join(med, on=keys, how="left")108    base_grade = (109        pl.when((pl.col("cell_median_tx") >= 150) & (pl.col("se_log") < 0.015)).then(pl.lit("A"))110        .when((pl.col("cell_median_tx") >= 75) & (pl.col("se_log") < 0.025)).then(pl.lit("B"))111        .when(pl.col("cell_median_tx") >= 40).then(pl.lit("C"))112        .when(pl.col("cell_median_tx") >= 15).then(pl.lit("D"))113        .otherwise(pl.lit("E"))114    )115    downgrade = {"A": "B", "B": "C", "C": "D", "D": "E", "E": "E"}116    canon = canon.with_columns(117        pl.when(pl.col("transactions").fill_nan(0) == 0)118        .then(base_grade.replace(downgrade)).otherwise(base_grade)119        .alias("reliability_grade"))120121    # ---------------- partial month (registration lag) ---------------- #122    max_date = cl["date"].max()123    last_period = max_date.strftime("%Y-%m")124    month_end = (dt.date(max_date.year + (max_date.month == 12),125                         max_date.month % 12 + 1, 1) - dt.timedelta(days=1))126    is_partial_last = max_date < month_end127    canon = canon.with_columns(128        ((pl.col("period") == last_period) & is_partial_last).alias("is_partial_month"))129130    # ---------------- effective sample size ---------------- #131    fp = pl.read_parquet(INTERIM_DIR / "rtd_residuals.parquet")132    sigma2 = {t: float(fp.filter(pl.col("propertyType") == t)["resid"].var())133              for t in PROPERTY_TYPES}134    sigma2["all"] = float(np.mean(list(sigma2.values())))135    canon = canon.with_columns(136        pl.col("property_type").replace_strict(sigma2, default=None).alias("_s2")137    ).with_columns(138        (pl.col("_s2") / (pl.col("se_log") ** 2).clip(1e-9, None))139        .round(1).alias("effective_sample_size")140    ).drop("_s2")141142    canon = canon.with_columns(143        pl.lit(MODEL_VERSION).alias("model_version"),144        pl.lit(data_vintage).alias("data_vintage"),145        pl.lit("monthly").alias("frequency"),146        pl.col("transactions").fill_nan(0).cast(pl.Int64),147    )148    cols = ["period", "geography_level", "geography_id", "geography_name",149            "property_type", "index", "index_smoothed", "representative_value",150            "transactions", "effective_sample_size",151            "monthly_pct", "three_month_pct", "six_month_pct", "yoy_pct",152            "lower_95", "upper_95", "reliability_grade", "shrinkage_weight",153            "is_partial_month", "model_version", "data_vintage"]154    canonical = canon.select(cols).sort(155        ["geography_level", "geography_id", "property_type", "period"])156    canonical.write_parquet(CANONICAL)157    print(f"Canonical: {canonical.height:,} rows -> {CANONICAL}")158159    # ---------------- liquidity + coverage ---------------- #160    liq = (tx_all.group_by([*keys, "period"]).len().rename({"len": "transactions"})161           .sort([*keys, "period"]))162    liq.write_parquet(PROCESSED_DIR / "monthly_liquidity.parquet")163164    months_total = cl["period"].n_unique()165    cov = (166        cl.group_by(["municipality", "geo_code", "region", "propertyType", "period"]).len()167        .group_by(["municipality", "geo_code", "region", "propertyType"])168        .agg(pl.col("len").median().alias("median_monthly_tx"),169             pl.col("len").sum().alias("total_tx"),170             (pl.len() / months_total * 100).round(1).alias("pct_months_active"))171        .with_columns(172            pl.when(pl.col("municipality").is_in(list(CITY_TARGETS)))173            .then(pl.lit("published"))174            .when((pl.col("median_monthly_tx") >= 15) & (pl.col("pct_months_active") >= 85))175            .then(pl.lit("conditional"))176            .otherwise(pl.lit("not_published"))177            .alias("coverage_status"))178        .sort("total_tx", descending=True)179    )180    cov.write_csv(TABLES_DIR / "coverage_monthly.csv")181    cov.write_parquet(PROCESSED_DIR / "coverage_matrix_monthly.parquet")182    print("Coverage:", dict(cov.group_by("coverage_status").len().iter_rows()))183184    # ---------------- assessment gap (monthly) ---------------- #185    gap = cl.filter(pl.col("totalArValue") > 0).with_columns(186        (pl.col("amount") / pl.col("totalArValue")).alias("ratio"))187    gap_cells = pl.concat([188        gap.select(pl.lit("quebec").alias("geography_id"),189                   pl.col("propertyType").alias("property_type"), "period", "ratio"),190        gap.select(pl.concat_str([pl.lit("region-"), pl.col("region_code")])191                   .alias("geography_id"),192                   pl.col("propertyType").alias("property_type"), "period", "ratio"),193    ])194    (gap_cells.group_by([*keys, "period"] if False else ["geography_id", "property_type", "period"])195     .agg(pl.col("ratio").median().round(4).alias("median_ratio"), pl.len().alias("n"))196     .sort(["geography_id", "property_type", "period"])197     .write_parquet(PROCESSED_DIR / "assessment_gap_monthly.parquet"))198199    # ---------------- vintages (monthly store) ---------------- #200    vd = VINTAGE_DIR / f"data_vintage={data_vintage}"201    vd.mkdir(parents=True, exist_ok=True)202    canonical.write_parquet(vd / "qhpi_monthly.parquet")203    fresh = canonical.select(204        *KEY,205        pl.col("index").alias("first_release_index"),206        pl.col("index_smoothed").alias("first_release_index_smoothed"),207        pl.lit(data_vintage).alias("first_release_vintage"))208    if FIRST_RELEASE.exists():209        existing = pl.read_parquet(FIRST_RELEASE)210        fresh = pl.concat([existing, fresh.join(existing.select(KEY), on=KEY, how="anti")])211    fresh.write_parquet(FIRST_RELEASE)212213    latest_period = canonical.filter(~pl.col("is_partial_month"))["period"].max()214    canonical.filter(pl.col("period") == latest_period) \215        .write_csv(TABLES_DIR / "index_latest_monthly.csv")216    print(f"Latest complete month: {latest_period}")217218219if __name__ == "__main__":220    main()221