""" replication_grid.py - the same 14 methods, re-run on three sports. Most tutorials on this site teach a method on one dataset. This script re-runs the core estimate of 14 of those methods on the three game-results files that ship with the site, so you can see which findings travel between sports and which do not: NBA nba_home_results.csv 2023-24 regular season, 1,230 games (the In-Season Tournament final is dropped, as official standings do) NFL nfl_games_lines.csv 1999-2025, 7,276 played games (only the teams, scores, season, game type and the Home/Neutral venue flag are read) EPL epl_results_2025_26.csv Premier League 2025-26, 380 matches It needs only numpy and pandas, makes no network request, and is deterministic: every random step uses numpy's default_rng with SEED below, so two runs print the same table. Run: python replication_grid.py # prints the table python replication_grid.py --out grid.csv # also writes it as CSV Methods (the tutorial each one comes from is named in the table's `tutorial` column): 1 home_points_share share of available points the home side took (draw/tie = half) 2 bootstrap mean home margin with a 95% percentile bootstrap interval 3 permutation sign-flip permutation p-value for "mean home margin = 0" 4 bayes Beta(1,1) prior + decisive home results -> posterior mean, 95% interval 5 power decisive games needed to detect that home edge (alpha .05, power .80) 6 regression_to_mean first-half vs second-half margin per game, team by team 7 luck_share share of the spread in team win% that coin-flip noise alone explains 8 gini Gini coefficient of team win% (wins + half of draws/ties) 9 pythagorean the exponent that best turns points for/against into win% 10 elo_k Elo K (and home edge) picked by log loss on a training block, scored on a later held-out block 11 ewma_half_life the half-life of an exponentially weighted form rating that best predicts the next game's margin (each half-life with its own fitted scale), picked and scored the same way 12 auc held-out AUC of the tuned Elo's pre-game home-win probability 13 outliers share of game margins outside Tukey's IQR fences / modified z > 3.5 14 rest home margin when the home side had more rest minus when it had less Training and held-out blocks: NBA and EPL are single seasons, so the first half of the games (by date) trains and the second half is held out; the NFL uses 1999 as a burn-in, trains on 2000-2012 and holds out 2013-2025. """ from __future__ import annotations import argparse import math import os import numpy as np import pandas as pd HERE = os.path.dirname(os.path.abspath(__file__)) SEED = 20261029 K_GRID = [4, 8, 12, 16, 20, 24, 28, 32, 40, 48, 60, 80] HFA_GRID = [0, 25, 50, 75, 100, 125, 150, 175, 200] HALF_LIVES = [1, 2, 3, 4, 5, 6, 8, 10, 12, 15, 20, 25, 30, 40, 50, 60, 80, 100] TUTORIAL_HALF_LIFE = math.log(0.5) / math.log(1 - 2 / 11) # the EWMA tutorial's span=10 CARRY = 1 / 3 # NFL Elo: each rating moves a third of the way back to 1500 between seasons MIN_PRIOR = 3 # EWMA forecasts are scored once both teams have played 3 games FRANCHISE = {"STL": "LA", "SD": "LAC", "OAK": "LV"} # NFL relocations: one rating per franchise # ----------------------------------------------------------------------------- data def load_nba(): d = pd.read_csv(os.path.join(HERE, "nba_home_results.csv")) d = d[d["date"] != "2023-12-09"] # In-Season Tournament final d = d.rename(columns={"home_team": "home", "away_team": "away", "home_pts": "hs", "away_pts": "as_"}) d["season"], d["regular"], d["neutral"] = 2024, True, False return _finish(d) def load_nfl(): d = pd.read_csv(os.path.join(HERE, "nfl_games_lines.csv"), usecols=["season", "game_type", "gameday", "away_team", "home_team", "away_score", "home_score", "location"]) d = d.dropna(subset=["home_score", "away_score"]) # unplayed 2026 schedule d = d.rename(columns={"gameday": "date", "home_team": "home", "away_team": "away", "home_score": "hs", "away_score": "as_"}) d["home"], d["away"] = d["home"].replace(FRANCHISE), d["away"].replace(FRANCHISE) d["regular"] = d["game_type"] == "REG" d["neutral"] = d["location"] == "Neutral" return _finish(d) def load_epl(): d = pd.read_csv(os.path.join(HERE, "epl_results_2025_26.csv")) d = d.rename(columns={"home_goals": "hs", "away_goals": "as_"}) d["season"], d["regular"], d["neutral"] = 2026, True, False return _finish(d) def _finish(d): d = d.copy() d["date"] = pd.to_datetime(d["date"]) d = d.sort_values("date", kind="mergesort").reset_index(drop=True) # stable: file order within a day d["hs"], d["as_"] = d["hs"].astype(int), d["as_"].astype(int) d["margin"] = d["hs"] - d["as_"] d["score"] = np.where(d["margin"] > 0, 1.0, np.where(d["margin"] < 0, 0.0, 0.5)) return d def split(sport, d): """(train mask, test mask) over the rows of d, in date order.""" if sport == "NFL": return ((d["season"] >= 2000) & (d["season"] <= 2012)).values, (d["season"] >= 2013).values half = len(d) // 2 idx = np.arange(len(d)) return idx < half, idx >= half # ----------------------------------------------------------------------------- models def run_elo(d, K, hfa, carry=CARRY): """Pre-game home-win probability for every game, in date order (online Elo).""" r, season, out = {}, None, np.empty(len(d)) for i, (s, h, a, sc, neu) in enumerate(zip(d["season"], d["home"], d["away"], d["score"], d["neutral"])): if s != season: r = {t: 1500 + (v - 1500) * (1 - carry) for t, v in r.items()} season = s rh, ra = r.get(h, 1500.0), r.get(a, 1500.0) p = 1 / (1 + 10 ** ((ra - rh - (0 if neu else hfa)) / 400)) out[i] = p r[h], r[a] = rh + K * (sc - p), ra - K * (sc - p) return out def log_loss(p, s): p = np.clip(p, 1e-12, 1 - 1e-12) return float(-np.mean(s * np.log(p) + (1 - s) * np.log(1 - p))) def tune_elo(d, train, test): s = d["score"].values _, K, hfa = min((log_loss(run_elo(d, k, h)[train], s[train]), k, h) for k in K_GRID for h in HFA_GRID) curve = {k: log_loss(run_elo(d, k, hfa)[test], s[test]) for k in K_GRID} base_p = s[train].mean() return { "K": K, "hfa": hfa, "test_ll": curve[K], "test_ll_tutorial": log_loss(run_elo(d, 20, 100)[test], s[test]), "test_ll_base": log_loss(np.full(test.sum(), base_p), s[test]), "test_best_K": min(curve, key=curve.get), "curve": curve, "p": run_elo(d, K, hfa), } def ewma_forecast(d, half_life, mu): """Forecast home margin = venue edge + home form - away form; form = EWMA of each team's own margins with the venue edge removed. Returns (forecast, both-teams-ready).""" alpha = 1 - 0.5 ** (1 / half_life) form, n = {}, {} pred, ready = np.empty(len(d)), np.zeros(len(d), bool) for i, (h, a, m, neu) in enumerate(zip(d["home"], d["away"], d["margin"], d["neutral"])): edge = 0.0 if neu else mu pred[i] = edge + form.get(h, 0.0) - form.get(a, 0.0) ready[i] = n.get(h, 0) >= MIN_PRIOR and n.get(a, 0) >= MIN_PRIOR for t, x in ((h, m - edge), (a, -(m - edge))): form[t] = x if t not in form else alpha * x + (1 - alpha) * form[t] n[t] = n.get(t, 0) + 1 return pred, ready def rmse(e): return float(np.sqrt(np.mean(e ** 2))) def tune_ewma(d, train, test): """Each half-life gets a scale factor (least squares on the training block) because a raw form gap over-reacts: part of any team's recent margin is noise (see regression_to_mean). The half-life with the lowest training RMSE is picked.""" m = d["margin"].values.astype(float) mu = float(m[train & ~d["neutral"].values].mean()) y = m - np.where(d["neutral"].values, 0.0, mu) # margin left after the venue edge fits = {} for hl in HALF_LIVES + [TUTORIAL_HALF_LIFE]: pred, ready = ewma_forecast(d, hl, mu) x = pred - (m - y) # home form minus away form a, b = train & ready, test & ready beta = float((x[a] * y[a]).sum() / (x[a] ** 2).sum()) fits[hl] = (rmse((y - beta * x)[a]), rmse((y - beta * x)[b]), beta, b) best = min(HALF_LIVES, key=lambda h: fits[h][0]) b = fits[best][3] return { "mu": mu, "half_life": best, "beta": fits[best][2], "test_rmse": fits[best][1], "test_rmse_tutorial": fits[TUTORIAL_HALF_LIFE][1], "test_rmse_base": rmse(y[b]), "test_best_half_life": min(HALF_LIVES, key=lambda h: fits[h][1]), "n_test": int(b.sum()), "curve": {h: fits[h][1] for h in HALF_LIVES}, } # ----------------------------------------------------------------------------- the other methods def team_rows(d): """One row per team per game, from that team's side (regular-season games only).""" g = d[d["regular"]] home = pd.DataFrame({"season": g["season"], "date": g["date"], "team": g["home"], "pf": g["hs"], "pa": g["as_"], "res": g["score"]}) away = pd.DataFrame({"season": g["season"], "date": g["date"], "team": g["away"], "pf": g["as_"], "pa": g["hs"], "res": 1 - g["score"]}) t = pd.concat([home, away]).sort_values(["season", "team", "date"], kind="mergesort") return t.reset_index(drop=True) def gini(x): x = np.sort(np.asarray(x, float)) n = len(x) return float((2 * np.arange(1, n + 1) - n - 1).dot(x) / (n * x.sum())) def auc(score, y): r = pd.Series(score).rank().values pos = y == 1 n1, n0 = pos.sum(), (~pos).sum() return float((r[pos].sum() - n1 * (n1 + 1) / 2) / (n1 * n0)) def methods(sport, d, rng): rows = [] def add(n, method, tutorial, metric, value, fmt): rows.append({"sport": sport, "n": n, "method": method, "tutorial": tutorial, "metric": metric, "value": fmt % value}) home = d[d["regular"] & ~d["neutral"]] # home-edge methods: regular season, real home games m = home["margin"].values.astype(float) s = home["score"].values n = len(home) add(n, "home_points_share", "quantifying-home-advantage-across-five-leagues", "home share of points (draw = half)", s.mean(), "%.4f") boots = rng.choice(m, size=(10_000, n), replace=True).mean(axis=1) lo, hi = np.percentile(boots, [2.5, 97.5]) add(n, "bootstrap", "bootstrap-confidence-intervals", "mean home margin", m.mean(), "%.3f") add(n, "bootstrap", "bootstrap-confidence-intervals", "95% CI low", lo, "%.3f") add(n, "bootstrap", "bootstrap-confidence-intervals", "95% CI high", hi, "%.3f") flips = rng.choice([-1.0, 1.0], size=(20_000, n)) null = np.abs((flips * m).mean(axis=1)) p = (np.sum(null >= abs(m.mean())) + 1) / (20_000 + 1) add(n, "permutation", "permutation-test-significance", "sign-flip p-value (20,000 flips)", p, "%.5f") w, l = int((s == 1).sum()), int((s == 0).sum()) draws = rng.beta(1 + w, 1 + l, size=200_000) add(w + l, "bayes", "bayesian-updating-beta-binomial", "posterior mean P(home win | decisive)", (1 + w) / (2 + w + l), "%.4f") add(w + l, "bayes", "bayesian-updating-beta-binomial", "95% interval low", np.percentile(draws, 2.5), "%.4f") add(w + l, "bayes", "bayesian-updating-beta-binomial", "95% interval high", np.percentile(draws, 97.5), "%.4f") p1 = w / (w + l) za, zb = 1.959963984540054, 0.8416212335729143 need = ((za * math.sqrt(0.25) + zb * math.sqrt(p1 * (1 - p1))) / (p1 - 0.5)) ** 2 add(w + l, "power", "statistical-power-sample-size", "decisive games to detect it", math.ceil(need), "%d") t = team_rows(d) t["k"] = t.groupby(["season", "team"]).cumcount() t["g"] = t.groupby(["season", "team"])["k"].transform("size") t["md"] = t["pf"] - t["pa"] first = t[t["k"] < t["g"] // 2].groupby(["season", "team"])["md"].mean() second = t[t["k"] >= t["g"] // 2].groupby(["season", "team"])["md"].mean() a, b = first.values, second.loc[first.index].values add(len(a), "regression_to_mean", "regression-to-the-mean", "r(first half, second half)", np.corrcoef(a, b)[0, 1], "%.3f") add(len(a), "regression_to_mean", "regression-to-the-mean", "slope of second half on first", np.polyfit(a, b, 1)[0], "%.3f") ts = t.groupby(["season", "team"]).agg(res=("res", "sum"), g=("res", "size"), pf=("pf", "sum"), pa=("pa", "sum")).reset_index() ts["wp"] = ts["res"] / ts["g"] obs, noise = [], [] for _, grp in ts.groupby("season"): games = d[d["regular"] & (d["season"] == grp["season"].iloc[0])] p_draw = (games["score"] == 0.5).mean() obs.append(grp["wp"].var(ddof=0)) noise.append((1 - p_draw) * 0.25 / grp["g"].mean()) add(len(ts), "luck_share", "team-luck-with-pdo", "share of win% variance from chance", np.mean(noise) / np.mean(obs), "%.3f") add(len(ts), "gini", "lorenz-curve-gini-coefficient", "Gini of team win% (mean over seasons)", np.mean([gini(grp["wp"]) for _, grp in ts.groupby("season")]), "%.3f") ks = np.round(np.arange(0.5, 20.0001, 0.01), 2) errs = [rmse(ts["wp"].values - ts["pf"] ** k / (ts["pf"] ** k + ts["pa"] ** k)) for k in ks] add(len(ts), "pythagorean", "pythagorean-wins-expected-standings", "best exponent", ks[int(np.argmin(errs))], "%.2f") train, test = split(sport, d) elo = tune_elo(d, train, test) nt = int(test.sum()) add(nt, "elo_k", "build-an-elo-rating-system", "K picked on training block", elo["K"], "%d") add(nt, "elo_k", "build-an-elo-rating-system", "home edge picked (Elo points)", elo["hfa"], "%d") add(nt, "elo_k", "build-an-elo-rating-system", "held-out log loss, picked K", elo["test_ll"], "%.4f") add(nt, "elo_k", "build-an-elo-rating-system", "held-out log loss, K=20 home=100", elo["test_ll_tutorial"], "%.4f") add(nt, "elo_k", "build-an-elo-rating-system", "held-out log loss, home-rate only", elo["test_ll_base"], "%.4f") add(nt, "elo_k", "build-an-elo-rating-system", "best K on held-out block (hindsight)", elo["test_best_K"], "%d") ew = tune_ewma(d, train, test) add(ew["n_test"], "ewma_half_life", "exponentially-weighted-moving-average", "half-life picked (games)", ew["half_life"], "%d") add(ew["n_test"], "ewma_half_life", "exponentially-weighted-moving-average", "scale on the form gap (fitted)", ew["beta"], "%.3f") add(ew["n_test"], "ewma_half_life", "exponentially-weighted-moving-average", "held-out RMSE, picked half-life", ew["test_rmse"], "%.3f") add(ew["n_test"], "ewma_half_life", "exponentially-weighted-moving-average", "held-out RMSE, span 10 (half-life 3.45), own scale", ew["test_rmse_tutorial"], "%.3f") add(ew["n_test"], "ewma_half_life", "exponentially-weighted-moving-average", "held-out RMSE, home edge only", ew["test_rmse_base"], "%.3f") add(ew["n_test"], "ewma_half_life", "exponentially-weighted-moving-average", "best half-life on held-out block (hindsight)", ew["test_best_half_life"], "%d") dec = test & (d["score"].values != 0.5) add(int(dec.sum()), "auc", "roc-curve-and-auc-from-scratch", "held-out AUC, tuned Elo", auc(elo["p"][dec], d["score"].values[dec]), "%.3f") q1, q3 = np.percentile(m, [25, 75]) iqr = q3 - q1 add(n, "outliers", "outlier-detection-iqr-mad", "share outside 1.5 IQR fences", np.mean((m < q1 - 1.5 * iqr) | (m > q3 + 1.5 * iqr)), "%.4f") med = np.median(m) mad = np.median(np.abs(m - med)) add(n, "outliers", "outlier-detection-iqr-mad", "share with modified z > 3.5", np.mean(np.abs(0.6745 * (m - med) / mad) > 3.5), "%.4f") t["rest"] = t.groupby(["season", "team"])["date"].diff().dt.days rest = t.set_index(["season", "date", "team"])["rest"] hr = rest.reindex(list(zip(home["season"], home["date"], home["home"]))).values ar = rest.reindex(list(zip(home["season"], home["date"], home["away"]))).values ok = ~np.isnan(hr) & ~np.isnan(ar) more, less = ok & (hr > ar), ok & (hr < ar) add(int(more.sum() + less.sum()), "rest", "rest-days-and-back-to-backs", "home margin, home more rested minus less rested", m[more].mean() - m[less].mean(), "%.3f") return rows def run(): rng = np.random.default_rng(SEED) rows = [] for sport, loader in (("NBA", load_nba), ("NFL", load_nfl), ("EPL", load_epl)): rows += methods(sport, loader(), rng) return pd.DataFrame(rows) def main(): ap = argparse.ArgumentParser(description=__doc__.split("\n")[1]) ap.add_argument("--out", help="also write the table to this CSV path") args = ap.parse_args() table = run() wide = table.pivot_table(index=["method", "metric"], columns="sport", values="value", aggfunc="first", sort=False)[["NBA", "NFL", "EPL"]] with pd.option_context("display.width", 140, "display.max_rows", 200, "display.max_colwidth", 52): print(wide.to_string()) if args.out: table.to_csv(args.out, index=False, lineterminator="\n") print("\nwrote", args.out) if __name__ == "__main__": main()