"""Reproducible calibration analysis of resolved Polymarket binary markets (Prediction Market Math).

Input:  polymarket-resolved-markets-2023-2026.csv, collected from the public Polymarket Gamma and CLOB APIs:
        resolved Yes/No markets with volume >= $50,000 that closed between 13 March 2023 and 5 October 2026; "Yes" price
        24 hours and 7 days before close (last hourly print at or before that time).
Output: printed tables and polymarket_calibration_results.json

Notes
- Markets from multi-outcome events ("neg_risk") are correlated within an event; results are shown with and without them.
- ROI for buying Yes at price p: (outcome - p) / p per market; for buying No: ((1 - outcome) - (1 - p)) / (1 - p). Fees ignored.
Run: python polymarket_calibration_analysis.py
"""
import csv, json, math, os, statistics

HERE = os.path.dirname(os.path.abspath(__file__))
DATA = os.path.join(HERE, "polymarket-resolved-markets-2023-2026.csv")
BINS = [(i / 10, (i + 1) / 10) for i in range(10)]

def load():
    with open(DATA, encoding="utf-8") as f:
        return list(csv.DictReader(f))

def ci95(values):
    if len(values) < 2: return [None, None]
    m = statistics.fmean(values); se = statistics.pstdev(values) / math.sqrt(len(values))
    return [round((m - 1.96 * se) * 100, 2), round((m + 1.96 * se) * 100, 2)]

def analyse(rows, col):
    pairs = [(float(r[col]), int(r["outcome_yes"])) for r in rows if r[col] not in ("", None)]
    pairs = [(min(max(p, 0.0), 1.0), y) for p, y in pairs]
    n = len(pairs); base = statistics.fmean(y for _, y in pairs)
    brier = statistics.fmean((p - y) ** 2 for p, y in pairs)
    logloss = statistics.fmean(-math.log(max(1e-6, p if y else 1 - p)) for p, y in pairs)
    table, rel, res = [], 0.0, 0.0
    for lo, hi in BINS:
        b = [(p, y) for p, y in pairs if lo <= p < hi or (hi == 1.0 and p == 1.0)]
        if not b: continue
        mp, of = statistics.fmean(p for p, _ in b), statistics.fmean(y for _, y in b)
        rel += len(b) / n * (mp - of) ** 2; res += len(b) / n * (of - base) ** 2
        yes_roi = [(y - p) / p for p, y in b if p > 0.005]
        no_roi = [((1 - y) - (1 - p)) / (1 - p) for p, y in b if p < 0.995]
        table.append({"price_band": f"{lo:.1f}-{hi:.1f}", "markets": len(b), "mean_price_pct": round(mp * 100, 2),
                      "resolved_yes_pct": round(of * 100, 2),
                      "buy_yes_roi_pct": round(statistics.fmean(yes_roi) * 100, 2) if yes_roi else None, "buy_yes_roi_ci95_pct": ci95(yes_roi),
                      "buy_no_roi_pct": round(statistics.fmean(no_roi) * 100, 2) if no_roi else None, "buy_no_roi_ci95_pct": ci95(no_roi)})
    return {"markets": n, "base_rate_yes_pct": round(base * 100, 2), "brier": round(brier, 4), "log_loss": round(logloss, 4),
            "murphy": {"reliability": round(rel, 4), "resolution": round(res, 4), "uncertainty": round(base * (1 - base), 4)},
            "calibration": table}

def main():
    rows = load()
    single = [r for r in rows if r["neg_risk"] in ("False", "false", "0", "")]
    res = {"source": "Polymarket Gamma and CLOB public APIs", "markets_total": len(rows),
           "close_dates": [min(r["close_date"] for r in rows), max(r["close_date"] for r in rows)],
           "all_markets": {"h24": analyse(rows, "yes_price_24h"), "d7": analyse(rows, "yes_price_7d")},
           "excluding_neg_risk": {"h24": analyse(single, "yes_price_24h"), "d7": analyse(single, "yes_price_7d")}}
    with open(os.path.join(HERE, "polymarket_calibration_results.json"), "w", encoding="utf-8") as f:
        json.dump(res, f, indent=2)
    print(json.dumps(res, indent=2))

if __name__ == "__main__":
    main()