""" Tutorial 98 - Regression discontinuity: what a close win is worth. Every Monday a team that won by one point is filed with the winners and a team that lost by one point with the losers, as if the W carried information the margin does not. Regression discontinuity (RD) is the method built for exactly that question. When a rule hands out a treatment at a cutoff - you win if your final margin is above zero - the cases just above and just below the cutoff are as alike as the data allows, so a jump in a later outcome right at the cutoff is the treatment's effect. The design: every regular-season NFL game since 1999, seen from both sides; the team's own final margin is the running variable, zero is the cutoff, and the outcome is the team's NEXT regular-season game - did it win, what line did the market post, did it beat that line. A straight line is fitted on each side of zero with triangular weights that fade with distance from the cutoff, and the jump is the gap between the two lines at zero. RD has a check that matters more than its estimate: nothing decided BEFORE the game is allowed to jump at the cutoff. Here something does. Narrow winners went into their games as slightly bigger favourites than narrow losers, so the closest NFL games are not quite coin flips. Written by hand: the running variable, triangular kernel weights, a weighted least-squares line on each side of the cutoff, the jump, a bootstrap that resamples whole games, a bandwidth sweep, covariate balance checks and placebo cutoffs. No scipy, no statsmodels, no rdrobust. Every number the tutorial page quotes is asserted below. Run: python downloads/98_regression_discontinuity_close_games.py Data: nflverse games table (github.com/nflverse/nfldata), June 2026 snapshot, bundled as nfl_games_lines.csv, plus the overtime flag from the same snapshot, bundled as nfl_games_overtime.csv. """ import math import os from fractions import Fraction import matplotlib.pyplot as plt import numpy as np import pandas as pd from matplotlib.lines import Line2D from matplotlib.patches import Patch import sdt_common as sdt sdt.init("regression-discontinuity-close-games") HERE = os.path.dirname(os.path.abspath(__file__)) CSV = os.path.join(HERE, "nfl_games_lines.csv") OT_CSV = os.path.join(HERE, "nfl_games_overtime.csv") SEED = 20261004 REPS = 2000 # bootstrap resamples of whole games H = 9 # bandwidth: margins 1-8 (every one-score game) get weight SWEEP = (4, 6, 9, 12, 16, 22) # bandwidths for the sensitivity table PLACEBOS = np.arange(9.5, 30.0, 1.0) # fake cutoffs where nothing changes hands def near(a, b, tol=5e-4): return abs(a - b) < tol # --- the method, written out in full ----------------------------------------------------- def wls_line(x, y, w): """Weighted least-squares straight line y = a + b*x, from five weighted sums. These are the normal equations for a line, solved directly: b = (S0*Sxy - Sx*Sy) / (S0*Sxx - Sx**2) a = (Sy - b*Sx) / S0 with S0 = sum(w), Sx = sum(w*x), Sxx = sum(w*x*x), Sy = sum(w*y), Sxy = sum(w*x*y). """ s0, sx, sxx = w.sum(), (w * x).sum(), (w * x * x).sum() sy, sxy = (w * y).sum(), (w * x * y).sum() b = (s0 * sxy - sx * sy) / (s0 * sxx - sx * sx) return (sy - b * sx) / s0, b def side_line(x, y, h, right, freq=None): """The local line on one side of the cutoff: |x| < h, triangular weights 1 - |x|/h. freq is an optional per-row count (the bootstrap uses it); rows with a missing outcome are skipped. """ keep = ((x > 0) if right else (x < 0)) & (np.abs(x) < h) & ~np.isnan(y) w = 1 - np.abs(x[keep]) / h if freq is not None: w = w * freq[keep] return wls_line(x[keep], y[keep], w) def rd_jump(x, y, h, freq=None): """Right-hand line at zero minus left-hand line at zero: the discontinuity.""" return side_line(x, y, h, True, freq)[0] - side_line(x, y, h, False, freq)[0] class Prepared: """One jump to bootstrap, with each side's rows and kernel weights cut out once.""" def __init__(self, x, y, code, h): self.sides = [] for right in (True, False): keep = ((x > 0) if right else (x < 0)) & (np.abs(x) < h) & ~np.isnan(y) self.sides.append((x[keep], y[keep], 1 - np.abs(x[keep]) / h, code[keep])) def jump(self, counts): (xr, yr, kr, cr), (xl, yl, kl, cl) = self.sides return wls_line(xr, yr, kr * counts[cr])[0] - wls_line(xl, yl, kl * counts[cl])[0] def game_bootstrap(prepared, n_games, reps, seed): """Resample whole games with replacement and recompute every jump each time. Drawing n_games games with replacement is the same as giving every game a count (0, 1, 2, ...) that sums to n_games and using that count as a row weight - so the two rows of one game, its winner and its loser, always travel together. """ rng = np.random.default_rng(seed) draws = {name: np.empty(reps) for name in prepared} for r in range(reps): counts = np.bincount(rng.integers(0, n_games, n_games), minlength=n_games).astype(float) for name, p in prepared.items(): draws[name][r] = p.jump(counts) return draws def interval(d): lo, hi = np.percentile(d, [2.5, 97.5]) return float(lo), float(hi) # --- 1. the design: running variable, cutoff, outcome ------------------------------------ games = pd.read_csv(CSV) ot = pd.read_csv(OT_CSV) games = games.merge(ot, on="game_id", how="left", validate="one_to_one") reg = games[(games.game_type == "REG") & games.result.notna()].copy() sides = [] for me, sign, home_venue in (("home", 1, 1), ("away", -1, -1)): sides.append(pd.DataFrame({ "game_id": reg.game_id, "season": reg.season, "gameday": reg.gameday, "team": reg[f"{me}_team"], "margin": sign * reg.result, # the running variable: own final margin "line": sign * reg.spread_line, # points this team was favoured by "venue": np.where(reg.location == "Home", home_venue, 0), "overtime": reg.overtime})) tg = (pd.concat(sides, ignore_index=True) .sort_values(["team", "season", "gameday"], kind="stable").reset_index(drop=True)) by = tg.groupby(["team", "season"], sort=False) tg["next_margin"] = by.margin.shift(-1) # the same team's next regular-season game tg["next_line"] = by.line.shift(-1) tg["next_win"] = np.where(tg.next_margin.isna(), np.nan, (tg.next_margin > 0) + 0.5 * (tg.next_margin == 0)) tg["next_cover"] = tg.next_margin - tg.next_line # margin against the posted line tg["pts"] = (tg.margin > 0) + 0.5 * (tg.margin == 0) played_before = by.cumcount() tg["winpct_before"] = (by.pts.cumsum() - tg.pts) / played_before # 0/0 -> NaN in week one tg["prev_margin"] = by.margin.shift(1) dec = tg[tg.margin != 0].reset_index(drop=True) # every decided game, from both sides out = dec[dec.next_margin.notna()].reset_index(drop=True) code_of = {g: i for i, g in enumerate(pd.unique(dec.game_id))} dec_code = dec.game_id.map(code_of).to_numpy() out_code = out.game_id.map(code_of).to_numpy() N_GAMES = len(code_of) counts_by_margin = dec.margin.value_counts() mirror = all(counts_by_margin.get(m, 0) == counts_by_margin.get(-m, 0) for m in counts_by_margin.index) abs_m = dec.margin.abs() share_3_7 = float(abs_m.isin([3, 7]).mean()) one_score = float((abs_m <= 8).mean()) rows_next = dec.assign(has=dec.next_margin.notna()).groupby("game_id").has.sum() finale_both, finale_one = int((rows_next == 0).sum()), int((rows_next == 1).sum()) with sdt.snippet("design"): print(f"regular-season games played {reg.season.min()}-{reg.season.max()}: {len(reg):,}") print(f" ties: {int((reg.result == 0).sum())}, every one after overtime; " f"overtime games in all: {int(reg.overtime.sum())}") print(f"rows, one per team per game: {len(tg):,}") print(f" decided games, both sides: {len(dec):,} rows from {N_GAMES:,} games") print(f" rows whose team plays again that regular season: {len(out):,}") print(f" rows with no next game: {len(dec) - len(out):,} " f"({finale_both} games where neither side plays again, {finale_one} where one side does not)") print("running variable = own final margin, cutoff = 0, treated = won") print(f" decided by 8 points or fewer: {one_score:.2%}; by exactly 3 or 7: {share_3_7:.2%}") print(f" rows at +m equal rows at -m for every m: {mirror}") assert len(games) == 7548 and len(reg) == 6967 and len(tg) == 13934 assert reg.season.min() == 1999 and reg.season.max() == 2025 assert int((reg.result == 0).sum()) == 15 and int(reg.overtime.sum()) == 415 assert (reg.loc[reg.result == 0, "overtime"] == 1).all() # a tie needs overtime assert reg.spread_line.notna().all() # every game has a line assert len(dec) == 13904 and N_GAMES == 6952 and len(dec) == 2 * N_GAMES assert len(out) == 13043 and len(dec) - len(out) == 861 assert finale_both == 429 and finale_one == 3 and 2 * finale_both + finale_one == 861 last_week = reg.groupby("season").week.transform("max") finale_ids = set(rows_next[rows_next == 0].index) assert set(reg.loc[reg.game_id.isin(finale_ids), "week"] - last_week[reg.game_id.isin(finale_ids)]) == {0} one_sided = reg[reg.game_id.isin(rows_next[rows_next == 1].index)] assert sorted(one_sided.season) == [1999, 2000, 2001] and set(one_sided.week) == {16} assert mirror # the mirror is exact assert (dec.groupby("game_id").margin.sum() == 0).all() # +m and -m, game by game assert (dec.groupby("game_id").line.sum() == 0).all() # and the line too assert set(dec.venue.unique()) == {-1, 0, 1} assert near(one_score, 0.5075, 5e-5) and near(share_3_7, 0.2417, 5e-5) assert int(counts_by_margin.idxmax()) in (3, -3) # 3 is the most common margin # --- 2. look before fitting: the cells next to the cutoff ---------------------------------- cells = out[out.margin.abs() <= 3].groupby("margin").agg( rows=("next_win", "size"), next_win=("next_win", "mean"), next_line=("next_line", "mean")) pre = dec[dec.margin.abs() <= 3].groupby("margin").line.agg(["size", "mean"]) cells["pre_line"] = pre["mean"] with sdt.snippet("cells"): print("margin rows next-game win % next-game line | pre-game line (all rows)") for m, r in cells.iterrows(): print(f"{int(m):+5d} {int(r.rows):5d} {100 * r.next_win:14.2f} {r.next_line:+13.2f}" f" | {r.pre_line:+.4f} ({int(pre.loc[m, 'size'])})") gap = {k: 100 * (cells.loc[k, "next_win"] - cells.loc[-k, "next_win"]) for k in (1, 2, 3)} assert [int(cells.loc[m, "rows"]) for m in (-3, -2, -1, 1, 2, 3)] == [993, 263, 278, 279, 263, 994] assert [int(pre.loc[m, "size"]) for m in (1, 2, 3)] == [293, 287, 1046] assert near(100 * cells.loc[1, "next_win"], 52.15, 5e-3) and near(100 * cells.loc[-1, "next_win"], 44.24, 5e-3) assert near(100 * cells.loc[2, "next_win"], 50.38, 5e-3) and near(100 * cells.loc[-2, "next_win"], 44.49, 5e-3) assert near(100 * cells.loc[3, "next_win"], 50.00, 5e-3) and near(100 * cells.loc[-3, "next_win"], 50.15, 5e-3) assert near(gap[1], 7.91, 5e-3) and near(gap[2], 5.89, 5e-3) and near(gap[3], -0.15, 5e-3) assert near(pre.loc[1, "mean"], 0.7594, 5e-5) and near(pre.loc[2, "mean"], 0.9321, 5e-5) assert near(pre.loc[3, "mean"], 1.0315, 5e-5) assert all(pre.loc[m, "mean"] == -pre.loc[-m, "mean"] for m in (1, 2, 3)) # exact mirror assert cells.loc[3, "rows"] > cells.loc[1, "rows"] + cells.loc[2, "rows"] assert cells.loc[-3, "rows"] > cells.loc[-1, "rows"] + cells.loc[-2, "rows"] assert min(cells.loc[[-2, -1, 1, 2], "rows"]) == 263 and max(cells.loc[[-2, -1, 1, 2], "rows"]) == 279 assert f"{gap[1]:.2f}" == "7.91" and f"{gap[2]:.2f}" == "5.89" and f"{gap[3]:.2f}" == "-0.15" # --- 3. the local line, by hand ---------------------------------------------------------- # Worked example: the pre-game line on the winning side at bandwidth 4. With whole-number # margins a row-level weighted fit equals a fit through the three cell means, each # weighted by (rows in the cell) x (kernel weight), so the arithmetic fits on a page. EX_H = 4 ex = dec[(dec.margin > 0) & (dec.margin < EX_H)].groupby("margin").line.agg(["size", "sum"]) ex_n = {int(m): int(r["size"]) for m, r in ex.iterrows()} ex_sum = {int(m): Fraction(str(r["sum"])) for m, r in ex.iterrows()} # half-point lines ex_k = {m: Fraction(EX_H - m, EX_H) for m in ex_n} # 3/4, 2/4, 1/4 ex_W = {m: ex_n[m] * ex_k[m] for m in ex_n} ex_mean = {m: ex_sum[m] / ex_n[m] for m in ex_n} W = sum(ex_W.values()) xbar = sum(ex_W[m] * m for m in ex_n) / W ybar = sum(ex_W[m] * ex_mean[m] for m in ex_n) / W slope = (sum(ex_W[m] * (m - xbar) * (ex_mean[m] - ybar) for m in ex_n) / sum(ex_W[m] * (m - xbar) ** 2 for m in ex_n)) at_zero = ybar - slope * xbar x_dec, x_out = dec.margin.to_numpy(float), out.margin.to_numpy(float) y_line = dec.line.to_numpy(float) row_fit = rd_jump(x_dec, y_line, EX_H) OUTCOMES = [("next_win", "next-game win rate", 100), ("next_margin", "next-game margin", 1), ("next_line", "next-game line", 1), ("next_cover", "next-game cover margin", 1)] fit = {} for col, _label, scale in OUTCOMES: y = out[col].to_numpy(float) fit[col] = (scale * side_line(x_out, y, H, True)[0], scale * side_line(x_out, y, H, False)[0]) fit["line"] = (side_line(x_dec, y_line, H, True)[0], side_line(x_dec, y_line, H, False)[0]) with sdt.snippet("fit"): print(f"worked example: pre-game line, winning side, bandwidth {EX_H}") print("margin rows sum of lines mean line kernel rows x kernel") for m in sorted(ex_n): print(f"{m:6d} {ex_n[m]:4d} {float(ex_sum[m]):12.1f} {float(ex_mean[m]):9.5f}" f" {str(ex_k[m]):>6} {float(ex_W[m]):13.2f}") print(f"weighted mean margin {xbar} = {float(xbar):.5f}, " f"weighted mean line {ybar} = {float(ybar):.5f}") print(f"slope {float(slope):.5f} per point -> line at zero " f"{float(ybar):.5f} - {float(slope):.5f} x {float(xbar):.5f} = {float(at_zero):+.5f}") print(f"mirror: the losing side ends at {float(-at_zero):+.5f}, jump {float(2 * at_zero):+.5f}" f" (row-by-row fit: {row_fit:+.5f})") print() print(f"bandwidth {H}: margins 1-{H - 1} each side, weights {H - 1}/{H} down to 1/{H}") print(f"{'':24s}{'just above 0':>13s}{'just below 0':>13s}{'jump':>9s}") for col, label, _scale in OUTCOMES + [("line", "pre-game line", 1)]: a, b = fit[col] print(f"{label:24s}{a:13.3f}{b:13.3f}{a - b:+9.3f}") jump = {col: fit[col][0] - fit[col][1] for col in fit} # the worked example, exactly assert ex_n == {1: 293, 2: 287, 3: 1046} assert ex_k == {1: Fraction(3, 4), 2: Fraction(1, 2), 3: Fraction(1, 4)} assert ex_W == {1: Fraction(879, 4), 2: Fraction(287, 2), 3: Fraction(523, 2)} # 219.75 143.5 261.5 assert W == Fraction(2499, 4) # 624.75 assert near(float(ex_W[3] / W), 0.4186, 5e-5) # the 3-point cell: 42% of the weight assert ex_sum == {1: Fraction(445, 2), 2: Fraction(535, 2), 3: Fraction(1079)} # 222.5 267.5 1079 assert near(float(ex_mean[1]), 0.75939, 5e-6) and near(float(ex_mean[2]), 0.93206, 5e-6) assert near(float(ex_mean[3]), 1.03155, 5e-6) assert xbar == Fraction(5165, 2499) and ybar == Fraction(1521, 1666) assert near(float(xbar), 2.06683, 5e-6) and near(float(ybar), 0.91297, 5e-6) assert near(float(slope), 0.13535, 5e-6) and near(float(at_zero), 0.63322, 5e-6) assert 2 * at_zero == Fraction(3028511, 2391343) # the jump, exactly assert near(float(2 * at_zero), 1.26645, 5e-6) assert near(0.91297 - 0.13535 * 2.06683, 0.63322, 5e-6) # the rounded hand sum holds assert abs(row_fit - float(2 * at_zero)) < 1e-9 # cell fit == row-by-row fit assert abs(side_line(x_dec, y_line, EX_H, False)[0] + float(at_zero)) < 1e-9 # the mirror # the hand-written solver against numpy's own least squares (polyfit weights are sqrt(w)) keep = (x_dec > 0) & (x_dec < H) kw = 1 - x_dec[keep] / H b_np, a_np = np.polyfit(x_dec[keep], y_line[keep], 1, w=np.sqrt(kw)) a_me, b_me = wls_line(x_dec[keep], y_line[keep], kw) assert abs(a_np - a_me) < 1e-9 and abs(b_np - b_me) < 1e-9 # the bandwidth-9 table assert near(fit["next_win"][0], 48.955) and near(fit["next_win"][1], 46.055) assert near(jump["next_win"], 2.900) assert near(fit["next_margin"][0], -1.052) and near(fit["next_margin"][1], -0.849) assert near(jump["next_margin"], -0.203) assert near(fit["next_line"][0], -0.217) and near(fit["next_line"][1], -0.448) assert near(jump["next_line"], 0.232) assert near(fit["next_cover"][0], -0.836) and near(fit["next_cover"][1], -0.401) assert near(jump["next_cover"], -0.435) assert near(fit["line"][0], 0.604) and abs(fit["line"][0] + fit["line"][1]) < 1e-9 assert near(jump["line"], 1.209) assert abs(jump["next_margin"] - (jump["next_line"] + jump["next_cover"])) < 1e-9 # margin = line + cover assert f"{fit['line'][0]:.3f}" == "0.604" and f"{jump['line']:.3f}" == "1.209" assert f"{jump['next_win']:.2f}" == "2.90" and f"{jump['next_win']:.1f}" == "2.9" assert f"{jump['next_line']:.2f}" == "0.23" and f"{jump['next_cover']:.2f}" == "-0.43" assert f"{jump['line']:.2f}" == "1.21" # --- 4. how sure: a bootstrap over whole games, and a bandwidth sweep ---------------------- prepared = {} for col, _label, scale in OUTCOMES: prepared[col] = Prepared(x_out, scale * out[col].to_numpy(float), out_code, H) prepared["line"] = Prepared(x_dec, y_line, dec_code, H) for h in SWEEP: prepared[f"win_h{h}"] = Prepared(x_out, 100 * out.next_win.to_numpy(float), out_code, h) prepared[f"cover_h{h}"] = Prepared(x_out, out.next_cover.to_numpy(float), out_code, h) prepared[f"line_h{h}"] = Prepared(x_dec, y_line, dec_code, h) # balance extras for step 5 prepared["venue"] = Prepared(x_dec, dec.venue.to_numpy(float), dec_code, H) prepared["winpct_before"] = Prepared(x_dec, 100 * dec.winpct_before.to_numpy(float), dec_code, H) prepared["prev_margin"] = Prepared(x_dec, dec.prev_margin.to_numpy(float), dec_code, H) no_ot = (dec.overtime == 0).to_numpy() prepared["line_no_ot"] = Prepared(x_dec[no_ot], y_line[no_ot], dec_code[no_ot], H) early = (dec.season <= 2011).to_numpy() prepared["line_1999_2011"] = Prepared(x_dec[early], y_line[early], dec_code[early], H) prepared["line_2012_2025"] = Prepared(x_dec[~early], y_line[~early], dec_code[~early], H) win_side = x_dec > 0 for c in PLACEBOS: prepared[f"placebo_{c}"] = Prepared(x_dec[win_side] - c, y_line[win_side], dec_code[win_side], H) draws = game_bootstrap(prepared, N_GAMES, REPS, SEED) est = {name: p.jump(np.ones(N_GAMES)) for name, p in prepared.items()} ci = {name: interval(d) for name, d in draws.items()} with sdt.snippet("boot"): print(f"bandwidth {H}, {REPS:,} bootstrap resamples of {N_GAMES:,} games") print(f"{'':24s}{'jump':>9s} 95% interval") for col, label, _scale in OUTCOMES + [("line", "pre-game line", 1)]: lo, hi = ci[col] print(f"{label:24s}{est[col]:+9.3f} {lo:+.3f} to {hi:+.3f}") print() print("bandwidth next-game win rate (pts) next-game cover (pts) pre-game line (pts)") for h in SWEEP: a, c, b = f"win_h{h}", f"cover_h{h}", f"line_h{h}" print(f"{h:9d} {est[a]:+6.2f} ({ci[a][0]:+6.2f} to {ci[a][1]:+6.2f})" f" {est[c]:+5.2f} ({ci[c][0]:+5.2f} to {ci[c][1]:+5.2f})" f" {est[b]:+5.2f} ({ci[b][0]:+5.2f} to {ci[b][1]:+5.2f})") # the shortcut equals the bootstrap it stands in for: one literal resample, rows repeated _rng = np.random.default_rng(SEED) _counts = np.bincount(_rng.integers(0, N_GAMES, N_GAMES), minlength=N_GAMES) _rep = np.repeat(np.arange(len(dec)), _counts[dec_code]) assert abs(rd_jump(x_dec[_rep], y_line[_rep], H) - prepared["line"].jump(_counts.astype(float))) < 1e-9 assert abs(draws["line"][0] - prepared["line"].jump(_counts.astype(float))) < 1e-12 for name in ("next_win", "next_margin", "next_line", "next_cover", "line"): assert abs(est[name] - jump[name]) < 1e-9 # prepared == plain fit assert abs(est["win_h9"] - est["next_win"]) < 1e-9 and abs(est["line_h9"] - est["line"]) < 1e-9 for name, lo_hi in {"next_win": (-2.947, 8.987), "next_margin": (-1.851, 1.538), "next_line": (-0.429, 0.928), "next_cover": (-1.920, 1.115), "line": (0.265, 2.149)}.items(): assert near(ci[name][0], lo_hi[0]) and near(ci[name][1], lo_hi[1]) for name in ("next_win", "next_margin", "next_line", "next_cover"): assert ci[name][0] < 0 < ci[name][1] # every outcome: no jump assert ci["line"][0] > 0 # the covariate: a jump win_sweep = {4: (12.47, 0.46, 25.28), 6: (7.84, -1.61, 17.74), 9: (2.90, -2.95, 8.99), 12: (1.71, -3.22, 6.67), 16: (1.11, -3.00, 5.32), 22: (1.35, -2.20, 4.99)} for h, (e_, lo_, hi_) in win_sweep.items(): assert near(est[f"win_h{h}"], e_, 5e-3) assert near(ci[f"win_h{h}"][0], lo_, 5e-3) and near(ci[f"win_h{h}"][1], hi_, 5e-3) assert ci["win_h4"][0] > 0 # the narrow-window mirage assert all(ci[f"win_h{h}"][0] < 0 < ci[f"win_h{h}"][1] for h in SWEEP if h > 4) # gone at every wider window assert all(est[f"win_h{h}"] < est["win_h4"] for h in SWEEP if h > 4) assert all(ci[f"cover_h{h}"][0] < 0 < ci[f"cover_h{h}"][1] for h in SWEEP) # cover never jumps line_sweep = {4: 1.266, 6: 1.142, 9: 1.209, 12: 1.165, 16: 1.139, 22: 1.125} for h, e_ in line_sweep.items(): assert near(est[f"line_h{h}"], e_) assert all(1.1 < est[f"line_h{h}"] < 1.3 for h in SWEEP) # the same size everywhere assert all(ci[f"line_h{h}"][0] > 0 for h in SWEEP if h >= 9) assert all(ci[f"line_h{h}"][0] < 0 for h in SWEEP if h < 9) assert f"{ci['next_win'][0]:.1f}" == "-2.9" and f"{ci['next_win'][1]:.1f}" == "9.0" assert f"{ci['line'][0]:.2f}" == "0.26" and f"{ci['line'][1]:.2f}" == "2.15" assert f"{ci['line'][0]:.3f}" == "0.265" and f"{ci['line'][1]:.3f}" == "2.149" assert f"{ci['win_h4'][1] - ci['win_h4'][0]:.1f}" == "24.8" # interval width at h = 4 assert f"{ci['win_h16'][1] - ci['win_h16'][0]:.1f}" == "8.3" # and at h = 16 assert f"{min(est[f'line_h{h}'] for h in SWEEP):.2f}" == "1.12" assert f"{max(est[f'line_h{h}'] for h in SWEEP):.2f}" == "1.27" with np.errstate(invalid="ignore", divide="ignore"): assert np.isnan(rd_jump(x_dec, y_line, 2)) # h = 2: one margin a side # --- 5. the check that matters: was anything decided before kickoff? ------------------------ one = reg[reg.result.abs() == 1] fav_won = int(((np.sign(one.result) == np.sign(one.spread_line)) & (one.spread_line != 0)).sum()) dog_won = int(((np.sign(one.result) == -np.sign(one.spread_line)) & (one.spread_line != 0)).sum()) pickem = int((one.spread_line == 0).sum()) n_fav = fav_won + dog_won coin_tail = Fraction(sum(math.comb(n_fav, j) for j in range(fav_won, n_fav + 1)), 2 ** n_fav) plac = np.array([est[f"placebo_{c}"] for c in PLACEBOS]) plac_excl0 = sum(1 for c in PLACEBOS if ci[f"placebo_{c}"][0] > 0 or ci[f"placebo_{c}"][1] < 0) era_gap = est["line_1999_2011"] - est["line_2012_2025"] era_se = float(np.std(draws["line_1999_2011"] - draws["line_2012_2025"], ddof=1)) with sdt.snippet("balance"): print(f"pre-determined variable at the cutoff (bandwidth {H}) jump 95% interval") for name, label in (("line", "pre-game line, points"), ("winpct_before", "season win % before the game"), ("prev_margin", "margin in the previous game"), ("venue", "venue (+1 home, -1 away)")): print(f" {label:44s}{est[name]:+7.3f} {ci[name][0]:+.3f} to {ci[name][1]:+.3f}") print(f"games decided by one point: {len(one)}; favourite won {fav_won}, " f"underdog won {dog_won}, pick'em {pickem}") print(f" a fair coin gives the favourite {fav_won}+ of {n_fav} with probability " f"{float(coin_tail):.4f}") print() print("is it overtime, or one era?") for name, label in (("line_no_ot", "without the overtime games"), ("line_1999_2011", "1999-2011 only"), ("line_2012_2025", "2012-2025 only")): print(f" {label:44s}{est[name]:+7.3f} {ci[name][0]:+.3f} to {ci[name][1]:+.3f}") print(f" difference between the eras {era_gap:+.3f}, bootstrap sd {era_se:.3f} " f"({era_gap / era_se:.2f} sd)") print() print(f"placebo cutoffs at {PLACEBOS[0]:g}, {PLACEBOS[1]:g}, ... {PLACEBOS[-1]:g} " f"({len(PLACEBOS)} of them), same bandwidth:") print(f" mean jump {plac.mean():+.3f}, sd {plac.std(ddof=1):.3f}, " f"largest |jump| {np.abs(plac).max():.3f}; intervals excluding zero: {plac_excl0}") print(f" the real cutoff: {est['line']:+.3f}") assert near(est["winpct_before"], 3.527) and near(ci["winpct_before"][0], 0.265) assert near(ci["winpct_before"][1], 6.697) and ci["winpct_before"][0] > 0 assert int(dec.winpct_before.isna().sum()) == int((played_before.loc[tg.margin != 0] == 0).sum()) assert near(est["prev_margin"], 0.962) and ci["prev_margin"][0] < 0 < ci["prev_margin"][1] assert near(ci["prev_margin"][0], -0.740) and near(ci["prev_margin"][1], 2.661) assert near(est["venue"], 0.079) and ci["venue"][0] < 0 < ci["venue"][1] assert near(ci["venue"][0], -0.085) and near(ci["venue"][1], 0.238) assert len(one) == 293 and (fav_won, dog_won, pickem) == (159, 133, 1) assert near(fav_won / n_fav, 0.5445, 5e-5) assert near(float(coin_tail), 0.0717, 5e-5) # weak on its own assert near(est["line_no_ot"], 1.124) and ci["line_no_ot"][0] > 0 # not overtime assert near(ci["line_no_ot"][0], 0.125) and near(ci["line_no_ot"][1], 2.126) assert int(no_ot.sum()) == len(dec) - 2 * (415 - 15) # 400 decided OT games dropped assert near(est["line_1999_2011"], 2.033) and near(est["line_2012_2025"], 0.548) assert ci["line_1999_2011"][0] > 0 and ci["line_2012_2025"][0] < 0 < ci["line_2012_2025"][1] assert near(ci["line_1999_2011"][0], 0.473) and near(ci["line_1999_2011"][1], 3.567) assert near(ci["line_2012_2025"][0], -0.626) and near(ci["line_2012_2025"][1], 1.750) assert near(era_gap, 1.485) and near(era_se, 0.994) and era_gap / era_se < 1.96 # not established assert len(PLACEBOS) == 21 and PLACEBOS[0] == 9.5 and PLACEBOS[-1] == 29.5 assert near(plac.mean(), -0.096) and near(plac.std(ddof=1), 0.427) assert near(np.abs(plac).max(), 0.965) and plac_excl0 == 0 assert np.abs(plac).max() < est["line"] # the real cutoff beats all 21 assert f"{era_gap / era_se:.2f}" == "1.49" # --- 6. the exhibit --------------------------------------------------------------------------- BROWN, BLUE, GREY, INK, BAND = sdt.sport_color("football"), "#2C5E8A", "#8A8577", "#20242B", "#E9DCC3" win_cells = out[out.margin.abs() <= 21].groupby("margin").agg(n=("next_win", "size"), v=("next_win", "mean")) line_cells = dec[dec.margin.abs() <= 21].groupby("margin").agg(n=("line", "size"), v=("line", "mean")) inside = win_cells[(win_cells.index.to_series().abs() < H)] assert len(inside) == 16 assert f"{100 * inside.v.min():.2f}" == "44.24" and f"{100 * inside.v.max():.2f}" == "58.70" fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12.4, 5.0)) for ax, cells_, scale, xs, yv, lab in ( (ax1, win_cells, 100, x_out, 100 * out.next_win.to_numpy(float), "next-game win rate (%)"), (ax2, line_cells, 1, x_dec, y_line, "pre-game line (points favoured by)")): ax.axvspan(-(H - 0.5), H - 0.5, color=BAND, zorder=0, lw=0) ax.axvline(0, color=GREY, lw=1.0, ls="--", zorder=1) ax.scatter(cells_.index, scale * cells_.v, s=8 + cells_.n / 9, color=BROWN, alpha=0.75, linewidths=0, zorder=2) for right in (True, False): a, b = side_line(xs, yv, H, right) grid = np.linspace(0, H - 1, 50) * (1 if right else -1) ax.plot(grid, a + b * grid, color=BLUE, lw=2.0, zorder=3) ax.scatter([0], [a], s=46, facecolor="#FBF7EE", edgecolor=BLUE, linewidths=1.8, zorder=4) ax.set_xlim(-21.8, 21.8) ax.set_xlabel("this game's final margin (team's points minus opponent's)") ax.set_ylabel(lab) ax1.set_title("Next week: a gap inside the noise", fontsize=12) ax2.set_title("Before kickoff: a jump that should not be there", fontsize=12) lo, hi = ci["next_win"] ax1.text(0.03, 0.95, f"jump {est['next_win']:+.1f} pts\n95% interval {lo:+.1f} to {hi:+.1f}", transform=ax1.transAxes, ha="left", va="top", fontsize=9, color=INK) lo, hi = ci["line"] ax2.text(0.03, 0.95, f"jump {est['line']:+.2f} pts\n95% interval {lo:+.2f} to {hi:+.2f}", transform=ax2.transAxes, ha="left", va="top", fontsize=9, color=INK) handles = [Line2D([], [], ls="", marker="o", markersize=7, color=BROWN, alpha=0.75, label="mean at each margin (size = rows)"), Line2D([], [], color=BLUE, lw=2.0, label=f"local line, bandwidth {H}"), Line2D([], [], ls="", marker="o", markersize=7, markerfacecolor="#FBF7EE", markeredgecolor=BLUE, markeredgewidth=1.8, label="each line's value at zero"), Patch(facecolor=BAND, label="games inside the bandwidth")] fig.legend(handles=handles, loc="upper center", ncol=4, fontsize=8.6, frameon=False, bbox_to_anchor=(0.5, 1.0)) # one key for both panels, above them fig.tight_layout(rect=(0, 0, 1, 0.94)) sdt.save_fig(fig, "rd_close_wins", source="nflverse games table (github.com/nflverse/nfldata), regular seasons 1999-2025", asof="June 2026") print("\nall asserts passed")