Learning ObjectivesBy the end of this chapter, you will be able to:

1. Build Monte Carlo simulations for football scenarios
2. Simulate games and seasons with realistic uncertainty
3. Calculate playoff probabilities using simulation
4. Analyze strategic scenarios through simulation
5. Quantify uncertainty through resampling methods
6. Apply sensitivity analysis to test model assumptions
7. Use bootstrapping to estimate confidence intervals
8. Design and interpret what-if scenario analyses

Introduction

Football is inherently uncertain. Will the quarterback complete the pass? Will the running back break a tackle? Will the kicker make the field goal? Traditional analytics often provides point estimates—the average EPA of a play, the expected win probability—but these single numbers mask the underlying uncertainty that defines every moment of football.

Consider a team with a 65% playoff probability heading into Week 17. What does this number actually mean? It doesn't mean the team will make 65% of the playoffs—they'll either make it or they won't. Rather, it means that if we could replay the remaining games thousands of times under identical conditions, the team would make the playoffs in about 65% of those replays. This is the fundamental insight of Monte Carlo simulation: we can't predict the future with certainty, but we can quantify the range of possible futures and their relative likelihoods.

Monte Carlo simulation offers a powerful framework for embracing and quantifying this uncertainty. Instead of asking "what will happen?", we ask "what could happen?" By simulating thousands or millions of possible outcomes, we can:

• Estimate playoff probabilities across uncertain futures, accounting for schedule strength and randomness
• Evaluate strategic decisions under various scenarios, from fourth-down calls to draft pick trades
• Build player projections that capture uncertainty, showing ranges rather than single numbers
• Construct optimal rosters using portfolio theory, balancing upside and consistency
• Stress-test assumptions through sensitivity analysis, identifying which inputs matter most
• Generate confidence intervals for any metric, from EPA to win probability

The beauty of simulation is its flexibility. Unlike analytical methods that require mathematical tractability, simulation can handle any level of complexity. Want to model player injuries, weather conditions, and referee tendencies all at once? Simulation can do it. Need to account for coaching adjustments, player fatigue, and momentum shifts? Simulation handles it naturally. This makes it an indispensable tool in modern football analytics, complementing the regression models and machine learning approaches we've explored in previous chapters.

In this chapter, we'll start with the foundational principles of Monte Carlo simulation, then progress through increasingly sophisticated applications: from simple field goal models to complex season simulators, from strategic decision analysis to portfolio-based roster construction. Along the way, we'll explore practical considerations like choosing the right number of simulations, validating model assumptions, and communicating uncertainty to stakeholders who prefer certainty.

Monte Carlo Basics

The Core Principle

At its heart, Monte Carlo simulation is deceptively simple. Rather than trying to calculate probabilities through complex mathematics, we simulate random outcomes many times and count what happens. If we want to know the probability a team makes the playoffs, we simulate the rest of the season 10,000 times and count how many times they make it. If they made it 6,742 times out of 10,000, our estimate is 67.42%.

This approach works because of the Law of Large Numbers, a fundamental theorem in probability theory. As the number of simulations increases, the average of the simulated outcomes converges to the true expected value. The more simulations we run, the closer we get to the "true" probability—though of course, we never know that true value with certainty.

Monte Carlo simulation follows a simple four-step recipe:

1. Define the uncertainty: Identify random variables and their distributions. What aspects of the problem are uncertain? For a field goal, this might be wind speed, kick angle, and kicker skill. For a game, it might be turnovers, third-down conversions, and big plays.

2. Generate random samples: Draw values from these distributions. This is where the "Monte" in Monte Carlo comes in—we're essentially rolling dice repeatedly. Modern computers make this trivial with pseudo-random number generators.

3. Calculate outcomes: Compute results for each sample. Given the random inputs, what happens? If the kicker's accuracy is 87% and we randomly generate a number, did they make the kick or miss it?

4. Aggregate results: Summarize across all simulations. After running thousands of iterations, what patterns emerge? What's the average outcome? What's the distribution? What percentile thresholds matter?

The power of this approach is that step 3 can be arbitrarily complex. We can model intricate dependencies, conditional probabilities, and nonlinear relationships that would be mathematically intractable. The computer doesn't care—it just runs the simulation over and over.

Random Number Generation

All simulation starts with random number generation. Modern computers use pseudo-random number generators (PRNGs) that produce sequences that appear random but are deterministic given a seed. This determinism is actually a feature, not a bug—it allows us to reproduce our results exactly by setting the same seed.

The most common distributions in football simulation are:

• Uniform: All values in a range equally likely. Used for coin flips, random game outcomes, and generating other distributions.
• Normal (Gaussian): Bell-curved distribution. Models many natural phenomena like play outcomes, scoring margins, and measurement errors.
• Binomial: Success/failure trials. Perfect for field goals, conversion attempts, and win/loss outcomes.
• Poisson: Count data with a given rate. Can model scoring events, turnovers, or sacks per game.

Let's visualize these fundamental distributions to build intuition:

#| label: setup-r
#| message: false
#| warning: false

library(tidyverse)
library(nflfastR)
library(nflplotR)
library(gt)
library(gridExtra)

# Set seed for reproducibility
# Using 2024 ensures we get the same "random" numbers every time
set.seed(2024)

# Generate random numbers from different distributions
# Uniform: flat distribution, all values equally likely
uniform_samples <- runif(1000, min = 0, max = 1)

# Normal: bell curve centered at 0 with spread of 1
normal_samples <- rnorm(1000, mean = 0, sd = 1)

# Poisson: count distribution, lambda=3 means avg of 3 events
poisson_samples <- rpois(1000, lambda = 3)

# Create comparison plot to visualize distribution shapes
tibble(
  uniform = uniform_samples,
  normal = normal_samples,
  poisson = poisson_samples
) %>%
  pivot_longer(everything(), names_to = "distribution", values_to = "value") %>%
  ggplot(aes(x = value, fill = distribution)) +
  geom_histogram(alpha = 0.7, bins = 30) +
  facet_wrap(~distribution, scales = "free") +
  labs(
    title = "Common Probability Distributions",
    subtitle = "1,000 random samples from each distribution",
    x = "Value",
    y = "Count"
  ) +
  theme_minimal() +
  theme(legend.position = "none")

#| label: setup-py
#| message: false
#| warning: false

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scipy import stats

# Set seed for reproducibility
np.random.seed(2024)

# Generate random numbers from different distributions
uniform_samples = np.random.uniform(0, 1, 1000)
normal_samples = np.random.normal(0, 1, 1000)
poisson_samples = np.random.poisson(3, 1000)

# Create comparison plot
fig, axes = plt.subplots(1, 3, figsize=(12, 4))

axes[0].hist(uniform_samples, bins=30, alpha=0.7, color='#F8766D')
axes[0].set_title('Uniform Distribution')
axes[0].set_xlabel('Value')
axes[0].set_ylabel('Count')

axes[1].hist(normal_samples, bins=30, alpha=0.7, color='#00BA38')
axes[1].set_title('Normal Distribution')
axes[1].set_xlabel('Value')

axes[2].hist(poisson_samples, bins=range(0, 12), alpha=0.7, color='#619CFF')
axes[2].set_title('Poisson Distribution')
axes[2].set_xlabel('Value')

plt.tight_layout()
plt.show()

These visualizations show how 1,000 random samples from each distribution create different patterns. Notice how the uniform distribution spreads evenly, the normal clusters around zero, and the Poisson produces discrete counts. Understanding these shapes helps you choose the right distribution for different football phenomena.

