TeX 55.1%
Python 44.2%
1#!/usr/bin/env python32# Author: Simon-Pierre Boucher — simon-pierre.boucher@uqo.ca3"""Phase 7 — LaTeX table fragments (paper/tables/*.tex) and number macros (paper/tables/numbers.tex) from results/.4Every number in the manuscript comes from these files; missing inputs are skipped with a note."""5from __future__ import annotations67import json8import re9import 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 C17from sharp_hpi.plotting import METHOD_LABELS1819R = C.RESULTS; T = C.ROOT / "paper" / "tables"; T.mkdir(parents=True, exist_ok=True)20MACROS: dict[str, str] = {}21WORDS = {"0": "Zero", "1": "One", "2": "Two", "3": "Three", "4": "Four", "5": "Five", "6": "Six", "7": "Seven", "8": "Eight", "9": "Nine"}222324def mname(name: str) -> str:25 s = re.sub(r"[^A-Za-z0-9]", "", name)26 return "".join(WORDS[ch] if ch.isdigit() else ch for ch in s)272829def macro(name: str, value, fmt: str = "{:,.1f}") -> None:30 if value is None or (isinstance(value, float) and np.isnan(value)):31 MACROS[mname(name)] = "n/a"32 else:33 MACROS[mname(name)] = fmt.format(value).replace(",", "{,}")343536def esc(s) -> str:37 return str(s).replace("&", "\\&").replace("%", "\\%").replace("_", "\\_").replace("−", "$-$").replace("→", "$\\to$").replace("#", "\\#").replace("≤", "$\\le$").replace("≥", "$\\ge$").replace("×", "$\\times$")383940def table(name: str, header: list[str], rows: list[list], caption: str, label: str, align: str | None = None, notes: str = "", size: str = "\\small", longtable: bool = False) -> None:41 align = align or "l" + "r" * (len(header) - 1)42 body = "\n".join(" & ".join(esc(c) for c in r) + " \\\\" for r in rows)43 if longtable:44 tex = f"""{{{size}45\\begin{{longtable}}{{{align}}}46\\caption{{{caption}}}\\label{{{label}}}\\\\47\\toprule48{' & '.join(esc(h) for h in header)} \\\\49\\midrule50\\endfirsthead51\\toprule52{' & '.join(esc(h) for h in header)} \\\\53\\midrule54\\endhead55\\bottomrule56\\endfoot57{body}58\\end{{longtable}}59}}60{('\\noindent{\\footnotesize ' + notes + '}') if notes else ''}61"""62 else:63 tex = f"""\\begin{{table}}[htbp]64\\centering{size}65\\caption{{{caption}}}\\label{{{label}}}66\\begin{{threeparttable}}67\\begin{{tabular}}{{{align}}}68\\toprule69{' & '.join(esc(h) for h in header)} \\\\70\\midrule71{body}72\\bottomrule73\\end{{tabular}}74{('\\begin{tablenotes}\\footnotesize\\item ' + notes + '\\end{tablenotes}') if notes else ''}75\\end{{threeparttable}}76\\end{{table}}77"""78 (T / f"{name}.tex").write_text(tex); print("table", name)798081def f1(x): return f"{x:.1f}"82def f2(x): return f"{x:.2f}"83def f3(x): return f"{x:.3f}"84def pct(x): return f"{100 * x:.1f}"85def n0(x): return f"{x:,.0f}".replace(",", "{,}")868788def ok(*p):89 miss = [str(x) for x in p if not Path(x).exists()]90 if miss:91 print("skip:", miss)92 return not miss939495# ---- data tables -----------------------------------------------------------96def t_data():97 man = json.loads((C.INTERIM / "raw_manifest.json").read_text())98 macro("NRawTx", man["tables"]["transactions"]["rows"], "{:,.0f}"); macro("NUnitsStock", man["tables"]["units_stock"]["rows"], "{:,.0f}"); macro("NUnitsMatched", man["tables"]["units_matched"]["rows"], "{:,.0f}")99 led = pl.read_csv(R / "cleaning" / "exclusion_ledger.csv")100 lab = {None: "Retained: headline sample", "class_not_in_headline": "Class not in headline (cottage, mobile, 6+ units, unclassified)", "assessment_ratio_outlier": "Price/assessment ratio outside ±log 4 of year–class–region median",101 "land": "Vacant land / road / forest (roll)", "non_residential": "Non-residential use code (roll)", "same_parcel_same_date_diff_price": "Same parcel, same day, different prices",102 "no_roll_characteristics": "No assessment-roll characteristics", "shared_parcel": "Parcel with more than 10 sales in 2021–2026", "built_after_sale": "Roll building year after the sale year",103 "multi_parcel_candidate": "Multi-parcel deed (same day, price, street or 1-km cell)", "exact_duplicate": "Exact duplicate row", "price_above_cap": "Price above \\$10 million"}104 rows = [[lab.get(r["exclude_reason_headline"], r["exclude_reason_headline"]), n0(r["n"]), pct(r["share_of_raw"])] for r in led.iter_rows(named=True)]105 n_head = led.filter(pl.col("exclude_reason_headline").is_null())["n"][0]; macro("NHeadline", n_head, "{:,.0f}"); macro("ShareHeadline", 100 * n_head / man["tables"]["transactions"]["rows"], "{:.1f}")106 table("tab_exclusions", ["Rule (first match, in priority order)", "Rows", "% of raw"], rows, "Population definition: exclusion ledger of the Vrai Prix snapshot (769,995 registered sales, 2021-01-04 to 2026-10-07)", "tab:exclusions",107 notes="Rules are applied in the order listed; a row is counted once, under the first rule it triggers. Independent flag counts are in results/cleaning/flag\\_counts.csv.")108 hs = pl.read_csv(R / "cleaning" / "headline_sample_by_year_class.csv").filter(pl.col("property_class").is_not_null())109 piv = hs.pivot(on="property_class", index="year", values="n").fill_null(0)110 rows = [[str(r["year"]), n0(r.get("house", 0)), n0(r.get("condo", 0)), n0(r.get("plex", 0)), n0(r.get("house", 0) + r.get("condo", 0) + r.get("plex", 0))] for r in piv.sort("year").iter_rows(named=True)]111 tot = [n0(piv[c].sum()) for c in ["house", "condo", "plex"]]112 rows.append(["Total", *tot, n0(sum(piv[c].sum() for c in ["house", "condo", "plex"]))])113 for c in ["house", "condo", "plex"]:114 macro(f"NHeadline{c}", piv[c].sum(), "{:,.0f}")115 table("tab_sample_by_year", ["Year", "Houses", "Condominiums", "Plexes", "Total"], rows, "Headline estimation sample by year and property class", "tab:sample_by_year", notes="2026 covers January 1 to October 7 (partial).")116 comp = pl.read_csv(R / "profile" / "completeness_by_year.csv")117 cols = [("floor_area", "Floor area"), ("year_built", "Year built"), ("land_area", "Lot area"), ("nb_etages", "Storeys"), ("lien_physique", "Physical link"), ("genre_construction", "Building genre"), ("building_type", "Provider style (banned)"), ("id_provinc", "Roll identifier")]118 rows = [[lab] + [pct(v) for v in comp[c].to_list()] for c, lab in cols]119 table("tab_completeness_supp", ["Variable (% non-missing)"] + [str(y) for y in comp["year"].to_list()], rows, "Completeness of characteristics by year of sale, all raw rows (including land, non-residential and unmatched sales)", "tab:completeness_supp",120 notes="The provider's building style is backfilled over time and is excluded from every model (time-correlated missingness). Lower shares in 2021 reflect the larger share of land and unclassified rows, not missing data in the headline sample (Table~\\ref{tab:completeness}).")121 hs2 = pl.read_csv(R / "cleaning" / "headline_sample_by_year_class.csv").filter(pl.col("property_class").is_not_null())122 rows = []123 for cls, lab in [("house", "Houses"), ("condo", "Condominiums"), ("plex", "Plexes")]:124 x = hs2.filter(pl.col("property_class") == cls).sort("year")125 rows.append([f"{lab}: floor area"] + [pct(1 - v) for v in x["share_floor_area_missing"].to_list()])126 rows.append([f"{lab}: year built"] + [pct(1 - v) for v in x["share_year_built_missing"].to_list()])127 rows.append([f"{lab}: new construction (\\% of sales)"] + [pct(v) for v in x["share_new_construction"].to_list()])128 table("tab_completeness", ["Headline sample (% of sales)"] + [str(y) for y in sorted(hs2["year"].unique().to_list())], rows, "Completeness of the two main characteristics and share of new constructions in the headline sample, by year and class", "tab:completeness",129 notes="Lot area, storeys, physical link and building genre are complete for more than 90\\% of headline sales in every year (results/profile/completeness\\_by\\_year.csv). The provider's building style (12\\% present in January 2021, 87\\% in January 2025) is excluded from every model. New constructions (sold in the roll's construction year) fall to zero in 2025--2026 because the roll vintage lags construction.")130 bands = pl.read_csv(R / "profile" / "municipality_size_bands.csv")131 table("tab_mun_bands", ["Sales 2021–2026", "Municipalities", "Sales"], [[r["band"], n0(r["municipalities"]), n0(r["transactions"])] for r in bands.iter_rows(named=True)],132 "Municipalities by number of roll-matched sales, 2021–2026", "tab:mun_bands")133 macro("NMunicipalities", bands["municipalities"].sum(), "{:,.0f}")134 v = pl.read_csv(R / "vintages" / "vintage_summary.csv")135 table("tab_vintages", ["Extraction", "Rows", "Last date", "Lag (days)", "Final rows ≤ last date", "Share present", "Deleted", "Price revised"],136 [[r["vintage"], n0(r["vintage_rows"]), r["vintage_last_date"], str(r["vintage_max_lag_days"]), n0(r["final_rows_dated_le_last"]), f"{100*r['share_present']:.2f}\\%", str(r["deleted_since_vintage"]), str(r["amount_revised"])] for r in v.iter_rows(named=True)],137 "Earlier extractions of the transaction table compared with the final snapshot (publication lag audit)", "tab:vintages")138 macro("PubLagMin", v["vintage_max_lag_days"].min(), "{:.0f}"); macro("PubLagMax", v["vintage_max_lag_days"].max(), "{:.0f}"); macro("SharePresentMin", 100 * v["share_present"].min(), "{:.2f}")139 geo = pl.read_csv(R / "geo" / "roll_vs_spatial_agreement.csv"); macro("RollSpatialAgree", 100 * geo["share_roll_code_agrees"][0], "{:.2f}")140 fs = pl.read_csv(R / "geo" / "fine_location_summary.csv"); d = dict(zip(fs["metric"], fs["value"]))141 macro("NGridOneKm", d["distinct_grid1k"], "{:,.0f}"); macro("NMrc", d["distinct_mrc"], "{:.0f}"); macro("NMunWithSales", d["distinct_municipalities"], "{:,.0f}")142 cov = pl.read_csv(R / "profile" / "coverage_month_type.csv").group_by("month").agg(pl.col("n").sum()).sort("month")143 macro("NJanTwentyOne", cov["n"][0], "{:,.0f}"); macro("NOctTwentySix", cov["n"][-1], "{:,.0f}"); macro("AvgMonthly", cov.filter(pl.col("month").str.slice(0, 4).is_in(["2022", "2023", "2024", "2025"]))["n"].mean(), "{:,.0f}")144145146def t_stage1():147 if not ok(R / "panel" / "stage1_meta.json"):148 return149 meta = json.loads((R / "panel" / "stage1_meta.json").read_text())150 keys = [("lfa", "log floor area (centred)"), ("lfa2", "log floor area$^2$"), ("lfa3", "log floor area$^3$"), ("lla", "log lot area (centred)"), ("lla2", "log lot area$^2$"),151 ("C(age_bin)[T.a01_2-5]", "Age 2–5"), ("C(age_bin)[T.a03_11-20]", "Age 11–20"), ("C(age_bin)[T.a04_21-40]", "Age 21–40"), ("C(age_bin)[T.a06_61-80]", "Age 61–80"), ("C(age_bin)[T.a08_120+]", "Age 120+"),152 ("C(storeys_c)[T.2]", "Two storeys"), ("C(subtype_c)[T.semi]", "Semi-detached"), ("C(subtype_c)[T.row]", "Row"), ("C(units_c)[T.3]", "Triplex"), ("C(units_c)[T.5]", "5 units"),153 ("lfa:tt", "log floor area × time"), ("lla:tt", "log lot area × time"), ("fa_missing", "Floor area missing"), ("age_missing", "Age missing")]154 rows = [[lab] + [f3(meta[c]["coef"][k]) if k in meta[c]["coef"] else "—" for c in ["house", "condo", "plex"]] for k, lab in keys]155 rows += [["Observations"] + [n0(meta[c]["n"]) for c in ["house", "condo", "plex"]], ["Robust residual s.d."] + [f3(meta[c]["sigma_robust"]) for c in ["house", "condo", "plex"]],156 ["Within $R^2$"] + [f3(meta[c]["r2_within"]) for c in ["house", "condo", "plex"]]]157 table("tab_stage1", ["Regressor", "Houses", "Condominiums", "Plexes"], rows, "Stage-1 robust hedonic regressions (selected coefficients; dependent variable log price)", "tab:stage1",158 notes="Huber-IRLS ($k=1.345$) with fine-location fixed effects (nested metric grid, finest cell with at least 15 sales) and month effects; time index $tt$ runs from $-0.5$ (January 2021) to $0.5$ (October 2026). Full coefficient lists in results/panel/stage1\\_meta.json.")159 for c in ["house", "condo", "plex"]:160 macro(f"SigmaRobust{c}", meta[c]["sigma_robust"], "{:.3f}"); macro(f"RtwoWithin{c}", meta[c]["r2_within"], "{:.2f}"); macro(f"BetaFa{c}", meta[c]["coef"]["lfa"], "{:.3f}"); macro(f"BetaFaTT{c}", meta[c]["coef"]["lfa:tt"], "{:.3f}")161 for lvl in ["region", "mrc", "municipality"]:162 m, v = meta[c]["priors"][lvl]; macro(f"PriorTau{c}{lvl}", 100 * np.exp(m / 2), "{:.2f}"); macro(f"PriorV{c}{lvl}", v, "{:.1f}")163 rows = [[f"{c} / {lvl}", f"{100*np.exp(meta[c]['priors'][lvl][0]/2):.2f}", f"{meta[c]['priors'][lvl][1]:.2f}"] for c in ["house", "condo", "plex"] for lvl in ["region", "mrc", "municipality"]]164 table("tab_priors", ["Class / level", "Prior median $\\tau$ (\\% per month)", "Prior variance of $\\log\\tau^2$"], rows, "Empirical-Bayes priors for the innovation standard deviation of cell deviations", "tab:priors",165 notes="Grid EM over all cells of a level with at least six observed months; $\\tau$ is the monthly standard deviation of the latent deviation from the parent path.")166167168def t_panel():169 if not ok(R / "panel" / "grade_distribution.csv", R / "panel" / "reliability_cells.csv"):170 return171 g = pl.read_csv(R / "panel" / "grade_distribution.csv")172 piv = g.pivot(on="reliability_grade", index=["geography_level", "property_type"], values="cell_months").fill_null(0)173 rows = []174 for r in piv.sort(["geography_level", "property_type"]).iter_rows(named=True):175 tot = sum(r.get(k, 0) for k in "ABCDE")176 rows.append([f"{r['geography_level']} / {r['property_type']}"] + [f"{100*r.get(k,0)/tot:.0f}" for k in "ABCDE"] + [n0(tot)])177 table("tab_grades", ["Level / class", "A", "B", "C", "D", "E", "Cell-months"], rows, "Distribution of reliability grades across cell-months (percent of cell-months)", "tab:grades",178 notes="A: effective sample size $\\ge 150$ and posterior s.e. of the log level $< 1.5$\\%; B: $\\ge 75$ and $< 2.5$\\%; C: $\\ge 30$ and $< 3.5$\\%; D: $\\ge 15$ and $< 5$\\%; E otherwise (published for transparency, flagged).")179 rc = pl.read_csv(R / "panel" / "reliability_cells.csv", schema_overrides={"geography_id": pl.Utf8})180 macro("NCells", rc.height, "{:,.0f}"); macro("NCellsMun", rc.filter(pl.col("geography_level") == "municipality").height, "{:,.0f}")181 gm = rc.filter(pl.col("geography_level") == "municipality").group_by("property_type").agg(pl.col("mean_gain").median().alias("g"), pl.col("median_monthly_n").median().alias("n"))182 for r in gm.iter_rows(named=True):183 macro(f"MedianGainMun{r['property_type']}", r["g"], "{:.2f}"); macro(f"MedianMonthlyNMun{r['property_type']}", r["n"], "{:.1f}")184 panel = pl.read_parquet(C.PROCESSED / "index_panel.parquet")185 for cls in ["house", "condo", "plex"]:186 p = panel.filter((pl.col("geography_level") == "province") & (pl.col("property_type") == cls)).sort("month")187 macro(f"ProvSep{cls}", p.filter(pl.col("month") == "2026-09")["index_value"][0], "{:.1f}"); macro(f"ProvYoySep{cls}", p.filter(pl.col("month") == "2026-09")["annual_pct_change"][0], "{:.1f}")188 macro(f"ProvPeak{cls}", p["index_value"].max(), "{:.1f}"); macro(f"ProvPeakMonth{cls}", p.filter(pl.col("index_value") == p["index_value"].max())["month"][0], "{}")189 # municipal dispersion of 5-year change (houses), September 2026, grades A–C190 m = panel.filter((pl.col("geography_level") == "municipality") & (pl.col("property_type") == "house") & (pl.col("month") == "2026-09") & pl.col("reliability_grade").is_in(["A", "B", "C"]))191 macro("NMunHouseABCSep", m.height, "{:,.0f}"); macro("MunHousePTen", m["index_value"].quantile(0.1), "{:.0f}"); macro("MunHousePFifty", m["index_value"].quantile(0.5), "{:.0f}"); macro("MunHousePNinety", m["index_value"].quantile(0.9), "{:.0f}")192 big = panel.filter((pl.col("geography_level") == "municipality") & (pl.col("month") == "2026-09")).join(193 panel.filter((pl.col("geography_level") == "municipality") & (pl.col("month") == "2026-09")).group_by("geography_id").agg(pl.col("actual_transaction_count").sum()), on="geography_id")194 tops = rc.filter(pl.col("geography_level") == "municipality").group_by("geography_id").agg(pl.col("n_total").sum().alias("n")).sort("n", descending=True).head(12)["geography_id"].to_list()195 rows = []196 for gid in tops:197 r = [panel.filter((pl.col("geography_id") == gid) & (pl.col("geography_level") == "municipality"))["geography_name"][0]]198 for cls in ["house", "condo", "plex"]:199 s = panel.filter((pl.col("geography_id") == gid) & (pl.col("geography_level") == "municipality") & (pl.col("property_type") == cls) & (pl.col("month") == "2026-09"))200 r += [f"{s['index_value'][0]:.1f} [{s['lower_95'][0]:.0f}, {s['upper_95'][0]:.0f}] {s['reliability_grade'][0]}" if s.height else "—"]201 rows.append(r)202 table("tab_top_municipalities", ["Municipality", "Houses", "Condominiums", "Plexes"], rows, "SHARP-HPI, September 2026 (2021 average = 100), twelve largest municipalities by sales: index [95\\% interval] grade", "tab:top_mun", size="\\footnotesize")203204205def t_benchmarks():206 if not ok(R / "benchmarks" / "series.parquet"):207 return208 s = pl.read_parquet(R / "benchmarks" / "series.parquet")209 if (R / "benchmarks" / "series_hmts.parquet").exists():210 s = pl.concat([s, pl.read_parquet(R / "benchmarks" / "series_hmts.parquet")])211 methods = ["A_median", "A_mean", "B_stratified_median", "C_pooled_td_ols", "C_pooled_td", "D_rtd_w7_mean", "D_rtd_w13_mean", "D_rtd_w13_movement", "D_rtd_w13_window", "D_rtd_w25_mean", "E_hi_laspeyres_chained", "E_hi_paasche_chained",212 "E_hi_fisher_chained", "E_hi_fisher_direct", "F_geks_w13_movement", "F_geks_w25_movement", "G_hmts_w4", "G_hmts_w13", "H_repeat_sales_bmn", "H_repeat_sales_cs", "J_gbm_oot_residual", "K_qhpi", "S_sharp_notv", "S_sharp"]213 labels = {**METHOD_LABELS, "A_mean": "Mean price", "C_pooled_td_ols": "Pooled time dummy (OLS)", "D_rtd_w7_mean": "Rolling TD, 7 m, mean splice", "D_rtd_w13_movement": "Rolling TD, 13 m, movement splice",214 "D_rtd_w13_window": "Rolling TD, 13 m, window splice", "D_rtd_w25_mean": "Rolling TD, 25 m, mean splice", "E_hi_laspeyres_chained": "Hedonic imputation Laspeyres (chained)", "E_hi_paasche_chained": "Hedonic imputation Paasche (chained)",215 "E_hi_fisher_direct": "Hedonic imputation Fisher (direct)", "F_geks_w13_movement": "GEKS-Fisher, 13 m", "G_hmts_w4": "HMTS, 4 m", "H_repeat_sales_bmn": "Repeat sales (BMN)", "S_sharp_notv": "SHARP stage 1 without time-varying prices"}216 rows = []; tim = pl.read_csv(R / "benchmarks" / "timing.csv") if (R / "benchmarks" / "timing.csv").exists() else None217 for m in methods:218 r = [labels.get(m, m)]219 for cls in ["house", "condo", "plex"]:220 sub = s.filter((pl.col("class") == cls) & (pl.col("level") == "province") & (pl.col("method") == m)).sort("month")221 if sub.height == 0:222 r += ["—", "—"]; continue223 li = sub["log_index"].to_numpy(); li = li[:len(li) - 1] # drop partial month224 okm = ~np.isnan(li); tot = li[okm][-1] - li[okm][0]; sd = np.nanstd(np.diff(li))225 r += [f"{100*tot:.1f}", f"{100*sd:.2f}"]226 macro(f"Prov{m}{cls}Total", 100 * tot, "{:.1f}")227 t = tim.filter((pl.col("method") == m) & (pl.col("level") == "province") & (pl.col("class") == "house"))["seconds"].mean() if tim is not None else None228 r += [f"{t:.0f}" if t is not None and t == t else "—"]229 rows.append(r)230 table("tab_benchmarks_province", ["Method", "H: total", "H: s.d. $\\Delta$", "C: total", "C: s.d. $\\Delta$", "P: total", "P: s.d. $\\Delta$", "sec."], rows,231 "Province-level benchmark indices, January 2021 to September 2026: cumulative log change (×100) and standard deviation of monthly log changes (×100)", "tab:bench_prov", size="\\footnotesize",232 notes="H = houses, C = condominiums, P = plexes. Rolling TD = rolling time dummy; GEKS on bilateral double-imputation Fisher links with movement splice; HMTS = hedonic multilateral time series (Ishaak et al., 2024); QHPI = re-implementation of the QWHPI v2.1 architecture; last column = wall-clock seconds for the house sample on one Apple M-series core.")233 # region-level agreement table: RMSE of monthly changes vs SHARP and vs repeat sales, by method, averaged over regions234 rows = []235 for m in ["A_median", "B_stratified_median", "C_pooled_td", "D_rtd_w13_mean", "E_hi_fisher_chained", "F_geks_w25_movement", "G_hmts_w13", "K_qhpi_direct_local_td", "S_sharp_filtered", "S_sharp"]:236 r = [labels.get(m, "QHPI v2.1 architecture" if m.startswith("K_") else m)]237 for cls in ["house", "condo", "plex"]:238 sub = s.filter((pl.col("class") == cls) & (pl.col("level") == "region") & pl.col("method").is_in([m, "H_repeat_sales_cs", "S_sharp"]))239 vals_rs, vals_sh, sds = [], [], []240 for reg, g in sub.group_by("geo_id"):241 a = g.filter(pl.col("method") == m).sort("month")["log_index"].to_numpy()242 b = g.filter(pl.col("method") == "H_repeat_sales_cs").sort("month")["log_index"].to_numpy()243 c = g.filter(pl.col("method") == "S_sharp").sort("month")["log_index"].to_numpy()244 if len(a) == 0:245 continue246 d12a = a[12:] - a[:-12]; sds.append(np.nanstd(np.diff(a)))247 if len(b):248 d12b = b[12:] - b[:-12]; okk = ~np.isnan(d12a) & ~np.isnan(d12b)249 if okk.sum() > 12:250 vals_rs.append(np.sqrt(np.mean((d12a[okk] - d12b[okk]) ** 2)))251 if len(c) and m != "S_sharp":252 d12c = c[12:] - c[:-12]; okk = ~np.isnan(d12a) & ~np.isnan(d12c)253 if okk.sum() > 12:254 vals_sh.append(np.sqrt(np.mean((d12a[okk] - d12c[okk]) ** 2)))255 r += [f"{100*np.nanmean(sds):.2f}" if sds else "—", f"{100*np.nanmean(vals_rs):.1f}" if vals_rs else "—"]256 macro(f"RegSd{m}{cls}", 100 * np.nanmean(sds) if sds else None, "{:.2f}"); macro(f"RegRmseRS{m}{cls}", 100 * np.nanmean(vals_rs) if vals_rs else None, "{:.1f}")257 rows.append(r)258 table("tab_benchmarks_region", ["Method", "H: s.d. $\\Delta$", "H: RMSE vs RS", "C: s.d. $\\Delta$", "C: RMSE vs RS", "P: s.d. $\\Delta$", "P: RMSE vs RS"], rows,259 "Region-level benchmarks (16 administrative regions): mean standard deviation of monthly log changes (×100) and RMSE of 12-month log changes against the Case–Shiller repeat-sales index (×100)", "tab:bench_region", size="\\footnotesize",260 notes="Repeat sales are an independent (pair-based) arbiter, not the truth: they carry renovation and selection biases. Lower s.d. is not better per se; the Monte Carlo (Section 7) separates noise from signal.")261262263def t_simulation():264 if not ok(R / "simulation" / "summary.csv"):265 return266 d = pl.read_csv(R / "simulation" / "summary.csv")267 methods = ["A_median", "B_stratified", "C_pooled_td", "D_rtd13", "H_repeat_sales", "K_qhpi", "S_sharp", "S_sharp_filtered", "S_noEB", "S_noTV", "S_LLT", "S_norobust", "S_coarseloc", "S_sharp_bottomup"]268 lab = {"A_median": "Median", "B_stratified": "Stratified median", "C_pooled_td": "Pooled TD (Huber)", "D_rtd13": "Rolling TD 13 m", "H_repeat_sales": "Repeat sales", "K_qhpi": "QHPI v2.1 arch.", "S_sharp": "SHARP-HPI (smoothed)",269 "S_sharp_filtered": "SHARP-HPI (filtered)", "S_noEB": "SHARP: no EB pooling of $\\tau^2$", "S_noTV": "SHARP: fixed quality prices", "S_LLT": "SHARP: local linear trend", "S_norobust": "SHARP: Gaussian", "S_coarseloc": "SHARP: municipality-level location FE", "S_sharp_bottomup": "SHARP: stock-weighted bottom-up"}270 scen = ["baseline", "turning_point", "flat_quality_up", "rising_quality_down", "drifting_hedonics", "heavy_tails", "sparse", "no_local_dynamics", "idiosyncratic_towns"]271 scen = [s for s in scen if s in d["scenario"].unique().to_list()]272 for lvl in ["municipality", "region", "province"]:273 rows = []274 for m in methods:275 r = [lab.get(m, m)]276 for sc in scen:277 x = d.filter((pl.col("scenario") == sc) & (pl.col("level") == lvl) & (pl.col("method") == m))278 r.append(f"{100*x['rmse_ret'][0]:.2f} ({100*x['rmse_ret_mcse'][0]:.2f})" if x.height and x["rmse_ret"][0] == x["rmse_ret"][0] else "—")279 if x.height:280 macro(f"MC{lvl}{m}{sc}Rmse", 100 * x["rmse_ret"][0], "{:.2f}"); macro(f"MC{lvl}{m}{sc}Cov", 100 * x["cov95"][0], "{:.0f}"); macro(f"MC{lvl}{m}{sc}Bias", 100 * x["bias_ret"][0], "{:.3f}")281 rows.append(r)282 table(f"tab_mc_rmse_{lvl}", ["Method"] + [s.replace("_", " ") for s in scen], rows, f"Monte Carlo: RMSE of monthly log changes (×100) at the {lvl} level, by scenario (Monte Carlo standard error in parentheses)", f"tab:mc_rmse_{lvl}", size="\\scriptsize", longtable=False,283 notes="Averages over cells and replications; see Table~\\ref{tab:mc_design} for the scenarios. Lower is better; the final two scenarios are the deliberate counter-examples.")284 rows = []285 for m in ["C_pooled_td", "K_qhpi", "S_sharp", "S_sharp_filtered", "S_noEB", "S_LLT", "S_norobust"]:286 r = [lab.get(m, m)]287 for sc in scen:288 x = d.filter((pl.col("scenario") == sc) & (pl.col("level") == "municipality") & (pl.col("method") == m))289 r.append(f"{100*x['cov95'][0]:.0f}" if x.height and x["cov95"][0] == x["cov95"][0] else "—")290 rows.append(r)291 table("tab_mc_coverage", ["Method"] + [s.replace("_", " ") for s in scen], rows, "Monte Carlo: empirical coverage (\\%) of nominal 95\\% intervals for municipal index levels", "tab:mc_coverage", size="\\scriptsize")292 rows = []293 for m in ["A_median", "C_pooled_td", "K_qhpi", "S_sharp", "S_sharp_filtered"]:294 r = [lab.get(m, m)]295 for sc in scen:296 x = d.filter((pl.col("scenario") == sc) & (pl.col("level") == "municipality") & (pl.col("method") == m))297 r.append(f"{100*x['bias_ret'][0]:+.3f}" if x.height else "—")298 rows.append(r)299 table("tab_mc_bias", ["Method"] + [s.replace("_", " ") for s in scen], rows, "Monte Carlo: mean bias of monthly log changes (×100) at the municipality level", "tab:mc_bias", size="\\scriptsize")300 rows = []301 for sc in scen:302 x = d.filter((pl.col("scenario") == sc) & (pl.col("level") == "municipality") & (pl.col("method") == "S_sharp"))303 y = d.filter((pl.col("scenario") == sc) & (pl.col("level") == "municipality") & (pl.col("method") == "K_qhpi"))304 rows.append([sc.replace("_", " "), str(int(x["reps"][0])), f"{x['lag_xcorr'][0]:+.2f}", f"{100*x['revision_rmse'][0]:.2f}", f"{y['lag_xcorr'][0]:+.2f}" if y.height else "—", f"{x['secs'][0]:.1f}" if x["secs"][0] == x["secs"][0] else "—"])305 table("tab_mc_lag", ["Scenario", "Reps", "SHARP lag (months)", "SHARP revision RMSE", "QHPI lag", "sec./rep"], rows, "Monte Carlo: timing diagnostics at the municipality level", "tab:mc_lag", size="\\small",306 notes="Lag = argmax of the cross-correlation between estimated and true 3-month changes (negative = estimate lags). Revision RMSE = RMSE between filtered (real-time) and smoothed monthly changes (×100).")307308309def t_vintages():310 if not ok(R / "vintages" / "revisions_by_release.csv", R / "vintages" / "revision_policies_province.csv"):311 return312 d = pl.read_csv(R / "vintages" / "revisions_by_release.csv")313 rows = []314 for cls in ["house", "condo", "plex"]:315 for lvl in ["province", "region", "municipality"]:316 r = [f"{cls} / {lvl}"]317 for k in [1, 2, 3, 6, 12]:318 x = d.filter((pl.col("class") == cls) & (pl.col("level") == lvl) & (pl.col("release") == k))319 r.append(f"{100*x['mae_rev_m1_vs_final'][0]:.2f}" if x.height else "—")320 x1 = d.filter((pl.col("class") == cls) & (pl.col("level") == lvl) & (pl.col("release") == 1))321 r.append(f"{100*x1['mae_rev_12m_vs_final'][0]:.2f}" if x1.height else "—"); r.append(f"{100*x1['share_m1_rev_gt_1pct'][0]:.0f}" if x1.height else "—")322 rows.append(r)323 if x1.height:324 macro(f"RevFirst{cls}{lvl}", 100 * x1["mae_rev_m1_vs_final"][0], "{:.2f}"); macro(f"RevFirstTwelve{cls}{lvl}", 100 * x1["mae_rev_12m_vs_final"][0], "{:.2f}")325 table("tab_revisions", ["Class / level", "Rel. 1", "Rel. 2", "Rel. 3", "Rel. 6", "Rel. 12", "12-m change, rel. 1", "% $|$rev$|>1$ pt"], rows,326 "Real-time vintages: mean absolute revision of the monthly log change (×100) between the $k$-th release and the final retrospective estimate", "tab:revisions", size="\\footnotesize",327 notes="Releases simulated monthly from July 2022 to October 2026 with a 10-day publication delay; municipalities with at least 1,500 sales. Final = smoothed full-sample panel.")328 p = pl.read_csv(R / "vintages" / "revision_policies_province.csv")329 table("tab_policies", ["Class", "Months", "MAE first release", "Bias first release", "MAE third release", "RMSE first", "RMSE third", "Mean $|\\Delta|$"],330 [[r["class"], str(r["months"]), f"{100*r['mae_first_release_vs_final']:.2f}", f"{100*r['bias_first_release_vs_final']:+.3f}", f"{100*r['mae_third_release_vs_final']:.2f}", f"{100*r['rmse_first']:.2f}", f"{100*r['rmse_third']:.2f}", f"{100*r['mean_abs_true_change']:.2f}"] for r in p.iter_rows(named=True)],331 "Provincial monthly change (×100): first and third release versus final", "tab:policies", notes="Full re-estimation policy (smoothed at each vintage). A revision-limited policy freezes the value at the third release; the frozen-first-release (movement-splice) policy is in Table~\\ref{tab:revisions_frozen}.")332 for r in p.iter_rows(named=True):333 macro(f"PolFirst{r['class']}", 100 * r["mae_first_release_vs_final"], "{:.2f}"); macro(f"PolThird{r['class']}", 100 * r["mae_third_release_vs_final"], "{:.2f}"); macro(f"PolMeanAbs{r['class']}", 100 * r["mean_abs_true_change"], "{:.2f}")334 if (R / "vintages" / "revisions_by_release_filtered.csv").exists():335 d = pl.read_csv(R / "vintages" / "revisions_by_release_filtered.csv")336 rows = []337 for cls in ["house", "condo", "plex"]:338 for lvl in ["province", "region", "municipality"]:339 x1 = d.filter((pl.col("class") == cls) & (pl.col("level") == lvl) & (pl.col("release") == 1)); x12 = d.filter((pl.col("class") == cls) & (pl.col("level") == lvl) & (pl.col("release") == 12))340 if x1.height:341 rows.append([f"{cls} / {lvl}", f"{100*x1['mae_rev_m1_vs_final'][0]:.2f}", f"{100*x12['mae_rev_m1_vs_final'][0]:.2f}" if x12.height else "—", f"{100*x1['mae_rev_12m_vs_final'][0]:.2f}", f"{100*x1['mae_rev_m1_vs_last_rt'][0]:.2f}"])342 macro(f"FrozenFirst{cls}{lvl}", 100 * x1["mae_rev_m1_vs_final"][0], "{:.2f}")343 table("tab_revisions_frozen", ["Class / level", "Rel. 1 vs final", "Rel. 12 vs final", "12-m change, rel. 1", "Rel. 1 vs last real-time"], rows,344 "Frozen-first-release (one-sided filtered) policy: mean absolute gap of the monthly log change (×100) between the real-time value and the final smoothed estimate", "tab:revisions_frozen", size="\\footnotesize",345 notes="Under this policy a month's value changes only through stage-1 re-estimation; the gap to the final series is the price of never revising (filter versus smoother), and it does not decline with later releases.")346347348def t_validation():349 V = R / "validation"350 if ok(V / "repeat_sales_agreement.csv"):351 d = pl.read_csv(V / "repeat_sales_agreement.csv")352 rows = []353 for lvl in ["province", "region", "municipality"]:354 for cls in ["house", "condo", "plex"]:355 x = d.filter((pl.col("level") == lvl) & (pl.col("class") == cls))356 if x.height == 0:357 continue358 rows.append([f"{lvl} / {cls}", str(x.height), n0(x["pairs"].sum()), f2(x["corr_level"].median()), f2(x["corr_12m_changes"].median()), f"{100*x['rmse_level'].median():.1f}", f"{100*(x['total_change_rs'] - x['total_change_sharp']).median():+.1f}"])359 macro(f"RSCorrLevel{lvl}{cls}", x["corr_level"].median(), "{:.2f}"); macro(f"RSGapTotal{lvl}{cls}", 100 * (x["total_change_rs"] - x["total_change_sharp"]).median(), "{:+.1f}")360 table("tab_repeat_sales", ["Level / class", "Cells", "Pairs", "Corr. levels (median)", "Corr. 12-m changes (median)", "RMSE levels (median, ×100)", "RS − SHARP total change (median, ×100)"], rows,361 "Agreement between SHARP-HPI and Case–Shiller repeat-sales indices estimated on the same cells (cells with at least 300 pairs)", "tab:repeat_sales", size="\\footnotesize",362 notes="Repeat-sales pairs: same assessment unit, holding period of at least six months, both sales in the headline sample.")363 if ok(V / "composition_shock.csv"):364 d = pl.read_csv(V / "composition_shock.csv")365 rows = [[f"Region {str(r['region']).zfill(2)}", METHOD_LABELS.get(r["method"], r["method"]), f"{100*r['original']:+.1f}", f"{100*r['shock_drop_small_from_2024']:+.1f}", f"{100*r['shock_effect']:+.1f}"] for r in d.sort(["region", "method"]).iter_rows(named=True)]366 table("tab_composition", ["Region", "Method", "2023→2024 change, original", "Same, shocked sample", "Shock effect"], rows,367 "Synthetic composition shock (houses): the smallest 40\\% of houses are removed from the sample from January 2024; change in the 2023-to-2024 mean log index (×100)", "tab:composition", size="\\footnotesize")368 for r in d.iter_rows(named=True):369 macro(f"CompShock{r['method']}{str(r['region']).zfill(2)}", 100 * r["shock_effect"], "{:+.1f}")370 if ok(V / "placebo_summary.csv"):371 d = pl.read_csv(V / "placebo_summary.csv")372 rows = [[f"{r['class']} / {r['level']}", str(r["cells"]), f"{100*r['sd_dev']:.2f}", f"{100*r['max_abs_12m_change']:.2f}", f"{100*r['share_12m_outside_95']:.1f}", f"{100*r['tau']:.2f}"] for r in d.iter_rows(named=True)]373 table("tab_placebo", ["Class / level", "Cells", "s.d. of deviation (×100)", "Median max $|$12-m change$|$ (×100)", "% of 12-m changes outside 95\\% band", "Median $\\tau$ (×100)"], rows,374 "Placebo: sale months permuted within municipality (three replications); estimated deviations should be flat", "tab:placebo", size="\\footnotesize")375 for r in d.iter_rows(named=True):376 macro(f"PlaceboFP{r['class']}{r['level']}", 100 * r["share_12m_outside_95"], "{:.1f}")377 if ok(V / "moran_residuals.csv"):378 d = pl.read_csv(V / "moran_residuals.csv").filter(pl.col("moran_I_knn10").is_not_nan())379 agg = d.group_by(["class", "scope"]).agg(pl.len().alias("months"), pl.col("moran_I_knn10").mean().alias("I"), pl.col("z").mean().alias("z"), (pl.col("z").abs() > 1.96).mean().alias("sig"))380 rows = [[f"{r['class']} / {r['scope']}", str(r["months"]), f3(r["I"]), f2(r["z"]), pct(r["sig"])] for r in agg.sort(["class", "scope"]).iter_rows(named=True)]381 table("tab_moran", ["Class / scope", "Months tested", "Mean Moran's $I$ (10-NN)", "Mean $z$", "% months significant"], rows, "Residual spatial autocorrelation of stage-1 residuals (June of each year)", "tab:moran", size="\\footnotesize")382 if ok(V / "bootstrap_vs_model_se.csv"):383 d = pl.read_csv(V / "bootstrap_vs_model_se.csv")384 rows = [[f"{r['level']} {r['geo_id']}", n0(r["n"]), f"{100*r['model_se_level_mean']:.2f}", f"{100*r['boot_se_level_mean']:.2f}", f2(r["ratio_level"]), f"{100*r['model_se_12m_mean']:.2f}", f"{100*r['boot_se_12m_mean']:.2f}", f2(r["ratio_12m"])] for r in d.iter_rows(named=True)]385 table("tab_bootstrap", ["Cell (condominiums)", "Sales", "Model s.e. level", "Bootstrap s.e. level", "Ratio", "Model s.e. 12-m", "Bootstrap s.e. 12-m", "Ratio"], rows,386 "Cluster bootstrap (40 draws, resampling fine-location cells) versus model-based standard errors, condominiums (×100)", "tab:bootstrap", size="\\footnotesize")387 macro("BootRatioLevelMedian", d["ratio_level"].median(), "{:.2f}"); macro("BootRatioTwelveMedian", d["ratio_12m"].median(), "{:.2f}")388 if ok(V / "aggregation_sensitivity.csv"):389 d = pl.read_csv(V / "aggregation_sensitivity.csv")390 lab = {"stock_units": "Housing stock (units)", "stock_value": "Housing stock (assessed value)", "tx": "Transactions 2021–2026", "eq": "Equal municipality weights"}391 rows = [[r["class"], lab[r["weights"]], str(r["cells"]), f"{100*r['total_change_bottomup']:.1f}", f"{100*r['total_change_direct']:.1f}", f"{100*r['rmse_vs_direct']:.2f}", f"{100*r['max_abs_gap']:.1f}"] for r in d.iter_rows(named=True)]392 table("tab_aggregation", ["Class", "Weights", "Municipalities", "Bottom-up total change", "Direct total change", "RMSE vs direct", "Max gap"], rows,393 "Provincial aggregation: stock-weighted, transaction-weighted and equal-weighted bottom-up aggregates of municipal indices versus the direct provincial path (log ×100, January 2021 to September 2026)", "tab:aggregation", size="\\footnotesize")394 for r in d.iter_rows(named=True):395 macro(f"Agg{r['class']}{r['weights']}Gap", 100 * (r["total_change_bottomup"] - r["total_change_direct"]), "{:+.1f}"); macro(f"AggRmse{r['class']}{r['weights']}", 100 * r["rmse_vs_direct"], "{:.1f}")396 if ok(V / "ablation_real.csv"):397 d = pl.read_csv(V / "ablation_real.csv")398 lab = {"full": "Full model", "noEB": "No EB pooling of $\\tau^2$", "noTV": "Fixed quality prices", "LLT": "Local linear trend", "norobust": "Gaussian innovations", "coarseloc": "Municipality-level location FE", "skipMRC": "No MRC level"}399 rows = []400 for v in ["full", "noEB", "noTV", "LLT", "norobust", "coarseloc", "skipMRC"]:401 r = [lab[v]]402 for cls in ["house", "condo", "plex"]:403 x = d.filter((pl.col("variant") == v) & (pl.col("class") == cls) & (pl.col("level") == "municipality"))404 r += [f"{100*x['rmse_monthly_change_vs_full'][0]:.2f}" if x.height else "—", f2(x["oos_innovation_rmse_z"][0]) if x.height else "—", f"{100*x['oos_innovation_share_within_95'][0]:.0f}" if x.height else "—"]405 if x.height:406 macro(f"Abl{v}{cls}municipality", 100 * x["rmse_monthly_change_vs_full"][0], "{:.2f}"); macro(f"AblZ{v}{cls}", x["oos_innovation_rmse_z"][0], "{:.2f}"); macro(f"AblBand{v}{cls}", 100 * x["oos_innovation_share_within_95"][0], "{:.0f}")407 rows.append(r)408 table("tab_ablation_real", ["Variant", "H: $\\Delta$ vs full", "H: OOS $z$", "H: \\% in band", "C: $\\Delta$ vs full", "C: OOS $z$", "C: \\% in band", "P: $\\Delta$ vs full", "P: OOS $z$", "P: \\% in band"], rows,409 "Real-data ablations at the municipality level: RMSE of monthly changes against the full model (×100) and one-step-ahead standardized innovation diagnostics", "tab:ablation_real", size="\\scriptsize",410 notes="OOS $z$ = root mean square of the standardized one-step-ahead prediction error of the cell's monthly mean deviation (1 if calibrated); \\% in band = share of innovations inside $\\pm 1.96$.")411 if ok(V / "trimming_sensitivity.csv"):412 d = pl.read_csv(V / "trimming_sensitivity.csv")413 lab = {"headline": "Headline sample", "with_ratio_outliers": "Assessment-ratio outliers kept", "drop_type_conflicts": "Type conflicts dropped", "drop_new_construction": "New construction dropped", "trim_1_99_monthly": "Monthly 1–99\\% price trim", "provider_only_classes": "Roll-classified sales only"}414 rows = [[r["class"], lab[r["variant"]], n0(r["n"]), f"{100*r['total_change_huber']:.1f}", f"{100*r['total_change_ols']:.1f}", f"{100*r['rmse_monthly_vs_headline']:.2f}", f"{100*r['max_abs_level_gap']:.1f}"] for r in d.iter_rows(named=True)]415 table("tab_trimming", ["Class", "Sample", "Sales", "Total change (Huber)", "Total change (OLS)", "RMSE monthly vs headline", "Max level gap"], rows, "Sensitivity of the provincial time-dummy path to the sample screens (log ×100)", "tab:trimming", size="\\footnotesize")416 if ok(V / "external_comparison.csv", V / "external_yearly.csv"):417 d = pl.read_csv(V / "external_comparison.csv"); y = pl.read_csv(V / "external_yearly.csv")418 def yr(sn, g):419 r = y.filter((pl.col("series_name") == sn) & (pl.col("geography") == g))420 return (r["2026"][0] - r["2021"][0]) if r.height else np.nan421 rows = []422 for r in d.sort(["sharp_geo", "sharp_series"]).iter_rows(named=True):423 gs = yr(r["sharp_series"], r["sharp_geo"]); ge = yr(r["external_series"], r["external_geo"])424 rows.append([r["sharp_series"].replace("SHARP ", "").replace(" stock-weighted", " (CMM/CMQ)"), r["external_source"] + ", " + r["external_series"].replace("CREA ", "").replace("Teranet ", ""), r["external_geo"], f"{100*gs:.1f}", f"{100*ge:.1f}", f"{100*(gs-ge):+.1f}", f2(r["corr_12m"]), f2(r["corr_1m"])])425 key = f"Ext{r['sharp_series']}{r['external_series']}{r['external_geo']}"426 macro(key + "Corr", r["corr_12m"], "{:.2f}"); macro(key + "Gap", 100 * (gs - ge), "{:+.1f}"); macro(key + "Sharp", 100 * gs, "{:.1f}"); macro(key + "Ext", 100 * ge, "{:.1f}")427 table("tab_external", ["SHARP series", "External series", "Geography", "SHARP 2026 vs 2021", "External 2026 vs 2021", "Gap", "Corr. 12-m", "Corr. 1-m"], rows,428 "Comparison with external indices: change of the annual mean log index from 2021 to 2026 (January--September), log ×100, and correlations of 12-month and 1-month log changes", "tab:external", size="\\scriptsize",429 notes="Teranet--National Bank: repeat sales on all registry transfers, CMA, all residential types. CREA MLS HPI: benchmark hedonic on brokered resales, CMA / provincial aggregate. SHARP metropolitan series = stock-weighted aggregates of municipal cells inside the Montréal (CMM) and Québec (CMQ) metropolitan communities; all series rebased to 2021 average = 100 before comparison.")430 rs = pl.read_csv(V / "external_repeat_sales_check.csv").filter(pl.col("min_hold_months") == 6)431 rows = [[r["sample"], n0(r["pairs"]), f"{100*r['rs_2026_vs_2021']:.1f}", f"{100*r['median_2026_vs_2021']:.1f}"] for r in rs.iter_rows(named=True)]432 table("tab_external_rs", ["Sample (this data)", "Pairs", "Repeat sales 2026 vs 2021", "Median price 2026 vs 2021"], rows,433 "Repeat-sales and median-price changes computed on the Vrai Prix data for the geographies of the external indices (annual means, log ×100)", "tab:external_rs", size="\\footnotesize",434 notes="Case--Shiller weighted repeat sales, holding period at least six months. Teranet Montréal CMA (all types) rises by " + f"{100*yr('Teranet NSA (repeat sales)', 'Montréal CMA'):.1f}" + " and Québec CMA by " + f"{100*yr('Teranet NSA (repeat sales)', 'Québec CMA'):.1f}" + " on the same basis.")435 for r in rs.iter_rows(named=True):436 macro("ExtRS" + r["sample"], 100 * r["rs_2026_vs_2021"], "{:.1f}"); macro("ExtMed" + r["sample"], 100 * r["median_2026_vs_2021"], "{:.1f}")437 macro("TeranetMtlYr", 100 * yr("Teranet NSA (repeat sales)", "Montréal CMA"), "{:.1f}"); macro("TeranetQcYr", 100 * yr("Teranet NSA (repeat sales)", "Québec CMA"), "{:.1f}")438 macro("CreaSFMtlYr", 100 * yr("CREA Single Family NSA", "Montréal CMA"), "{:.1f}"); macro("CreaSFQcYr", 100 * yr("CREA Single Family NSA", "Québec CMA"), "{:.1f}"); macro("CreaSFProvYr", 100 * yr("CREA Single Family NSA", "Quebec (province)"), "{:.1f}")439 macro("CreaAptMtlYr", 100 * yr("CREA Apartment NSA", "Montréal CMA"), "{:.1f}"); macro("CreaAptQcYr", 100 * yr("CREA Apartment NSA", "Québec CMA"), "{:.1f}")440441442def t_mc_design():443 from sharp_hpi.simulate import SCENARIOS444 rows = [[s.name.replace("_", " "), f"{s.drift*100:.1f}" + (f" → {s.break_drift*100:.1f} at m{s.break_month}" if s.break_month else ""), f"{s.tau_prov*100:.1f}/{s.tau_region*100:.1f}/{s.tau_mrc*100:.1f}/{s.tau_mun*100:.2f}",445 f"{s.quality_trend*100:+.1f}", f"{s.beta_fa_drift:+.2f}", f"{s.t_df:.1f} / {s.contamination*100:.0f}\\%", f"{s.scale_sales:.2f}"] for s in SCENARIOS.values()]446 table("tab_mc_design", ["Scenario", "Drift (\\%/m)", "$\\tau$ prov/reg/MRC/mun (\\%)", "Quality trend (\\%/m)", "$\\beta_{fa}$ drift", "$t$ d.f. / contamination", "Sales scale"], rows,447 "Monte Carlo scenarios (60 months, 4 regions × 3 MRCs × 5 municipalities, log-normal municipality sizes)", "tab:mc_design", size="\\footnotesize",448 notes="Quality trend = monthly drift of the mean log floor area of sold homes (composition); $\\beta_{fa}$ drift = change of the floor-area elasticity over the sample; sales scale multiplies expected monthly sales (baseline median $\\approx$ 3.3 per municipality).")449450451if __name__ == "__main__":452 for fn in [t_data, t_stage1, t_panel, t_benchmarks, t_simulation, t_vintages, t_validation, t_mc_design]:453 try:454 fn()455 except Exception as ex:456 import traceback; traceback.print_exc(); print("FAILED", fn.__name__, repr(ex))457 (T / "numbers.tex").write_text("% generated by scripts/61_tables.py — do not edit\n" + "\n".join(f"\\newcommand{{\\{k}}}{{{v}}}" for k, v in sorted(MACROS.items())) + "\n")458 used = set()459 for f in list((C.ROOT / "paper" / "sections").glob("*.tex")) + list((C.ROOT / "paper" / "appendix").glob("*.tex")):460 used |= set(re.findall(r"\\([A-Z][A-Za-z]+)(?=[^A-Za-z]|$)", f.read_text()))461 latex_cmds = {"E", "Var", "WPnumber", "WPtitle", "WPsubtitle", "WPdate", "WPversion", "WPabstract", "WPkeywords", "WPjel", "WPauthor", "WPaffiliation", "WPemail", "WPaddress", "DeclareMathOperator", "LARGE", "Large", "Delta", "WPabstract"}462 missing = sorted(u for u in used - set(MACROS) - latex_cmds if not u.startswith("WP"))463 (T / "numbers_missing.tex").write_text("% macros used in the text but not yet produced by results (must be empty in the final version)\n" + "\n".join(f"\\providecommand{{\\{k}}}{{\\textbf{{??{k}??}}}}" for k in missing) + "\n")464 print(len(MACROS), "macros;", len(missing), "missing:", missing)465