"""
Tutorial 92 - From point spread to win probability: calibrate a market line on
26 seasons of NFL games.

A point spread is a prediction about the MARGIN; a win probability is a
prediction about the SIGN of the margin. This tutorial converts one into the
other three ways and checks each against what actually happened in 7,261
decided NFL games with a closing line (1999-2025, from the nflverse games
table): the raw empirical win rate at each spread, the classic normal-margin
model (one number - the standard deviation of margin around the spread - does
all the work), and a one-feature logistic regression fit by Newton's method.
It ends by pricing the opening week of 2026 from the lines bundled in the file
and cross-checking those probabilities against the moneylines the same market
posted, with the bookmaker's margin removed.

Runs entirely offline from the bundled nfl_games_lines.csv next to this script
(no API, no key). Run:  python downloads/92_point_spread_to_win_probability.py

Data: nflverse games table (github.com/nflverse/nfldata), snapshot June 2026.
"""
import os
from math import erf, sqrt

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

import sdt_common as sdt

sdt.init("point-spread-to-win-probability")
HERE = os.path.dirname(os.path.abspath(__file__))
CSV = os.path.join(HERE, "nfl_games_lines.csv")

# --- 1. load and learn the sign convention -----------------------------------------
games = pd.read_csv(CSV)
played = games.dropna(subset=["result", "spread_line"]).copy()
played["result"] = played["result"].astype(int)
ties = played[played.result == 0]
dec = played[played.result != 0].copy()          # decided games only
dec["home_win"] = (dec.result > 0).astype(int)

with sdt.snippet("load"):
    print(f"{len(games)} rows in the file, seasons {games.season.min()}-{games.season.max()}")
    print(f"{games.result.notna().sum()} played, {games.spread_line.notna().sum()} carry a spread, "
          f"{len(played)} both")
    print(f"{len(ties)} ties dropped -> {len(dec)} decided games with a closing line\n")
    sdt.show_df(dec[["game_id", "away_team", "home_team", "result", "spread_line"]], n=3)
    r = np.corrcoef(dec.spread_line, dec.result)[0, 1]
    print(f"\ncorr(spread_line, result) = {r:+.3f}  -> positive spread = home team favored")
    print(f"home teams won {100 * dec.home_win.mean():.1f}% and were favored in "
          f"{100 * (dec.spread_line > 0).mean():.1f}% of games "
          f"({(dec.spread_line == 0).sum()} pick'ems)")

assert len(games) == 7548 and len(played) == 7276 and len(ties) == 15 and len(dec) == 7261
assert 0.35 < r < 0.45, r
assert 56.0 < 100 * dec.home_win.mean() < 56.6
assert 64.5 < 100 * (dec.spread_line > 0).mean() < 65.5

# --- 2. the empirical answer: fold to the favorite and bin ------------------------
fav = dec[dec.spread_line != 0].copy()
fav["fav_spread"] = fav.spread_line.abs()
fav["fav_win"] = (np.sign(fav.result) == np.sign(fav.spread_line)).astype(int)
edges = [0, 2.5, 3, 3.5, 6.5, 7, 7.5, 9.5, 13.5, 30]
labels = ["0.5-2.5", "3", "3.5", "4-6.5", "7", "7.5", "8-9.5", "10-13.5", "14+"]
fav["fav_margin"] = fav.result * np.sign(fav.spread_line)      # margin from the favorite's side
fav["covered"] = (fav.fav_margin > fav.fav_spread).astype(int)   # beat the spread (pushes count as no)
fav["bucket"] = pd.cut(fav.fav_spread, bins=edges, labels=labels, right=True)
table = fav.groupby("bucket", observed=True).agg(
    games=("fav_win", "size"), fav_win_pct=("fav_win", "mean"), covered_pct=("covered", "mean"))
table[["fav_win_pct", "covered_pct"]] *= 100

