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/07_canonical.py7# Purpose : Step 7 — assemble data/processed/qwhpi_weekly.parquet (canonical8#           dataset), coverage matrix, assessment gap, vintages, liquidity.9# =============================================================================10"""Canonical dataset build (Execution Order step 7, schema §4.9)."""1112from __future__ import annotations1314import datetime as dt15import sys16import unicodedata17from pathlib import Path1819sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src"))2021import numpy as np22import polars as pl2324from qwhpi.clean import CLEAN_PARQUET, research_sample25from qwhpi.config import (26    CITY_TARGETS, INTERIM_DIR, PROCESSED_DIR, PROPERTY_TYPES, TABLES_DIR,27    ensure_dirs,28)29from qwhpi.index import (30    add_growth, add_reliability, add_representative_value, compute_basket,31)32from qwhpi.nowcast import partial_week_flags33from qwhpi.seasonal import seasonality_test34from qwhpi.vintages import snapshot, update_first_release3536MODEL_VERSION = "0.1.0"37CANONICAL = PROCESSED_DIR / "qwhpi_weekly.parquet"383940def slug(name: str) -> str:41    s = unicodedata.normalize("NFKD", name).encode("ascii", "ignore").decode()42    return s.lower().replace(" ", "-").replace("'", "")434445def city_slug(name: str) -> str:46    return "quebec-city" if name == "Québec" else slug(name)474849def main() -> None:50    ensure_dirs()51    data_vintage = dt.date.today().isoformat()52    hier = pl.read_parquet(INTERIM_DIR / "hierarchical_indexes.parquet")53    cl = research_sample(pl.read_parquet(CLEAN_PARQUET))54    cl = cl.with_columns(pl.col("week").dt.strftime("%Y-%m-%d").alias("week_str"))5556    # ------------------------------------------------------------------ #57    # Map transactions to geography_ids (for baskets & liquidity)58    # ------------------------------------------------------------------ #59    cell_frames = []60    for gid_expr, filt in [61        (pl.lit("quebec"), pl.lit(True)),62        (pl.concat_str([pl.lit("region-"), pl.col("region_code")]),63         pl.col("region_code").is_not_null()),64    ]:65        cell_frames.append(66            cl.filter(filt).select(67                gid_expr.alias("geography_id"),68                pl.col("propertyType").alias("property_type"),69                pl.col("week_str").alias("week"),70                "log_amount", "amount",71            )72        )73    for city in CITY_TARGETS:74        sub = cl.filter(pl.col("municipality") == city)75        if sub.height:76            cell_frames.append(sub.select(77                pl.lit(city_slug(city)).alias("geography_id"),78                pl.col("propertyType").alias("property_type"),79                pl.col("week_str").alias("week"),80                "log_amount", "amount",81            ))82    tx_cells = pl.concat(cell_frames)83    tx_cells_all = pl.concat([84        tx_cells,85        tx_cells.with_columns(pl.lit("all").alias("property_type")),86    ])8788    # ------------------------------------------------------------------ #89    # Canonical columns90    # ------------------------------------------------------------------ #91    basket = compute_basket(tx_cells_all, hier)92    canon = add_representative_value(hier, basket)93    canon = add_growth(canon)94    canon = add_reliability(canon)9596    prov_vol = (97        cl.group_by("week_str").len()98        .rename({"week_str": "week", "len": "volume"})99    )100    canon = canon.join(partial_week_flags(prov_vol), on="week", how="left") \101        .with_columns(pl.col("is_partial_week").fill_null(False))102103    # Effective sample size: information units of the filtered estimate.104    sigma2_map = {105        t: float(pl.read_parquet(INTERIM_DIR / "stage1_residuals.parquet")106                 .filter(pl.col("propertyType") == t)["resid"].var())107        for t in PROPERTY_TYPES108    }109    sigma2_map["all"] = float(np.mean(list(sigma2_map.values())))110    canon = canon.with_columns(111        pl.col("property_type").replace_strict(sigma2_map, default=None)112        .alias("_sigma2")113    ).with_columns(114        (pl.col("_sigma2") / (pl.col("se_log") ** 2).clip(1e-9, None))115        .round(1).alias("effective_sample_size")116    ).drop("_sigma2")117118    canon = canon.with_columns(119        pl.lit(MODEL_VERSION).alias("model_version"),120        pl.lit(data_vintage).alias("data_vintage"),121        pl.col("transactions").fill_nan(0).cast(pl.Int64),122    )123124    cols = ["week", "geography_level", "geography_id", "geography_name",125            "property_type", "index", "index_smoothed", "representative_value",126            "transactions", "effective_sample_size",127            "weekly_pct", "four_week_pct", "thirteen_week_pct",128            "twenty_six_week_pct", "yoy_pct",129            "lower_95", "upper_95", "reliability_grade", "shrinkage_weight",130            "is_partial_week", "model_version", "data_vintage"]131    canonical = canon.select(cols).sort(132        ["geography_level", "geography_id", "property_type", "week"])133    canonical.write_parquet(CANONICAL)134    print(f"Canonical table: {canonical.height:,} rows -> {CANONICAL}")135136    # ------------------------------------------------------------------ #137    # Coverage matrix (all municipalities × types)138    # ------------------------------------------------------------------ #139    weeks_total = cl["week_str"].n_unique()140    liq = (141        cl.group_by(["municipality", "geo_code", "region", "propertyType", "week_str"])142        .len()143        .group_by(["municipality", "geo_code", "region", "propertyType"])144        .agg(pl.col("len").median().alias("median_weekly_tx"),145             pl.col("len").sum().alias("total_tx"),146             (pl.len() / weeks_total * 100).round(1).alias("pct_weeks_active"))147    )148    published_cities = {city_slug(c) for c in CITY_TARGETS}149    liq = liq.with_columns(150        pl.when(pl.col("municipality").is_in(list(CITY_TARGETS)))151        .then(pl.lit("published"))152        .when((pl.col("median_weekly_tx") >= 5) & (pl.col("pct_weeks_active") >= 80))153        .then(pl.lit("conditional"))154        .otherwise(pl.lit("not_published"))155        .alias("coverage_status")156    ).sort("total_tx", descending=True)157    liq.write_csv(TABLES_DIR / "coverage.csv")158    liq.write_parquet(PROCESSED_DIR / "coverage_matrix.parquet")159    counts = liq.group_by("coverage_status").len()160    print("Coverage matrix:", dict(counts.iter_rows()))161162    # Weekly liquidity per published cell.163    weekly_liq = (164        tx_cells_all.group_by(["geography_id", "property_type", "week"]).len()165        .rename({"len": "transactions"}).sort(["geography_id", "property_type", "week"])166    )167    weekly_liq.write_parquet(PROCESSED_DIR / "weekly_liquidity.parquet")168    (169        weekly_liq.group_by(["geography_id", "property_type"])170        .agg(pl.col("transactions").median().alias("median"),171             pl.col("transactions").mean().round(1).alias("mean"))172        .sort(["geography_id", "property_type"])173        .write_csv(TABLES_DIR / "weekly_liquidity.csv")174    )175176    # ------------------------------------------------------------------ #177    # Assessment Gap Index (secondary module — never mixed with QWHPI)178    # ------------------------------------------------------------------ #179    gap = (180        cl.filter(pl.col("totalArValue") > 0)181        .with_columns((pl.col("amount") / pl.col("totalArValue")).alias("ratio"))182    )183    gap_cells = pl.concat([184        gap.select(pl.lit("quebec").alias("geography_id"),185                   pl.col("propertyType").alias("property_type"),186                   pl.col("week_str").alias("week"), "ratio"),187        gap.select(pl.concat_str([pl.lit("region-"), pl.col("region_code")])188                   .alias("geography_id"),189                   pl.col("propertyType").alias("property_type"),190                   pl.col("week_str").alias("week"), "ratio"),191    ])192    assessment = (193        gap_cells.group_by(["geography_id", "property_type", "week"])194        .agg(pl.col("ratio").median().round(4).alias("median_ratio"),195             pl.len().alias("n"))196        .sort(["geography_id", "property_type", "week"])197    )198    assessment.write_parquet(PROCESSED_DIR / "assessment_gap.parquet")199    print(f"Assessment gap: {assessment.height:,} rows")200201    # ------------------------------------------------------------------ #202    # Seasonality check (documents the NSA decision)203    # ------------------------------------------------------------------ #204    seas_rows = []205    for gid, ptype in [("quebec", "all"), ("quebec", "unifamilial"),206                       ("montreal", "condo")]:207        s = canonical.filter((pl.col("geography_id") == gid)208                             & (pl.col("property_type") == ptype)).sort("week")209        res = seasonality_test(np.log(s["index_smoothed"].to_numpy()))210        seas_rows.append({"cell": f"{gid}/{ptype}", **res})211    seas = pl.DataFrame(seas_rows)212    seas.write_csv(TABLES_DIR / "seasonality_tests.csv")213    print(seas)214215    # ------------------------------------------------------------------ #216    # Vintages + latest snapshot table217    # ------------------------------------------------------------------ #218    snapshot(canonical, data_vintage)219    update_first_release(canonical, data_vintage)220221    latest_week = canonical.filter(~pl.col("is_partial_week"))["week"].max()222    latest = canonical.filter(pl.col("week") == latest_week)223    latest.write_csv(TABLES_DIR / "index_latest.csv")224    print(f"Latest complete week: {latest_week} | {latest.height} series")225226227if __name__ == "__main__":228    main()229