A Simple Example: Field Goal Accuracy

Let's start with a simple Monte Carlo simulation to build intuition: estimating field goal percentage. Suppose we know a kicker has an 85% success rate. We could calculate that after 1,000 kicks, they'd expect to make 850 of them. But simulation lets us see the natural variation around that expectation.

In real games, a kicker might get "hot" or "cold" due to random chance alone. How much variation should we expect? Simulation answers this by playing out the randomness explicitly. Each simulated kick is a Bernoulli trial (success or failure), and we aggregate the results to see what percentage of kicks the kicker makes across our simulation.

This simple example illustrates all four steps of Monte Carlo simulation: we define uncertainty (85% success rate), generate random samples (1,000 kicks), calculate outcomes (make or miss for each), and aggregate results (overall percentage). Let's see it in action:

#| label: fg-simple-r
#| message: false
#| warning: false

# Simulate 1,000 field goal attempts
# Assume 85% success rate (typical for NFL kickers in moderate conditions)
n_kicks <- 1000
make_probability <- 0.85

# Simulate outcomes (1 = make, 0 = miss)
# rbinom generates random binomial trials: n_kicks trials, each with make_probability
kick_results <- rbinom(n_kicks, size = 1, prob = make_probability)

# Calculate make percentage from our simulation
make_pct <- mean(kick_results)

cat("Simulated make percentage:", round(make_pct * 100, 1), "%\n")
cat("True probability:", make_probability * 100, "%\n")
cat("Difference:", round(abs(make_pct - make_probability) * 100, 2), "%\n")

#| label: fg-simple-py
#| message: false
#| warning: false

# Simulate 1,000 field goal attempts
n_kicks = 1000
make_probability = 0.85

# Simulate outcomes using binomial distribution
# Each kick is an independent trial with 85% success rate
kick_results = np.random.binomial(1, make_probability, n_kicks)

# Calculate make percentage
make_pct = kick_results.mean()

print(f"Simulated make percentage: {make_pct*100:.1f}%")
print(f"True probability: {make_probability*100:.0f}%")
print(f"Difference: {abs(make_pct - make_probability)*100:.2f}%")

Notice that our simulated make percentage isn't exactly 85%—it might be 84.7% or 85.3% due to random chance. This sampling error decreases as we increase the number of simulations. With 1,000 kicks, we typically get within 1-2% of the true value. With 10,000 kicks, we'd get much closer. This relationship between sample size and accuracy is fundamental to understanding Monte Carlo methods.

Convergence with Sample Size

A key principle of Monte Carlo simulation is that accuracy improves with more samples. This isn't just theoretical—we can demonstrate it directly by running simulations with increasing sample sizes and watching our estimates converge to the true value.

The relationship follows a mathematical law: estimation error decreases proportional to 1/√n, where n is the sample size. Double your simulations, cut your error by about 1.4×. Quadruple them, cut it by 2×. This square root relationship means we hit diminishing returns quickly—going from 1,000 to 10,000 simulations helps a lot, but going from 100,000 to 1,000,000 helps much less.

Let's visualize this convergence to build intuition about how many simulations we need for reliable results:

#| label: fig-convergence-r
#| fig-cap: "Monte Carlo convergence with sample size"
#| fig-width: 10
#| fig-height: 6

# Run simulations with increasing sample sizes
# Start at 10, go up to 10,000, checking every 10 kicks
sample_sizes <- seq(10, 10000, by = 10)
true_probability <- 0.85

# For each sample size, simulate that many kicks and calculate make %
estimates <- map_dbl(sample_sizes, function(n) {
  mean(rbinom(n, size = 1, prob = true_probability))
})

# Plot convergence
tibble(
  n = sample_sizes,
  estimate = estimates,
  error = abs(estimates - true_probability)
) %>%
  ggplot(aes(x = n, y = estimate)) +
  geom_line(color = "#619CFF", linewidth = 1) +
  # Red dashed line shows true value
  geom_hline(yintercept = true_probability,
             linetype = "dashed", color = "red", linewidth = 1) +
  # Shaded band shows ±2% margin of error
  geom_ribbon(aes(ymin = true_probability - 0.02,
                  ymax = true_probability + 0.02),
              alpha = 0.2, fill = "red") +
  scale_x_continuous(labels = scales::comma) +
  scale_y_continuous(labels = scales::percent, limits = c(0.75, 0.95)) +
  labs(
    title = "Monte Carlo Convergence to True Value",
    subtitle = "Field goal make probability (true value = 85%)",
    x = "Number of Simulations",
    y = "Estimated Make Probability",
    caption = "Red band shows ±2% margin | Estimates stabilize around 1,000 simulations"
  ) +
  theme_minimal() +
  theme(plot.title = element_text(face = "bold"))

#| label: fig-convergence-py
#| fig-cap: "Monte Carlo convergence with sample size - Python"
#| fig-width: 10
#| fig-height: 6

# Run simulations with increasing sample sizes
sample_sizes = range(10, 10001, 10)
true_probability = 0.85

# Simulate field goals for each sample size
estimates = [np.random.binomial(1, true_probability, n).mean()
             for n in sample_sizes]

# Plot convergence
plt.figure(figsize=(10, 6))
plt.plot(sample_sizes, estimates, color='#619CFF', linewidth=2)
plt.axhline(y=true_probability, color='red', linestyle='--', linewidth=2,
            label='True Probability')
plt.fill_between(sample_sizes,
                 true_probability - 0.02,
                 true_probability + 0.02,
                 alpha=0.2, color='red')
plt.xlabel('Number of Simulations', fontsize=12)
plt.ylabel('Estimated Make Probability', fontsize=12)
plt.title('Monte Carlo Convergence to True Value\nField goal make probability (true value = 85%)',
          fontsize=14, fontweight='bold')
plt.text(0.98, 0.02, 'Red band shows ±2% margin | Estimates stabilize around 1,000 simulations',
         transform=plt.gca().transAxes, ha='right', fontsize=9)
plt.ylim(0.75, 0.95)
plt.legend()
plt.tight_layout()
plt.show()

The plot shows our estimate bouncing around when sample sizes are small, then settling down as we run more simulations. By 1,000 simulations, we're consistently within ±2% of the true value. By 10,000, we're even tighter. This visual proof demonstrates why we typically use 10,000 simulations for production analytics—it's where we get stable, trustworthy results without excessive computation time.

Game Simulation Models

Modeling Game Outcomes

Moving from simple field goal simulations to full game simulations requires modeling the complex interaction between two teams competing against each other. The simplest approach uses team strength ratings—single numbers representing each team's quality—combined with random variation to generate realistic game outcomes.

This approach builds on the ratings systems we explored in Chapter 33 (ELO and Ratings). The key insight is that game outcomes aren't deterministic—even when a strong team plays a weak team, the underdog sometimes wins due to turnovers, injuries, referee calls, and pure luck. We model this by adding random noise to the expected outcome.

The typical model has three components:

1. Team ratings: Numbers representing offensive/defensive strength, typically calibrated so 0 = average team, positive = above average, negative = below average. A +7 rating means the team beats an average team by about 7 points on average.

2. Home field advantage: The home team typically has a 2-3 point edge due to crowd noise, travel fatigue, and referee bias. We add this constant to the home team's expected performance.

3. Random variation: Even if we expect a team to win by 10 points, the actual margin follows a distribution centered at 10 with some spread (standard deviation). Research shows NFL game margins have a standard deviation around 13-14 points, reflecting the high variance in game outcomes.