with sdt.snippet("bins"):
    print(f"{len(fav)} decided games with a favorite\n")
    print("favorite's win rate, and cover rate, by size of the spread:")
    print(table.round(1).to_string())
    print(f"\nall favorites: won {100 * fav.fav_win.mean():.1f}%, "
          f"covered {100 * fav.covered.mean():.1f}%")
    at3 = fav[fav.fav_spread == 3]
    push3 = 100 * (at3.fav_margin == 3).mean()
    print(f"pushes at exactly 3: {push3:.1f}% of {len(at3)} games (counted as not covered)")
    cover3 = 100 * at3[at3.fav_margin != 3].covered.mean()
    print(f"cover rate at 3 among the {100 - push3:.1f}% of games that were decided: {cover3:.1f}%")

assert len(fav) == 7230
assert 8.5 < push3 < 9.5 and len(at3) == int(table.loc["3", "games"])
assert 47 < cover3 < 48.5 and abs(cover3 - table.loc["3", "covered_pct"] / (1 - push3 / 100)) < 0.05
assert table.loc["3", "games"] > 700 and 58 < table.loc["3", "fav_win_pct"] < 62
assert 74 < table.loc["7", "fav_win_pct"] < 78 and 45 < table.loc["7", "covered_pct"] < 50
assert 52 < table.loc["0.5-2.5", "fav_win_pct"] < 54
assert 82 < table.loc["10-13.5", "fav_win_pct"] < 84
assert (table.covered_pct < 50.5).all() and 45 < 100 * fav.covered.mean() < 50
assert table.loc["14+", "fav_win_pct"] > 88
assert 66 < 100 * fav.fav_win.mean() < 68

# --- 3. the normal-margin model: one number does all the work ----------------------
resid = dec.result - dec.spread_line
sigma = resid.std(ddof=1)


def phi(z):
    """Standard normal CDF, vectorised, no scipy needed."""
    return np.vectorize(lambda t: 0.5 * (1 + erf(t / sqrt(2))))(np.asarray(z, dtype=float))


def p_home_normal(spread, s=sigma):
    return phi(np.asarray(spread, dtype=float) / s)


with sdt.snippet("normal"):
    print(f"margin minus spread: mean {resid.mean():+.2f}, sd {sigma:.2f} points "
          f"over {len(dec)} games")
    print("(a mean near zero is the market being unbiased; the sd is the whole model)\n")
    print("P(favorite wins) = Phi(spread / sd):")
    for s in [1, 2.5, 3, 4, 6.5, 7, 10, 14]:
        print(f"  {s:>4} -point favorite  {100 * p_home_normal(s):5.1f}%")
    sd_early = (dec[dec.season <= 2012].result - dec[dec.season <= 2012].spread_line).std(ddof=1)
    sd_late = (dec[dec.season >= 2013].result - dec[dec.season >= 2013].spread_line).std(ddof=1)
    print(f"\nsd by era: 1999-2012 {sd_early:.2f}, 2013-2025 {sd_late:.2f}")

assert 13.0 < sigma < 13.4, sigma
assert abs(resid.mean()) < 0.25
assert 58.5 < 100 * p_home_normal(3) < 59.5
assert 69.5 < 100 * p_home_normal(7) < 70.5
assert 85 < 100 * p_home_normal(14) < 86.5
assert abs(sd_early - sd_late) < 1.0

# --- 4. a one-feature logistic regression, fit by Newton's method ------------------
X = np.column_stack([np.ones(len(dec)), dec.spread_line.to_numpy(float)])
y = dec.home_win.to_numpy(float)
beta = np.zeros(2)
for it in range(25):
    p = 1 / (1 + np.exp(-X @ beta))
    W = p * (1 - p)
    grad = X.T @ (y - p)
    hess = X.T @ (X * W[:, None])
    step = np.linalg.solve(hess, grad)
    beta += step
    if np.abs(step).max() < 1e-10:
        break
a, b = beta


def p_home_logit(spread):
    return 1 / (1 + np.exp(-(a + b * np.asarray(spread, dtype=float))))


