TeX 55.1%
Python 44.2%
1#!/usr/bin/env python32# Author: Simon-Pierre Boucher — simon-pierre.boucher@uqo.ca3"""Phase 4 — Monte Carlo experiments on synthetic markets with known latent paths (protocol §12).45Usage: python3 scripts/20_simulation.py [--reps 100] [--workers 14] [--scenarios a,b,c]6Writes results/simulation/metrics.parquet (one row per scenario × rep × method × level × cell) and summary.csv7(Monte Carlo means with Monte Carlo standard errors).8"""9from __future__ import annotations1011import os1213os.environ.setdefault("OMP_NUM_THREADS", "1"); os.environ.setdefault("OPENBLAS_NUM_THREADS", "1"); os.environ.setdefault("MKL_NUM_THREADS", "1")14import argparse15import sys16import time17import warnings18from multiprocessing import Pool19from pathlib import Path2021sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src"))22warnings.filterwarnings("ignore")23import numpy as np24import polars as pl2526from sharp_hpi import config as C27from sharp_hpi import estimators as E28from sharp_hpi.hedonic import assign_loc_fe, build_features29from sharp_hpi.qhpi import qhpi_cell, qhpi_fullpool_residuals, qhpi_province30from sharp_hpi.sharp import cell_log_index, run_sharp31from sharp_hpi.simulate import SCENARIOS, Scenario, simulate32from sharp_hpi.statespace import fit_cell3334OUT = C.RESULTS / "simulation"35OUT.mkdir(parents=True, exist_ok=True)363738def _metrics(est: np.ndarray, se: np.ndarray | None, truth: np.ndarray, est_rt: np.ndarray | None = None) -> dict:39 ok = ~np.isnan(est)40 if ok.sum() < 13:41 return {"rmse_ret": np.nan, "bias_ret": np.nan, "rmse_lvl": np.nan, "cov95": np.nan, "share_est": float(ok.mean()), "mae_ret": np.nan, "corr_ret": np.nan, "lag_xcorr": np.nan, "revision_rmse": np.nan, "cov95_ret": np.nan}42 e = est - np.nanmean(est[ok]); t = truth - np.mean(truth[ok])43 de = np.diff(est); dt = np.diff(truth); okd = ~np.isnan(de)44 out = {"rmse_ret": float(np.sqrt(np.mean((de[okd] - dt[okd]) ** 2))), "mae_ret": float(np.mean(np.abs(de[okd] - dt[okd]))),45 "bias_ret": float(np.mean(de[okd] - dt[okd])), "rmse_lvl": float(np.sqrt(np.mean((e[ok] - t[ok]) ** 2))), "share_est": float(ok.mean())}46 out["corr_ret"] = float(np.corrcoef(de[okd], dt[okd])[0, 1]) if okd.sum() > 5 and np.std(de[okd]) > 0 else np.nan47 # lag: argmax of cross-correlation of 3-month returns within [-3, 3]48 r3e = est[3:] - est[:-3]; r3t = truth[3:] - truth[:-3]; best, bl = -2, np.nan49 for L in range(-3, 4):50 a = r3e[max(0, L):len(r3e) + min(0, L)]; b = r3t[max(0, -L):len(r3t) + min(0, -L)]; m = ~np.isnan(a)51 if m.sum() > 8 and np.std(a[m]) > 0 and np.std(b[m]) > 0:52 c = np.corrcoef(a[m], b[m])[0, 1]53 if c > best:54 best, bl = c, L55 out["lag_xcorr"] = bl56 if se is not None:57 s = np.asarray(se, float); oks = ok & ~np.isnan(s)58 out["cov95"] = float(np.mean(np.abs(e[oks] - t[oks]) <= 1.96 * s[oks])) if oks.sum() > 5 else np.nan59 ds = np.sqrt(s[1:] ** 2 + s[:-1] ** 2); okr = okd & ~np.isnan(ds) # conservative (ignores covariance)60 out["cov95_ret"] = float(np.mean(np.abs(de[okr] - dt[okr]) <= 1.96 * ds[okr])) if okr.sum() > 5 else np.nan61 else:62 out["cov95"] = np.nan; out["cov95_ret"] = np.nan63 out["revision_rmse"] = float(np.sqrt(np.nanmean((np.diff(est_rt) - de) ** 2))) if est_rt is not None else np.nan64 return out656667def one_rep(args):68 name, rep = args69 sc = SCENARIOS[name]; sc = Scenario(**{**sc.__dict__, "seed": 1000 + rep})70 t0 = time.time()71 sim = simulate(sc)72 grid = sim.months73 df = build_features(sim.df)74 rows = []75 truth = {("municipality", c): sim.truth_mun[i] for i, c in enumerate(sim.mun_codes)}76 truth.update({("region", f"{r:02d}"): sim.truth_region[r] for r in range(sim.truth_region.shape[0])})77 truth.update({("mrc", f"R{m:02d}"): sim.truth_mrc[m] for m in range(sim.truth_mrc.shape[0])})78 truth[("province", "QC")] = sim.truth_prov7980 def add(method, level, gid, est, se=None, est_rt=None, secs=np.nan):81 m = _metrics(np.asarray(est, float), se, truth[(level, gid)], est_rt)82 rows.append({"scenario": name, "rep": rep, "method": method, "level": level, "geo_id": gid, "secs": secs, **m})8384 # --- province ---85 prov = assign_loc_fe(df, min_n=15).with_columns(pl.col("mrc_code").alias("loc_coarse"))86 for fn, kw, nm in [(E.median_index, {}, "A_median"), (E.stratified_median, {}, "B_stratified"), (E.pooled_td, {}, "C_pooled_td"),87 (E.rolling_td, {"window": 13, "splice": "mean"}, "D_rtd13"), (E.hedonic_imputation, {"kind": "fisher"}, "E_fisher_chained"),88 (E.repeat_sales, {}, "H_repeat_sales")]:89 t = time.time(); s = fn(prov, grid, **kw); add(nm, "province", "QC", s.log_index, s.se, secs=time.time() - t)90 t = time.time(); sK = qhpi_province(prov, grid); add("K_qhpi", "province", "QC", sK.log_index, secs=time.time() - t)91 qres, qsig = qhpi_fullpool_residuals(prov, grid)92 variants = {93 "S_sharp": dict(), "S_noEB": dict(eb=False), "S_noTV": dict(time_varying=False), "S_LLT": dict(model="llt"),94 "S_norobust": dict(robust=False), "S_coarseloc": dict(min_loc_n=10**9),95 }96 sharp_runs = {}97 for v, kw in variants.items():98 t = time.time(); r = run_sharp(df, grid, **kw); sharp_runs[v] = r99 add(v, "province", "QC", r.stage1.mu, r.stage1.mu_se, secs=time.time() - t)100 # --- sub-provincial cells ---101 for level, key, codes in [("region", "region_code", [f"{r:02d}" for r in range(sim.truth_region.shape[0])]),102 ("mrc", "mrc_code", [f"R{m:02d}" for m in range(sim.truth_mrc.shape[0])]),103 ("municipality", "mun_code", sim.mun_codes)]:104 for gid in codes:105 dc = df.filter(pl.col(key) == gid)106 if dc.height == 0:107 continue108 n_by_m = E._counts(dc, grid); med_n = np.median(n_by_m[n_by_m > 0]) if (n_by_m > 0).any() else 0109 s = E.median_index(dc, grid); add("A_median", level, gid, s.log_index)110 if dc.height >= 60:111 s = E.stratified_median(dc, grid); add("B_stratified", level, gid, s.log_index)112 dcc = assign_loc_fe(dc, min_n=15)113 if dc.height >= 120 and dc["month"].n_unique() >= 24:114 try:115 s = E.pooled_td(dcc, grid); add("C_pooled_td", level, gid, s.log_index, s.se)116 except Exception:117 pass118 if med_n >= 10:119 try:120 s = E.rolling_td(dcc, grid, 13, "mean"); add("D_rtd13", level, gid, s.log_index)121 except Exception:122 pass123 s = E.repeat_sales(dc, grid); add("H_repeat_sales", level, gid, s.log_index)124 try:125 path, var, meth, gain = qhpi_cell(dcc, qres.filter(pl.col(key) == gid), sK.log_index, grid, qsig); add("K_qhpi", level, gid, path, np.sqrt(var))126 except Exception:127 pass128 for v, r in sharp_runs.items():129 if (level, gid) in r.cells:130 p, var = cell_log_index(r, level, gid, "smoothed"); pf, _ = cell_log_index(r, level, gid, "filtered")131 add(v, level, gid, p, np.sqrt(var), est_rt=pf)132 if v == "S_sharp":133 add("S_sharp_filtered", level, gid, pf, np.sqrt(_))134 # aggregate coherence: stock-weighted bottom-up of municipal SHARP vs province path135 r = sharp_runs["S_sharp"]; w = sim.stock / sim.stock.sum()136 bu = np.zeros(len(grid))137 for i, c in enumerate(sim.mun_codes):138 if ("municipality", c) in r.cells:139 bu += w[i] * cell_log_index(r, "municipality", c, "smoothed")[0]140 add("S_sharp_bottomup", "province", "QC", bu)141 for r_ in rows:142 for k_, v_ in r_.items():143 if k_ not in ("scenario", "method", "level", "geo_id"):144 r_[k_] = float(v_) if v_ is not None else float("nan")145 out = pl.DataFrame(rows, infer_schema_length=None).with_columns(pl.lit(time.time() - t0).alias("rep_secs"), pl.lit(sim.df.height).alias("n_sales"))146 return out147148149def main() -> None:150 ap = argparse.ArgumentParser(); ap.add_argument("--reps", type=int, default=100); ap.add_argument("--workers", type=int, default=14)151 ap.add_argument("--scenarios", type=str, default=",".join(SCENARIOS)); ap.add_argument("--tag", type=str, default="")152 a = ap.parse_args()153 jobs = [(s, r) for s in a.scenarios.split(",") for r in range(a.reps)]154 t0 = time.time()155 with Pool(a.workers) as pool:156 parts = []157 for i, part in enumerate(pool.imap_unordered(one_rep, jobs)):158 parts.append(part)159 if (i + 1) % 10 == 0:160 print(f"{i + 1}/{len(jobs)} reps, {time.time() - t0:.0f}s", flush=True)161 pl.concat(parts).write_parquet(OUT / f"metrics{a.tag}.parquet")162 m = pl.concat(parts); m.write_parquet(OUT / f"metrics{a.tag}.parquet")163 summ = m.group_by(["scenario", "method", "level"]).agg(164 pl.len().alias("cells"), pl.col("rep").n_unique().alias("reps"),165 *[pl.col(c).mean().alias(c) for c in ["rmse_ret", "mae_ret", "bias_ret", "rmse_lvl", "cov95", "cov95_ret", "corr_ret", "lag_xcorr", "revision_rmse", "share_est", "secs"]],166 *[(pl.col(c).std() / pl.col(c).count().sqrt()).alias(c + "_mcse") for c in ["rmse_ret", "bias_ret", "rmse_lvl", "cov95"]],167 ).sort(["scenario", "level", "rmse_ret"])168 summ.write_csv(OUT / f"summary{a.tag}.csv")169 print("done", time.time() - t0)170171172if __name__ == "__main__":173 main()174