TeX 62.3%
Python 37.6%
1#!/usr/bin/env python32# Author: Simon-Pierre Boucher — simon-pierre.boucher@uqo.ca3"""LaTeX tables and number macros for WP23 from results/ (missing inputs skipped)."""4from __future__ import annotations56import json7import re8import sys9from pathlib import Path1011sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src"))12import numpy as np13import polars as pl1415from wp23 import config as C1617T = C.ROOT / "paper" / "tables"; T.mkdir(parents=True, exist_ok=True); R = C.RESULTS18MACROS: dict[str, str] = {}19WORDS = {str(i): w for i, w in enumerate(["Zero", "One", "Two", "Three", "Four", "Five", "Six", "Seven", "Eight", "Nine"])}20LABELS = {"naive_mun_median": "Naive: municipal median", "hedonic_ols": "Hedonic OLS", "ridge": "Ridge", "elastic_net": "Elastic net", "random_forest": "Random forest", "xgboost": "XGBoost", "lightgbm": "LightGBM",21 "svr": "SVR (subsample)", "mlp": "Quantile MLP", "hshn_partial": "HSHN (partial pooling)", "hshn_full": "HSHN (full pooling)", "hshn_none": "HSHN (no pooling)",22 "no_change": "No change", "mean_growth": "Mean growth", "momentum": "Momentum", "rw_drift": "Random walk with drift", "seasonal_naive": "Seasonal naive", "ar2_cell": "AR(2), cell", "local_linear_trend": "Local linear trend",23 "ridge_panel": "Ridge (panel)", "lightgbm_panel": "LightGBM (panel)", "mlp_panel": "MLP (panel)", "lstm": "LSTM", "gru": "GRU", "tcn": "TCN", "tft_lite": "TFT-lite", "dsgfn": "DSGFN", "raref": "RAREF", "bhpf": "BHPF", "qanss": "QANSS", "qrefa": "QREFA"}24ORDER_A = ["naive_mun_median", "hedonic_ols", "ridge", "elastic_net", "svr", "random_forest", "xgboost", "lightgbm", "mlp", "hshn_none", "hshn_full", "hshn_partial"]25ORDER_B = ["no_change", "mean_growth", "momentum", "rw_drift", "seasonal_naive", "ar2_cell", "local_linear_trend", "ridge_panel", "lightgbm_panel", "mlp_panel", "lstm", "gru", "tcn", "tft_lite", "bhpf", "raref", "qanss", "dsgfn", "qrefa"]262728def mname(n): return "".join(WORDS[c] if c.isdigit() else c for c in re.sub(r"[^A-Za-z0-9]", "", n))29def macro(n, v, fmt="{:,.1f}"):30 MACROS[mname(n)] = "n/a" if v is None or (isinstance(v, float) and np.isnan(v)) else fmt.format(v).replace(",", "{,}")31def esc(s): return str(s).replace("&", "\\&").replace("%", "\\%").replace("_", "\\_").replace("≤", "$\\le$").replace("≥", "$\\ge$").replace("×", "$\\times$")32def ok(*p):33 miss = [str(x) for x in p if not Path(x).exists()]34 if miss: print("skip:", miss)35 return not miss363738def table(name, header, rows, caption, label, align=None, notes="", size="\\small", fit=None):39 align = align or "l" + "r" * (len(header) - 1); body = "\n".join(" & ".join(esc(c) for c in r) + " \\\\" for r in rows)40 fit = (len(header) > 7) if fit is None else fit41 tab = f"""\\begin{{tabular}}{{{align}}}42\\toprule43{' & '.join(esc(h) for h in header)} \\\\44\\midrule45{body}46\\bottomrule47\\end{{tabular}}"""48 if fit: # wide table: scale the tabular to the text width; notes in a full-width parbox (threeparttable mis-measures scaled boxes)49 inner = "\\resizebox{\\textwidth}{!}{%\n" + tab + "}" + (("\n\\par\\vspace{3pt}\\parbox{\\textwidth}{\\footnotesize " + notes + "}") if notes else "")50 else:51 inner = "\\begin{threeparttable}\n" + tab + "\n" + (("\\begin{tablenotes}\\footnotesize\\item " + notes + "\\end{tablenotes}\n") if notes else "") + "\\end{threeparttable}"52 (T / f"{name}.tex").write_text(f"""\\begin{{table}}[htbp]53\\centering{size}54\\caption{{{caption}}}\\label{{{label}}}55{inner}56\\end{{table}}57"""); print("table", name)585960def t_data():61 macro("NBenchB", 14, "{:.0f}"); macro("NBenchA", 8, "{:.0f}")62 s = pl.read_parquet(C.WP22_ANALYTIC, columns=["tx_id", "property_class", "region_code", "mun_code", "month", "price", "floor_area_m2", "land_area_m2", "age", "storeys"])63 macro("NSales", s.height, "{:,.0f}"); macro("NMun", s["mun_code"].n_unique(), "{:,.0f}")64 rows = [["Floor area of storeys (m²)", "roll", "continuous; cubic in centred log", "all tasks"], ["Lot area (m²)", "roll", "quadratic in centred log", "all tasks"], ["Year built / age", "roll", "nine age bins; missingness indicator", "all tasks"],65 ["Storeys, building genre, physical link, dwellings (plex)", "roll", "categorical", "Task A"], ["Coordinates, municipality, MRC, region", "feed / SDA V2026-07", "location; nested metric cells (WP22)", "all tasks"],66 ["Assessment in force at sale", "feed", "not used as a predictor (post-dated for Task A test sales would be fine, but excluded for comparability)", "—"],67 ["Municipal / regional market features at origin", "computed from sales ≤ origin", "median log price and quality-stripped level (3/6/12/24 m), momentum, volumes, dispersion", "Task A"],68 ["10 nearest past sales (≤ 12 m)", "computed", "mean quality-stripped log price, mean distance, dispersion", "Task A"],69 ["Real-time SHARP-HPI vintage (13-month history)", "WP22 index_vintages (smoothed mode)", "1/3/6/12-m changes, volatilities, acceleration, counts, s.e.; parent (region, province) momentum", "Task B"],70 ["Bedrooms, bathrooms, parking, basement, condition, renovations, quality, listing, images", "—", "not available", "—"]]71 table("tab_variables", ["Variable", "Source", "Treatment", "Used in"], rows, "Variables available and their use", "tab:variables", align="p{4.2cm}p{3cm}p{5.2cm}p{1.6cm}", size="\\footnotesize")72 cov = s.group_by(["region_code", "property_class"]).agg(pl.len().alias("n")).pivot(on="property_class", index="region_code", values="n").fill_null(0).sort("region_code")73 rows = [[r["region_code"], f"{r.get('house',0):,}".replace(",", "{,}"), f"{r.get('condo',0):,}".replace(",", "{,}"), f"{r.get('plex',0):,}".replace(",", "{,}")] for r in cov.iter_rows(named=True)]74 if ok(R / "task_a" / "timing_trees.csv"):75 ta = pl.concat([pl.read_csv(f) for f in (R / "task_a").glob("timing_*.csv")]); macro("NOriginsA", ta["origin"].n_unique(), "{:.0f}"); macro("NTestA", ta.group_by("origin").agg(pl.col("n_test").first())["n_test"].sum(), "{:,.0f}")76 if ok(C.PROCESSED / "task_b_panel.parquet"):77 pb = pl.read_parquet(C.PROCESSED / "task_b_panel.parquet"); macro("NCellsB", pb.select(["class", "level", "geo_id"]).unique().height, "{:.0f}"); macro("NRowsB", pb.height, "{:,.0f}"); macro("NOriginsB", pb["origin"].n_unique(), "{:.0f}")78 for h in C.HORIZONS: macro(f"NOriginsBh{h}", pb.filter(pl.col("h") == h)["origin"].n_unique(), "{:.0f}")79 table("tab_coverage", ["Region", "Houses", "Condominiums", "Plexes"], rows, "Headline sales by administrative region and class (WP22 sample reused)", "tab:coverage", size="\\footnotesize")80 # selection: sold vs stock81 if ok(C.WP22_STOCK):82 st = pl.read_parquet(C.WP22_STOCK, columns=["cubf", "nb_logements", "lien_physique", "annee_construction", "aire_etages_m2", "superficie_terrain_m2"]).filter(pl.col("cubf") == 1000)83 st = st.with_columns(pl.when(pl.col("nb_logements") >= 6).then(None).when(pl.col("nb_logements").is_between(2, 5)).then(pl.lit("plex")).when((pl.col("nb_logements") <= 1) & (pl.col("lien_physique") == "Intégré")).then(pl.lit("condo")).when(pl.col("nb_logements") <= 1).then(pl.lit("house")).otherwise(None).alias("cls")).filter(pl.col("cls").is_not_null())84 rows = []85 for cls in ["house", "condo", "plex"]:86 a = s.filter(pl.col("property_class") == cls); b = st.filter(pl.col("cls") == cls)87 rows.append([cls, f"{a['floor_area_m2'].median():.0f}", f"{b['aire_etages_m2'].median():.0f}", f"{a['age'].median():.0f}", f"{2026 - b['annee_construction'].median():.0f}", f"{a['land_area_m2'].median():.0f}", f"{b['superficie_terrain_m2'].median():.0f}", f"{a.height:,}".replace(",", "{,}"), f"{b.height:,}".replace(",", "{,}")])88 macro(f"SelFa{cls}Sold", a["floor_area_m2"].median(), "{:.0f}"); macro(f"SelFa{cls}Stock", b["aire_etages_m2"].median(), "{:.0f}"); macro(f"SelAge{cls}Sold", a["age"].median(), "{:.0f}"); macro(f"SelAge{cls}Stock", 2026 - b["annee_construction"].median(), "{:.0f}")89 table("tab_selection", ["Class", "Floor area sold", "Floor area stock", "Age sold", "Age stock", "Lot sold", "Lot stock", "Sales", "Stock units"], rows, "Transaction selection: medians of sold dwellings (2021--2026) versus the assessment-roll stock", "tab:selection", size="\\footnotesize",90 notes="Stock = residential use code 1000 in the 2026 roll, classified with the WP22 rules; age of the stock as of 2026.")91 rows = [["HSHN", "embeddings 8 (geo) / 4 (h); encoder 256-128 GELU, dropout 0.1; AdamW, one-cycle lr $2\\times10^{-3}$, 8 epochs, batch 4,096; municipality weight decay $5\\times10^{-3}$"],92 ["MLP (A)", "256-128 GELU, dropout 0.1; pinball loss on 5 quantiles; same optimiser"], ["LightGBM (A)", "800 trees, lr 0.04, 63 leaves, min child 30, subsample 0.8; quantile models 400 trees"], ["XGBoost (A)", "600 trees, lr 0.05, depth 8, subsample 0.8"], ["Random forest (A)", "150 trees, depth 18, min leaf 5, max features 0.5"], ["SVR (A)", "RBF, C=3, $\\epsilon$=0.05, 15,000-sale subsample"],93 ["Panel MLP / LSTM / GRU / TCN / TFT-lite (B)", "hidden 64 (MLP) / 32 (RNN, attention) / 16 channels (TCN); 60 epochs, batch 512, AdamW lr $2\\times10^{-3}$, weight decay $10^{-4}$; pinball loss"],94 ["DSGFN (B)", "node encoder 32; one attention layer; prior adjacency: Gaussian distance kernel (60 km) $\\times$ co-movement; 60 epochs"], ["RAREF (B)", "3-state Gaussian HMM (60 EM iterations) on provincial real-time growth; 3 ridge experts ($\\alpha$=3)"],95 ["BHPF (B)", "conjugate Gaussian hierarchy; 15 variance-component iterations; priors $\\tau_b^2$=1 (slopes)"], ["QANSS (B)", "32 tanh units; 200 Adam steps, lr $5\\times10^{-3}$; Gaussian likelihood"], ["QREFA (B)", "NNLS stacking of QANSS, BHPF, DSGFN (40 epochs), RAREF on training origins"],96 ["LightGBM (B)", "300 trees, lr 0.03, 15 leaves, min child 40; quantile models"], ["Ridge (B)", "$\\alpha$=3 on standardised features"]]97 table("tab_hyper", ["Model", "Hyper-parameters (fixed in advance, never tuned on test origins)"], rows, "Hyper-parameters", "tab:hyper", align="p{3.4cm}p{11.5cm}", size="\\footnotesize")98 rows = [["Hedonic coefficients", "WP22 stage 1 re-estimated at every vintage; Task A features never use coefficients"], ["Geographic clustering / location cells", "WP22 nested grid rebuilt from data ≤ origin at every vintage; Task A uses raw coordinates"],99 ["Normalisation", "standardisation constants from training rows only; index levels enter as changes"], ["Feature selection", "fixed feature lists, no selection on test data"], ["Imputation", "WP22 conditional medians (no price information); indicators time-stable"],100 ["Target encoding", "none; one-hot with training categories"], ["Spatial graph", "DSGFN adjacency recomputed at each origin from centroids and the real-time history ≤ origin"], ["Hyper-parameter tuning", "fixed a priori (Table~\\ref{tab:hyper})"],101 ["Index revisions", "predictors = vintage published for the origin; target = final series (revision component reported)"], ["Partial last month", "October 2026 excluded; last target September 2026"]]102 table("tab_leakage", ["Item", "Treatment"], rows, "Leakage audit", "tab:leakage", align="p{4cm}p{11cm}", size="\\footnotesize")103104105def scores_table(tag, order, caption, label, hs, ycol_desc):106 f = R / f"eval_{tag}" / "scores.csv"107 if not ok(f): return108 d = pl.read_csv(f, infer_schema_length=None).filter(pl.col("group") == "all"); mcs = pl.read_csv(R / f"eval_{tag}" / "mcs.csv", infer_schema_length=None) if (R / f"eval_{tag}" / "mcs.csv").exists() else None109 tests = pl.read_csv(R / f"eval_{tag}" / "tests_vs_reference.csv", infer_schema_length=None) if (R / f"eval_{tag}" / "tests_vs_reference.csv").exists() else None110 rows = []111 for m in order:112 if m not in d["model"].to_list(): continue113 r = [LABELS.get(m, m)]114 for h in hs:115 x = d.filter((pl.col("model") == m) & (pl.col("h").cast(pl.Utf8) == str(h)))116 if x.height == 0: r += ["—", "—"]; continue117 mae = 100 * x["mae"][0]; star = ""118 if tests is not None:119 tt = tests.filter((pl.col("model") == m) & (pl.col("h").cast(pl.Utf8) == str(h)))120 if tt.height and tt["dm_p"][0] == tt["dm_p"][0]:121 star = "***" if tt["dm_p"][0] < 0.01 else "**" if tt["dm_p"][0] < 0.05 else "*" if tt["dm_p"][0] < 0.10 else ""122 inm = ""123 if mcs is not None:124 mm = mcs.filter((pl.col("model") == m) & (pl.col("h").cast(pl.Utf8) == str(h)))125 if mm.height and mm["in_mcs_90"][0]: inm = "$^{\\dagger}$"126 r += [f"{mae:.2f}{star}{inm}", f"{x['mase'][0]:.2f}" if x["mase"][0] == x["mase"][0] else "—"]127 macro(f"{tag}MAE{m}h{h}", mae, "{:.2f}"); macro(f"{tag}MASE{m}h{h}", x["mase"][0], "{:.2f}"); macro(f"{tag}CRPS{m}h{h}", 100 * x["crps"][0] if x["crps"][0] == x["crps"][0] else None, "{:.2f}"); macro(f"{tag}Cov{m}h{h}", 100 * x["cov90"][0] if x["cov90"][0] == x["cov90"][0] else None, "{:.0f}")128 rows.append(r)129 header = ["Model"] + [f"{h}: MAE" if i % 2 == 0 else f"{h}: MASE" for h in hs for i in range(2)]130 table(f"tab_scores_{tag}", header, rows, caption, label, size="\\scriptsize",131 notes=f"{ycol_desc} MAE in log points ×100; MASE relative to the naive reference. Stars: Diebold--Mariano test of equal squared-error loss against the reference on per-origin mean losses (* 10\\%, ** 5\\%, *** 1\\%); $\\dagger$ = in the 90\\% model confidence set of the horizon.")132 rows = []133 for m in order:134 if m not in d["model"].to_list(): continue135 r = [LABELS.get(m, m)]136 for h in hs:137 x = d.filter((pl.col("model") == m) & (pl.col("h").cast(pl.Utf8) == str(h)))138 r += [f"{100*x['crps'][0]:.2f}" if x.height and x["crps"][0] == x["crps"][0] else "—", f"{100*x['cov90'][0]:.0f}" if x.height and x["cov90"][0] == x["cov90"][0] else "—", f"{100*x['width90'][0]:.1f}" if x.height and x["width90"][0] == x["width90"][0] else "—"]139 rows.append(r)140 header = ["Model"] + [f"{h}: {k}" for h in hs for k in ("CRPS", "cov90", "width90")]141 table(f"tab_prob_{tag}", header, rows, caption.replace("Point", "Probabilistic").replace("accuracy", "accuracy (CRPS ×100, 90\\% interval coverage \\% and width ×100)"), label.replace("scores", "prob"), size="\\scriptsize")142143144def t_results():145 f = R / "eval_a" / "scores.csv"146 if f.exists():147 d = pl.read_csv(f, infer_schema_length=None).filter(pl.col("group") == "all")148 for m in ["lightgbm", "hshn_partial", "hedonic_ols"]:149 a = d.filter((pl.col("model") == m) & (pl.col("h") == "1-3")); b = d.filter((pl.col("model") == m) & (pl.col("h") == "7-12"))150 if a.height and b.height: macro(f"aMAEDelta{m}", 100 * (b["mae"][0] - a["mae"][0]), "{:.1f}")151 scores_table("a", ORDER_A, "Task A — Point accuracy of future sale-price forecasts by horizon band (months after the origin)", "tab:scores_a", ["1-3", "4-6", "7-12"], "Target: log sale price.")152 scores_table("b", ORDER_B, "Task B — Point accuracy of quality-adjusted index-change forecasts by horizon (months)", "tab:scores_b", [1, 3, 6, 12], "Target: final SHARP-HPI log change.")153 for tag, grp, lab in [("b", "level", "level"), ("b", "class", "class"), ("b", "size_level", "cell size"), ("a", "mun_size", "municipality size"), ("a", "property_class", "class"), ("a", "segment", "price segment")]:154 f = R / f"eval_{tag}" / "scores.csv"155 if not ok(f): continue156 d = pl.read_csv(f, infer_schema_length=None).filter(pl.col("group") == grp)157 if d.height == 0: continue158 groups = sorted(d[grp].unique().to_list()); order = ORDER_B if tag == "b" else ORDER_A; hs = [3] if tag == "b" else ["4-6"]159 rows = []160 for m in order:161 if m not in d["model"].to_list(): continue162 r = [LABELS.get(m, m)]163 for g in groups:164 x = d.filter((pl.col("model") == m) & (pl.col(grp) == g) & (pl.col("h").cast(pl.Utf8) == str(hs[0])))165 r.append(f"{100*x['mae'][0]:.2f}" if x.height else "—")166 rows.append(r)167 table(f"tab_het_{tag}_{grp}", ["Model"] + [str(g) for g in groups], rows, f"Task {tag.upper()} — MAE (×100) by {lab}, horizon {hs[0]}", f"tab:het_{tag}_{grp}", size="\\scriptsize")168169170def t_ablation():171 f = R / "ablation_b" / "ablation_summary.csv"172 if not ok(f): return173 d = pl.read_csv(f, infer_schema_length=None).filter(pl.col("level") == "municipality")174 lab = {"dsgfn": "DSGFN (full)", "dsgfn_no_distance": "DSGFN without distance kernel", "dsgfn_no_comovement": "DSGFN without co-movement", "dsgfn_no_graph": "DSGFN with uninformative prior graph", "qanss": "QANSS (full)", "qanss_no_parent": "QANSS without parent drifts",175 "qrefa": "QREFA (NNLS stack)", "equal_weight_combo": "Equal-weight combination of the four", "qrefa_without_qanss": "QREFA without QANSS", "qrefa_without_bhpf": "QREFA without BHPF", "qrefa_without_dsgfn": "QREFA without DSGFN", "qrefa_without_raref": "QREFA without RAREF",176 "ridge_panel": "Ridge panel (full features)", "ridge_no_momentum": "Ridge without own momentum", "ridge_no_parent": "Ridge without parent (region/province) features", "ridge_no_volume": "Ridge without volume features"}177 rows = []178 for m in lab:179 r = [lab[m]]180 for h in (1, 3, 6):181 x = d.filter((pl.col("model") == m) & (pl.col("h") == h))182 if x.height == 0: r += ["—", "—"]; continue183 diff = x["mae_diff_vs_ref"][0]; lo, hi = x["diff_lo"][0], x["diff_hi"][0]184 sig = "" if diff is None or diff != diff else ("*" if (lo > 0 or hi < 0) else "")185 r += [f"{100*x['mae'][0]:.2f}", (f"{100*diff:+.2f}{sig}" if diff == diff and diff is not None else "ref")]186 macro(f"Abl{m}h{h}", 100 * x["mae"][0], "{:.2f}")187 if diff == diff and diff is not None: macro(f"AblDiff{m}h{h}", 100 * diff, "{:+.2f}")188 rows.append(r)189 table("tab_ablation_b", ["Variant", "h=1: MAE", "vs ref", "h=3: MAE", "vs ref", "h=6: MAE", "vs ref"], rows, "Task B ablations at the municipality level: MAE (×100) and difference to the reference variant (per-origin mean, * = 95\\% stationary-bootstrap interval excludes zero)", "tab:ablation_b", size="\\scriptsize",190 notes="Reference variants: DSGFN, QANSS, QREFA and the full ridge panel. Positive differences mean the ablated variant is worse.", fit=True)191 a = R / "ablation_b" / "dsgfn_attention_summary.csv"192 if a.exists():193 x = pl.read_csv(a, infer_schema_length=None)194 for r in x.iter_rows(named=True):195 macro(f"AttTop{r['h']}", r["top_weight"], "{:.2f}"); macro(f"AttDist{r['h']}", r["median_dist_top_km"], "{:.0f}"); macro(f"AttSameReg{r['h']}", 100 * r["same_region"], "{:.0f}"); macro(f"AttSameMrc{r['h']}", 100 * r["same_mrc"], "{:.0f}")196197198def t_robust():199 f = R / "robustness_b" / "summary.csv"200 if not ok(f): return201 d = pl.read_csv(f, infer_schema_length=None).filter(pl.col("level") == "municipality")202 variants = ["expanding", "rolling_24", "min_n_10", "winsorised", "seed_1", "seed_2", "bandwidth"]; labv = {"expanding": "Expanding (baseline)", "rolling_24": "Rolling 24 origins", "min_n_10": "Cells with $n \\ge 10$", "winsorised": "Winsorised target", "seed_1": "Seed 1", "seed_2": "Seed 2", "bandwidth": "Bandwidth 30 / 120 km"}203 models = ["no_change", "momentum", "ar2_cell", "ridge_panel", "qanss", "tft_lite", "dsgfn", "dsgfn_bw30", "dsgfn_bw120"]204 rows = []205 for m in models:206 r = [LABELS.get(m, m.replace("dsgfn_bw", "DSGFN, bandwidth ").replace("30", "30 km").replace("120", "120 km") if m.startswith("dsgfn_bw") else LABELS.get(m, m))]207 for v in variants:208 x = d.filter((pl.col("model") == m) & (pl.col("variant") == v)); r.append(f"{100*x['mae'][0]:.2f}" if x.height else "—")209 if x.height: macro(f"Rob{m}{v}", 100 * x["mae"][0], "{:.2f}"); macro(f"RobCov{m}{v}", 100 * x["cov90"][0], "{:.0f}")210 rows.append(r)211 table("tab_robust_b", ["Model"] + [labv[v] for v in variants], rows, "Task B robustness at the municipality level, horizon 3: MAE (×100) under alternative designs", "tab:robust_b", size="\\scriptsize",212 notes="Rolling: models estimated on the 24 most recent origins only. Seeds: neural models re-estimated with different random seeds (the baseline is seed 0). Winsorised: targets clipped at the 1st and 99th percentiles in training and test.", fit=True)213214215def t_interpret():216 f = R / "interpret" / "shap_importance_wide.csv"217 if not ok(f): return218 d = pl.read_csv(f, infer_schema_length=None).head(15); cols = [c for c in d.columns if c not in ("feature", "mean")]219 lab = {"knn_q": "10 nearest past sales: quality-stripped log price", "lfa": "log floor area (centred)", "age_f": "age", "origin_t": "origin month (time index)", "lla": "log lot area (centred)", "mun_q_12": "municipality: quality-stripped level, 12 m",220 "y_m": "northing", "x_m": "easting", "mun_medlp_3": "municipality: median log price, 3 m", "knn_dist_km": "distance to the 10 nearest sales", "knn_q_sd": "dispersion of the 10 nearest sales", "mun_n_3": "municipality: sales, 3 m", "mun_n_12": "municipality: sales, 12 m",221 "reg_q_12": "region: quality-stripped level, 12 m", "reg_q_3": "region: level, 3 m", "lfa2": "log floor area squared", "mun_mom_3_12": "municipality: momentum 3 vs 12 m", "reg_mom_3_12": "region: momentum", "mun_vol_ratio": "municipality: volume ratio"}222 rows = [[lab.get(r["feature"], r["feature"])] + [f"{100*r[c]:.1f}" for c in cols] + [f"{100*r['mean']:.1f}"] for r in d.iter_rows(named=True)]223 table("tab_shap", ["Feature"] + cols + ["Mean"], rows, "Task A — mean absolute SHAP value (×100, log-price units) of the LightGBM model at three origins, 4,000 test sales each", "tab:shap", size="\\footnotesize",224 notes="Importance is predictive, not causal. Origin-month importance declines as the training history lengthens.", fit=True)225 for r in d.iter_rows(named=True):226 macro(f"Shap{r['feature']}", 100 * r["mean"], "{:.1f}")227 for c in cols: macro(f"Shap{r['feature']}{c}", 100 * r[c], "{:.1f}")228 a = pl.read_csv(R / "interpret" / "ale.csv", infer_schema_length=None).filter(pl.col("origin") == "2024-03")229 for feat in ["mun_mom_3_12", "knn_q", "mun_vol_ratio"]:230 x = a.filter(pl.col("feature") == feat).sort("bin"); macro(f"ALE{feat}Range", 100 * (x["ale"].max() - x["ale"].min()), "{:.1f}"); macro(f"ALE{feat}Top", 100 * x["ale"][-1], "{:+.1f}"); macro(f"ALE{feat}Bottom", 100 * x["ale"][0], "{:+.1f}")231232233def t_extra():234 X = R / "extra"235 if not ok(X / "tests_holm_b.csv"): return236 t = pl.read_csv(X / "tests_holm_b.csv", infer_schema_length=None)237 rows, rows2 = [], []238 for m in ORDER_B:239 if m == "no_change": continue240 r, r2 = [LABELS.get(m, m)], [LABELS.get(m, m)]241 for h in (1, 3, 6, 12):242 x = t.filter((pl.col("model") == m) & (pl.col("h") == h))243 has = x.height and x["dm_stat"][0] is not None244 r += [f"{x['dm_stat'][0]:+.2f}" if has else "—", f"{x['dm_p_holm'][0]:.2f}" if has else "—"]245 r2 += [f"{x['dm_p'][0]:.3f}" if has else "—", f"{x['gw_p'][0]:.3f}" if has and x["gw_p"][0] is not None else "—"]246 if has: macro(f"HolmP{m}h{h}", x["dm_p_holm"][0], "{:.2f}"); macro(f"DMP{m}h{h}", x["dm_p"][0], "{:.3f}"); macro(f"GWP{m}h{h}", x["gw_p"][0], "{:.3f}")247 rows.append(r); rows2.append(r2)248 table("tab_tests_b", ["Model"] + [f"h={h}: {k}" for h in (1, 3, 6, 12) for k in ("DM", "Holm $p$")], rows, "Task B — Diebold--Mariano statistics against the no-change forecast (negative = better than no change) and Holm-adjusted $p$-values across the 18 models of each horizon", "tab:tests_b", size="\\scriptsize", fit=True,249 notes="HAC variance with $h-1$ Bartlett lags and the Harvey--Leybourne--Newbold correction; $t_{n-1}$ reference distribution; Holm step-down over the 18 two-sided $p$-values of each horizon.")250 table("tab_gw_b", ["Model"] + [f"h={h}: {k}" for h in (1, 3, 6, 12) for k in ("DM $p$", "GW $p$")], rows2, "Task B — unadjusted Diebold--Mariano $p$-values and Giacomini--White conditional-predictive-ability $p$-values against the no-change forecast", "tab:gw_b", size="\\scriptsize", fit=True,251 notes="The Giacomini--White test regresses the loss differential on a constant and its own lag; a rejection means the relative accuracy of the two forecasts is predictable from the past.")252 c = pl.read_csv(X / "computational_cost.csv", infer_schema_length=None)253 rows = [[r["task"], LABELS.get(r["model"], r["model"]), f"{r['mean_seconds_per_fit']:.1f}", f"{r['total_seconds']:.0f}", str(r["fits"]), f"{r['mean_train_rows']:,.0f}".replace(",", "{,}")] for r in c.iter_rows(named=True)]254 table("tab_cost", ["Task", "Model", "Seconds per fit", "Total seconds", "Fits", "Mean training rows"], rows, "Computational cost (wall-clock on a 32-core Mac Studio M3 Ultra; tree models 10 threads, neural models 8 threads)", "tab:cost", size="\\scriptsize")255 for r in c.iter_rows(named=True): macro(f"Cost{r['task']}{r['model']}", r["mean_seconds_per_fit"], "{:.1f}")256 y = pl.read_csv(X / "error_by_year_b.csv", infer_schema_length=None).filter(pl.col("h") == 3)257 yrs = sorted(y["year"].unique().to_list()); rows = []258 for m in KEYS_B:259 x = y.filter(pl.col("model") == m)260 if x.height == 0: continue261 rows.append([LABELS.get(m, m)] + [f"{100*x.filter(pl.col('year') == yr)['mae'][0]:.2f} ({x.filter(pl.col('year') == yr)['mase'][0]:.2f})" if x.filter(pl.col("year") == yr).height else "—" for yr in yrs])262 for yr in yrs:263 z = x.filter(pl.col("year") == yr)264 if z.height: macro(f"Year{m}{yr}", z["mase"][0], "{:.2f}")265 table("tab_year_b", ["Model"] + [f"Origins in {yr}" for yr in yrs], rows, "Task B, horizon 3 — MAE (×100) and MASE (in parentheses) by year of the forecast origin", "tab:year_b", size="\\scriptsize", notes="2023 covers origins from October to December only; 2026 covers January to June.")266 g = pl.read_csv(X / "error_by_regime_b.csv", infer_schema_length=None).filter(pl.col("h") == 3); regs = ["rising", "flat", "falling"]; rows = []267 for m in KEYS_B:268 x = g.filter(pl.col("model") == m)269 if x.height == 0: continue270 rows.append([LABELS.get(m, m)] + [f"{100*x.filter(pl.col('regime') == rg)['mae'][0]:.2f} ({x.filter(pl.col('regime') == rg)['mase'][0]:.2f})" if x.filter(pl.col("regime") == rg).height else "—" for rg in regs])271 for rg in regs:272 z = x.filter(pl.col("regime") == rg)273 if z.height: macro(f"Regime{m}{rg}", z["mase"][0], "{:.2f}")274 ns = {rg: (int(g.filter((pl.col("regime") == rg) & (pl.col("model") == "no_change"))["origins"][0]) if g.filter((pl.col("regime") == rg) & (pl.col("model") == "no_change")).height else 0) for rg in regs}275 table("tab_regime_b", ["Model"] + [f"{rg} ({ns[rg]} origins)" for rg in regs], rows, "Task B, horizon 3 — MAE (×100) and MASE (in parentheses) by ex-post provincial regime of the target window (rising: provincial house index change above +2 log points over the horizon; falling: below $-2$; flat otherwise)", "tab:regime_b", size="\\scriptsize")276 for rg in regs: macro(f"RegimeOrigins{rg}", ns[rg], "{:.0f}")277 d = pl.read_csv(X / "directional_b.csv", infer_schema_length=None); rows = []278 for m in KEYS_B:279 if m == "no_change": continue280 r = [LABELS.get(m, m)]281 for h in (1, 3, 6, 12):282 x = d.filter((pl.col("model") == m) & (pl.col("h") == h))283 r += [f"{100*x['directional_accuracy'][0]:.0f}" if x.height else "—", f"{100*x['directional_accuracy_at_turning_points'][0]:.0f}" if x.height and x["directional_accuracy_at_turning_points"][0] is not None else "—"]284 if x.height: macro(f"Dir{m}h{h}", 100 * x["directional_accuracy"][0], "{:.0f}"); macro(f"DirTP{m}h{h}", 100 * x["directional_accuracy_at_turning_points"][0] if x["directional_accuracy_at_turning_points"][0] is not None else None, "{:.0f}")285 rows.append(r)286 table("tab_direction_b", ["Model"] + [f"h={h}: {k}" for h in (1, 3, 6, 12) for k in ("all", "turning pts")], rows, "Task B — directional accuracy (\\% of forecasts with the sign of the realised change, changes above 0.5 log points) overall and at turning points (origins where the realised change reverses sign)", "tab:direction_b", size="\\scriptsize",287 notes="The no-change forecast has no sign; momentum has the sign of the last three-month change.", fit=True)288 dc = pl.read_csv(X / "decline_b.csv", infer_schema_length=None); rows = []289 for m in KEYS_B:290 r = [LABELS.get(m, m)]291 for h in (1, 3, 6, 12):292 x = dc.filter((pl.col("model") == m) & (pl.col("h") == h))293 r += [f"{x['brier_decline'][0]:.3f}" if x.height else "—", f"{x['brier_skill_decline'][0]:+.2f}" if x.height else "—"]294 if x.height: macro(f"BrierSkill{m}h{h}", x["brier_skill_decline"][0], "{:+.2f}"); macro(f"ShareDeclinesh{h}", 100 * x["share_declines"][0], "{:.0f}")295 rows.append(r)296 table("tab_decline_b", ["Model"] + [f"h={h}: {k}" for h in (1, 3, 6, 12) for k in ("Brier", "skill")], rows, "Task C — probability of an index decline over the horizon: Brier score and Brier skill score relative to the climatological frequency of declines", "tab:decline_b", size="\\scriptsize",297 notes="Probabilities are read off the five-quantile forecasts by linear interpolation of the quantile function. Skill = $1 - $ Brier / Brier of the constant climatological probability.")298 cb = pl.read_csv(X / "combinations_b.csv", infer_schema_length=None); lab = {"eq_lgbm_qanss": "Equal weights: LightGBM + QANSS", "eq_lgbm_qanss_ridge_ar2": "Equal weights: LightGBM + QANSS + ridge + AR(2)", "eq_lgbm_ar2": "Equal weights: LightGBM + AR(2)", "eq_all_with_nochange": "Equal weights: four models + no change"}299 rows = []300 for name in ["lightgbm_panel", "qanss", "ridge_panel", "ar2_cell", "qrefa", "eq_lgbm_qanss", "eq_lgbm_ar2", "eq_lgbm_qanss_ridge_ar2", "eq_all_with_nochange"]:301 x = cb.filter(pl.col("combination") == name); r = [lab.get(name, LABELS.get(name, name))]302 for h in (1, 3, 6, 12):303 z = x.filter(pl.col("h") == h); r.append(f"{z['mase'][0]:.3f}" if z.height else "—")304 if z.height: macro(f"Combo{name}h{h}", z["mase"][0], "{:.3f}")305 rows.append(r)306 table("tab_combo_b", ["Forecast", "h=1", "h=3", "h=6", "h=12"], rows, "Task B — MASE of equal-weight combinations versus their components (same cells and origins)", "tab:combo_b", size="\\footnotesize")307 md = pl.read_csv(X / "municipal_mase_summary_b.csv", infer_schema_length=None)308 rows = [[LABELS.get(r["model"], r["model"]), f"{r['median_mase']:.2f}", f"{r['p10']:.2f}", f"{r['p90']:.2f}", f"{100*r['share_beating_no_change']:.0f}", str(r["cells"])] for r in md.iter_rows(named=True)]309 table("tab_mun_dist_b", ["Model", "Median MASE", "P10", "P90", "% of cells beating no change", "Cells"], rows, "Task B, horizon 3 — distribution of the MASE across the 192 municipal cells (64 municipalities × 3 classes)", "tab:mun_dist_b", size="\\footnotesize")310 for r in md.iter_rows(named=True): macro(f"MunShare{r['model']}", 100 * r["share_beating_no_change"], "{:.0f}"); macro(f"MunMedMASE{r['model']}", r["median_mase"], "{:.2f}"); macro(f"MunPninety{r['model']}", r["p90"], "{:.2f}")311 de = pl.read_csv(X / "dollar_errors_a.csv", infer_schema_length=None)312 rows = [[r["property_class"], LABELS.get(r["model"], r["model"]), f"{100*r['mae_log']:.1f}", f"{r['median_price']:,.0f}".replace(",", "{,}"), f"{r['median_abs_dollar_error']:,.0f}".replace(",", "{,}")] for r in de.iter_rows(named=True)]313 table("tab_dollar_a", ["Class", "Model", "MAE (log ×100)", "Median price (\\$)", "Median absolute error (\\$)"], rows, "Task A, horizons 1--3 months — economic magnitude of forecast errors", "tab:dollar_a", size="\\footnotesize")314 for r in de.iter_rows(named=True): macro(f"Dollar{r['model']}{r['property_class']}", r["median_abs_dollar_error"], "{:,.0f}")315 sr = pl.read_csv(X / "selection_rates.csv", infer_schema_length=None).filter(pl.col("quintile").is_not_null()); rows = []316 for cls in ["house", "condo", "plex"]:317 for dim, lab2 in [("fa_q", "floor area"), ("age_q", "age")]:318 x = sr.filter((pl.col("class") == cls) & (pl.col("dimension") == dim)).sort("quintile"); rows.append([f"{cls}, by {lab2} quintile (1 = smallest / newest)"] + [f"{v:.1f}" for v in x["annual_sale_rate_pct"].to_list()])319 for q, v in zip(x["quintile"].to_list(), x["annual_sale_rate_pct"].to_list()): macro(f"SaleRate{cls}{q}", v, "{:.1f}")320 table("tab_selection_rates", ["Class and dimension", "Q1", "Q2", "Q3", "Q4", "Q5"], rows, "Transaction selection — annual sale rate of the roll stock (\\% of units sold per year, 2021--2026) by quintile of floor area and of age", "tab:selection_rates", size="\\footnotesize")321322323KEYS_B = ["no_change", "mean_growth", "momentum", "ar2_cell", "ridge_panel", "lightgbm_panel", "tft_lite", "bhpf", "raref", "qanss", "dsgfn", "qrefa"]324325326def t_simulation():327 f = R / "simulation" / "summary.csv"328 if not ok(f): return329 d = pl.read_csv(f, infer_schema_length=None); scen = ["rw", "momentum", "regimes", "spillover", "sparse"]; scen = [s for s in scen if s in d["scenario"].to_list()]330 for lvl in ["municipality", "region", "province"]:331 rows = []332 for m in ORDER_B:333 x = d.filter((pl.col("model") == m) & (pl.col("level") == lvl) & (pl.col("h") == 3))334 if x.height == 0: continue335 r = [LABELS.get(m, m)]336 for s in scen:337 y = x.filter(pl.col("scenario") == s); r.append(f"{y['rel_mae'][0]:.2f} ({y['rel_mae_mcse'][0]:.2f})" if y.height else "—")338 if y.height: macro(f"Sim{lvl}{m}{s}", 100 * y["mae"][0], "{:.2f}"); macro(f"SimRel{lvl}{m}{s}", y["rel_mae"][0], "{:.2f}"); macro(f"SimCov{lvl}{m}{s}", 100 * y["cov90"][0], "{:.0f}")339 rows.append(r)340 table(f"tab_sim_{lvl}", ["Model"] + scen, rows, f"Simulation — MAE of 3-month index-change forecasts relative to the no-change forecast at the {lvl} level, by scenario (Monte Carlo s.e. over replications)", f"tab:sim_{lvl}", size="\\scriptsize",341 notes="rw: random walk with a constant drift; momentum: AR(1) drifts at every level; regimes: Markov-switching provincial drift; spillover: municipal drift follows the lagged regional drift; sparse: momentum with 35\\% of the sales. Values below 1 beat the no-change forecast. Classic and tree models: 10 replications; neural models: 6.")342343344if __name__ == "__main__":345 for fn in [t_data, t_results, t_ablation, t_robust, t_interpret, t_extra, t_simulation]:346 try: fn()347 except Exception as ex:348 import traceback; traceback.print_exc(); print("FAILED", fn.__name__)349 (T / "numbers.tex").write_text("% generated by scripts/50_tables.py\n" + "\n".join(f"\\newcommand{{\\{k}}}{{{v}}}" for k, v in sorted(MACROS.items())) + "\n")350 used = set()351 for f in list((C.ROOT / "paper" / "sections").glob("*.tex")) + list((C.ROOT / "paper" / "appendix").glob("*.tex")):352 used |= set(re.findall(r"\\([A-Z][A-Za-z]+)(?=[^A-Za-z]|$)", f.read_text()))353 skip = {"E", "Var", "DeclareMathOperator", "LARGE", "Large", "Delta", "WPabstract", "N"}354 missing = sorted(u for u in used - set(MACROS) - skip if not u.startswith("WP"))355 (T / "numbers_missing.tex").write_text("% must be empty in the final version\n" + "\n".join(f"\\providecommand{{\\{k}}}{{\\textbf{{??{k}??}}}}" for k in missing) + "\n")356 print(len(MACROS), "macros;", len(missing), "missing:", missing)357