Let's implement this simple but effective model:

#| label: game-sim-basic-r
#| message: false
#| warning: false

# Function to simulate a single game between two teams
simulate_game <- function(team_a_rating, team_b_rating,
                         home_advantage = 2.5,
                         sd_points = 13.5) {

  # Expected margin (positive = team A wins)
  # This is the "most likely" outcome before adding randomness
  expected_margin <- team_a_rating - team_b_rating + home_advantage

  # Add random variation using normal distribution
  # Actual game might deviate from expectation due to luck, turnovers, etc.
  actual_margin <- rnorm(1, mean = expected_margin, sd = sd_points)

  # Determine winner (positive margin = team A wins)
  team_a_wins <- actual_margin > 0

  # Estimate scores assuming both teams score around 24 points on average
  # Then adjust based on margin (winner gets more, loser gets less)
  list(
    margin = actual_margin,
    team_a_wins = team_a_wins,
    team_a_score = round(24 + actual_margin / 2),
    team_b_score = round(24 - actual_margin / 2)
  )
}

# Simulate Chiefs vs Bills (both strong teams)
chiefs_rating <- 6.5  # Chiefs are 6.5 points above average team
bills_rating <- 5.8   # Bills are 5.8 points above average team

set.seed(2024)
game_result <- simulate_game(chiefs_rating, bills_rating, home_advantage = 2.5)

cat("Simulated Game Result:\n")
cat("Chiefs:", game_result$team_a_score, "\n")
cat("Bills:", game_result$team_b_score, "\n")
cat("Margin:", round(game_result$margin, 1), "\n")
cat("Winner:", ifelse(game_result$team_a_wins, "Chiefs", "Bills"), "\n")

#| label: game-sim-basic-py
#| message: false
#| warning: false

def simulate_game(team_a_rating, team_b_rating,
                  home_advantage=2.5, sd_points=13.5):
    """
    Simulate a single game between two teams.

    Parameters:
    - team_a_rating: Team A strength (points above average)
    - team_b_rating: Team B strength (points above average)
    - home_advantage: Points added to home team (default 2.5)
    - sd_points: Standard deviation of game outcomes (default 13.5)

    Returns:
    - Dictionary with margin, winner, and scores
    """

    # Expected margin (positive = team A wins)
    expected_margin = team_a_rating - team_b_rating + home_advantage

    # Add random variation
    actual_margin = np.random.normal(expected_margin, sd_points)

    # Determine winner
    team_a_wins = actual_margin > 0

    return {
        'margin': actual_margin,
        'team_a_wins': team_a_wins,
        'team_a_score': round(24 + actual_margin / 2),
        'team_b_score': round(24 - actual_margin / 2)
    }

# Simulate Chiefs vs Bills (both strong teams)
chiefs_rating = 6.5
bills_rating = 5.8

np.random.seed(2024)
game_result = simulate_game(chiefs_rating, bills_rating, home_advantage=2.5)

print("Simulated Game Result:")
print(f"Chiefs: {game_result['team_a_score']}")
print(f"Bills: {game_result['team_b_score']}")
print(f"Margin: {game_result['margin']:.1f}")
print(f"Winner: {'Chiefs' if game_result['team_a_wins'] else 'Bills'}")

This single simulation shows one possible outcome: the Chiefs might win by 3 points, or the Bills might pull an upset. The beauty of Monte Carlo methods is that we don't stop at one simulation—we run thousands to see the full range of possibilities.

Monte Carlo Game Probabilities

A single game simulation tells us one possible outcome, but it doesn't tell us probabilities. To estimate win probability, we simulate the game thousands of times and count how often each team wins. If the Chiefs win 6,850 out of 10,000 simulations, their estimated win probability is 68.5%.

This approach has several advantages over analytical methods. First, it handles asymmetry naturally—win probability doesn't have to be a smooth function of ratings. Second, it captures tail risk—we see not just the average outcome but the full distribution including blowouts and close games. Third, it's easily extended—want to account for weather? Just modify the simulation logic.

The key parameter is how many simulations to run. Too few and our estimates will be noisy; too many and we waste computation time. For game-level probabilities, 10,000 simulations provides estimates accurate to about ±1%, which is more than sufficient given the uncertainty in our team ratings themselves.

#| label: game-monte-carlo-r
#| message: false
#| warning: false

# Simulate 10,000 games between Chiefs and Bills
# This lets us estimate win probability with ~1% accuracy
n_sims <- 10000

# Run simulations and store all results
results <- map_dfr(1:n_sims, function(i) {
  # Simulate one game
  game <- simulate_game(chiefs_rating, bills_rating, home_advantage = 2.5)

  # Store result in a data frame row
  tibble(
    sim = i,
    margin = game$margin,
    chiefs_win = game$team_a_wins,
    chiefs_score = game$team_a_score,
    bills_score = game$team_b_score
  )
})

# Calculate win probability by counting wins
chiefs_win_prob <- mean(results$chiefs_win)

cat("Based on", scales::comma(n_sims), "simulations:\n")
cat("Chiefs win probability:", round(chiefs_win_prob * 100, 1), "%\n")
cat("Bills win probability:", round((1 - chiefs_win_prob) * 100, 1), "%\n")
cat("Average margin:", round(mean(results$margin), 1), "points\n")
cat("Margin std dev:", round(sd(results$margin), 1), "points\n")

#| label: game-monte-carlo-py
#| message: false
#| warning: false

# Simulate 10,000 games between Chiefs and Bills
n_sims = 10000

# Store results in a list, then convert to DataFrame
results = []
for i in range(n_sims):
    game = simulate_game(chiefs_rating, bills_rating, home_advantage=2.5)
    results.append({
        'sim': i,
        'margin': game['margin'],
        'chiefs_win': game['team_a_wins'],
        'chiefs_score': game['team_a_score'],
        'bills_score': game['team_b_score']
    })

results_df = pd.DataFrame(results)

# Calculate win probability
chiefs_win_prob = results_df['chiefs_win'].mean()

print(f"Based on {n_sims:,} simulations:")
print(f"Chiefs win probability: {chiefs_win_prob*100:.1f}%")
print(f"Bills win probability: {(1-chiefs_win_prob)*100:.1f}%")
print(f"Average margin: {results_df['margin'].mean():.1f} points")
print(f"Margin std dev: {results_df['margin'].std():.1f} points")

The results show the Chiefs win about 68% of simulations, giving them a roughly 2:1 edge despite playing another elite team. This demonstrates how small rating differences (0.7 points) plus home field advantage (2.5 points) compound into meaningful win probability gaps. The Bills still win about 1 in 3 games—which is why we play the games rather than just calculating probabilities.

Visualizing Game Outcomes

Numbers tell part of the story, but visualizations reveal patterns that summary statistics miss. By plotting the distribution of simulated game outcomes, we can see the full range of possibilities: not just that the Chiefs win 68% of games, but how they win and by how much.

The histogram of margins will show several interesting features. First, it should be roughly bell-shaped (normal), centered around our expected 3-point margin. Second, it will be bimodal or at least lumpy due to discrete scoring (touchdowns are worth 7 points). Third, we'll see significant overlap around zero—many games are close, representing the Bills' 32% win probability.

Color-coding the histogram by winner makes it easy to see the probability mass on each side. The red area (Chiefs wins) should be about twice the size of the blue area (Bills wins), visualizing the 68-32 split. Let's create this visualization:

