Premise and approach

I wanted to do a clean statistical build-out of a real-world situation during this Utah chapter, but my original plan to try and predict snowfall in LCC proved to be much more complex than expected. As such, given the upcoming World Cup and my personal interests in soccer, I decided to think through and build out a model for the tournament, assisted, as usual, by the beast Claude Code.

The work contained here has strong contribution from Claude - I would not have been able to get anywhere near this far without its help. At the same time, I've contributed key insights and advancements where Claude would have otherwise fell into a loop, lost track of the overall objective, or simply given up from data limitation excuses. It goes to show what the current state of Human-AI collaboration looks like.

Philosophy

Sports analytics is a rich field with some very smart people involved, and I don't necessarily see this beating out commercial models or highly-specialized pundits. What I do see this as is a clean bottom-up statistical framework for thinking about a team-sport tournament, and I hope that the statistical methodology framed below does the Yale Department of Statistics justice.

In our model, we're not just picking a "winner" for each match-up - instead, we're modeling the underlying expected number of goals each team will score, then use that to feed a prediction for each matchup. We'll build this model off the (https://www.ajbuckeconbikesail.net/wkpapers/Airports/MVPoisson/soccer_betting.pdf)[Dixon-Coles model], a well-established model for soccer that modifies an underlying Poisson distribution model for goals scored. The work below adds a couple additional pieces to simply fitting parameters of the Dixon-Coles model to historical international matchups, summarized in the schematic below.

---

The Poisson foundation

We begin with all of our modeling by understanding the Poisson distribution, which models the occurrence of rare, independent events (already an approximation, since goals are, by their very nature, very dependent on previous goals). Suspending our disbelief, if we model a specific team $i$ as expected to score $\mu$ goals against team $j$, then the probability of actually scoring exactly $x$ goals is:

The parameter $\mu$ is what we are trying to estimate. It is not a fixed constant — it depends on who is playing and their relative strengths. In our naive model, we define $\mu$ as a product of three factors:

The product form models how these factors compound - if $\alpha_i = 1.5$ (a strong attack, 50% above average) and $\beta_j = 1.4$ (a very leaky defense, 40% above average), the expected goals rate is $1.5 \times 1.4 = 2.1$ (110% more likely to score), not $1.5 + 1.4 - 1 = 1.9$ (90% more likely to score). Previous work with football data shows that the multiplicative model works better than additive alternatives.

The expected goals for the other team are definitionally symmetric:

With these two parameters, we can apply the Poisson model to set up an independent joint bivariate distribution. Again, this is a simplification - obviously, a team down a goal will play more conservatively, changing the math, but we're coming up with a approximate model, not a perfect description.

Under independence, the joint probability of a specific scoreline $(x, y)$ is just the product of the two marginal Poisson probabilities:

To get outcome probabilities, we sum over all scorelines in the appropriate region. For all practical purposes we can truncate at ten goals per team, which captures over 99.9% of the probability mass for any realistic $\mu, \nu < 4$:

For World Cup matches played at neutral sites — which describes the entire 2026 tournament — we set $\gamma = 1$.

---

The Dixon-Coles correction

The independent Poisson model has a known failure mode: it systematically underestimates the frequency of 0-0 draws, overestimates 1-1 draws, and slightly misestimates 1-0 and 0-1 results. The reason is that real matches have game states. When a match is 0-0, both teams are often pushing harder to break the deadlock — the game opens up, but also tightens tactically. When one team goes ahead, they may sit back and protect, suppressing further scoring. These dynamics create a mild negative correlation between the two teams' goals precisely at low-scoring outcomes, which the independence assumption cannot capture.

Dixon and Coles addressed this with an elegant correction. Rather than replacing the Poisson model entirely, they multiply the joint probability by a correction factor $\tau$ that only modifies the four cells where $x + y \leq 1$:

The single parameter $\rho \leq 0$ controls the magnitude of the adjustment. With $\rho < 0$, the factor $\tau(0,0) = 1 - \rho\mu\nu > 1$ inflates the 0-0 probability, and probability mass is redistributed away from the 1-1 cell. The correction is valid — that is, the corrected probabilities still sum to one across all scorelines — by construction.

Empirically, fits on European league data give $\rho \approx -0.13$. International football, which features fewer high-scoring matches and more conservative defensive tactics, tends to produce slightly more negative values. Our fit on 1,908 WC-team matches yielded $\rho \approx -0.11$, consistent with this range. We treat $\rho$ as a free parameter and estimate it jointly with the team strength parameters.

