Exploring the difficulty of estimating win probability: a simulation study

2.1 Introducing random walk football

We begin by describing the rules of random walk football. Random walk football begins at midfield (yardline L/2, where L is an even integer). Each play, the ball moves left or right by one yardline with equal probability. If the ball reaches the left end of the field (yardline 0), team one scores a touchdown, worth +1 point. If the ball reaches the right end of the field (yardline L), team two scores a touchdown, worth −1 point. The ball resets to midfield after each touchdown. After T plays, the game ends. If the game is still tied after T plays, a fair coin is flipped to determine the winner.

Formally, the outcome of the tth play of the gth game is

The game starts at midfield, X_g0 = L/2, and the game begins tied, S_g0 = 0. The field position at the start of play t is

and the score differential at the start of play t is

The binary win/loss response column is

The true win probability

of random walk football is computed explicitly using dynamic programming,

and

2.2 Random walk football play-by-play data

A simulated observational dataset of random walk football plays consists of G games, each with T plays per game. Such a dataset has the form

For each play of game g, we record the timestep t, the field position X _gt , the score differential S _gt , and a binary variable y _gt indicating whether the team with possession wins the game. The dependence structure of random walk football (and real football) play-by-play data manifests in the outcome vector {y _gt }: plays from the same game share the same draw of the winner of the game. Formally, y _gt ≡ y_g,T+1 as in Equation (4). Importantly, we also know the underlying true win probability W P g t at each play, which we use to evaluate win probability estimates.

Throughout this study, we let T = 56, the average number of first-down plays per game in the dataset of American football plays from Brill et al. (2025). We use L = 4 yardlines so that the average number of plays between each score is similar to the average number of first-down plays in a game of American football.

We want to assess the impact of the dependence structure of the response variable of observational football data on the accuracy of a statistical win probability estimator. To do so, we compare the accuracy of a W P estimator fit from datasets of varying degrees of dependence. We introduce a parameter K that controls the degree of dependence of the win/loss outcome variable in a generated dataset: we keep a random subsample of K plays per generated game. K = 1 – keep just 1 randomly sampled play in each of the G simulated games – reflects independence across all plays in the filtered dataset because each game is generated independently. When K = 1, the outcome variable y _gt reflects an independent draw of the win/loss outcome of the game since each play is filtered from a separate game. K = T – keep all T plays in each of the G simulated games – reflects full dependence within each game and is equivalent to the original dataset. When K = T, the outcome variable y _gt for each play t in game g reflects the same draw of the win/loss outcome of the game. Integer values of K between 1 and T reflect intermediate degrees of dependence because just 1 < K < T plays per game share the same draw of the response variable. For concreteness, we visualize example datasets for K = 1, K = 3, and K = T in Table 4.

Table 4:

Visualizing example generated datasets with varying degrees of dependence K. Each dataset has the same nominal sample size, ζ ⋅ T plays (rows). The variables are: field position x, time t, score differential s, binary win/loss y. In the K = 1 dataset (a), all plays are independent because they are generated from separate independent games. In the K = 3 dataset (b), groups of 3 plays are generated from the same game. Each group of 3 plays shares the same draw of the response y and plays from different games are independent. The K = T dataset (c) includes all T plays from each generated game. Within each game, all plays share the same draw of the response y, and plays from different games are independent.

Generating a dataset with G games and T = 56 plays per game, keeping just K randomly sampled plays per game, has a nominal sample size (number of rows) of G ⋅ K plays. We want to compare the performance of a win probability estimator fit from datasets of varying degrees of dependence K that have the same nominal sample size. To do so, given K, we generate G = round(ζ ⋅ T/K) games and keep K plays per game. This yields a dataset consisting of (approximately) ζ ⋅ T plays, which is independent of K. As T = 56 throughout this paper, ζ parameterizes the sample size.

2.3 Ways to estimate win probability

From a dataset of football plays, we want to estimate win probability. Win probability is not a counting statistic or an observable quantity. Rather, it is defined by a model that is estimated from data. Estimating win probability is one of the classic modeling tasks of sports analytics. Win probability estimates arise broadly from one of three classes of models: simple mathematical models, probabilistic state-space models, and statistical models. We give an overview of these classes below. For a more thorough review we refer to the reader to Baumer et al. (2023).

