spb/wp23_uqo

1#!/usr/bin/env python32# Author: Simon-Pierre Boucher — simon-pierre.boucher@uqo.ca3"""Summarise the forecasting Monte Carlo: MAE / RMSE / 90% coverage by scenario × level × horizon × model, with Monte Carlo s.e. over replications."""4from __future__ import annotations56import sys7from pathlib import Path89sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src"))10import numpy as np11import polars as pl1213from wp23 import config as C1415OUT = C.RESULTS / "simulation"16files = sorted(OUT.glob("predictions_*.parquet"))17p = pl.concat([pl.read_parquet(f) for f in files], how="diagonal")18p = p.with_columns((pl.col("pred") - pl.col("target")).abs().alias("ae"), ((pl.col("pred") - pl.col("target")) ** 2).alias("se"), ((pl.col("target") >= pl.col("q05")) & (pl.col("target") <= pl.col("q95"))).cast(pl.Float64).alias("cov90"))19rep = p.group_by(["scenario", "rep", "level", "h", "model"]).agg(pl.col("ae").mean().alias("mae"), pl.col("se").mean().sqrt().alias("rmse"), pl.col("cov90").mean().alias("cov90"), pl.len().alias("n"))20nc = rep.filter(pl.col("model") == "no_change").select(["scenario", "rep", "level", "h", pl.col("mae").alias("mae_nc")])21rep = rep.join(nc, on=["scenario", "rep", "level", "h"], how="left").with_columns((pl.col("mae") / pl.col("mae_nc")).alias("rel_mae"))22summ = rep.group_by(["scenario", "level", "h", "model"]).agg(pl.col("rep").n_unique().alias("reps"), pl.col("mae").mean(), (pl.col("mae").std() / pl.col("mae").count().sqrt()).alias("mae_mcse"), pl.col("rmse").mean(), pl.col("cov90").mean(), pl.col("rel_mae").mean().alias("rel_mae"), (pl.col("rel_mae").std() / pl.col("rel_mae").count().sqrt()).alias("rel_mae_mcse")).sort(["scenario", "level", "h", "mae"])23summ.write_csv(OUT / "summary.csv"); pl.Config.set_tbl_rows(100); print(summ.filter((pl.col("level") == "municipality") & (pl.col("h") == 3)).select(["scenario", "model", "reps", "mae", "mae_mcse", "rel_mae", "cov90"]))24