#| label: fig-game-outcomes-r
#| fig-cap: "Distribution of simulated game outcomes"
#| fig-width: 10
#| fig-height: 6

results %>%
  ggplot(aes(x = margin)) +
  # Histogram colored by winner
  geom_histogram(aes(fill = chiefs_win),
                 bins = 50, alpha = 0.8, color = "white") +
  # Vertical line at even game (0 margin)
  geom_vline(xintercept = 0, linetype = "dashed",
             color = "black", linewidth = 1) +
  # Vertical line at average margin
  geom_vline(xintercept = mean(results$margin),
             linetype = "solid", color = "darkgreen", linewidth = 1) +
  # Team colors for Chiefs (red) and Bills (blue)
  scale_fill_manual(
    values = c("TRUE" = "#E31837", "FALSE" = "#00338D"),
    labels = c("Chiefs Win", "Bills Win")
  ) +
  labs(
    title = "Simulated Game Outcomes: Chiefs vs Bills",
    subtitle = paste0("Chiefs win ", round(chiefs_win_prob * 100, 1),
                     "% of simulations (n=", scales::comma(n_sims), ")"),
    x = "Margin of Victory (positive = Chiefs win)",
    y = "Count",
    fill = "Winner",
    caption = "Green line shows average margin | Black dashed line shows even game"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(face = "bold"),
    legend.position = "top"
  )

#| label: fig-game-outcomes-py
#| fig-cap: "Distribution of simulated game outcomes - Python"
#| fig-width: 10
#| fig-height: 6

plt.figure(figsize=(10, 6))

# Separate data by winner for colored histogram
chiefs_wins = results_df[results_df['chiefs_win']]['margin']
bills_wins = results_df[~results_df['chiefs_win']]['margin']

# Plot histograms
plt.hist(chiefs_wins, bins=50, alpha=0.8, color='#E31837',
         label='Chiefs Win', edgecolor='white', linewidth=0.5)
plt.hist(bills_wins, bins=50, alpha=0.8, color='#00338D',
         label='Bills Win', edgecolor='white', linewidth=0.5)

# Add reference lines
plt.axvline(x=0, color='black', linestyle='--', linewidth=2,
            label='Even Game')
plt.axvline(x=results_df['margin'].mean(), color='darkgreen', linewidth=2,
            label='Average Margin')

plt.xlabel('Margin of Victory (positive = Chiefs win)', fontsize=12)
plt.ylabel('Count', fontsize=12)
plt.title(f"Simulated Game Outcomes: Chiefs vs Bills\nChiefs win {chiefs_win_prob*100:.1f}% of simulations (n={n_sims:,})",
          fontsize=14, fontweight='bold')
plt.legend(loc='upper right')
plt.tight_layout()
plt.show()

The visualization makes clear what the 68% win probability means in practice: the distribution is centered on a small Chiefs advantage, but with enough variance that the Bills win a substantial minority of games. The overlapping distributions show why football outcomes remain uncertain even when one team is clearly favored.

Season Simulation and Playoff Odds

Building a Season Simulator

Calculating playoff probabilities requires simulating entire seasons, not just individual games. This adds complexity: we need to model 17 games for each team, account for opponent strength varying by schedule, and aggregate results to determine playoff seeding. The 7-team playoff format means we need to rank all teams by record and apply tiebreakers.

Season simulation builds on our game simulation function but adds scheduling logic. Each team plays 17 games against different opponents, with some at home and some away. The sequence of opponents matters less than their overall strength, so we can simplify by having teams play random opponents from the pool. In reality, NFL schedules have structure (division games, same-finish games), but for playoff probability estimation, random schedules work reasonably well.

The key insight is that season outcomes involve compounding uncertainty. A team's final record depends on 17 game outcomes, each uncertain. Small differences in team quality (1-2 points of rating) compound over 17 games into meaningful differences in expected wins. Yet randomness remains large enough that a team with expected 10 wins might finish anywhere from 7-13 wins with reasonable probability.

Let's build a simplified but functional season simulator for the AFC:

#| label: season-sim-setup-r
#| message: false
#| warning: false
#| cache: true

# Create simplified team ratings for AFC
# Ratings based on point differential expectations vs average team
afc_teams <- tibble(
  team = c("KC", "BUF", "BAL", "MIA", "CIN", "JAX", "CLE", "HOU",
           "PIT", "LAC", "IND", "TEN", "LV", "NYJ", "NE", "DEN"),
  rating = c(6.5, 5.8, 5.2, 4.8, 4.5, 3.2, 2.8, 2.5,
             2.0, 1.5, 0.5, -1.0, -2.0, -2.5, -3.5, -4.0),
  division = c("West", "East", "North", "East", "North", "South", "North", "South",
               "North", "West", "South", "South", "West", "East", "East", "West")
)

# Function to simulate a complete season
# Returns final win totals for all teams
simulate_season <- function(teams, n_games = 17) {

  # Initialize results with zero wins
  results <- teams %>%
    select(team, rating) %>%
    mutate(wins = 0)

  # Simulate each team's season
  for(i in 1:nrow(teams)) {
    team_rating <- teams$rating[i]

    # Play against random opponents (simplified scheduling)
    # In reality, NFL schedule has structure, but random works for our purposes
    opponents <- sample(setdiff(teams$team, teams$team[i]), n_games)
    opponent_ratings <- teams$rating[teams$team %in% opponents]

    # Simulate each of the 17 games
    for(j in 1:n_games) {
      # Randomly assign home/away (50-50 split)
      is_home <- runif(1) > 0.5
      home_adv <- ifelse(is_home, 2.5, -2.5)

      # Simulate the game
      game <- simulate_game(team_rating, opponent_ratings[j],
                           home_advantage = home_adv)

      # Add win if team won
      if(game$team_a_wins) {
        results$wins[i] <- results$wins[i] + 1
      }
    }
  }

  # Return results sorted by wins (playoff seeding order)
  results %>%
    mutate(losses = n_games - wins) %>%
    arrange(desc(wins))
}

# Run one simulation to see sample output
set.seed(2024)
season_result <- simulate_season(afc_teams)

season_result %>%
  left_join(afc_teams %>% select(team, division), by = "team") %>%
  select(team, division, wins, losses, rating) %>%
  head(10) %>%
  gt() %>%
  cols_label(
    team = "Team",
    division = "Division",
    wins = "Wins",
    losses = "Losses",
    rating = "Rating"
  ) %>%
  fmt_number(
    columns = rating,
    decimals = 1
  ) %>%
  tab_header(
    title = "Simulated Season Results",
    subtitle = "Top 10 teams (one simulation)"
  ) %>%
  tab_style(
    style = cell_fill(color = "#E8F4EA"),
    locations = cells_body(rows = 1:7)
  )

#| label: season-sim-setup-py
#| message: false
#| warning: false
#| cache: true

# Create simplified team ratings for AFC
afc_teams = pd.DataFrame({
    'team': ["KC", "BUF", "BAL", "MIA", "CIN", "JAX", "CLE", "HOU",
             "PIT", "LAC", "IND", "TEN", "LV", "NYJ", "NE", "DEN"],
    'rating': [6.5, 5.8, 5.2, 4.8, 4.5, 3.2, 2.8, 2.5,
               2.0, 1.5, 0.5, -1.0, -2.0, -2.5, -3.5, -4.0],
    'division': ["West", "East", "North", "East", "North", "South", "North", "South",
                 "North", "West", "South", "South", "West", "East", "East", "West"]
})