Stern (1994), for example, uses a simple mathematical model to estimate win probability in basketball. Supposing possession-level score differential outcomes are independent and approximately Gaussian, he models the score differential process by a Brownian motion with drift μ points advantage for the home team and variance σ². He uses probit regression to estimate p(l, t), the probability the home team wins if they are leading by l points after t seconds of game time, p ( l , t ) = Φ ( l + ( 1 − t ) μ ) / ( 1 − t ) σ 2 . The benefit of such mathematical models is their simplicity: they have closed-form solutions. The drawback is they rely on unreasonable assumptions (e.g., normality), which fail in particular towards the end of a game and work for a limited set of game-state variables (e.g., just score differential and time).

State-space models simplify a game into a series of transitions between game-states. Transition probabilities are estimated from play-level data and are then propagated into win probability by simulating games. Win probability in baseball is commonly estimated using state-space models going back to Lindsey (1961). These models work well in baseball because the game consists of discrete events (i.e., the game is divided into 9 innings, each of which feature a sequence of individual pitcher-batter matchups) and there are just a few important game-state variables (e.g., base-state, outs, runs, and inning). When implemented correctly, state-space models are sensible ways to estimate W P . However, they are difficult in practice, as they require: a careful encoding of the convoluted rules of a sport into a set of states and the actions between those states, careful estimation of transition probabilities, and enough computing power to run enough simulated games to achieve desired granularity. Each of these can be nontrivial depending on the complexity of the sport.

Finally, statistical models are fit entirely from historical data. Given the results of a set of observed plays, statistical models fit the relationship between a binary win/loss outcome variable and certain game-state variables using regression or machine learning approaches. Notably, these models are widely used today by analysts of American football (Baldwin 2021; Lock and Nettleton 2014). They form the foundation of open-source fourth-down recommendations (Brill et al. 2025). These models are popular due to the accessibility of rich publicly available data (e.g., play-by-play data from nflFastR Carl and Baldwin 2022) and off-the-shelf machine learning tools (e.g., X G B o o s t Chen and Guestrin 2016). They are widely assumed to be reasonable and trustworthy because of the modeling flexibility of machine learning algorithms. For these reasons, in this paper we focus on statistical win probability estimators. The binary win/loss outcome variable, however, is noisy and features a strong dependence structure–each play in the same game shares the same draw of the win/loss outcome. Accordingly, Brill et al. (2025) found that these estimators are subject to substantial uncertainty and produce wide confidence intervals, even when fit from a large dataset featuring 229,635 first-down plays across 4,101 games and 16 years. We suspect statistical W P estimators have high variance, exacerbated by the dependence structure of observational play-by-play data.

Continuing the tradition of Baldwin (2021), throughout this paper we estimate win probability using X G B o o s t . The covariates are x = (t, x, s), where t denotes time, x denotes field position, and s denotes score differential. The outcome variable is binary win/loss y. We use half of the games from the training set as a validation set to tune X G B o o s t models.

2.4 Bias-variance decomposition

In this section, we analyze the bias-variance decomposition of an X G B o o s t win probability estimator. We compare a W P estimator fit from datasets having the same nominal sample size ζ ⋅ T but generated with varying degrees of dependence K. Given a combination of data generating parameters, we generate M = 100 training datasets. We fit a win probability estimator from each dataset, { W P ̂ ( m ) } m = 1 M . We also generate M = 100 out-of-sample testing datasets D test ( m ) m = 1 M using (G = 10,000, T = 56, K = 1). Then, we calculate the squared bias of the mth estimator by

and the variance by

The root mean squared error is R M S E m = b i a s m 2 + v a r m . We then calculate the average squared bias, variance, and RMSE across the M simulations and their standard errors.

First, let ζ = 4,101 to mimic the dataset of real American football plays. We compare the accuracy of a W P estimator fit from a (G = round(ζ ⋅ T/K), T = 56, K) dataset as the degree of dependence K varies. The sample size in each of these datasets is (approximately) the same, G ⋅ K = ζ ⋅ T. We visualize this bias-variance decomposition as K varies in Figure 1. As the strength K of the correlation increases, accuracy decreases linearly. Fixing the number of observations in the dataset but increasing the degree of dependence across outcomes reduces model accuracy. Intuitively, this makes sense because the degree of dependence K is inversely proportional to the amount of available independent data – there are (approximately) ζ ⋅ T/K independent draws of the response variable. Less independent data produces less accurate estimators.

Figure 1:

Squared bias (left), variance (middle), and RMSE of a win probability estimator fit from a (G = round(ζ ⋅ T/K), T = 56, K) dataset as a function of K, where ζ = 4,101. The dots denote the average values across M = 100 simulations and the bars denote plus/minus twice the standard errors. The gray line is the regression line.

