""" Tutorial 96 - Rank correlation from scratch: Kendall's tau and Spearman's rho. Pearson's r asks whether two columns move together on a straight line. That is a strong assumption and an awkward one for standings, where the quantity everyone argues about is an ORDER: who finished above whom. Rank correlation drops the line and counts orderings instead. The test case: how much of an NFL season carries over to the next one. Every pair of teams in season t either kept its order in season t+1 (concordant) or flipped it (discordant), and the whole method is the difference between those two counts. Win totals tie constantly, and when they do the naive coefficient cannot reach 1 even against itself. Both versions are written out here so the gap is visible. Written by hand: the concordant/discordant pair count, Kendall's tau-a and tau-b, Spearman's rho on midranks, the textbook sum-of-d-squared shortcut (shown failing), the normal approximation, an EXACT null by enumerating every permutation at n = 6, and a permutation null where enumeration is impossible. No scipy, no statsmodels. Run: python downloads/96_rank_correlation_kendall_and_spearman.py Data: nflverse games table (github.com/nflverse/nfldata), June 2026 snapshot, bundled as nfl_games_lines.csv; NHL public API (api-web.nhle.com), retrieved August 2026, bundled as nhl_team_pdo_two_seasons.csv. """ import itertools import math import os import matplotlib.pyplot as plt import numpy as np import pandas as pd import sdt_common as sdt sdt.init("rank-correlation-kendall-and-spearman") HERE = os.path.dirname(os.path.abspath(__file__)) CSV = os.path.join(HERE, "nfl_games_lines.csv") NHL_CSV = os.path.join(HERE, "nhl_team_pdo_two_seasons.csv") # Three franchises changed city (and therefore code) inside this file. A move is not # a new team, so the codes are folded before anything is counted. FRANCHISE = {"STL": "LA", "SD": "LAC", "OAK": "LV"} FIRST, LAST = 1999, 2025 # --- the method, written out in full ------------------------------------------------- def pair_counts(x, y): """Classify every pair of observations. This is the whole of rank correlation. For each pair (i, j) compare the two variables' directions. Agree -> concordant, disagree -> discordant. A pair tied on x only, on y only, or on both is none of those, and the three tie buckets are kept separate because tau-b needs them apart. """ x, y = np.asarray(x, dtype=float), np.asarray(y, dtype=float) C = D = Tx = Ty = Txy = 0 n = len(x) for i in range(n): for j in range(i + 1, n): a, b = np.sign(x[i] - x[j]), np.sign(y[i] - y[j]) if a == 0 and b == 0: Txy += 1 elif a == 0: Tx += 1 elif b == 0: Ty += 1 elif a * b > 0: C += 1 else: D += 1 return C, D, Tx, Ty, Txy def tau_a(x, y): """Kendall's tau-a: (C - D) over every pair there is. Ignores ties, so with ties present it cannot reach 1 even when the two orders agree as far as they can.""" C, D, _tx, _ty, _txy = pair_counts(x, y) n = len(x) return (C - D) / (n * (n - 1) / 2) def tau_b(x, y): """Kendall's tau-b: the same numerator over a denominator that drops, from each variable separately, the pairs that variable cannot order (Kendall 1945).""" C, D, Tx, Ty, _txy = pair_counts(x, y) return (C - D) / math.sqrt((C + D + Tx) * (C + D + Ty)) def midranks(v): """Ranks, with tied values sharing the average of the ranks they span.""" v = np.asarray(v, dtype=float) order = v.argsort() r = np.empty(len(v), dtype=float) r[order] = np.arange(1, len(v) + 1) for u in np.unique(v): m = v == u if m.sum() > 1: r[m] = r[m].mean() return r def spearman_rho(x, y): """Spearman's rho, by its definition: Pearson's r computed on the midranks.""" return float(np.corrcoef(midranks(x), midranks(y))[0, 1]) def spearman_shortcut(x, y): """The 1 - 6*sum(d^2)/(n(n^2-1)) formula from every textbook. It is algebraically identical to Pearson-on-ranks ONLY when no value is tied. With ties it is not the same number, and this script measures the gap. """ n = len(x) d = midranks(x) - midranks(y) return 1 - 6 * float((d ** 2).sum()) / (n * (n * n - 1)) def tau_normal_z(x, y): """The large-sample z for tau, from the null variance of (C - D) with NO ties.""" C, D, _tx, _ty, _txy = pair_counts(x, y) n = len(x) return 3 * (C - D) / math.sqrt(n * (n - 1) * (2 * n + 5) / 2) def two_sided_p(z): """Normal tail, from the standard library: erfc does the whole job.""" return math.erfc(abs(z) / math.sqrt(2)) def tau_a_null_sd(n): """Standard deviation of tau-a under independence, assuming no ties.""" return math.sqrt(2 * (2 * n + 5) / (9 * n * (n - 1))) def exact_null_tau_a(x, y): """Every permutation of y against a fixed x: the null distribution, complete. Feasible while n! is small. At n = 6 that is 720 relabelings and the answer is exact - no approximation, no random number generator. """ y = np.asarray(y, dtype=float) return np.array([tau_a(x, y[list(p)]) for p in itertools.permutations(range(len(y)))]) def sign_matrix(v): """sign(v_i - v_j) for every pair at once - the vectorised half of the fast path.""" v = np.asarray(v, dtype=float) return np.sign(v[:, None] - v[None, :]) def tau_b_fast(sx, y): """tau-b from a precomputed sign matrix. Same answer as tau_b, many times quicker, which is what makes a 10,000-shuffle permutation null tolerable to run.""" sy = sign_matrix(y) num = float((sx * sy).sum()) / 2.0 # (C - D) nx = float((sx != 0).sum()) / 2.0 # pairs x can order ny = float((sy != 0).sum()) / 2.0 # pairs y can order return num / math.sqrt(nx * ny) # --- 1. build one row per team-season -------------------------------------------------- games = pd.read_csv(CSV) reg = games[(games.game_type == "REG") & games.result.notna()].copy() sides = [] for me, opp in (("home", "away"), ("away", "home")): s = reg[["season", f"{me}_team", f"{me}_score", f"{opp}_score"]].copy() s.columns = ["season", "team", "pf", "pa"] sides.append(s) tg = pd.concat(sides, ignore_index=True) tg["team"] = tg.team.replace(FRANCHISE) tg["w"] = (tg.pf > tg.pa) + 0.5 * (tg.pf == tg.pa) # a tie is half a win ts = (tg.groupby(["season", "team"]) .agg(gp=("w", "size"), wins=("w", "sum"), pf=("pf", "sum"), pa=("pa", "sum")) .reset_index()) ts["pdiff"] = ts.pf - ts.pa def season(y): return ts[ts.season == y].set_index("team") # how often does each column refuse to order a pair of teams? tied_wins = tied_pdiff = all_pairs = 0 for yr in range(FIRST, LAST + 1): a = season(yr) vw, vp = a.wins.to_numpy(float), a.pdiff.to_numpy(float) n_y = len(a) all_pairs += n_y * (n_y - 1) // 2 tied_wins += sum(1 for i in range(n_y) for j in range(i + 1, n_y) if vw[i] == vw[j]) tied_pdiff += sum(1 for i in range(n_y) for j in range(i + 1, n_y) if vp[i] == vp[j]) with sdt.snippet("build"): print(f"{len(games)} rows in the file -> {len(reg)} played regular-season games " f"-> {len(ts)} team-seasons, {ts.season.min()}-{ts.season.max()}") print(f" tie games kept as half a win: {int((reg.result == 0).sum())}") sizes = ts.groupby("season").size() print(f" teams per season: {sizes.min()} (1999-2001) to {sizes.max()} " f"(Houston arrives in {int(sizes[sizes == sizes.max()].index.min())})") gp = ts.groupby("season").gp.agg(["min", "max"]) odd = gp[gp["min"] != gp["max"]] print(f" games per team: 16 through 2020, 17 from 2021 - except " f"{', '.join(str(i) for i in odd.index)}, where " f"{', '.join(sorted(season(int(odd.index[0])).query('gp != 17').index))} played 16") print(f"\nwithin-season pairs of teams, {FIRST}-{LAST}: {all_pairs}") print(f" tied on wins: {tied_wins:>5} ({tied_wins / all_pairs:.2%})") print(f" tied on point differential: {tied_pdiff:>3} ({tied_pdiff / all_pairs:.2%})") print("ties are the whole reason this page has two versions of Kendall's tau.") assert len(games) == 7548 and len(reg) == 6967 and len(ts) == 861 assert int((reg.result == 0).sum()) == 15 assert ts.season.min() == FIRST and ts.season.max() == LAST assert set(ts.groupby("season").size()) == {31, 32} assert all(len(season(y)) == 31 for y in (1999, 2000, 2001)) assert all(len(season(y)) == 32 for y in range(2002, LAST + 1)) assert "HOU" not in season(2001).index and "HOU" in season(2002).index assert not {"STL", "SD", "OAK"} & set(ts.team) # relocations folded assert set(ts.query("season <= 2020").gp) == {16} assert set(ts.query("season >= 2023").gp) == {17} assert sorted(season(2022).query("gp != 17").index) == ["BUF", "CIN"] assert all_pairs == 13299 and tied_wins == 1039 and tied_pdiff == 40 assert tied_wins / all_pairs > 25 * (tied_pdiff / all_pairs) # wins tie 26x more often assert abs(tied_wins / all_pairs - 0.0781) < 5e-4 assert abs(tied_pdiff / all_pairs - 0.0030) < 5e-4 assert float(ts.wins.sum()) == float(len(reg)) # every game produced one win # --- 2. count the pairs by hand, on a group with no ties ------------------------------- EX_YEAR = 2005 ex_prev = season(EX_YEAR).pdiff.sort_values(ascending=False).head(6) EX_TEAMS = list(ex_prev.index) ex_x = ex_prev.to_numpy(float) ex_y = season(EX_YEAR + 1).loc[EX_TEAMS, "pdiff"].to_numpy(float) exC, exD, exTx, exTy, exTxy = pair_counts(ex_x, ex_y) ex_tau = tau_a(ex_x, ex_y) ex_null = exact_null_tau_a(ex_x, ex_y) ex_one = float((ex_null >= ex_tau - 1e-12).mean()) ex_two = float((np.abs(ex_null) >= abs(ex_tau) - 1e-12).mean()) with sdt.snippet("hand"): print(f"the six best point differentials of {EX_YEAR}, and the same teams in " f"{EX_YEAR + 1} (point differential almost never ties):\n") print(" team " + f"{EX_YEAR} " + f"{EX_YEAR + 1}") for t, a_v, b_v in zip(EX_TEAMS, ex_x, ex_y): print(f" {t:<5} {a_v:+6.0f} {b_v:+6.0f}") print(f"\nevery one of the {len(EX_TEAMS) * (len(EX_TEAMS) - 1) // 2} pairs, " f"classified:") for i in range(len(EX_TEAMS)): for j in range(i + 1, len(EX_TEAMS)): kind = ("concordant" if (ex_x[i] - ex_x[j]) * (ex_y[i] - ex_y[j]) > 0 else "discordant") print(f" {EX_TEAMS[i]:<4} vs {EX_TEAMS[j]:<4} " f"{EX_YEAR} {ex_x[i] - ex_x[j]:+5.0f} " f"{EX_YEAR + 1} {ex_y[i] - ex_y[j]:+5.0f} {kind}") print(f"\n concordant {exC}, discordant {exD}, tied {exTx + exTy + exTxy}") print(f" tau-a = ({exC} - {exD}) / {len(EX_TEAMS) * (len(EX_TEAMS) - 1) // 2} " f"= {ex_tau:.6f}") print(f"\nis a third of a tau worth anything at n = {len(EX_TEAMS)}? Enumerate the " f"whole null:") print(f" {len(ex_null)} relabelings of the second column, every one of them") print(f" distinct values of tau-a: {len(np.unique(np.round(ex_null, 10)))}") print(f" null mean {ex_null.mean():.12f}, null sd {ex_null.std():.6f}") print(f" the no-ties formula says the sd should be " f"{tau_a_null_sd(len(EX_TEAMS)):.6f}") print(f" exact P(tau >= {ex_tau:.4f}) = {ex_one:.6f}, two-sided {ex_two:.6f}") print("six teams cannot tell you anything. The arithmetic is still the arithmetic.") assert EX_TEAMS == ["IND", "SEA", "DEN", "CAR", "PIT", "NYG"] assert [int(v) for v in ex_x] == [192, 181, 137, 132, 131, 108] assert [int(v) for v in ex_y] == [67, -6, 14, -35, 38, -7] assert (exTx, exTy, exTxy) == (0, 0, 0) # the point of picking this group assert (exC, exD) == (10, 5) and exC + exD == 15 assert abs(ex_tau - 1 / 3) < 1e-12 # (10 - 5) / 15, exactly assert len(ex_null) == 720 == math.factorial(6) assert abs(ex_null.mean()) < 1e-12 # the null is centred on zero assert len(np.unique(np.round(ex_null, 10))) == 16 # C - D is odd, from -15 to +15 # with no ties anywhere, the textbook null sd is not an approximation at all assert abs(ex_null.std() - tau_a_null_sd(6)) < 1e-6 assert abs(ex_null.std() - 0.35486) < 5e-6 assert abs(ex_one - 0.234722) < 1e-5 and abs(ex_two - 0.469444) < 1e-5 assert ex_two > 0.05 # nothing is established here assert abs(two_sided_p(tau_normal_z(ex_x, ex_y)) - 0.469444) > 0.01 # n=6 is too small assert tau_a(ex_x, ex_x) == 1.0 and tau_a(ex_x, -ex_x) == -1.0 assert abs(tau_b(ex_x, ex_y) - ex_tau) < 1e-12 # no ties -> tau-b IS tau-a # --- 3. what ties do, and the correction that answers them ----------------------------- CON_YEAR = 2023 playoffs = games[games.game_type.isin(["WC", "DIV", "CON", "SB"]) & games.result.notna()] con = playoffs[(playoffs.season == CON_YEAR) & (playoffs.game_type == "CON")] CON4 = sorted(set(con.home_team.replace(FRANCHISE)) | set(con.away_team.replace(FRANCHISE))) c4_x = season(CON_YEAR).loc[CON4, "wins"].to_numpy(float) c4_y = season(CON_YEAR + 1).loc[CON4, "wins"].to_numpy(float) c4C, c4D, c4Tx, c4Ty, c4Txy = pair_counts(c4_x, c4_y) c4_null = exact_null_tau_a(c4_x, c4_y) prev, nxt = season(CON_YEAR), season(CON_YEAR + 1) both = sorted(set(prev.index) & set(nxt.index)) X = prev.loc[both, "wins"].to_numpy(float) Y = nxt.loc[both, "wins"].to_numpy(float) C, D, Tx, Ty, Txy = pair_counts(X, Y) n32 = len(both) total_pairs = n32 * (n32 - 1) // 2 with sdt.snippet("ties"): print(f"the four teams that reached the {CON_YEAR} conference championship games, " f"by regular-season wins:\n") print(f" team {CON_YEAR} {CON_YEAR + 1}") for t, a_v, b_v in zip(CON4, c4_x, c4_y): print(f" {t:<5} {a_v:5.0f} {b_v:5.0f}") print(f"\n concordant {c4C}, discordant {c4D}, tied on {CON_YEAR} only {c4Tx}, " f"tied on {CON_YEAR + 1} only {c4Ty}, tied on both {c4Txy}") print(f" tau-a = ({c4C} - {c4D}) / 6 = {tau_a(c4_x, c4_y):+.6f}") print(f" tau-b = ({c4C} - {c4D}) / sqrt({c4C + c4D + c4Tx} x {c4C + c4D + c4Ty}) " f"= {tau_b(c4_x, c4_y):+.6f}") print(f"\nthe whole league, {CON_YEAR} wins against {CON_YEAR + 1} wins " f"({n32} teams, {total_pairs} pairs):") print(f" concordant {C}, discordant {D}") print(f" tied on {CON_YEAR} only {Tx}, on {CON_YEAR + 1} only {Ty}, on both {Txy} " f"-> {Tx + Ty + Txy} pairs no order at all") print(f" of the {C + D} pairs both seasons DO order, {C / (C + D):.2%} kept their order") print(f"\n tau-a {tau_a(X, Y):.6f} tau-b {tau_b(X, Y):.6f}") print() print(f" the clearest way to see why tau-b exists: score the {CON_YEAR} column " f"against ITSELF.") print(f" perfect agreement by construction, and tau-a still reads " f"{tau_a(X, X):.6f},") print(f" because {Tx + Txy} of the {total_pairs} pairs are tied inside " f"that one column") print(f" tau-b reads {tau_b(X, X):.6f}") print("tau-b divides by what each column can actually order.") print("That is the only difference, and it is why tau-b is the one to quote.") assert CON4 == ["BAL", "DET", "KC", "SF"] assert [int(v) for v in c4_x] == [13, 12, 11, 12] assert [int(v) for v in c4_y] == [12, 15, 15, 6] assert (c4C, c4D, c4Tx, c4Ty, c4Txy) == (1, 3, 1, 1, 0) assert c4C + c4D + c4Tx + c4Ty + c4Txy == 6 assert abs(tau_a(c4_x, c4_y) - (-1 / 3)) < 1e-12 # (1 - 3) / 6 assert abs(tau_b(c4_x, c4_y) - (-0.4)) < 1e-12 # -2 / sqrt(5 * 5) assert tau_b(c4_x, c4_y) < tau_a(c4_x, c4_y) # correction can go either way assert len(c4_null) == 24 == math.factorial(4) assert c4_null.std() < tau_a_null_sd(4) # ties shrink the real null assert n32 == 32 and total_pairs == 496 assert (C, D, Tx, Ty, Txy) == (269, 148, 42, 31, 6) assert C + D + Tx + Ty + Txy == total_pairs assert C + D == 417 and abs(C / (C + D) - 0.645084) < 1e-5 assert abs(tau_a(X, Y) - 0.243952) < 1e-5 assert abs(tau_b(X, Y) - 0.266833) < 1e-5 assert tau_b(X, Y) > tau_a(X, Y) # here the correction lifts it assert Tx + Txy == 48 # pairs the 2023 column cannot order assert abs(tau_a(X, X) - 448 / 496) < 1e-12 # perfect agreement is NOT 1 assert abs(tau_a(X, X) - 0.903226) < 1e-5 assert tau_b(X, X) == 1.0 # tau-b is exactly 1 assert tau_a(X, X) < tau_b(X, X) # --- 4. Spearman's rho, and the shortcut that does not survive ties -------------------- rho_exact = spearman_rho(X, Y) rho_fast = spearman_shortcut(X, Y) pearson = float(np.corrcoef(X, Y)[0, 1]) shortcut_rows = [] for yr in range(FIRST, LAST): a, b = season(yr), season(yr + 1) com = sorted(set(a.index) & set(b.index)) xv, yv = a.loc[com, "wins"].to_numpy(float), b.loc[com, "wins"].to_numpy(float) shortcut_rows.append((yr, spearman_rho(xv, yv), spearman_shortcut(xv, yv))) SC = pd.DataFrame(shortcut_rows, columns=["season", "rho", "shortcut"]) SC["gap"] = SC.shortcut - SC.rho with sdt.snippet("spearman"): print(f"the same {CON_YEAR} -> {CON_YEAR + 1} table, three coefficients:\n") print(f" Kendall tau-b (pairs that kept their order) {tau_b(X, Y):.6f}") print(f" Spearman rho (Pearson on the midranks) {rho_exact:.6f}") print(f" Pearson r (on the raw win totals) {pearson:.6f}") print("\nrho and tau-b are not rivals and not on the same scale: tau-b is a share of " "pairs,\nrho is a correlation of ranks. rho is reliably the larger number.") print(f"\nthe textbook shortcut 1 - 6*sum(d^2)/(n(n^2-1)) on the same column: " f"{rho_fast:.6f}") print(f" that is {rho_fast - rho_exact:+.6f} away from rho, because it is only " f"equal to Pearson-on-ranks when nothing is tied.") print(f"\nacross all {len(SC)} consecutive-season pairs:") print(" season rho shortcut gap") for r in SC.sort_values("gap", ascending=False).head(4).itertuples(): print(f" {r.season} {r.rho:.4f} {r.shortcut:.4f} {r.gap:+.4f}") print(f" ...") print(f" the shortcut is ABOVE the real value in {int((SC.gap > 0).sum())} of " f"{len(SC)} pairs, never below, by at most {SC.gap.max():.6f}") print("a one-directional error is a bias, not a rounding difference.") print("Compute rho as Pearson on the midranks.") assert abs(rho_exact - 0.349467) < 1e-5 assert abs(rho_fast - 0.359421) < 1e-5 assert abs(pearson - 0.345234) < 1e-5 assert rho_exact > tau_b(X, Y) # rho is the bigger scale assert rho_fast > rho_exact # the shortcut overstates assert int((SC.gap > 0).sum()) == 26 == len(SC) # in EVERY pair, not most assert (SC.gap > 0).all() and SC.gap.min() > 0.004 assert abs(SC.gap.max() - 0.015791) < 1e-5 assert int(SC.loc[SC.gap.idxmax(), "season"]) == 2006 # not 2016, which is the tau-b MIN # with no ties the shortcut and the definition agree to machine precision assert abs(spearman_shortcut(ex_x, ex_y) - spearman_rho(ex_x, ex_y)) < 1e-12 assert abs(spearman_rho(ex_x, ex_x) - 1.0) < 1e-12 # --- 5. every consecutive pair of seasons, and whether any of it is real ---------------- rows = [] for yr in range(FIRST, LAST): a, b = season(yr), season(yr + 1) com = sorted(set(a.index) & set(b.index)) xw = a.loc[com, "wins"].to_numpy(float) xp = a.loc[com, "pdiff"].to_numpy(float) yw = b.loc[com, "wins"].to_numpy(float) rows.append((yr, len(com), tau_a(xw, yw), tau_b(xw, yw), spearman_rho(xw, yw), float(np.corrcoef(xw, yw)[0, 1]), tau_b(xp, yw), tau_b(xw, xp))) PAIRS = pd.DataFrame(rows, columns=["season", "n", "tau_a", "tau_b", "rho", "pearson", "tb_pdiff", "tb_wins_pdiff"]) z32 = tau_normal_z(X, Y) p_normal = two_sided_p(z32) SHUFFLES = 10000 rng = np.random.default_rng(96) sx32 = sign_matrix(X) null_b = np.empty(SHUFFLES) null_a = np.empty(SHUFFLES) for i in range(SHUFFLES): perm = rng.permutation(Y) null_b[i] = tau_b_fast(sx32, perm) sp = sign_matrix(perm) null_a[i] = float((sx32 * sp).sum()) / 2.0 / total_pairs p_perm = float((np.abs(null_b) >= abs(tau_b(X, Y)) - 1e-12).mean()) with sdt.snippet("pairs"): print(f"every consecutive pair of seasons, {FIRST}-{LAST}:\n") print(" season n tau-a tau-b rho Pearson") for r in PAIRS.itertuples(): print(f" {r.season} {r.n:>3} {r.tau_a:.4f} {r.tau_b:.4f} {r.rho:.4f} " f"{r.pearson:.4f}") print(f"\n median tau-a {PAIRS.tau_a.median():.4f} tau-b " f"{PAIRS.tau_b.median():.4f} rho {PAIRS.rho.median():.4f} " f"Pearson {PAIRS.pearson.median():.4f}") print(f" tau-b is positive in {int((PAIRS.tau_b > 0).sum())} of {len(PAIRS)} pairs, " f"from {PAIRS.tau_b.min():.4f} ({int(PAIRS.loc[PAIRS.tau_b.idxmin(), 'season'])}) " f"to {PAIRS.tau_b.max():.4f} " f"({int(PAIRS.loc[PAIRS.tau_b.idxmax(), 'season'])})") print(f"\nis one season's tau-b distinguishable from zero? Take {CON_YEAR}:") print(f" normal approximation z = {z32:.4f}, two-sided p = {p_normal:.5f}") print(f" {SHUFFLES} shuffles two-sided p = {p_perm:.5f}") print(f" null sd of tau-a: {null_a.std():.6f} shuffled, " f"{tau_a_null_sd(n32):.6f} from the no-ties formula") print(f"\nat n = {n32} the tie correction to the null barely matters " f"({Tx + Ty + Txy} tied pairs of {total_pairs}).") print(f"at n = 4, where {c4Tx + c4Ty + c4Txy} of 6 pairs are tied, the same formula " f"claims sd {tau_a_null_sd(4):.4f} and the exact enumeration says " f"{c4_null.std():.4f}.") assert len(PAIRS) == 26 and PAIRS.season.min() == 1999 and PAIRS.season.max() == 2024 assert set(PAIRS.n) == {31, 32} assert int((PAIRS.tau_b > 0).sum()) == 26 # every single pair, positive assert (PAIRS.tau_b > PAIRS.tau_a).all() # ties always cost tau-a here assert (PAIRS.rho > PAIRS.tau_b).all() assert abs(PAIRS.tau_b.median() - 0.2049) < 5e-4 assert abs(PAIRS.tau_a.median() - 0.1895) < 5e-4 assert abs(PAIRS.rho.median() - 0.3046) < 5e-4 assert abs(PAIRS.pearson.median() - 0.2850) < 5e-4 assert abs(PAIRS.tau_b.min() - 0.0889) < 5e-4 assert int(PAIRS.loc[PAIRS.tau_b.idxmin(), "season"]) == 2016 assert abs(PAIRS.tau_b.max() - 0.4493) < 5e-4 assert int(PAIRS.loc[PAIRS.tau_b.idxmax(), "season"]) == 2013 assert abs(z32 - 1.962191) < 1e-5 and abs(p_normal - 0.04974) < 1e-4 assert 0.04 < p_perm < 0.055 # the two agree at this n assert abs(p_perm - p_normal) < 0.01 assert abs(null_a.std() - tau_a_null_sd(n32)) < 0.002 # ties are 16% of pairs, so close assert abs(null_a.mean()) < 0.01 assert c4_null.std() < 0.9 * tau_a_null_sd(4) # at n=4 the formula is far off assert p_normal < 0.05 and p_perm < 0.05 # one season, barely # --- 6. where the persistence goes ------------------------------------------------------ field_rows = [] for yr in range(FIRST, LAST): d = playoffs[playoffs.season == yr] field = sorted(set(d.home_team.replace(FRANCHISE)) | set(d.away_team.replace(FRANCHISE))) a, b = season(yr), season(yr + 1) f = [t for t in field if t in a.index and t in b.index] field_rows.append((yr, len(f), tau_b(a.loc[f, "wins"].to_numpy(float), b.loc[f, "wins"].to_numpy(float)))) FIELD = pd.DataFrame(field_rows, columns=["season", "size", "tau_b"]) FIELD["full"] = PAIRS.tau_b.to_numpy() f23 = sorted(set(playoffs[playoffs.season == CON_YEAR].home_team.replace(FRANCHISE)) | set(playoffs[playoffs.season == CON_YEAR].away_team.replace(FRANCHISE))) f23_x = season(CON_YEAR).loc[f23, "wins"].to_numpy(float) f23_y = season(CON_YEAR + 1).loc[f23, "wins"].to_numpy(float) f23C, f23D, f23Tx, f23Ty, f23Txy = pair_counts(f23_x, f23_y) beats = int((PAIRS.tb_pdiff > PAIRS.tau_b).sum()) sign_p = 2 * sum(math.comb(len(PAIRS), i) for i in range(beats, len(PAIRS) + 1)) / 2 ** len(PAIRS) nhl = pd.read_csv(NHL_CSV) nhl["team"] = nhl.team.replace({"Utah Hockey Club": "Utah Mammoth"}) s_old, s_new = sorted(nhl.season.unique()) NA = nhl[nhl.season == s_old].set_index("team") NB = nhl[nhl.season == s_new].set_index("team") ncom = sorted(set(NA.index) & set(NB.index)) nx = NA.loc[ncom, "points"].to_numpy(float) ny = NB.loc[ncom, "points"].to_numpy(float) with sdt.snippet("restrict"): print("restrict the same measurement to the teams that made the playoffs that " "season:\n") print(f" full league median tau-b {FIELD.full.median():.4f} " f"positive in {int((FIELD.full > 0).sum())} of {len(FIELD)}") print(f" playoff field median tau-b {FIELD.tau_b.median():.4f} " f"positive in {int((FIELD.tau_b > 0).sum())} of {len(FIELD)}") print(f" the field is the lower of the two in " f"{int((FIELD.tau_b < FIELD.full).sum())} of {len(FIELD)} pairs") print(f"\n {CON_YEAR}'s {len(f23)} playoff teams, into {CON_YEAR + 1}: " f"concordant {f23C}, discordant {f23D}, tied {f23Tx + f23Ty + f23Txy} " f"-> tau-b {tau_b(f23_x, f23_y):+.4f}") print("\ndoes last season's point differential order next season better than its " "wins do?\n") print(f" median tau-b from wins {PAIRS.tau_b.median():.4f}") print(f" median tau-b from point differential {PAIRS.tb_pdiff.median():.4f}") print(f" point differential wins {beats} of the {len(PAIRS)} head-to-heads, " f"sign-test two-sided p = {sign_p:.4f}") print(f" (the two orders agree within a season at a median tau-b of " f"{PAIRS.tb_wins_pdiff.median():.4f}, so they are mostly the same ranking)") print(f"\none season pair from another sport, for scale - NHL points, " f"{s_old} into {s_new}:") print(f" tau-b {tau_b(nx, ny):.4f} rho {spearman_rho(nx, ny):.4f} " f"Pearson {float(np.corrcoef(nx, ny)[0, 1]):.4f} " f"two-sided p {two_sided_p(tau_normal_z(nx, ny)):.4f}") assert list(FIELD["size"].unique()) == [12, 14] assert (FIELD[FIELD.season <= 2019]["size"] == 12).all() assert (FIELD[FIELD.season >= 2020]["size"] == 14).all() assert abs(FIELD.tau_b.median() - 0.0963) < 5e-4 assert abs(FIELD.full.median() - 0.2049) < 5e-4 assert int((FIELD.tau_b > 0).sum()) == 16 # against 26 of 26 league-wide assert int((FIELD.tau_b < FIELD.full).sum()) == 20 assert FIELD.tau_b.median() < 0.5 * FIELD.full.median() assert len(f23) == 14 and (f23C, f23D) == (35, 35) assert abs(tau_b(f23_x, f23_y)) < 1e-12 # exactly zero: 35 against 35 assert beats == 17 and abs(sign_p - 0.1686) < 5e-4 assert sign_p > 0.05 # so: not established assert abs(PAIRS.tb_pdiff.median() - 0.2227) < 5e-4 assert PAIRS.tb_pdiff.median() > PAIRS.tau_b.median() assert abs(PAIRS.tb_wins_pdiff.median() - 0.7741) < 5e-4 assert len(ncom) == 32 and s_old == "2024-25" and s_new == "2025-26" assert abs(tau_b(nx, ny) - 0.1423) < 5e-4 assert abs(spearman_rho(nx, ny) - 0.2096) < 5e-4 assert two_sided_p(tau_normal_z(nx, ny)) > 0.05 # one pair settles nothing # the fast path and the readable path are the same estimator assert abs(tau_b_fast(sign_matrix(X), Y) - tau_b(X, Y)) < 1e-12 assert abs(tau_b_fast(sign_matrix(c4_x), c4_y) - tau_b(c4_x, c4_y)) < 1e-12 # --- 7. the exhibit --------------------------------------------------------------------- BROWN, GREY = sdt.sport_color("football"), "#8A8577" fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10.6, 4.9), gridspec_kw={"width_ratios": [1.5, 1]}) ax1.bar(PAIRS.season, PAIRS.tau_b, color=BROWN, width=0.72, zorder=2, label="Kendall tau-b (wins)") ax1.plot(PAIRS.season, PAIRS.rho, color="#2C5E8A", lw=1.5, marker="o", markersize=3.2, zorder=3, label="Spearman rho") # the median gets a legend entry rather than a floating label: a free-standing # annotation here lands on the bars and on its own dashed line. ax1.axhline(float(PAIRS.tau_b.median()), color=GREY, lw=1.2, ls="--", zorder=4, label=f"median tau-b {PAIRS.tau_b.median():.3f}") ax1.set_ylim(0, 0.78) ax1.set_xlabel("season, measured against the season after it") ax1.set_ylabel("rank correlation with next season's wins") ax1.set_title("Every pair of seasons since 1999", fontsize=11.5) ax1.legend(loc="upper right", fontsize=8.2, frameon=False) ax1.grid(True, axis="y") ax2.scatter(FIELD.full, FIELD.tau_b, s=34, color=BROWN, alpha=0.85, linewidths=0, zorder=3) lim = [-0.28, 0.52] ax2.plot(lim, lim, color="#2C5E8A", lw=1.2, ls="--", zorder=2) ax2.axhline(0, color=GREY, lw=0.9, zorder=1) # every point has x >= 0.088, so the upper-left quadrant is empty ax2.text(-0.255, 0.42, f"{int((FIELD.tau_b < FIELD.full).sum())} of {len(FIELD)} below the line:" + chr(10) + "order holds worse", fontsize=8.2, color="#2C5E8A", ha="left", va="center") ax2.set_xlim(lim) ax2.set_ylim(lim) ax2.set_xlabel("tau-b, all 32 teams") ax2.set_ylabel("tau-b, that season's playoff field only") ax2.set_title("The signal is in the bad teams", fontsize=11.5) ax2.grid(True, axis="both") fig.tight_layout() sdt.save_fig(fig, "rank_persistence", source="nflverse games table (github.com/nflverse/nfldata)", asof="June 2026") print("\nall asserts passed")