---

Composite data decay

This is one of the main ways that we refine the naive model that simply trains the alpha and beta parameters of the Dixon-Coles setup with historical data. The most naive approach, of course, would simply be to treat each game as contributing the same amount of signal. Obviously, that doesn't make sense, so the slightly less naive approach is to weight each match result's contribution to the model fitting by some decay factor. Instead of just using that, our model integrates three additional signals to weight how much each match informs the training, combining it into a composite factor:

Each factor captures a distinct mechanism by which a historical result becomes more or less representative of the current team.

Calendar decay

The first factor, $e^{-\phi d_k}$, is standard exponential decay over elapsed days. Even with the same manager and squad, teams evolve: tactics develop, individual form rises and falls, set-piece routines change. We set the half-life to 600 days (~20 months), meaning a result from 600 days ago receives half the weight of a result from today. This is calibrated for international football, where squad turnover is slower and tactical systems are more stable than at club level.

Competition quality

The second factor, $s_k$, reflects that not all matches are equally diagnostic. A win over a top-10 opponent at a major tournament tells you much more about a team's current ceiling than a friendly against a mid-table qualifier. We weight matches by competition tier:

Competition	Weight
FIFA World Cup	1.0
Continental championships (Euro, Copa América, AFCON, Asian Cup)	0.9
Nations Leagues (UEFA, CONCACAF)	0.65
World Cup qualifiers	0.6
International friendlies	0.3

For the major footballing nations with abundant high-quality data, friendly results end up contributing very little total weight. For smaller nations whose few World Cup appearances are separated by decades, qualifier and friendly matches remain important training signal.

Manager epoch discount

The third factor, $\delta^{n_k}$, penalizes results from previous managerial regimes. Here $n_k$ is the number of manager changes that have occurred on either team between match $k$ and today, and $\delta \in (0, 1)$ is a per-transition discount.

A new manager represents a partial reset of tactical identity — the same eleven players will press differently, organize differently, and execute set pieces differently. Results from a prior regime are still informative because individual player quality persists, but they are weaker guides to what this team will do today. We use $\delta = 0.85$, meaning each manager change reduces a match's effective weight by ~15%. This is intentionally moderate: the calendar and competition-quality components already do most of the down-weighting work, so the epoch factor handles the residual discontinuity that smooth decay misses without over-penalizing teams with recent transitions.

As an example: a 2-1 win under a manager who left three seasons ago, after two subsequent regime changes, receives $0.85^2 = 0.72$ of the weight it would have had under the current manager — on top of whatever calendar and competition discounts already apply.

Player continuity (designed, pending data)

The full specification includes a fourth factor, $c_k^{\,\alpha}$, measuring the importance-weighted overlap between the squad that played in match $k$ and the current national squad:

A match where nine of eleven starters are still active has $c_k \approx 0.9$; a result from six years ago featuring players who have since retired might have $c_k \approx 0.2$. This factor is designed and its weighting is defined, but not yet incorporated into the current fit — it is pending historical squad roster data collection. The current model runs on calendar decay, competition quality, and manager epochs only.

---

Model architecture

The estimation pipeline runs in three sequential stages, each a residual layer on top of the prior.

Stage 1 fits one scalar attack and one scalar defence parameter per team from historical match outcomes alone — the classic Dixon-Coles MLE. It captures how many goals each team typically scores and concedes across all opponents.

Stage 2 adds a bilinear style interaction between every team pair, also trained from match outcomes. It captures matchup-specific edges: whether the way team $i$ attacks creates a systematic advantage against the way team $j$ defends, beyond what raw scalar strength predicts.

Stage 3 introduces a player attribute correction derived from current squad ratings. It captures a complementary signal: given who is actually in the squad today, does the current roster composition suggest an adjustment to what historical match data alone predicts? Stages 1 and 2 are retrospective — they answer how a team has performed. Stage 3 is prospective — it answers how strong this particular group of players looks right now.

Each stage is fit with prior-stage parameters frozen, learning only the residual that earlier stages cannot explain. Stages 1 and 2 are fully estimated and active. Stage 3 is designed and implemented in parallel — its integration into the live simulation is pending confirmed 23-man squad announcements ahead of the tournament.

---

Stage 1 — DC base layer