Next, as a function of sample size ζ ⋅ T (with T = 56), we compare the accuracy of a W P estimator fit from three types of datasets. First, we consider a (G = ζ, K = T) dataset, which keeps each generated play per game. This dataset mimics the historical dataset of real American football plays. Then, we consider a (G = ζ, K = 1) dataset, which keeps just one randomly sampled play per game. This dataset consists entirely of independent outcomes, and can be formed wholly from a (G = ζ, K = T) dataset, but its sample size is much smaller (ζ rather than ζ ⋅ T). Finally, we consider a (G = ζ ⋅ T, K = 1) dataset, which has the same sample size (number of rows) ζ ⋅ T as the first dataset, but consists entirely of independent outcomes.

In Figure 2 we visualize the bias-variance decomposition of a win probability estimator fit from these three types of datasets as a function of ζ. The x-axis is log₄(ζ) because 4⁶ = 4,096 ≈ 4,101 is the sample size of the historical dataset of American football plays. We see that it is much better to use all plays per game rather than one independent play per game. Despite the strong dependence structure, keeping all the plays elucidates information about the structure of the covariate space. We also see that it would be much better if the plays had independent outcomes. This suggests that the dependence structure reduces the effective sample size of the dataset. We explore the extent of this reduction in the next section.

Figure 2:

Squared bias (left), variance (middle), and RMSE (right) of a win probability estimator fit from three datasets as a function of ζ. The (G = ζ, K = T) dataset (blue) involves keeping each generated play per game, which mimics the historical dataset of real American football plays. The (G = ζ, K = 1) dataset (red) is formed by keeping just 1 play per game. The (G = ζ ⋅ T, K = 1) dataset (orange) has the same sample size (number of rows) ζ ⋅ T as the first dataset but consists entirely of independent outcomes.

Interestingly, as shown in Figure 3, we see that the bias is much worse at the beginning of the game. Game-states with larger score differentials in the early game occur less frequently than other game-states. Also, intuitively it is easier to tie later game-states to the ultimate win/loss outcome, which is determined at the end of the game.

Figure 3:

Bias (y-axis) versus time n (x-axis) and score differential s (color) at midfield (field position x = 2) of win probability estimated from a (G = 4,101, T = 56, K = T) dataset. The lines denote the average values across M = 100 simulations and the shaded regions denote plus-minus twice the standard errors.

2.5 Effective sample size

As discussed in the previous section, we can calculate the accuracy of a win probability estimator fit from a (G = ζ, T = 56, K) dataset, denoted RMSE(ζ, K). The sample size (number of rows) of such a dataset with ζ = 4,101 and K = T, which mimics the historical dataset of American football plays, is ζ ⋅ T. We saw that the dependence structure of this dataset reduces the accuracy of our estimator, but we would like to understand the extent of this reduction. In particular, we are interested in the effective sample size (ESS) of that dataset. The ESS is the sample size ζ′ ⋅ T of a (G = ζ′ ⋅ T, T = 56, K = 1) dataset consisting of independent outcomes, which produces an estimator having the same accuracy as one fit from the original dataset. For brevity, we refer to the sample size as ζ and the effective sample size as ζ′, dropping the T since we use T = 56 throughout this study.

To estimate this effective sample size, we begin by fitting the K = 1 and K = T accuracy curves ζ ↦ RMSE(ζ, K) from Figure 2b. For each curve, we fit a biexponential model using nonlinear least squares. Then, as a function of ζ, the ESS is the value ζ′ satisfying RMSE(ζ′, K = 1) = RMSE(ζ, K = T). In Figure 4 we visualize the ESS ζ′ as a function of ζ. The ESS of a (ζ = 4,101, K = T) dataset is ζ′ = 2,291, or 56 % of the nominal sample size. This result is striking: we estimate that the historical dataset of American football plays (where ζ = 4,101) consists of about half as much data as suggested by the number of plays. In other words, we are effectively fitting win probability models from 8 years, not 16 years, worth of independent win/loss outcomes. Real American football is exponentially more complex than random walk football. Its game-state space is much larger, so we expect the ESS to be even smaller in real life.

Figure 4:

The effective sample size ζ′ of a (G = ζ, T = 56, K = T) dataset as a function of ζ. The red dot denotes ζ = 4,101, the number of games in the historical dataset of real American football plays.

If we halved the size of our K = T dataset, fitting a win probability model from ζ = 2,050 games (8 seasons), we estimate the effective sample size is ζ′ = 645, or just 31 % of the nominal sample size. This mimics fitting a win probability model from just recent data. If we doubled the size of our K = T dataset, fitting a win probability model from ζ = 8,202 games (32 seasons), we estimate the effective sample size is ζ′ = 6,911, or 84 % of the nominal sample size.

2.6 Coverage of bootstrapped confidence intervals

