"""Reproducible analysis of real English Premier League closing lines (Sports Betting Math). Input: epl-closing-lines-2019-2026.csv (built from football-data.co.uk E0 files, see README in this folder) Output: printed tables and epl_closing_lines_results.json Methods - Overround: sum of implied probabilities minus 1, per bookmaker set (Pinnacle close, market average close). - De-vig: multiplicative, additive, power and Shin, fitted per match on Pinnacle closing odds. - Calibration: Brier score and log loss of each de-vig method against full-time results. - Favourite-longshot bias: realised flat-stake ROI at closing odds by odds band, with 95% intervals. - Closing line value (CLV): bets at Pinnacle opening odds; CLV = opening odds x fair closing probability - 1. Run: python epl_closing_lines_analysis.py """ import csv, json, math, os, random, statistics HERE = os.path.dirname(os.path.abspath(__file__)) DATA = os.path.join(HERE, "epl-closing-lines-2019-2026.csv") OUTCOMES = ("H", "D", "A") def implied(odds): return [1 / o for o in odds] def devig_multiplicative(odds): p = implied(odds); s = sum(p) return [x / s for x in p] def devig_additive(odds): p = implied(odds); s = sum(p) return [x - (s - 1) / len(p) for x in p] def devig_power(odds): p = implied(odds); lo, hi = 1.0, 3.0 for _ in range(80): k = (lo + hi) / 2 if sum(x ** k for x in p) > 1: lo = k else: hi = k return [x ** ((lo + hi) / 2) for x in p] def shin_fit(odds): """Return (fair probabilities, z) for Shin's model.""" p = implied(odds); s = sum(p) def probs(z): return [(math.sqrt(z * z + 4 * (1 - z) * x * x / s) - z) / (2 * (1 - z)) for x in p] lo, hi = 0.0, 0.5 for _ in range(80): z = (lo + hi) / 2 if sum(probs(z)) > 1: lo = z else: hi = z z = (lo + hi) / 2 return probs(z), z def devig_shin(odds): return shin_fit(odds)[0] def ci95(values): m = statistics.fmean(values) se = statistics.pstdev(values) / math.sqrt(len(values)) if len(values) > 1 else 0.0 return m, m - 1.96 * se, m + 1.96 * se def load(): with open(DATA, encoding="utf-8") as f: rows = list(csv.DictReader(f)) out = [] for r in rows: try: close = [float(r["PinnacleClose" + k]) for k in OUTCOMES] avg = [float(r["AvgClose" + k]) for k in OUTCOMES] open_ = [float(r["PinnacleOpen" + k]) for k in OUTCOMES] except ValueError: continue if r["FTR"] not in OUTCOMES or min(close + avg + open_) <= 1: continue out.append({"season": r["Season"], "ftr": r["FTR"], "close": close, "avg": avg, "open": open_}) return out def main(): matches = load() res = {"matches": len(matches), "seasons": sorted({m["season"] for m in matches})} # 1. Overround by season by_season = {} for m in matches: d = by_season.setdefault(m["season"], {"pin": [], "avg": []}) d["pin"].append((sum(implied(m["close"])) - 1) * 100) d["avg"].append((sum(implied(m["avg"])) - 1) * 100) res["overround_by_season"] = {s: {"matches": len(d["pin"]), "pinnacle_close_pct": round(statistics.fmean(d["pin"]), 2), "market_average_close_pct": round(statistics.fmean(d["avg"]), 2)} for s, d in sorted(by_season.items())} res["overround_overall"] = {"pinnacle_close_pct": round(statistics.fmean(x for d in by_season.values() for x in d["pin"]), 2), "market_average_close_pct": round(statistics.fmean(x for d in by_season.values() for x in d["avg"]), 2)} # 2. Calibration of de-vig methods (Pinnacle close) methods = {"multiplicative": devig_multiplicative, "additive": devig_additive, "power": devig_power, "shin": devig_shin} cal = {} for name, fn in methods.items(): brier, logloss = [], [] for m in matches: p = fn(m["close"]) y = [1.0 if m["ftr"] == k else 0.0 for k in OUTCOMES] brier.append(sum((pi - yi) ** 2 for pi, yi in zip(p, y))) logloss.append(-math.log(max(1e-12, p[OUTCOMES.index(m["ftr"])]))) cal[name] = {"brier": round(statistics.fmean(brier), 5), "log_loss": round(statistics.fmean(logloss), 5)} res["devig_calibration_pinnacle_close"] = cal # 3. Shin minus multiplicative probability by odds band bands = [("favourite (<1.70)", 0, 1.70), ("mid (1.70-3.50)", 1.70, 3.50), ("longshot (>3.50)", 3.50, 1e9)] delta = {b[0]: [] for b in bands} for m in matches: mult, shin = devig_multiplicative(m["close"]), devig_shin(m["close"]) for i, o in enumerate(m["close"]): for name, lo, hi in bands: if lo <= o < hi: delta[name].append((shin[i] - mult[i]) * 100) res["shin_minus_multiplicative_pp"] = {k: {"selections": len(v), "mean_pp": round(statistics.fmean(v), 3)} for k, v in delta.items()} # 4. Favourite-longshot bias: realised ROI by odds band (flat 1-unit stakes on every selection) roi_bands = [(1.0, 1.5), (1.5, 2.0), (2.0, 3.0), (3.0, 5.0), (5.0, 10.0), (10.0, 1e9)] flb = {} for src in ("close", "avg"): table = [] for lo, hi in roi_bands: rets, hits, imp = [], 0, [] for m in matches: for i, o in enumerate(m[src]): if lo <= o < hi: won = m["ftr"] == OUTCOMES[i] rets.append(o - 1 if won else -1.0); hits += won imp.append(devig_multiplicative(m[src])[i]) if not rets: continue mean, low, high = ci95(rets) table.append({"band": f"{lo:.2f}-{hi:.2f}" if hi < 1e9 else f">{lo:.2f}", "selections": len(rets), "hit_rate_pct": round(hits / len(rets) * 100, 2), "fair_implied_pct": round(statistics.fmean(imp) * 100, 2), "roi_pct": round(mean * 100, 2), "roi_ci95_pct": [round(low * 100, 2), round(high * 100, 2)]}) flb["pinnacle_close" if src == "close" else "market_average_close"] = table res["favourite_longshot_bias_roi"] = flb # 5. Closing line value: bets at Pinnacle opening odds clv, ret = [], [] for m in matches: fair_close = devig_shin(m["close"]) for i in range(3): clv.append(m["open"][i] * fair_close[i] - 1) ret.append(m["open"][i] - 1 if m["ftr"] == OUTCOMES[i] else -1.0) mx, my = statistics.fmean(clv), statistics.fmean(ret) beta = sum((x - mx) * (y - my) for x, y in zip(clv, ret)) / sum((x - mx) ** 2 for x in clv) pos = [r for c, r in zip(clv, ret) if c > 0]; neg = [r for c, r in zip(clv, ret) if c <= 0] res["clv_pinnacle_open_vs_close"] = { "selections": len(clv), "mean_clv_pct": round(mx * 100, 3), "beta_return_on_clv": round(beta, 3), "positive_clv": {"selections": len(pos), "roi_pct": round(statistics.fmean(pos) * 100, 2), "roi_ci95_pct": [round(v * 100, 2) for v in ci95(pos)[1:]]}, "non_positive_clv": {"selections": len(neg), "roi_pct": round(statistics.fmean(neg) * 100, 2), "roi_ci95_pct": [round(v * 100, 2) for v in ci95(neg)[1:]]}, } # Bootstrap over matches (the three outcomes of one match stay together): ROI gap between # positive and non-positive CLV bets, and the slope of return on CLV. def gap_and_beta(rows): c = [r[0] for r in rows]; y = [r[1] for r in rows] cx, cy = statistics.fmean(c), statistics.fmean(y) b = sum((u - cx) * (v - cy) for u, v in zip(c, y)) / sum((u - cx) ** 2 for u in c) return statistics.fmean([v for u, v in rows if u > 0]) - statistics.fmean([v for u, v in rows if u <= 0]), b per_match = [list(zip(clv[k:k + 3], ret[k:k + 3])) for k in range(0, len(clv), 3)] rng = random.Random(20261016) gaps, betas = [], [] for _ in range(2000): sample = [row for _ in per_match for row in per_match[rng.randrange(len(per_match))]] g, b = gap_and_beta(sample); gaps.append(g); betas.append(b) gaps.sort(); betas.sort() res["clv_pinnacle_open_vs_close"]["roi_gap_positive_minus_non_positive_pp"] = round((statistics.fmean(pos) - statistics.fmean(neg)) * 100, 2) res["clv_pinnacle_open_vs_close"]["roi_gap_ci95_pp_bootstrap_by_match"] = [round(gaps[49] * 100, 2), round(gaps[1950] * 100, 2)] res["clv_pinnacle_open_vs_close"]["beta_ci95_bootstrap_by_match"] = [round(betas[49], 2), round(betas[1950], 2)] # 6. Shin insider parameter z fitted per match def pct(v, q): v = sorted(v); return v[min(len(v) - 1, int(q * len(v)))] res["shin_z"] = {} for src, key in (("close", "pinnacle_close"), ("avg", "market_average_close")): zs = [shin_fit(m[src])[1] for m in matches] res["shin_z"][key] = {"median": round(statistics.median(zs), 4), "p10": round(pct(zs, 0.10), 4), "p90": round(pct(zs, 0.90), 4), "max": round(max(zs), 4)} with open(os.path.join(HERE, "epl_closing_lines_results.json"), "w", encoding="utf-8") as f: json.dump(res, f, indent=2) print(json.dumps(res, indent=2)) if __name__ == "__main__": main()