def simulate_season(teams, n_games=17):
    """Simulate a full season for all teams."""

    # Initialize results
    results = teams[['team', 'rating']].copy()
    results['wins'] = 0

    # Each team plays n_games
    for i, row in teams.iterrows():
        team_rating = row['rating']

        # Play against random opponents
        opponents = teams[teams['team'] != row['team']].sample(n_games)

        # Simulate each game
        for _, opp in opponents.iterrows():
            # Random home/away
            is_home = np.random.random() > 0.5
            home_adv = 2.5 if is_home else -2.5

            game = simulate_game(team_rating, opp['rating'],
                               home_advantage=home_adv)

            if game['team_a_wins']:
                results.loc[results['team'] == row['team'], 'wins'] += 1

    results['losses'] = n_games - results['wins']
    return results.sort_values('wins', ascending=False)

# Run one simulation
np.random.seed(2024)
season_result = simulate_season(afc_teams)

season_result_display = season_result.merge(
    afc_teams[['team', 'division']],
    on='team'
)[['team', 'division', 'wins', 'losses', 'rating']].head(10)

print("\nSimulated Season Results (Top 10 teams - one simulation):")
print(season_result_display.to_string(index=False))
print("\nTop 7 teams make playoffs (highlighted)")

One simulation shows a plausible season outcome, with top teams winning 11-13 games and weak teams winning 4-7 games. Notice some randomness: a team with a mediocre rating might sneak into the top 7, or a strong team might underperform. This is football—variance matters. The question is: over thousands of simulated seasons, how often does each team make the playoffs?

Monte Carlo Season Simulations

Now comes the power of Monte Carlo simulation: we run thousands of seasons to estimate playoff probabilities. One simulated season showed the Chiefs winning 13 games, but what if they got unlucky with injuries or close games? What if their opponents over-performed? By simulating 10,000 seasons, we capture all these possibilities and estimate each team's probability of making the playoffs.

The process is straightforward: run simulate_season() 10,000 times, storing the results each time. Then for each team, calculate what percentage of simulations they finished in the top 7. This percentage is our playoff probability estimate. Teams with high ratings will make it most of the time; teams with low ratings rarely; teams in the middle have uncertain fates.