We have seen that machine learning win probability estimators fit from noisy and highly correlated observational data have high variance. Due to the dependence structure of historical football data, the effective sample size is much smaller than the nominal sample size. Therefore, we want to quantify uncertainty in win probability point estimates. The point estimates alone may not be trustworthy. The bootstrap is a natural choice to capture such uncertainty since it is non-parametric and does not make strong assumptions. Hence, in this section we explore the efficacy of bootstrapped win probability confidence intervals.

We begin with the standard (i.i.d.) bootstrap, which assumes each row (play) of the dataset is independently drawn. In the standard bootstrap, each of B bootstrapped datasets are formed by re-sampling G ⋅ T plays uniformly with replacement (recall G is the number of games, T is the number of plays per game, and G ⋅ T is the total number of plays in a random walk football observational dataset). The assumptions of the standard bootstrap do not apply to observational play-by-play data due to its dependence structure. Hence, we also try the cluster bootstrap, in which each of B bootstrapped datasets are formed by re-sampling G games uniformly with replacement, keeping each observed play within each re-sampled game. Finally, in the randomized cluster bootstrap, each of B bootstrapped datasets are formed by re-sampling G games uniformly with replacement, and within each game re-sampling T plays uniformly with replacement.

Each type of bootstrap produces B bootstrapped datasets D train ( b ) b = 1 B from the training dataset D train . We then fit a win probability model to each bootstrapped dataset using X G B o o s t , { W P ̂ b } b = 1 B . From these, we form a 90 % confidence interval for W P ( x ) at game-state x by the 5th and 95th quantiles of { W P ̂ b ( x ) } b = 1 B . Letting B = 101 in this section, we form a 90 % confidence interval by [ W P ̂ ( 6 ) ( x ) , W P ̂ ( 96 ) ( x ) ] . To avoid substantially low coverage near the extremes ( W P ( x ) ≈ 0 or W P ( x ) ≈ 1 ), we widen our confidence intervals when W P ̂ ( x ) < 0.025 to have a lower bound of 0 and when W P ̂ ( x ) > 0.975 to have an upper bound of 1.

We evaluate the efficacy of these intervals by their coverage and width. For each type of bootstrap (standard bootstrap, cluster bootstrap, and randomized cluster bootstrap) and each simulation m ∈ {1, …, M = 100}, we compute the pointwise marginal coverage of bootstrapped confidence intervals,

This is the proportion of plays in mth held-out dataset whose true win probability lies inside the confidence interval. We also compute the mean width,

For each play in the mth held-out dataset, we calculate the width of the confidence interval, and then take the average across all the plays. We report the average and the standard error of these values { c o v e r a g e m } m = 1 M and { w i d t h m } m = 1 M in Table 5. To mimic the historical dataset of American football plays, each simulated dataset here consists of G = 4,101 games, T = 56 plays per game, and K = T plays per game that share the same outcome.

Even in our simplified setting of random walk football, each of these bootstrapped win probability confidence intervals are undercovered. Intuitively, naive bootstraps produce intervals that are too narrow because they involve resampling from observed data – which entails re-using observations without generating new ones – thus exploring a strictly smaller subspace of the (x, y) space than the true sampling distribution. In other words, the bootstrapped resampling distribution is a rough approximation of the true sampling distribution. The naive standard bootstrap in particular achieves dismally low marginal coverage. Even the randomized cluster bootstrap that accounts for the dependence structure does not achieve high enough coverage. We suspect this coverage issue would be even worse for real American football, which is exponentially more complex than random walk football.

2.7 The fractional bootstrap

Naive bootstrapped confidence intervals are not wide enough. A natural question arises: how wide do confidence intervals need to be so that nominally 90 % intervals actually achieve 90 % marginal coverage? Hence, in this section we explore alternative forms of the bootstrap to increase coverage.

The traditional method of tuning non-parametric bootstrapped confidence intervals is to calibrate the bootstrapped quantiles (DiCiccio and Efron 1996). For instance, instead of using the α/2th and (1 − α/2)th quantiles of { W P ̂ b ( x ) } b = 1 B to form a 1 − α confidence interval, use the β/2th and (1 − β/2)th quantiles for some β < α. In order for this traditional calibration method to work, B would have to be much larger than 101, likely an order of magnitude larger (e.g., B = 1,001). We, however, prefer to use lower values of B (e.g., closer to 101) for several reasons. It is much better to keep B small for applications that require evaluating bootstrapped predictions in real time. For example, a bootstrapped fourth-down decision recommendation from Brill et al. (2025) takes about 15 s when B = 101 and about 2.5 minutes when B = 1,001. The former can be run before a fourth-down play begins and the latter takes far too long. Additionally, storing 1,001 machine learning models is much more expensive than storing 101 of them. For these reasons, in this study we stray away from the traditional bootstrap calibration method.

