TeX 62.3%
Python 37.6%
1#!/usr/bin/env python32# Author: Simon-Pierre Boucher — simon-pierre.boucher@uqo.ca3"""Score Task A and Task B predictions: point and probabilistic metrics by model × horizon (× level/class/region/size),4Diebold–Mariano and Giacomini–White tests against the reference benchmark, Model Confidence Sets per horizon, bootstrap CIs.5Usage: python3 scripts/30_evaluate.py a|b"""6from __future__ import annotations78import 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 C16from wp23 import evaluation as E1718TAUS = [0.05, 0.25, 0.5, 0.75, 0.95]192021def prob_block(d: pl.DataFrame, ycol: str) -> dict:22 if "q50" not in d.columns or d["q50"].null_count() == d.height:23 return {"crps": np.nan, "pinball_mean": np.nan, "cov90": np.nan, "width90": np.nan, "wis90": np.nan, "cov50": np.nan}24 y = d[ycol].to_numpy(); Q = np.column_stack([d[f"q{int(t*100):02d}"].to_numpy() for t in TAUS]); ok = np.all(np.isfinite(Q), axis=1) & np.isfinite(y)25 y, Q = y[ok], Q[ok]26 if len(y) == 0:27 return {"crps": np.nan, "pinball_mean": np.nan, "cov90": np.nan, "width90": np.nan, "wis90": np.nan, "cov50": np.nan}28 crps = float(np.mean(E.crps_from_quantiles(y, Q, np.array(TAUS)))); pb = float(np.mean([np.mean(E.pinball(y, Q[:, j], t)) for j, t in enumerate(TAUS)]))29 i90 = E.interval_scores(y, Q[:, 0], Q[:, 4], 0.10); i50 = E.interval_scores(y, Q[:, 1], Q[:, 3], 0.50)30 return {"crps": crps, "pinball_mean": pb, "cov90": i90["coverage"], "width90": i90["width"], "wis90": i90["interval_score"], "cov50": i50["coverage"]}313233def evaluate(pred: pl.DataFrame, ycol: str, keys: list[str], ref_model: str, group_extra: list[str] | None = None, tag: str = "") -> None:34 OUT = C.RESULTS / f"eval_{tag}"; OUT.mkdir(parents=True, exist_ok=True)35 pred = pred.filter(pl.col("pred").is_not_null() & pl.col("pred").is_finite())36 # naive reference for MASE: the model named ref_model37 ref = pred.filter(pl.col("model") == ref_model).select(keys + [pl.col("pred").alias("pred_ref")])38 j = pred.join(ref, on=keys, how="left")39 rows = []40 for gkeys in ([["model", "h"]] + ([["model", "h"] + [g] for g in group_extra] if group_extra else [])):41 for g, d in j.group_by(gkeys):42 y = d[ycol].to_numpy(); f = d["pred"].to_numpy(); fr = d["pred_ref"].to_numpy()43 sc = E.point_scores(y, f, fr if np.isfinite(fr).any() else None); sc.update(prob_block(d, ycol))44 rows.append({**dict(zip(gkeys, g)), "group": "|".join(gkeys[2:]) or "all", **sc})45 scores = pl.DataFrame(rows, infer_schema_length=None); scores.write_csv(OUT / "scores.csv")46 # tests per horizon on per-origin mean losses (time-series of forecast-origin averages → serial dependence handled by HAC/bootstrap)47 tests, mcs_rows = [], []48 HBAND = {"1-3": 2, "4-6": 5, "7-12": 9}49 for h, dh in j.group_by("h"):50 h = h[0]; hi = HBAND.get(str(h), None); hi = int(h) if hi is None else hi51 per = dh.with_columns(((pl.col("pred") - pl.col(ycol)) ** 2).alias("sq"), (pl.col("pred") - pl.col(ycol)).abs().alias("ab")).group_by(["model", "origin"]).agg(pl.col("sq").mean(), pl.col("ab").mean()).sort("origin")52 wide_sq = per.pivot(on="model", index="origin", values="sq").sort("origin"); wide_ab = per.pivot(on="model", index="origin", values="ab").sort("origin")53 models = [c for c in wide_sq.columns if c != "origin"]54 if ref_model not in models:55 continue56 for m in models:57 if m == ref_model:58 continue59 la = wide_sq[m].to_numpy(); lb = wide_sq[ref_model].to_numpy(); ok = np.isfinite(la) & np.isfinite(lb)60 dm = E.diebold_mariano(la[ok], lb[ok], h=hi); gw = E.giacomini_white(la[ok], lb[ok], h=hi)61 lo_, hi_ = E.block_bootstrap_ci(la[ok] - lb[ok], B=500, mean_block=max(2, hi))62 tests.append({"h": h, "model": m, "reference": ref_model, "n_origins": int(ok.sum()), "mean_sq_diff": dm["mean_diff"], "dm_stat": dm["dm"], "dm_p": dm["p"], "gw_stat": gw["stat"], "gw_p": gw["p"], "boot_ci_lo": lo_, "boot_ci_hi": hi_})63 L = np.column_stack([wide_sq[m].to_numpy() for m in models]); okr = np.all(np.isfinite(L), axis=1)64 if okr.sum() >= 8:65 mcs = E.model_confidence_set(L[okr], models, alpha=0.10, B=500, mean_block=max(2, hi), seed=0)66 for m in models:67 mcs_rows.append({"h": h, "model": m, "mcs_p": mcs["mcs_pvalues"].get(m, np.nan), "in_mcs_90": m in mcs["in_mcs"]})68 pl.DataFrame(tests, infer_schema_length=None).write_csv(OUT / "tests_vs_reference.csv"); pl.DataFrame(mcs_rows, infer_schema_length=None).write_csv(OUT / "mcs.csv")69 # calibration (PIT-like): share of outcomes below each quantile70 cal = []71 for (m, h), d in j.group_by(["model", "h"]):72 if "q50" in d.columns and d["q50"].null_count() < d.height:73 y = d[ycol].to_numpy()74 for t in TAUS:75 q = d[f"q{int(t*100):02d}"].to_numpy(); ok = np.isfinite(q) & np.isfinite(y)76 cal.append({"model": m, "h": h, "tau": t, "share_below": float(np.mean(y[ok] < q[ok])) if ok.any() else np.nan})77 pl.DataFrame(cal, infer_schema_length=None).write_csv(OUT / "calibration.csv")78 print(scores.filter(pl.col("group") == "all").sort(["h", "mae"]).select(["h", "model", "n", "mae", "rmse", "mase", "bias", "crps", "cov90"]))798081def task_a() -> None:82 pred = pl.concat([pl.read_parquet(f) for f in sorted((C.RESULTS / "task_a").glob("predictions_*.parquet"))], how="diagonal")83 tg = pl.concat([pl.read_parquet(p) for p in sorted((C.RESULTS / "task_a").glob("targets_*.parquet"))])84 sales = pl.read_parquet(C.WP22_ANALYTIC, columns=["tx_id", "property_class", "region_code", "mun_code", "price"])85 size = sales.group_by("mun_code").agg(pl.len().alias("mun_sales")).with_columns(pl.when(pl.col("mun_sales") >= 5000).then(pl.lit("large")).when(pl.col("mun_sales") >= 500).then(pl.lit("medium")).otherwise(pl.lit("small")).alias("mun_size"))86 j = pred.join(tg, on=["tx_id", "h", "origin"]).join(sales, on="tx_id").join(size, on="mun_code")87 j = j.with_columns(pl.when(pl.col("h") <= 3).then(pl.lit("1-3")).when(pl.col("h") <= 6).then(pl.lit("4-6")).otherwise(pl.lit("7-12")).alias("hband"),88 pl.when(pl.col("price") < 250000).then(pl.lit("<250k")).when(pl.col("price") < 500000).then(pl.lit("250-500k")).otherwise(pl.lit(">500k")).alias("segment"))89 j = j.with_columns(pl.col("hband").alias("h")) # evaluate by horizon band90 evaluate(j, "y", ["tx_id", "h", "origin"], "naive_mun_median", ["property_class", "region_code", "mun_size", "segment"], tag="a")919293def task_b() -> None:94 pred = pl.concat([pl.read_parquet(f) for f in sorted((C.RESULTS / "task_b").glob("predictions_*.parquet"))], how="diagonal")95 pred = pred.with_columns(pl.when(pl.col("level") == "municipality").then(pl.when(pl.col("n_last") >= 50).then(pl.lit("mun_large")).otherwise(pl.lit("mun_small"))).otherwise(pl.col("level")).alias("size_level"))96 evaluate(pred, "target", ["class", "level", "geo_id", "origin", "h"], "no_change", ["level", "class", "size_level", "region_code"], tag="b")979899if __name__ == "__main__":100 {"a": task_a, "b": task_b}[sys.argv[1]]()101