This approach naturally accounts for strength of schedule variations (some teams get luckier with random opponents), win distribution (teams can have the same expected wins but different variance), and interaction effects (if Team A beats Team B, both teams' playoff odds change). All the complexity is handled implicitly through simulation.

#| label: season-monte-carlo-r
#| message: false
#| warning: false
#| cache: true

# Run 10,000 season simulations
# This will take 20-30 seconds - we're simulating 272,000 games!
n_season_sims <- 10000

playoff_results <- map_dfr(1:n_season_sims, function(i) {

  # Simulate one complete season
  season <- simulate_season(afc_teams)

  # Top 7 teams make playoffs
  season %>%
    mutate(sim = i) %>%
    mutate(playoff = row_number() <= 7) %>%
    select(sim, team, wins, playoff)

})

# Calculate playoff probabilities by team
# Count how often each team made playoffs across all simulations
playoff_odds <- playoff_results %>%
  group_by(team) %>%
  summarise(
    playoff_pct = mean(playoff) * 100,  # % of sims where team made playoffs
    avg_wins = mean(wins),               # Average wins across all sims
    .groups = "drop"
  ) %>%
  left_join(afc_teams %>% select(team, rating, division), by = "team") %>%
  arrange(desc(playoff_pct))

playoff_odds %>%
  head(10) %>%
  gt() %>%
  cols_label(
    team = "Team",
    division = "Division",
    rating = "Rating",
    avg_wins = "Avg Wins",
    playoff_pct = "Playoff %"
  ) %>%
  fmt_number(
    columns = c(rating, avg_wins),
    decimals = 1
  ) %>%
  fmt_number(
    columns = playoff_pct,
    decimals = 1
  ) %>%
  tab_header(
    title = "Playoff Probabilities",
    subtitle = paste0("Based on ", scales::comma(n_season_sims), " season simulations")
  ) %>%
  data_color(
    columns = playoff_pct,
    colors = scales::col_numeric(
      palette = c("white", "#00BFC4"),
      domain = c(0, 100)
    )
  )

#| label: season-monte-carlo-py
#| message: false
#| warning: false
#| cache: true

# Run 10,000 season simulations
n_season_sims = 10000

playoff_results = []

for i in range(n_season_sims):
    # Simulate season
    season = simulate_season(afc_teams)

    # Top 7 teams make playoffs
    season['sim'] = i
    season['playoff'] = season.index < 7

    playoff_results.append(season[['sim', 'team', 'wins', 'playoff']])

playoff_results_df = pd.concat(playoff_results, ignore_index=True)

# Calculate playoff probabilities
playoff_odds = (playoff_results_df
    .groupby('team')
    .agg(
        playoff_pct=('playoff', lambda x: x.mean() * 100),
        avg_wins=('wins', 'mean')
    )
    .reset_index()
    .merge(afc_teams[['team', 'rating', 'division']], on='team')
    .sort_values('playoff_pct', ascending=False)
)

print(f"\nPlayoff Probabilities (Based on {n_season_sims:,} season simulations):")
print(playoff_odds.head(10).to_string(index=False))

The table shows playoff probabilities ranging from near-100% (Chiefs) to near-0% (bottom teams), with interesting uncertainty in the middle. This matches our intuition: great teams almost always make it, terrible teams almost never do, and average teams face genuine uncertainty. Monte Carlo simulation quantifies this uncertainty precisely.

Visualizing Playoff Probabilities

Numbers in a table are informative, but a visualization reveals patterns more clearly. By plotting playoff probability against team rating, we can see the relationship between team quality and playoff chances. This visualization serves multiple purposes: it validates our model (better teams should have higher odds), reveals the strength of the relationship (how steep is the curve?), and identifies outliers (teams over- or under-performing their rating).

We expect an S-shaped (logistic) curve: very strong teams approach 100% playoff probability, very weak teams approach 0%, and teams around average are most uncertain. The steepness of this curve tells us how deterministic NFL seasons are—a steep curve means ratings are very predictive, while a shallow curve means lots of randomness.

Let's visualize this relationship with playoff probability on the y-axis and team rating on the x-axis:

#| label: fig-playoff-odds-r
#| fig-cap: "Playoff probabilities by team rating"
#| fig-width: 10
#| fig-height: 8

playoff_odds %>%
  ggplot(aes(x = rating, y = playoff_pct)) +
  # Points colored by division
  geom_point(aes(color = division), size = 4, alpha = 0.7) +
  # Smooth trend line showing relationship
  geom_smooth(method = "loess", se = TRUE, color = "black",
              linetype = "dashed", linewidth = 0.8) +
  # Team labels offset to avoid overlap
  geom_text(aes(label = team), hjust = -0.2, size = 3) +
  # 50% line showing playoff cutoff
  geom_hline(yintercept = 50, linetype = "dotted",
             color = "gray40", alpha = 0.5) +
  scale_y_continuous(labels = scales::percent_format(scale = 1),
                    limits = c(0, 100)) +
  labs(
    title = "Playoff Probability vs Team Rating",
    subtitle = paste0("Based on ", scales::comma(n_season_sims),
                     " season simulations"),
    x = "Team Rating (points above/below average)",
    y = "Playoff Probability",
    color = "Division",
    caption = "Dashed line shows fitted relationship | Better teams have higher playoff odds"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(face = "bold"),
    legend.position = "right"
  )

#| label: fig-playoff-odds-py
#| fig-cap: "Playoff probabilities by team rating - Python"
#| fig-width: 10
#| fig-height: 8

from scipy.interpolate import make_interp_spline

plt.figure(figsize=(10, 8))

# Color by division
divisions = playoff_odds['division'].unique()
colors = plt.cm.Set2(range(len(divisions)))
division_colors = dict(zip(divisions, colors))

for division in divisions:
    div_data = playoff_odds[playoff_odds['division'] == division]
    plt.scatter(div_data['rating'], div_data['playoff_pct'],
               label=division, s=100, alpha=0.7,
               color=division_colors[division])

# Add team labels
for _, row in playoff_odds.iterrows():
    plt.annotate(row['team'],
                (row['rating'], row['playoff_pct']),
                xytext=(5, 0), textcoords='offset points',
                fontsize=8)

# Add smooth trend line
x_smooth = np.linspace(playoff_odds['rating'].min(),
                       playoff_odds['rating'].max(), 100)
from scipy.interpolate import UnivariateSpline
spl = UnivariateSpline(playoff_odds['rating'], playoff_odds['playoff_pct'], s=500)
plt.plot(x_smooth, spl(x_smooth), 'k--', alpha=0.5, linewidth=2,
         label='Fitted Trend')

# Add 50% reference line
plt.axhline(y=50, color='gray', linestyle=':', alpha=0.5)

plt.xlabel('Team Rating (points above/below average)', fontsize=12)
plt.ylabel('Playoff Probability (%)', fontsize=12)
plt.title(f'Playoff Probability vs Team Rating\nBased on {n_season_sims:,} season simulations',
          fontsize=14, fontweight='bold')
plt.ylim(0, 100)
plt.legend(title='Division', bbox_to_anchor=(1.05, 1), loc='upper left')
plt.tight_layout()
plt.show()

The plot makes clear that team quality matters enormously for playoff probability, but randomness keeps most teams in the 20-80% range rather than at the extremes. Only truly elite or terrible teams have playoff destinies that approach certainty.

Strategic Scenario Analysis

Fourth Down Decision Simulation

Monte Carlo simulation excels at evaluating strategic decisions with uncertain outcomes. Consider a fourth-and-2 at the opponent's 35-yard line. The coach faces three options: go for it, kick a field goal, or punt. Each decision has different probability distributions of outcomes, and the "right" choice depends on these distributions, not just expected values.

Traditional expected value analysis might say "going for it is worth +1.2 points," but that masks the variance. Going for it might yield 7 points (touchdown), 3 points (field goal after conversion), or -2 points (turnover on downs). The field goal attempt might yield 3 points (make) or -1.5 points (miss, opponent gets good field position). Punting might yield -0.5 points (opponent starts backed up but gets the ball).

Monte Carlo simulation lets us see the full distribution of outcomes, not just the expected value. This is crucial because risk matters. A conservative coach might prefer the 3-point field goal with low variance over the higher-EV but volatile go-for-it decision. A coach down 10 points in the fourth quarter needs the touchdown upside. Simulation informs these context-dependent choices.

Let's model this decision with realistic probabilities:

#| label: fourth-down-sim-r
#| message: false
#| warning: false

# Scenario: 4th and 2 at opponent's 35-yard line
# Options: Go for it, Kick FG, Punt

simulate_fourth_down <- function(decision, n_sims = 10000) {

  results <- numeric(n_sims)

  for(i in 1:n_sims) {

    if(decision == "go") {
      # Going for it: 50% conversion rate on 4th and 2
      if(runif(1) < 0.50) {
        # Converted! Now simulate rest of drive
        # 40% TD, 20% FG, 40% turnover/punt from deeper in redzone
        outcome <- sample(c("td", "fg", "turnover"), 1,
                         prob = c(0.4, 0.2, 0.4))
        results[i] <- case_when(
          outcome == "td" ~ 7,        # Touchdown
          outcome == "fg" ~ 3,        # Field goal
          outcome == "turnover" ~ -2, # Turnover, opponent gets ball
          TRUE ~ 0
        )
      } else {
        # Failed conversion - turnover on downs at 35 yard line
        # Opponent has good field position, expected to score ~2 points
        results[i] <- -2
      }

    } else if(decision == "kick_fg") {
      # 52-yard field goal attempt (~70% success from 35 yard line)
      if(runif(1) < 0.70) {
        results[i] <- 3  # Made field goal
      } else {
        # Missed FG, opponent gets ball at 35 or touchback
        results[i] <- -1.5
      }

    } else if(decision == "punt") {
      # Punt from 35 typically nets to ~opponent 10 yard line
      # Opponent starts backed up, expected points ~-0.5
      results[i] <- -0.5
    }
  }

  tibble(
    decision = decision,
    expected_points = mean(results),
    sd_points = sd(results),
    p10 = quantile(results, 0.1),
    p90 = quantile(results, 0.9),
    values = list(results)
  )
}

# Simulate all three decisions
fourth_down_results <- bind_rows(
  simulate_fourth_down("go"),
  simulate_fourth_down("kick_fg"),
  simulate_fourth_down("punt")
)

fourth_down_results %>%
  select(decision, expected_points, sd_points, p10, p90) %>%
  gt() %>%
  cols_label(
    decision = "Decision",
    expected_points = "Expected Points",
    sd_points = "Std Dev",
    p10 = "10th %ile",
    p90 = "90th %ile"
  ) %>%
  fmt_number(
    columns = c(expected_points, sd_points, p10, p90),
    decimals = 2
  ) %>%
  tab_header(
    title = "4th and 2 at Opponent's 35",
    subtitle = "Expected value of each decision"
  ) %>%
  data_color(
    columns = expected_points,
    colors = scales::col_numeric(
      palette = c("#F8766D", "white", "#00BFC4"),
      domain = c(-1, 3)
    )
  )

#| label: fourth-down-sim-py
#| message: false
#| warning: false

def simulate_fourth_down(decision, n_sims=10000):
    """Simulate fourth down decision outcomes."""

    results = np.zeros(n_sims)

    for i in range(n_sims):

        if decision == "go":
            # 50% conversion rate
            if np.random.random() < 0.50:
                # Converted - simulate drive outcome
                outcome = np.random.choice(
                    ["td", "fg", "turnover"],
                    p=[0.4, 0.2, 0.4]
                )
                if outcome == "td":
                    results[i] = 7
                elif outcome == "fg":
                    results[i] = 3
                else:  # turnover
                    results[i] = -2
            else:
                # Failed conversion
                results[i] = -2

        elif decision == "kick_fg":
            # 70% success rate from 52 yards
            if np.random.random() < 0.70:
                results[i] = 3
            else:
                results[i] = -1.5

        elif decision == "punt":
            # Punt outcome
            results[i] = -0.5

    return {
        'decision': decision,
        'expected_points': results.mean(),
        'sd_points': results.std(),
        'p10': np.percentile(results, 10),
        'p90': np.percentile(results, 90),
        'values': results
    }

# Simulate all three decisions
fourth_down_results = pd.DataFrame([
    simulate_fourth_down("go"),
    simulate_fourth_down("kick_fg"),
    simulate_fourth_down("punt")
])

print("\n4th and 2 at Opponent's 35 - Expected value of each decision:")
print(fourth_down_results[['decision', 'expected_points', 'sd_points', 'p10', 'p90']]
      .to_string(index=False))

The analysis shows that while going for it and kicking the field goal have similar expected values, they have very different risk profiles. The field goal is safer but capped at 3 points; going for it is riskier but offers touchdown upside. Your choice depends on risk tolerance, game situation, and team strengths.

Visualizing Decision Distributions

The summary statistics tell us expected values and variance, but visualizing the full distribution reveals patterns that numbers alone miss. By plotting probability densities for each decision, we can see how often each outcome occurs, where the mass of the distribution sits, and how much overlap exists between options.

The visualization will show three overlapping distributions: one for each fourth-down decision. We expect "go for it" to be bimodal (peaks at -2 for failures and +7 for touchdowns), "kick FG" to have most mass at +3 with a smaller peak at -1.5, and "punt" to be concentrated around -0.5. The overlap tells us when the choice matters—if all three distributions were identical, the decision wouldn't matter at all.

#| label: fig-fourth-down-dist-r
#| fig-cap: "Distribution of outcomes for 4th down decisions"
#| fig-width: 10
#| fig-height: 6

fourth_down_results %>%
  select(decision, values) %>%
  unnest(values) %>%
  ggplot(aes(x = values, fill = decision)) +
  # Density curves show probability of each outcome
  geom_density(alpha = 0.6, linewidth = 1) +
  # Zero line shows breakeven point
  geom_vline(xintercept = 0, linetype = "dashed", color = "black") +
  scale_fill_manual(
    values = c("go" = "#00BFC4", "kick_fg" = "#F8766D", "punt" = "#7CAE00"),
    labels = c("Go For It", "Kick FG", "Punt")
  ) +
  labs(
    title = "4th and 2 at Opponent's 35 Yard Line",
    subtitle = "Distribution of expected points for each decision",
    x = "Expected Points Added",
    y = "Probability Density",
    fill = "Decision"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(face = "bold"),
    legend.position = "top"
  )

#| label: fig-fourth-down-dist-py
#| fig-cap: "Distribution of outcomes for 4th down decisions - Python"
#| fig-width: 10
#| fig-height: 6

plt.figure(figsize=(10, 6))

colors = {'go': '#00BFC4', 'kick_fg': '#F8766D', 'punt': '#7CAE00'}
labels = {'go': 'Go For It', 'kick_fg': 'Kick FG', 'punt': 'Punt'}

for _, row in fourth_down_results.iterrows():
    decision = row['decision']
    # Use KDE for smooth density
    from scipy.stats import gaussian_kde
    density = gaussian_kde(row['values'])
    xs = np.linspace(row['values'].min(), row['values'].max(), 200)
    plt.plot(xs, density(xs), linewidth=2, alpha=0.8,
             color=colors[decision], label=labels[decision])
    plt.fill_between(xs, density(xs), alpha=0.3, color=colors[decision])

plt.axvline(x=0, color='black', linestyle='--', linewidth=2)
plt.xlabel('Expected Points Added', fontsize=12)
plt.ylabel('Probability Density', fontsize=12)
plt.title('4th and 2 at Opponent\'s 35 Yard Line\nDistribution of expected points for each decision',
          fontsize=14, fontweight='bold')
plt.legend(title='Decision', loc='upper right')
plt.tight_layout()
plt.show()

The visualization makes the risk-reward tradeoff crystal clear: going for it offers the highest upside (touchdown peak near +7) but also more downside risk, while the field goal is safer with less variance but capped upside. The punt is the most predictable but gives up points on average.

Drive Outcome Simulation

Modeling Drive Results

Moving from single decisions to complete drives requires simulating sequences of plays. A drive simulation tracks field position, down-and-distance, and play outcomes over multiple plays until the drive ends (touchdown, field goal, turnover, or turnover on downs). This is more complex than game simulation but reveals insights about offensive efficiency, play-calling strategy, and scoring probability by field position.

The key components of a drive simulator are:

1. Play-by-play modeling: Each play advances the ball some yardage (positive or negative) based on play type and randomness
2. Down-and-distance tracking: After each play, update downs and yards to go for first down
3. Drive-ending events: Model touchdowns, turnovers, turnovers on downs, and field position for field goal attempts
4. Play-type selection: Balance between running and passing, potentially adjusting based on down and distance

We'll build a simplified drive simulator that captures the essential dynamics without play-by-play granularity. The model will randomly choose run or pass, simulate yardage gained from distributions, and track drive progress until an ending event occurs:

#| label: drive-sim-r
#| message: false
#| warning: false

# Simulate a complete drive starting at given yard line
simulate_drive <- function(start_yard = 25, offense_rating = 0) {

  yard_line <- start_yard
  down <- 1
  yards_to_go <- 10
  plays <- 0
  max_plays <- 20  # Safety valve to prevent infinite loops

  while(plays < max_plays) {
    plays <- plays + 1

    # Determine play type (60% pass, 40% run in modern NFL)
    play_type <- sample(c("pass", "run"), 1, prob = c(0.6, 0.4))

    # Simulate yards gained
    # Pass: higher mean, higher variance (more boom/bust)
    # Run: lower mean, lower variance (more consistent)
    if(play_type == "pass") {
      yards <- round(rnorm(1, mean = 7 + offense_rating, sd = 10))
    } else {
      yards <- round(rnorm(1, mean = 4 + offense_rating * 0.5, sd = 5))
    }

    # Cap yards at reasonable bounds (-15 to +75)
    yards <- max(-15, min(yards, 75))

    # 5% chance of turnover (interception or fumble)
    if(runif(1) < 0.05) {
      return(list(result = "turnover", plays = plays, yards = yard_line))
    }

    # Update field position and yards to go
    yard_line <- yard_line + yards
    yards_to_go <- yards_to_go - yards

    # Check for touchdown
    if(yard_line >= 100) {
      return(list(result = "touchdown", plays = plays, yards = 100))
    }

    # Check for safety (backed into own endzone)
    if(yard_line <= 0) {
      return(list(result = "safety", plays = plays, yards = 0))
    }

    # Check for first down
    if(yards_to_go <= 0) {
      down <- 1
      yards_to_go <- 10
    } else {
      down <- down + 1

      # Turnover on downs after 4th down
      if(down > 4) {
        return(list(result = "turnover_downs", plays = plays,
                   yards = yard_line))
      }
    }
  }

  # Max plays reached (shouldn't happen often)
  return(list(result = "timeout", plays = plays, yards = yard_line))
}

# Simulate 10,000 drives to see distribution of outcomes
n_drive_sims <- 10000

drive_results <- map_dfr(1:n_drive_sims, function(i) {
  drive <- simulate_drive(start_yard = 25, offense_rating = 1)
  tibble(
    sim = i,
    result = drive$result,
    plays = drive$plays,
    yards = drive$yards
  )
})

# Summarize results
drive_summary <- drive_results %>%
  group_by(result) %>%
  summarise(
    count = n(),
    pct = n() / n_drive_sims * 100,
    avg_plays = mean(plays),
    .groups = "drop"
  ) %>%
  arrange(desc(count))

drive_summary %>%
  gt() %>%
  cols_label(
    result = "Drive Result",
    count = "Count",
    pct = "Percentage",
    avg_plays = "Avg Plays"
  ) %>%
  fmt_number(
    columns = pct,
    decimals = 1
  ) %>%
  fmt_number(
    columns = avg_plays,
    decimals = 1
  ) %>%
  tab_header(
    title = "Drive Outcome Simulation",
    subtitle = paste0("Based on ", scales::comma(n_drive_sims), " simulated drives")
  ) %>%
  data_color(
    columns = pct,
    colors = scales::col_numeric(
      palette = c("white", "#619CFF"),
      domain = c(0, max(drive_summary$pct))
    )
  )

#| label: drive-sim-py
#| message: false
#| warning: false

def simulate_drive(start_yard=25, offense_rating=0):
    """Simulate a complete drive."""

    yard_line = start_yard
    down = 1
    yards_to_go = 10
    plays = 0
    max_plays = 20

    while plays < max_plays:
        plays += 1

        # Determine play type (60% pass, 40% run)
        play_type = np.random.choice(['pass', 'run'], p=[0.6, 0.4])

        # Simulate yards gained
        if play_type == 'pass':
            yards = round(np.random.normal(7 + offense_rating, 10))
        else:
            yards = round(np.random.normal(4 + offense_rating * 0.5, 5))

        # Cap yards
        yards = max(-15, min(yards, 75))

        # Check for turnover
        if np.random.random() < 0.05:
            return {'result': 'turnover', 'plays': plays, 'yards': yard_line}

        # Update field position
        yard_line += yards
        yards_to_go -= yards

        # Check for touchdown
        if yard_line >= 100:
            return {'result': 'touchdown', 'plays': plays, 'yards': 100}

        # Check for safety
        if yard_line <= 0:
            return {'result': 'safety', 'plays': plays, 'yards': 0}

        # Check for first down
        if yards_to_go <= 0:
            down = 1
            yards_to_go = 10
        else:
            down += 1
            if down > 4:
                return {'result': 'turnover_downs', 'plays': plays,
                       'yards': yard_line}

    return {'result': 'timeout', 'plays': plays, 'yards': yard_line}

# Simulate 10,000 drives
n_drive_sims = 10000

drive_results = pd.DataFrame([
    simulate_drive(start_yard=25, offense_rating=1)
    for _ in range(n_drive_sims)
])

# Summarize results
drive_summary = (drive_results
    .groupby('result')
    .agg(
        count=('result', 'count'),
        avg_plays=('plays', 'mean')
    )
    .reset_index()
)
drive_summary['pct'] = drive_summary['count'] / n_drive_sims * 100
drive_summary = drive_summary.sort_values('count', ascending=False)

print(f"\nDrive Outcome Simulation (Based on {n_drive_sims:,} simulated drives):")
print(drive_summary.to_string(index=False))

The results show a realistic distribution of drive outcomes, with scoring rates that match NFL data for above-average offenses. This type of simulation can be extended to model entire games (alternating drives for each team), inform fourth-down decisions (by simulating drive continuation probability), or evaluate play-calling strategies.

Bootstrapping and Resampling Methods

Understanding Bootstrap Confidence Intervals

Monte Carlo simulation generates data from a model, but what if we want to quantify uncertainty in real data? Bootstrapping provides the answer. This resampling technique estimates confidence intervals and standard errors without making distributional assumptions, making it perfect for football analytics where data often violates normality assumptions.

The bootstrap works by repeatedly resampling your observed data with replacement, calculating your statistic of interest on each resample, and using the distribution of those statistics to estimate uncertainty. If you want a confidence interval for average EPA per play, you resample your plays thousands of times and see how the average varies across resamples.

Why does this work? The bootstrap principle states that the relationship between your sample and the resampled data approximates the relationship between the true population and your sample. In other words, the variability you see in bootstrap resamples mimics the variability you'd see if you could collect new samples from the population.

Let's use bootstrapping to estimate confidence intervals for team offensive efficiency:

#| label: bootstrap-r
#| message: false
#| warning: false
#| cache: true

# Load real play-by-play data
pbp_2023 <- load_pbp(2023)

# Calculate team EPA per play
team_epa <- pbp_2023 %>%
  filter(season_type == "REG", !is.na(epa), play_type %in% c("pass", "run")) %>%
  group_by(posteam) %>%
  summarise(
    n_plays = n(),
    observed_epa = mean(epa),
    .groups = "drop"
  )

# Bootstrap function: resample plays and recalculate EPA
bootstrap_team_epa <- function(team_abbr, pbp_data, n_boot = 1000) {

  # Get plays for this team
  team_plays <- pbp_data %>%
    filter(posteam == team_abbr, !is.na(epa), play_type %in% c("pass", "run"))

  # Bootstrap: resample plays with replacement n_boot times
  boot_estimates <- map_dbl(1:n_boot, function(i) {
    # Sample n plays with replacement
    boot_sample <- team_plays %>%
      slice_sample(n = nrow(team_plays), replace = TRUE)

    # Calculate EPA for this bootstrap sample
    mean(boot_sample$epa)
  })

  # Calculate confidence intervals from bootstrap distribution
  tibble(
    team = team_abbr,
    observed_epa = mean(team_plays$epa),
    boot_mean = mean(boot_estimates),
    boot_se = sd(boot_estimates),
    ci_lower = quantile(boot_estimates, 0.025),  # 2.5th percentile
    ci_upper = quantile(boot_estimates, 0.975)   # 97.5th percentile
  )
}

# Run bootstrap for top 6 teams
top_teams <- team_epa %>%
  arrange(desc(observed_epa)) %>%
  head(6) %>%
  pull(posteam)

boot_results <- map_dfr(top_teams,
                        ~bootstrap_team_epa(.x, pbp_2023, n_boot = 1000))

boot_results %>%
  gt() %>%
  cols_label(
    team = "Team",
    observed_epa = "Observed EPA",
    boot_mean = "Bootstrap Mean",
    boot_se = "Std Error",
    ci_lower = "95% CI Lower",
    ci_upper = "95% CI Upper"
  ) %>%
  fmt_number(
    columns = c(observed_epa, boot_mean, boot_se, ci_lower, ci_upper),
    decimals = 3
  ) %>%
  tab_header(
    title = "Bootstrap Confidence Intervals for Team EPA",
    subtitle = "Based on 1,000 bootstrap resamples of 2023 regular season"
  )

#| label: bootstrap-py
#| message: false
#| warning: false
#| cache: true

# Would load real data here - using simulated for demonstration
# pbp_2023 = nfl.import_pbp_data([2023])

# Simulated team data for demonstration
np.random.seed(2024)

def bootstrap_team_epa(team_epa_values, n_boot=1000):
    """Bootstrap confidence intervals for team EPA."""

    boot_estimates = []

    for i in range(n_boot):
        # Resample with replacement
        boot_sample = np.random.choice(team_epa_values,
                                      size=len(team_epa_values),
                                      replace=True)
        boot_estimates.append(boot_sample.mean())

    boot_estimates = np.array(boot_estimates)

    return {
        'observed_epa': team_epa_values.mean(),
        'boot_mean': boot_estimates.mean(),
        'boot_se': boot_estimates.std(),
        'ci_lower': np.percentile(boot_estimates, 2.5),
        'ci_upper': np.percentile(boot_estimates, 97.5)
    }

# Simulate example data for demonstration
print("\nBootstrap Confidence Intervals for Team EPA")
print("(Based on 1,000 bootstrap resamples)")
print("\nExample: Team with 1000 plays, mean EPA = 0.15")

# Simulated data
team_plays = np.random.normal(0.15, 2.0, 1000)
result = bootstrap_team_epa(team_plays, n_boot=1000)

print(f"Observed EPA: {result['observed_epa']:.3f}")
print(f"Bootstrap Mean: {result['boot_mean']:.3f}")
print(f"Standard Error: {result['boot_se']:.3f}")
print(f"95% CI: [{result['ci_lower']:.3f}, {result['ci_upper']:.3f}]")