Is the market's probability any good?

panel = cu.make_prediction_markets(n_markets=240, steps=48, seed=11)
for name, table in panel.items():
    db.create_table(name, table.schema, time_column="ts_init")
    db.append(name, table, note="panel load")
db.snapshot("panel-v1", tables=list(panel), note="calibration study input")
resolutions = panel["resolutions"]
print(f"resolutions: {resolutions.num_rows:,} rows x {resolutions.num_columns} columns")
resolutions.to_pandas().head()

stamps = sorted({value.as_py() for value in panel["book_deltas"].column("ts_init")})
decision = stamps[len(stamps) // 2]
decision_ns = int(pd.Timestamp(decision, tz="UTC").value)
print(f"session runs {stamps[0]} .. {stamps[-1]}")
print(f"decision instant: {decision} ({decision_ns} ns)")

forecasts = db.sql(
    f"""
    SELECT instrument_id,
           (max(CASE WHEN side = 'buy' THEN price END)
          + max(CASE WHEN side = 'sell' THEN price END)) / 2 AS p_yes
    FROM h5i('book_deltas', 'panel-v1')
    WHERE outcome = 0 AND ts_init = to_timestamp_nanos({decision_ns})
    GROUP BY instrument_id
    ORDER BY instrument_id
    """
).to_pandas()
print(f"\n{len(forecasts)} forecasts")
print(f"quoted probability spans {forecasts.p_yes.min():.3f} .. {forecasts.p_yes.max():.3f}")

knowable = db.read("resolutions", time_end=decision_ns)
print(f"resolutions knowable at the decision instant: {knowable.num_rows}")
print(f"resolutions available now:                    {db.read('resolutions').num_rows}")
assert knowable.num_rows == 0

truth = cu.market_truth(panel).to_pandas()
scored = forecasts.merge(truth, on="instrument_id", validate="one_to_one")
scored["y"] = scored.yes_won.astype(float)
print(f"{len(scored)} scored markets, base rate {scored.y.mean():.3f}")
scored.head()

EDGES = [0.0, 0.10, 0.20, 0.35, 0.50, 0.65, 0.80, 0.90, 1.0]
scored["bucket"] = pd.cut(scored.p_yes, EDGES)
reliability = (
    scored.groupby("bucket", observed=True)
    .agg(n=("y", "size"), mean_quote=("p_yes", "mean"), realized=("y", "mean"))
    .assign(gap_pp=lambda f: (f.realized - f.mean_quote) * 100)
)
print(reliability.round(3).to_string())

low = reliability.head(3)
high = reliability.tail(3)
print(f"three cheapest buckets, mean gap: {np.average(low.gap_pp, weights=low.n):+.2f} pp")
print(f"three richest  buckets, mean gap: {np.average(high.gap_pp, weights=high.n):+.2f} pp")

fig, ax = plt.subplots(figsize=(6, 6))
ax.plot([0, 1], [0, 1], color="grey", lw=1, ls="--", label="perfect calibration")
ax.scatter(reliability.mean_quote, reliability.realized, s=reliability.n * 3,
           color="#c0392b", zorder=3, label="realized frequency")
for row in reliability.itertuples():
    ax.annotate(f"n={row.n}", (row.mean_quote, row.realized),
                textcoords="offset points", xytext=(6, -10), fontsize=7)
ax.set_title("Reliability curve, 240 binary markets")
ax.set_xlabel("quoted probability at the decision instant")
ax.set_ylabel("realized YES frequency")
ax.set_xlim(0, 1)
ax.set_ylim(0, 1)
ax.legend(loc="upper left", fontsize=8)
fig.tight_layout()

def brier_decomposition(p: pd.Series, y: pd.Series, edges: list[float]) -> dict[str, float]:
    """Murphy's three-term decomposition of the Brier score."""
    frame = pd.DataFrame({"p": p, "y": y})
    frame["bin"] = pd.cut(frame.p, edges)
    n = len(frame)
    base = frame.y.mean()
    groups = frame.groupby("bin", observed=True)
    counts = groups.size()
    reliability = float((counts * (groups.p.mean() - groups.y.mean()) ** 2).sum() / n)
    resolution = float((counts * (groups.y.mean() - base) ** 2).sum() / n)
    uncertainty = float(base * (1.0 - base))
    return {
        "brier": float(((frame.p - frame.y) ** 2).mean()),
        "reliability": reliability,
        "resolution": resolution,
        "uncertainty": uncertainty,
        "identity": reliability - resolution + uncertainty,
    }


def log_loss(p: pd.Series, y: pd.Series, floor: float = 1e-6) -> float:
    clipped = p.clip(floor, 1 - floor)
    return float(-(y * np.log(clipped) + (1 - y) * np.log(1 - clipped)).mean())


parts = brier_decomposition(scored.p_yes, scored.y, EDGES)
for key, value in parts.items():
    print(f"  {key:12} {value:.4f}")
print(f"  {'log loss':12} {log_loss(scored.p_yes, scored.y):.4f}")
assert abs(parts["brier"] - parts["identity"]) < 0.02

BIAS = 0.12
scored["p_recal"] = (0.5 + (scored.p_yes - 0.5) / (1.0 - BIAS)).clip(0.01, 0.99)
rows = []
for label, forecast in (
    ("market quote", scored.p_yes),
    ("base rate", pd.Series(scored.y.mean(), index=scored.index)),
    ("recalibrated", scored.p_recal),
):
    part = brier_decomposition(forecast, scored.y, EDGES)
    rows.append(
        {
            "forecast": label,
            "brier": part["brier"],
            "reliability": part["reliability"],
            "resolution": part["resolution"],
            "log_loss": log_loss(forecast, scored.y),
        }
    )
print(pd.DataFrame(rows).set_index("forecast").round(4).to_string())

horizon = []
for index in (4, 16, 28, 40):
    at = stamps[index]
    at_ns = int(pd.Timestamp(at, tz="UTC").value)
    frame = db.sql(
        f"""
        SELECT instrument_id,
               (max(CASE WHEN side = 'buy' THEN price END)
              + max(CASE WHEN side = 'sell' THEN price END)) / 2 AS p_yes
        FROM h5i('book_deltas', 'panel-v1')
        WHERE outcome = 0 AND ts_init = to_timestamp_nanos({at_ns})
        GROUP BY instrument_id
        """
    ).to_pandas().merge(truth, on="instrument_id")
    frame["y"] = frame.yes_won.astype(float)
    part = brier_decomposition(frame.p_yes, frame.y, EDGES)
    horizon.append(
        {
            "instant": str(at.time()),
            "steps_to_expiry": len(stamps) - 1 - index,
            "brier": part["brier"],
            "reliability": part["reliability"],
            "resolution": part["resolution"],
        }
    )
print(pd.DataFrame(horizon).set_index("instant").round(4).to_string())