The first stage fits one scalar attack parameter $\alpha_i$ and one scalar defence parameter $\beta_i$ per team. These are estimated from historical match scorelines via maximum likelihood — no shot-level data required. The composite decay weights ensure that stale, low-quality, or pre-regime-change matches contribute minimally to the final values.

The fit runs on 1,908 matches between the 48 WC-qualified teams from January 2010 onward. Fitted globals: $\hat{\rho} \approx -0.11$, home advantage $\hat{\gamma} \approx 1.17$ (neutralized at the 2026 tournament). The scalar parameters are then frozen before Stage 2 begins.

---

Stage 2 — Bilinear style layer

Scalar $\alpha_i, \beta_j$ summarize team quality as single numbers. They cannot capture matchup-specific dynamics: whether France's defensive block is particularly effective against teams that prefer wide diagonal balls, or less effective against high-tempo pressing sides. That information is invisible to a scalar model.

Each team is assigned a $K$-dimensional attack style vector $\mathbf{u}_i$ and defence style vector $\mathbf{v}_i$. Their dot product adds a signed scalar residual on top of the DC log-expected-goals:

The dot product is positive when team $i$'s attacking style creates a systematic advantage against team $j$'s defensive shape, and negative when it does not. With $K = 4$ dimensions the model can represent up to four independent axes of stylistic variation. In practice, the average matchup adjustment is about 2%; the largest is approximately +57% (Spain attacking France), a pattern visible across multiple recent encounters.

Style vectors are initialized with small Gaussian noise ($\sigma = 0.05$) to break the zero-gradient saddle at the origin, then estimated by maximizing the DC likelihood with L2 regularization ($\lambda = 2.0$) on $\mathbf{u}_i, \mathbf{v}_j$. Stage 1 parameters remain frozen throughout. Stage 2 converged in 116 iterations; the refined correlation parameter is $\hat{\rho} \approx -0.14$.

---

Stage 3 — Player attribute correction

Stages 1 and 2 draw on a single data source: historical match outcomes. They answer how a team has performed. Player attribute data answers a different question: how strong is this team right now, given who is actually in the squad?

These two signals are complementary but not redundant. A team that surged in form two years ago but has since lost key players will have inflated DC parameters relative to its current ability. A team that quietly upgraded its squad during a flat qualifying campaign will be underrated. A young nation with sparse international history but a generation of technically developed players currently dominating European club football has almost no training signal in the DC fit — its $\alpha, \beta$ will be conservative by default. Neither the scalar parameters nor the style vectors are designed to catch any of this: their training signal is entirely retrospective.

Stage 3 introduces a correction term $\gamma_{ij}$ derived from a separate linear model trained on current player ratings, then centered to remove systematic bias before being applied to the DC output.

Squad embeddings

Player ratings are sourced from the EA FC 26 database, which covers 18,400+ players and includes all 48 WC-qualified nations. For each team, the squad is divided into four positional groups — goalkeeper, defenders, midfielders, forwards — and the mean attribute vector is computed within each group. The full team embedding concatenates these four group means:

where $M$ is the number of attributes per group. This preserves positional structure: a team with elite defenders and average midfielders produces a different embedding from one with the reverse imbalance, even if the squad-wide averages are identical. Data-sparse nations (Kenya, Honduras, Jordan) where DC parameters are least reliable are fully covered, since every FIFA-registered player has ratings regardless of league visibility.

The linear model

For each historical match $k$, the attribute model independently predicts log goals for team $i$ as a linear function of both teams' positional embeddings:

where $\mathbf{e}_{i}^{\,\text{atk}}$ selects team $i$'s midfield and forward embeddings and $\mathbf{e}_{j}^{\,\text{def}}$ selects team $j$'s goalkeeper and defender embeddings. The training target is $\log(X_k + 0.5)$, where $X_k$ is the observed goals scored. The log offset maps zero-goal outcomes to a finite value and keeps the target in the same scale as the DC log-expected-goals outputs. The model is fit by ridge regression using the same composite decay weights as Stage 1 — meaning recent, high-quality matches dominate the fit here too.

Centering and application

The linear model is not trusted for its absolute level: it lacks Poisson structure, has no temporal decay within matches, and is trained using current ratings projected back onto historical matches where squads may have differed substantially. What it is trusted for is its relative signal — which matchups look better or worse given current squad composition, compared to the average.

For each training match $k$, compute the residual between the linear prediction and the Stage 1+2 DC output:

