TeX 55.1%
Python 44.2%
1#!/usr/bin/env python32# Author: Simon-Pierre Boucher — simon-pierre.boucher@uqo.ca3"""Phase 6 — comparison with external indices (Teranet–National Bank, CREA MLS HPI), harmonised to 2021 average = 100.4SHARP metropolitan aggregates are stock-weighted bottom-up combinations of municipal cells inside the CMM (Montréal) and5CMQ (Québec) metropolitan communities (SDA `comet_s`), the closest available match to the CMA used by Teranet/CREA.6Writes results/validation/external_comparison.csv and results/validation/external_series.csv."""7from __future__ import annotations89import sys10from pathlib import Path1112sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src"))13import numpy as np14import polars as pl1516from sharp_hpi import config as C1718OUT = C.RESULTS / "validation"19EXT = C.EXTERNAL / "external_indices"202122def rebase(df: pl.DataFrame, col: str = "value") -> pl.DataFrame:23 b = df.filter(pl.col("period").str.starts_with("2021")).select(pl.col(col).log().mean().alias("b"))["b"][0]24 return df.with_columns((pl.col(col).log() - b).alias("li"))252627def main() -> None:28 panel = pl.read_parquet(C.PROCESSED / "index_panel.parquet")29 geo = pl.read_parquet(C.GEO_TX, columns=["mun_code", "comet_code"]).unique().filter(pl.col("comet_code").is_not_null())30 stock = pl.read_csv(OUT / "stock_weights_by_municipality_class.csv")31 series = []32 for cls in ["house", "condo"]:33 mun = panel.filter((pl.col("geography_level") == "municipality") & (pl.col("property_type") == cls) & ~pl.col("partial_month_flag")).with_columns(np.log(pl.col("index_value") / 100).alias("li"))34 mun = mun.join(geo, left_on="geography_id", right_on="mun_code", how="inner").join(stock.filter(pl.col("cls") == cls).select(["geography_id", "stock_units"]), on="geography_id", how="inner")35 for com, nm in [("663", "Montréal (CMM)"), ("235", "Québec (CMQ)")]:36 a = mun.filter(pl.col("comet_code") == com).group_by("month").agg(((pl.col("li") * pl.col("stock_units")).sum() / pl.col("stock_units").sum()).alias("li")).sort("month")37 series.append(a.select(pl.lit("SHARP-HPI").alias("source"), pl.lit(f"SHARP {cls} stock-weighted").alias("series_name"), pl.lit(nm).alias("geography"), pl.col("month").alias("period"), pl.col("li")))38 p = panel.filter((pl.col("geography_level") == "province") & (pl.col("property_type") == cls) & ~pl.col("partial_month_flag")).sort("month")39 series.append(p.select(pl.lit("SHARP-HPI").alias("source"), pl.lit(f"SHARP {cls} province").alias("series_name"), pl.lit("Quebec (province)").alias("geography"), pl.col("month").alias("period"), np.log(pl.col("index_value") / 100).alias("li")))40 ext = []41 t = pl.read_csv(EXT / "teranet_nb_hpi_monthly.csv").filter(pl.col("series_name").str.contains("not seasonally")).filter(pl.col("geography").is_in(["Montréal CMA", "Québec CMA"]))42 for g, sub in t.group_by("geography"):43 ext.append(rebase(sub.sort("period")).select(pl.lit("Teranet-NB").alias("source"), pl.lit("Teranet NSA (repeat sales)").alias("series_name"), pl.col("geography"), "period", "li"))44 c = pl.read_csv(EXT / "crea_mls_hpi_monthly.csv").filter(pl.col("series_name").str.contains("NSA") & pl.col("series_name").str.ends_with("_HPI"))45 for g, sn in [("Montréal CMA", "Single_Family"), ("Montréal CMA", "Apartment"), ("Montréal CMA", "Composite"), ("Québec CMA", "Single_Family"), ("Québec CMA", "Apartment"), ("Québec CMA", "Composite")]:46 sub = c.filter((pl.col("geography") == g) & pl.col("series_name").str.contains(sn))47 if sub.height:48 ext.append(rebase(sub.sort("period")).select(pl.lit("CREA MLS HPI").alias("source"), pl.lit(f"CREA {sn.replace('_', ' ')} NSA").alias("series_name"), pl.col("geography"), "period", "li"))49 qc = c.filter(pl.col("geography").str.contains("Quebec") & ~pl.col("geography").str.contains("CMA"))50 for sn in ["Single_Family", "Apartment", "Composite"]:51 sub = qc.filter(pl.col("series_name").str.contains(sn))52 if sub.height:53 ext.append(rebase(sub.sort("period")).select(pl.lit("CREA MLS HPI").alias("source"), pl.lit(f"CREA {sn.replace('_', ' ')} NSA").alias("series_name"), pl.lit("Quebec (province)").alias("geography"), "period", "li"))54 allser = pl.concat(series + ext); allser.write_csv(OUT / "external_series.csv")55 # pairwise comparison: SHARP vs each external in the same metro / province56 rows = []57 pairs = [("SHARP house stock-weighted", "Montréal (CMM)", "Montréal CMA"), ("SHARP condo stock-weighted", "Montréal (CMM)", "Montréal CMA"),58 ("SHARP house stock-weighted", "Québec (CMQ)", "Québec CMA"), ("SHARP condo stock-weighted", "Québec (CMQ)", "Québec CMA"),59 ("SHARP house province", "Quebec (province)", "Quebec (province)"), ("SHARP condo province", "Quebec (province)", "Quebec (province)")]60 for sname, sgeo, egeo in pairs:61 s = allser.filter((pl.col("series_name") == sname) & (pl.col("geography") == sgeo)).select("period", pl.col("li").alias("li_s"))62 for (src, ename), e in allser.filter((pl.col("geography") == egeo) & (pl.col("source") != "SHARP-HPI")).group_by(["source", "series_name"]):63 j = s.join(e.select("period", pl.col("li").alias("li_e")), on="period").sort("period")64 if j.height < 24:65 continue66 a, b = j["li_s"].to_numpy(), j["li_e"].to_numpy()67 d12a, d12b = a[12:] - a[:-12], b[12:] - b[:-12]; d1a, d1b = np.diff(a), np.diff(b)68 rows.append({"sharp_series": sname, "sharp_geo": sgeo, "external_source": src, "external_series": ename, "external_geo": egeo, "months": j.height, "last_period": j["period"][-1],69 "rmse_level": float(np.sqrt(np.mean((a - b) ** 2))), "max_abs_gap": float(np.max(np.abs(a - b))), "total_change_sharp": float(a[-1] - a[0]), "total_change_external": float(b[-1] - b[0]),70 "corr_12m": float(np.corrcoef(d12a, d12b)[0, 1]), "corr_1m": float(np.corrcoef(d1a, d1b)[0, 1]), "sd_1m_sharp": float(np.std(d1a)), "sd_1m_external": float(np.std(d1b))})71 out = pl.DataFrame(rows); out.write_csv(OUT / "external_comparison.csv"); print(out)727374if __name__ == "__main__" and len(sys.argv) == 1:75 main()767778def yearly_and_repeat_sales() -> None:79 """Year-by-year gap decomposition and our own Case–Shiller repeat sales on the CMM / island of Montréal, for comparison with Teranet."""80 from sharp_hpi import estimators as E81 from sharp_hpi.hedonic import build_features82 f = build_features(pl.read_parquet(C.ANALYTIC_TX))83 geo = pl.read_parquet(C.GEO_TX, columns=["id", "comet_code"]).rename({"id": "tx_id"})84 f = f.join(geo, on="tx_id", how="left")85 grid = E.month_grid("2021-01", "2026-10")86 rows = []87 for nm, sub in [("Montréal (CMM) houses", f.filter((pl.col("comet_code") == "663") & (pl.col("property_class") == "house"))),88 ("Montréal (CMM) condos", f.filter((pl.col("comet_code") == "663") & (pl.col("property_class") == "condo"))),89 ("Montréal (CMM) all three classes", f.filter(pl.col("comet_code") == "663")),90 ("Island of Montréal (region 06) all", f.filter(pl.col("region_code") == "06")),91 ("Québec (CMQ) all three classes", f.filter(pl.col("comet_code") == "235")),92 ("Province houses", f.filter(pl.col("property_class") == "house")), ("Province all three classes", f)]:93 for hold in (6, 12):94 rs = E.repeat_sales(sub, grid, weighted=True, min_hold_months=hold); li = rs.log_index[:-1]; okm = ~np.isnan(li)95 med = E.median_index(sub, grid).log_index[:-1]96 yr = {}97 for y in range(2021, 2027):98 idx = [i for i, m in enumerate(grid[:-1]) if m.startswith(str(y))]99 yr[y] = (float(np.nanmean(li[idx])), float(np.nanmean(med[idx])))100 rows.append({"sample": nm, "min_hold_months": hold, "pairs": rs.meta.get("pairs", 0), "rs_total_change": float(li[okm][-1] - li[okm][0]), "median_total_change": float(med[-1] - med[0]),101 **{f"rs_{y}_vs_2021": yr[y][0] - yr[2021][0] for y in range(2022, 2027)}, **{f"median_{y}_vs_2021": yr[y][1] - yr[2021][1] for y in range(2022, 2027)}})102 pl.DataFrame(rows).write_csv(OUT / "external_repeat_sales_check.csv"); print(pl.DataFrame(rows))103 # yearly averages of every series for the gap decomposition104 s = pl.read_csv(OUT / "external_series.csv").with_columns(pl.col("period").str.slice(0, 4).alias("year")).filter(pl.col("period") <= "2026-09")105 y = s.group_by(["source", "series_name", "geography", "year"]).agg(pl.col("li").mean()).sort(["source", "series_name", "geography", "year"])106 y.pivot(on="year", index=["source", "series_name", "geography"], values="li").write_csv(OUT / "external_yearly.csv")107108109if __name__ == "__main__" and len(sys.argv) > 1 and sys.argv[1] == "yearly":110 yearly_and_repeat_sales()111