TeX 62.3%
Python 37.6%
1#!/usr/bin/env python32# Author: Simon-Pierre Boucher — simon-pierre.boucher@uqo.ca3"""Forecasting Monte Carlo on synthetic markets (protocol §14). Markets are generated with the WP22 simulator (known latent4quality-constant paths per municipality, nested hierarchy, heterogeneous sizes) under scenarios that differ in *forecastability*:5 rw pure random-walk latent paths (no predictable component: the no-change forecast is optimal)6 momentum latent drift follows an AR(1) (persistent momentum)7 regimes drift switches between two regimes by a Markov chain (regime adaptation should help)8 spillover each municipality's drift partly follows its region's lagged drift (spatial information is predictive)9 sparse momentum with 35 % of the sales (real-time estimates noisier)10The "real-time" index at an origin is the WP22 SHARP filtered estimate on the synthetic sales (so estimation noise and11pooling are realistic); the target is the true latent change. Each forecaster is evaluated across origins and horizons.12Usage: python3 scripts/40_simulation.py --reps 20 --family classic|torch|trees"""13from __future__ import annotations1415import argparse16import sys17import time18from pathlib import Path1920sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src"))21import numpy as np22import polars as pl2324from wp23 import config as C25from wp23.data import month_add2627sys.path.insert(0, str(C.WP22_SRC))28from sharp_hpi.hedonic import build_features # noqa: E40229from sharp_hpi.sharp import cell_log_index, run_sharp # noqa: E40230from sharp_hpi.simulate import Scenario, simulate # noqa: E4023132OUT = C.RESULTS / "simulation"; OUT.mkdir(parents=True, exist_ok=True)33T = 72; LAGS = 12; HORIZONS = (1, 3, 6, 12)343536def latent_paths(sc: Scenario, kind: str, rng: np.random.Generator, n_regions: int, n_mun: int, mun_region: np.ndarray):37 """Override the WP22 random-walk latent paths with forecastable dynamics. Returns (mu, dR, dU) arrays."""38 T_ = sc.T39 if kind == "rw":40 drift_p = np.zeros(T_) + 0.003; drift_r = np.zeros((n_regions, T_)); drift_u = np.zeros((n_mun, T_))41 elif kind in ("momentum", "sparse"):42 drift_p = np.zeros(T_); drift_r = np.zeros((n_regions, T_)); drift_u = np.zeros((n_mun, T_))43 for t in range(1, T_):44 drift_p[t] = 0.85 * drift_p[t - 1] + rng.normal(0, 0.003)45 drift_r[:, t] = 0.8 * drift_r[:, t - 1] + rng.normal(0, 0.002, n_regions)46 drift_u[:, t] = 0.7 * drift_u[:, t - 1] + rng.normal(0, 0.003, n_mun)47 elif kind == "regimes":48 P = np.array([[0.93, 0.07], [0.10, 0.90]]); s = 0; drift_p = np.zeros(T_); states = np.zeros(T_, int)49 for t in range(T_):50 s = rng.choice(2, p=P[s]); states[t] = s; drift_p[t] = (0.012 if s == 0 else -0.006) + rng.normal(0, 0.002)51 drift_r = np.zeros((n_regions, T_)); drift_u = np.zeros((n_mun, T_))52 for t in range(1, T_):53 drift_r[:, t] = 0.6 * drift_r[:, t - 1] + rng.normal(0, 0.002, n_regions); drift_u[:, t] = 0.5 * drift_u[:, t - 1] + rng.normal(0, 0.003, n_mun)54 elif kind == "spillover":55 drift_p = np.full(T_, 0.004); drift_r = np.zeros((n_regions, T_)); drift_u = np.zeros((n_mun, T_))56 for t in range(1, T_):57 drift_r[:, t] = 0.8 * drift_r[:, t - 1] + rng.normal(0, 0.003, n_regions)58 drift_u[:, t] = 0.3 * drift_u[:, t - 1] + 0.6 * drift_r[mun_region, t - 1] + rng.normal(0, 0.002, n_mun) # municipalities follow their region with a lag59 mu = np.cumsum(drift_p + rng.normal(0, 0.004, T_)); mu -= mu[0]60 dR = np.cumsum(drift_r + rng.normal(0, 0.002, (n_regions, T_)), axis=1); dR -= dR[:, [0]]61 dU = np.cumsum(drift_u + rng.normal(0, 0.003, (n_mun, T_)), axis=1); dU -= dU[:, [0]]62 return mu, dR, dU636465def one_market(kind: str, rep: int):66 rng = np.random.default_rng(5000 + rep)67 sc = Scenario(name=kind, T=T, seed=5000 + rep, scale_sales=0.35 if kind == "sparse" else 1.0, tau_prov=1e-4, tau_region=1e-4, tau_mrc=1e-4, tau_mun=1e-4, tau_mun_hetero=False)68 sim = simulate(sc)69 n_mun = len(sim.mun_codes); n_regions = sim.truth_region.shape[0]70 mu, dR, dU = latent_paths(sc, kind, rng, n_regions, n_mun, sim.mun_region)71 L_new = mu[None, :] + dR[sim.mun_region] + dU72 # re-price the sales: replace the simulator's latent component by the new one (log_price += L_new − L_old)73 old = sim.truth_mun; mi = {c: i for i, c in enumerate(sim.mun_codes)}; months = sim.months74 df = sim.df.with_columns(pl.col("mun_code").replace_strict(mi, default=None).alias("_u"), pl.col("month").replace_strict({m: i for i, m in enumerate(months)}, default=None).alias("_t"))75 u = df["_u"].to_numpy(); t = df["_t"].to_numpy()76 df = df.with_columns((pl.col("log_price") + pl.Series(L_new[u, t] - old[u, t])).alias("log_price")).with_columns(pl.col("log_price").exp().alias("price")).drop(["_u", "_t"])77 w = sim.stock / sim.stock.sum(); truth_reg = np.vstack([(w[sim.mun_region == r][:, None] * L_new[sim.mun_region == r]).sum(0) / w[sim.mun_region == r].sum() for r in range(n_regions)])78 truth_prov = (w[:, None] * L_new).sum(0)79 return df, L_new, truth_reg, truth_prov, sim808182def build_sim_panel(kind: str, rep: int) -> pl.DataFrame:83 df, L, truth_reg, truth_prov, sim = one_market(kind, rep)84 f = build_features(df); months = sim.months85 truth = {("municipality", c): L[i] for i, c in enumerate(sim.mun_codes)}; truth.update({("region", f"{r:02d}"): truth_reg[r] for r in range(truth_reg.shape[0])}); truth[("province", "QC")] = truth_prov86 rows = []87 origins = [months[i] for i in range(24, T - 1)]88 for o in origins:89 oi = months.index(o); avail = f.filter(pl.col("month") <= o)90 r = run_sharp(avail, months[: oi + 1], realtime_parent=True)91 cells = [("province", "QC", r.stage1.mu)] + [(lvl, gid, cell_log_index(r, lvl, gid, "filtered")[0]) for (lvl, gid) in r.cells if lvl in ("region", "municipality")]92 for lvl, gid, path in cells:93 if len(path) < LAGS + 1 or (lvl, gid) not in truth:94 continue95 H = path[-(LAGS + 1):]; n = np.array([avail.filter((pl.col("month") == m) & ((pl.col("mun_code") == gid) if lvl == "municipality" else (pl.col("region_code") == gid) if lvl == "region" else True)).height for m in months[oi - LAGS: oi + 1]], float)96 se = np.zeros(LAGS + 1) + 0.0197 feats = {"rt_level": H[-1], "d1": H[-1] - H[-2], "d3": H[-1] - H[-4], "d6": H[-1] - H[-7], "d12": H[-1] - H[0], "vol6": np.std(np.diff(H[-7:])), "vol12": np.std(np.diff(H)),98 "n_last": n[-1], "n_mean12": n.mean(), "n_ratio": (n[-3:].mean() + 1) / (n.mean() + 1), "se_last": 0.01, "accel": (H[-1] - H[-4]) - (H[-4] - H[-7])}99 tr = truth[(lvl, gid)]100 for h in HORIZONS:101 if oi + h >= T:102 continue103 rows.append({"class": "house", "level": lvl, "geo_id": gid, "origin": o, "h": h, "month_of_year": int(o[5:7]), "origin_t": oi, "hist": list(H), **feats,104 "target": float(tr[oi + h] - tr[oi]), "y_final_o": float(tr[oi]), "region_code": gid if lvl == "region" else (f"{sim.mun_region[sim.mun_codes.index(gid)]:02d}" if lvl == "municipality" else None),105 "cx": float(df.filter(pl.col("mun_code") == gid)["x_m"].mean()) if lvl == "municipality" else 0.0, "cy": float(df.filter(pl.col("mun_code") == gid)["y_m"].mean()) if lvl == "municipality" else 0.0,106 "grade_target": "A", "n_target": 0, "scenario": kind, "rep": rep})107 p = pl.DataFrame(rows, infer_schema_length=None)108 reg = p.filter(pl.col("level") == "region").select(["origin", "h", pl.col("geo_id").alias("region_code"), pl.col("d3").alias("reg_d3"), pl.col("d12").alias("reg_d12"), pl.col("vol12").alias("reg_vol12")]).unique()109 prov = p.filter(pl.col("level") == "province").select(["origin", "h", pl.col("d3").alias("prov_d3"), pl.col("d12").alias("prov_d12"), pl.col("vol12").alias("prov_vol12")]).unique()110 return p.join(reg, on=["origin", "h", "region_code"], how="left").join(prov, on=["origin", "h"], how="left")111112113def main() -> None:114 ap = argparse.ArgumentParser(); ap.add_argument("--reps", type=int, default=20); ap.add_argument("--family", type=str, default="classic"); ap.add_argument("--build", action="store_true"); ap.add_argument("--scenarios", type=str, default="rw,momentum,regimes,spillover,sparse"); ap.add_argument("--horizons", type=str, default="1,3,6,12")115 a = ap.parse_args()116 if a.build:117 for kind in a.scenarios.split(","):118 for rep in range(a.reps):119 t0 = time.time(); p = build_sim_panel(kind, rep); p.write_parquet(OUT / f"panel_{kind}_{rep}.parquet"); print("built", kind, rep, p.shape, f"{time.time() - t0:.0f}s", flush=True)120 return121 if a.family == "torch":122 from wp23.models_b_torch import registry_torch_fast as registry123 else:124 from wp23.models_b import registry125 rows = []126 for kind in a.scenarios.split(","):127 for rep in range(a.reps):128 f = OUT / f"panel_{kind}_{rep}.parquet"129 if not f.exists():130 continue131 p = pl.read_parquet(f)132 for h in [int(x) for x in a.horizons.split(",")]:133 ph = p.filter(pl.col("h") == h); origins = sorted(ph["origin"].unique().to_list())134 for k in range(0, len(origins), 6):135 o = origins[k]; test_o = origins[k:k + 6]136 tr = ph.filter(pl.col("origin").map_elements(lambda x: month_add(x, h), return_dtype=pl.Utf8) <= o)137 if tr["origin"].n_unique() < 8:138 continue139 te = ph.filter(pl.col("origin").is_in(test_o))140 for m in registry(a.family, fast=True):141 try:142 m.fit(tr); pr = m.predict(te); Q = m.predict_quantiles(te)143 except Exception as ex:144 print("FAILED", kind, rep, h, m.name, repr(ex), flush=True); continue145 y = te["target"].to_numpy()146 rows.append(te.select(["level", "geo_id", "origin", "h", "target"]).with_columns(pl.lit(kind).alias("scenario"), pl.lit(rep).alias("rep"), pl.lit(m.name).alias("model"), pl.Series("pred", np.asarray(pr, dtype=np.float64)), pl.Series("q05", np.asarray(Q[:, 0], dtype=np.float64)), pl.Series("q95", np.asarray(Q[:, 4], dtype=np.float64))))147 print(kind, rep, "done", flush=True)148 pl.concat(rows).write_parquet(OUT / f"predictions_{a.family}.parquet")149 print("done")150151152if __name__ == "__main__":153 main()154