with sdt.snippet("logistic"):
    print(f"logit P(home win) = {a:+.4f} + {b:.4f} x spread   (Newton, {it + 1} iterations)")
    print(f"a pick'em at home: {100 * p_home_logit(0):.1f}%   "
          f"one point of spread multiplies the odds by {np.exp(b):.3f}\n")
    print("favorite by   empirical (n)     normal   logistic")
    for s, lab in [(1.5, "0.5-2.5"), (3, "3"), (3.5, "3.5"), (5, "4-6.5"), (7, "7"),
                   (11, "10-13.5"), (16, "14+")]:
        emp = table.loc[lab]
        print(f"  {lab:>8}     {emp.fav_win_pct:5.1f}% ({int(emp.games):>4})    "
              f"{100 * p_home_normal(s):5.1f}%    {100 * p_home_logit(s):5.1f}%")
    print("(the normal and logistic columns are evaluated at a typical spread in the bucket)")

assert abs(a) < 0.06 and 0.13 < b < 0.15, (a, b)
assert 49 < 100 * p_home_logit(0) < 51.5
assert abs(100 * p_home_logit(3) - 100 * p_home_normal(3)) < 2.0
assert 0 < 100 * p_home_logit(7) - 100 * p_home_normal(7) < 3.0   # logistic is the steeper curve
assert 0 < 100 * p_home_logit(14) - 100 * p_home_normal(14) < 5.0

# --- 5. hold out the last six seasons and score all three --------------------------
train = dec[dec.season <= 2019]
test = dec[dec.season >= 2020]
sig_tr = (train.result - train.spread_line).std(ddof=1)
Xtr = np.column_stack([np.ones(len(train)), train.spread_line.to_numpy(float)])
ytr = train.home_win.to_numpy(float)
bt = np.zeros(2)
for _ in range(25):
    pt = 1 / (1 + np.exp(-Xtr @ bt))
    bt += np.linalg.solve(Xtr.T @ (Xtr * (pt * (1 - pt))[:, None]), Xtr.T @ (ytr - pt))
yt = test.home_win.to_numpy(float)
preds = {
    "always the base rate": np.full(len(test), ytr.mean()),
    "normal model": p_home_normal(test.spread_line, sig_tr),
    "logistic": 1 / (1 + np.exp(-(bt[0] + bt[1] * test.spread_line.to_numpy(float)))),
}
brier = {k: np.mean((v - yt) ** 2) for k, v in preds.items()}
ll = {k: -np.mean(yt * np.log(v) + (1 - yt) * np.log(1 - v)) for k, v in preds.items()}

with sdt.snippet("holdout"):
    print(f"fit on 1999-2019 ({len(train)} games), scored on 2020-2025 ({len(test)} games)")
    print(f"training sd {sig_tr:.2f}; training home-win rate {100 * ytr.mean():.1f}%\n")
    print(f"{'model':<22}{'Brier':>8}{'log loss':>10}")
    for k in preds:
        print(f"{k:<22}{brier[k]:8.4f}{ll[k]:10.4f}")
    skill = 100 * (1 - brier["normal model"] / brier["always the base rate"])
    print(f"\nnormal model beats the base rate by {skill:.1f}% of Brier score")
    hit = 100 * np.mean((preds["normal model"] > 0.5) == (yt == 1))
    print(f"picking the favorite was right {hit:.1f}% of the time in the hold-out")

assert len(train) == 5573 and len(test) == 1688
assert brier["normal model"] < brier["always the base rate"]
assert abs(brier["normal model"] - brier["logistic"]) < 0.002
assert 10 < skill < 16, skill
assert 64 < hit < 70, hit

# --- 6. the exhibit: every spread's real win rate against both curves ---------------
by_spread = dec.groupby("spread_line").agg(n=("home_win", "size"), rate=("home_win", "mean"))
by_spread = by_spread[by_spread.n >= 30]
grid = np.linspace(-17, 17, 200)
BROWN = sdt.sport_color("football")
fig, ax = plt.subplots(figsize=(8.8, 5.2))
ax.axhline(50, color="#C2B7A1", lw=0.8, zorder=1)
ax.axvline(0, color="#C2B7A1", lw=0.8, zorder=1)
ax.scatter(by_spread.index, 100 * by_spread.rate, s=by_spread.n / 4, color=BROWN,
           alpha=0.55, edgecolor="none", zorder=3,
           label="actual home win rate at that spread (dot area = games)")
ax.plot(grid, 100 * p_home_normal(grid), color="#20242B", lw=1.8, zorder=4,
        label=f"normal model, sd = {sigma:.1f}")