The mean residual $\bar{\delta}$ is the systematic offset between the two models — reflecting both scale differences and the fact that current ratings are an imperfect proxy for historical squad quality. Subtracting it from each prediction leaves only the matchup-relative signal:

By construction, $\gamma$ has zero mean across training data. Applying it does not inflate or deflate aggregate expected goals — it only shifts individual matchups relative to the DC baseline. The final expected goals for any WC 2026 match are:

The shrinkage factor $\lambda \in (0,\,1]$ governs how much weight to place on the attribute correction relative to historical outcomes. It is calibrated on held-out matches. A smaller $\lambda$ defers more to the historical record; a larger one gives more weight to current roster composition. The Dixon-Coles $\tau$ correction is applied to the final $\mu_{ij}^{\text{final}}$, not the intermediate DC output, so the low-score cell adjustment remains consistent with the corrected expected goals.

---

Tournament simulation

With the full scoreline distribution $\tilde{P}(x, y)$ available for any match between any two teams, we simulate the entire tournament structure by Monte Carlo. Each iteration proceeds as follows:

1. Group stage. For each of the 104 group-stage matches, sample a scoreline from $\tilde{P}(x, y)$. Accumulate group points (3 for a win, 1 for a draw, 0 for a loss) and rank teams by points, then goal difference, then goals scored, then head-to-head result.

2. Knockout rounds. Draws are resolved by extra time (modeled as a second independent Poisson draw at half the regular rate) and if still level, a penalty shootout modeled as a coin flip.

3. Repeat. Run $N = 100{,}000$ full simulations. The reported probability for any event is simply its frequency across all simulations. One hundred thousand draws is sufficient for stable estimates down to roughly 0.1% probability.

Before the simulation loop, the full $48 \times 48$ matrix of expected goals is precomputed once and reused across all draws, keeping per-simulation cost efficient.

---

Parameter estimation

Stage 1 — DC base. Attack and defence parameters $\alpha_i, \beta_j$ and correlation $\rho$ are estimated by maximizing the composite-weighted DC log-likelihood:

Log-link reparameterization ($\alpha_i = e^{a_i}$, $\beta_i = e^{b_i}$) ensures positivity without box constraints. Scale is fixed by constraining $\mathrm{mean}(\{a_i\}) = 0$, so attack ratings are relative to the 48-team geometric mean.

Stage 2 — Style vectors. With $\alpha_i, \beta_j$ frozen, style vectors and a refined $\rho$ are estimated by maximizing the same DC likelihood with a bilinear expected goals term, plus an L2 penalty:

Analytical gradients for both the Poisson terms and the DC correction cells are passed directly to L-BFGS-B, avoiding finite-difference approximation over the 385-parameter space. Vectors are initialized with small Gaussian noise ($\sigma = 0.05$) to break the zero-gradient saddle at the origin.

Stage 3 — Attribute correction. The ridge coefficient vector $\mathbf{w}$ is fit by minimizing composite-weighted squared error between predicted and observed log goals. The ridge penalty is cross-validated on held-out matches. The centering constant $\bar{\delta}$ and shrinkage factor $\lambda$ are computed on the training split.

---

Limitations

Temporal alignment in Stage 3. The current implementation trains the attribute correction using today's FC 26 ratings projected back onto historical matches, where the actual squads may have differed significantly. The correct approach aligns each match with the FIFA ratings edition closest to its date. Historical editions (FIFA 15–23) are available; this alignment is in progress. The centering step absorbs the systematic component of the mismatch, but matchup-specific noise remains.

Player continuity not yet active. The $c_k^{\,\alpha}$ squad-overlap component of the composite weight is designed but pending historical squad roster data. It would specifically correct cases where a team's historical record is dominated by a generation of players since retired — calendar decay handles this only loosely.

Style vectors are team-level, not player-level. The bilinear layer captures aggregate stylistic matchup effects but cannot represent individual player interactions: whether a particular striker's movement exploits a specific centre-back's positional tendencies, or whether a midfield pair creates unusual press resistance against a specific pressing shape.

No within-match dynamics. The Poisson model assumes a constant goal rate throughout the match. In reality, the rate shifts with game state — a team trailing by one in the 75th minute plays very differently from a team level at half-time. The Dixon-Robinson extension addresses this but requires modelling the full score path rather than just the final result.

Penalty shootouts as coin flips. The model treats shootouts as 50-50. Historical shootout records, specialist penalty takers, and goalkeeper save rates are all partially predictive and worth incorporating as the knockout rounds approach.