Instead, we introduce an alternative method to calibrate bootstrapped confidence intervals, the fractional bootstrap. It has the same time and storage complexity as traditional bootstrap methods. Specifically, we introduce a parameter ϕ ∈ (0, 1] denoting the fraction of data to be re-sampled in generating a bootstrapped dataset. By re-sampling less data than in the original dataset, we widen bootstrapped confidence intervals and increase coverage. In the fractional standard bootstrap, we re-sample T ⋅ G ⋅ ϕ plays (rows) uniformly with replacement. In the fractional cluster bootstrap, we re-sample G ⋅ ϕ games uniformly with replacement, keeping each observed play within each re-sampled game. Finally, in the fractional randomized cluster bootstrap, we re-sample G ⋅ ϕ games uniformly with replacement, and within each game re-sample T plays uniformly with replacement.

In Table 6 we report the results of applying the randomized cluster bootstrap to our simulation study for various values of ϕ. To mimic the historical dataset of American football plays, each simulated dataset consists of G = 4,101 games, T = 56 plays per game, and K = T plays per game that share the same outcome. As expected, lowering ϕ widens the confidence intervals and increases marginal coverage. In order to achieve 90 % marginal coverage, ϕ needs to be as small as 0.35. Those intervals have a mean width of 6.3 %. This result is striking: in our simplified setting of random walk football, win probability confidence intervals need to be extremely wide to achieve approximately valid coverage. This exemplifies the difficulty of accurately estimating win probability by fitting a machine learning model from noisy and highly correlated football game outcomes. These estimators are subject to large uncertainty.

Marginal coverage is a sufficient condition for confidence intervals to be “good,” but it is not a necessary requirement for decision making. Even with 90 % marginal coverage, it could be that C I ( x ) always covers W P ( x ) for 90 % of game-states x and never covers for the other 10 %. It may be disastrous to make decisions at the game-states for which intervals never cover. To check that these intervals achieve reasonable coverage across the space of game-states, we bin game-states x by their true win probability W P ( x ) and consider coverage in each bin. In Figure 5 we visualize coverage and its standard error across such bins. For bins near the middle ( W P ≈ 0.5 ) or the extremes ( W P ≈ 0 and W P ≈ 1 ), coverage is high enough. For other bins ( W P ≈ 0.3 and W P ≈ 0.7 ), the intervals remain undercovered, albeit slightly (coverage hovers around 85 %). The game is still competitive in those regions and strategic decisions matter. In future work, we recommend exploring more refined confidence intervals that achieve higher conditional coverage.

Figure 5:

Coverage of 90 % bootstrapped confidence intervals, via the fractional randomized cluster bootstrap for ϕ = 0.35, across the space of game-states x binned by W P ( x ) .

In practice, it is impossible to tune the bootstrap at all, using either the traditional method from DiCiccio and Efron (1996) or the fractional bootstrap. This is because we need to know true win probability in order to tune the β or ϕ values in these alternative bootstraps to achieve adequate coverage. While we know true win probability in our simulation setting, it is an unobservable quantity in real life. Coverage is not calculable for real data because the win/loss outcome is either 0 or 1 and win probability estimates lie in (0,1).

This severe limitation of the tuned bootstrap, nonetheless, does not render its development in this paper worthless – it helped us further illustrate that win probability estimates are subject to substantial uncertainty. Furthermore, though imperfect and difficult, we propose a way to tune the bootstrap with real data as follows. Succinctly, we suggest tuning the fractional bootstrap using a hyper-realistic generative win probability model. To do so, first, given the real historical play-by-play dataset, fit as realistic and granular a probabilistic state-space win probability model as possible. As discussed in Section 2.3, this entails fitting a play-level transition probability model, which is then propagated into win probability by simulating games. Then, generate M synthetic play-by-play datasets from the simulator, apply the ϕ-fractional randomized cluster bootstrap for various values of ϕ to each of them, and select the value of ϕ that achieves desired marginal coverage.

As stressed in Section 2.3, fitting a hyper-realistic football simulator is a delicate and extremely difficult task. This is because it requires a careful encoding of the convoluted rules of football into a set of states and the actions between those states and careful estimation of transition probabilities. Those we know who have created such simulators – including professional sports bettors, football analysts, and hedge fund analysts – do not make them publicly available because they are proprietary and because they use them to make money on the betting markets. One contact said it took him eight months to construct such a simulator.