ax.plot(grid, 100 * p_home_logit(grid), color="#2C5E8A", lw=1.4, ls="--", zorder=4,
        label="logistic regression")
for s in (3, 7):
    ax.annotate(f"home by {s}: {100 * p_home_normal(s):.0f}%", (s, 100 * p_home_normal(s)),
                xytext=(s + 2.2, 100 * p_home_normal(s) - 9), fontsize=9,
                arrowprops=dict(arrowstyle="-", color="#6C7079", lw=0.8))
ax.set_xlim(-17.5, 17.5)
ax.set_ylim(0, 100)
ax.set_xlabel("closing spread (positive = home team favored by that many points)")
ax.set_ylabel("home team win rate, %")
ax.set_title("The spread is already a win probability - once you know the sd (1999-2025)")
ax.legend(loc="upper left", fontsize=8.5, frameon=False)
sdt.save_fig(fig, "spread_curve", source="nflverse games table", asof="June 2026 snapshot")

assert (by_spread.n >= 30).all() and 25 <= len(by_spread) <= 45

# --- 7. price the opening week of 2026 and cross-check the moneyline ---------------
w1 = games[(games.season == 2026) & (games.week == 1)].copy()


def implied(ml):
    """American moneyline -> the probability it implies (still carrying the vig)."""
    ml = np.asarray(ml, dtype=float)
    return np.where(ml < 0, -ml / (-ml + 100), 100 / (ml + 100))


w1["p_home_spread"] = 100 * p_home_normal(w1.spread_line)
raw_h, raw_a = implied(w1.home_moneyline), implied(w1.away_moneyline)
w1["overround"] = 100 * (raw_h + raw_a - 1)
w1["p_home_ml"] = 100 * raw_h / (raw_h + raw_a)
w1["gap"] = w1.p_home_spread - w1.p_home_ml
w1 = w1.sort_values(["gameday", "game_id"])

with sdt.snippet("week1"):
    print(f"2026 week 1: {len(w1)} games, lines as bundled in the June 2026 nflverse snapshot\n")
    show = w1[["gameday", "away_team", "home_team", "location", "spread_line",
               "p_home_spread", "home_moneyline", "p_home_ml", "gap"]].copy()
    show.columns = ["date", "away", "home", "venue", "spread", "p(home) spread",
                    "home ML", "p(home) ML", "gap"]
    print(show.round(1).to_string(index=False))
    print(f"\nbookmaker overround on the moneylines: {w1.overround.mean():.1f}% average")
    print(f"spread-model vs vig-free moneyline: mean gap {w1.gap.mean():+.1f} pts, "
          f"largest {w1.gap.abs().max():.1f} pts")
    gap_logit = 100 * p_home_logit(w1.spread_line) - w1.p_home_ml
    print(f"the steeper logistic curve instead: mean gap {gap_logit.mean():+.1f} pts, "
          f"largest {gap_logit.abs().max():.1f} pts")
    biggest = w1.loc[w1.spread_line.abs().idxmax()]
    print(f"biggest favorite: {biggest.home_team} by {biggest.spread_line} over "
          f"{biggest.away_team} -> {biggest.p_home_spread:.1f}% by the spread, "
          f"{biggest.p_home_ml:.1f}% by the moneyline")

assert len(w1) == 16 and w1.spread_line.notna().all() and w1.home_moneyline.notna().all()
assert (w1.game_id == "2026_01_NE_SEA").any()
assert float(w1.loc[w1.game_id == "2026_01_NE_SEA", "spread_line"].iloc[0]) == 3.5
assert (w1.loc[w1.game_id == "2026_01_SF_LA", "location"] == "Neutral").all()
assert biggest.home_team == "LAC" and biggest.spread_line == 11.5
assert 3 < w1.overround.mean() < 6, w1.overround.mean()
assert w1.gap.abs().max() < 5.0, w1.gap.abs().max()
assert abs(gap_logit.mean()) < abs(w1.gap.mean())   # the steeper curve sits closer to the moneyline
assert abs(w1.gap.mean()) < 2.0
print("\nall asserts passed")