Dixon-Coles Model Explained: Soccer Prediction Math

In 1997, Mark Dixon and Stuart Coles published "Modelling Association Football Scores and Inefficiencies in the Football Betting Market" in the journal Applied Statistics. Twenty-eight years later, this single paper is still the gold standard for academic soccer prediction. It is referenced in nearly every modern soccer ML paper, used as the baseline that machine learning models must beat, and implemented in production at sportsbooks, prediction sites, and academic research groups.

What makes Dixon-Coles remarkable is its parsimony. The entire model uses just 2N + 2 parameters: an attack rating and a defense rating for each of the N teams in the league, plus a global home advantage parameter (γ) and a low-score correction parameter (ρ). For a 20-team league like the EPL, that is just 42 parameters total. For comparison, modern XGBoost models for soccer prediction typically use 40-100 features with hundreds of internal tree splits — and rarely beat Dixon-Coles by more than 1 percentage point of accuracy.

The reason Dixon-Coles still works comes down to a fundamental property of soccer: goals follow a Poisson distribution. They are discrete count events that occur at a roughly constant rate within a match. The model captures this structure directly. There is no feature-engineering layer to mess up, no overfitting risk from noisy auxiliary stats, no hyperparameter tuning that could drift away from the true data-generating process. Dixon-Coles is what statisticians call a "structurally correct" model — its mathematical form matches the underlying physics of the problem.

For any match between home team h and away team a, Dixon-Coles assumes home goals X and away goals Y are jointly distributed as:

X ~ Poisson(λ)
Y ~ Poisson(μ)

where:
  λ = exp(α_h + β_a + γ)    # home rate
  μ = exp(α_a + β_h)        # away rate

α_i = team i's attack strength (positive = good)
β_i = team i's defense strength (positive = bad — concedes more)
γ   = home advantage (positive — home teams score more)

Notice the asymmetry: home rate uses home team's attack plus away team's defense plus home advantage; away rate uses only away team's attack plus home team's defense (no home boost). This captures the empirical reality that home teams score about 0.3 more goals per match than away teams across the top European leagues.

The exponential function ensures rates are always positive. The team strengths α and β are constrained so that the average attack and average defense are zero (an identifiability constraint — without it, you could shift all α values by a constant and shift γ to compensate). This means α and β are interpretable as deviations from league average.

λ = exp(0.42 + 0.24 + 0.21) = exp(0.87) ≈ 2.39 goals
μ = exp(-0.31 + (-0.28)) = exp(-0.59) ≈ 0.55 goals

If you stop at independent Poisson, you have what statisticians call "the simple Poisson model" — it works adequately but consistently underestimates the frequency of low-scoring results. Real soccer matches end 0-0, 0-1, 1-0, or 1-1 more often than independent Poisson predicts. Across 13,495 European matches in our analysis, the empirical 0-0 rate was 8.4% while pure Poisson predicted 7.1%. The 1-1 rate was 11.2% empirically versus 9.8% modeled.

Dixon and Coles introduced a single correction parameter ρ (rho) that adjusts these four cells:

τ(x, y, λ, μ, ρ) =
    1 - λμρ        if (x, y) = (0, 0)
    1 + λρ         if (x, y) = (0, 1)
    1 + μρ         if (x, y) = (1, 0)
    1 - ρ          if (x, y) = (1, 1)
    1              otherwise

Adjusted joint probability:
P(X=x, Y=y) = Pois(x; λ) × Pois(y; μ) × τ(x, y, λ, μ, ρ)

With ρ > 0 (typical fitted values are 0.05 to 0.15), the correction inflates 0-0 and 1-1 probabilities (drawn low-scoring outcomes) while slightly deflating 0-1 and 1-0 (decisive low-scoring outcomes). After applying the correction, you renormalize the joint distribution to sum to 1.

The original Dixon-Coles paper introduced time-decay weighting to handle the fact that team strength changes over time. A match from two seasons ago should not contribute as much information about a team's current ability as a match from last week. The decay function:

weight_i = exp(-ξ × days_ago_i)

where ξ controls the half-life:
  ξ = 0.00231 → 300-day half-life
  ξ = 0.00650 → 107-day half-life  ← we use this
  ξ = 0.01155 → 60-day half-life

We use ξ = 0.0065, which means a match from 107 days ago contributes half as much information as a match from yesterday. This handles squad turnover, manager firings, January transfer windows, and form swings naturally — without requiring explicit feature engineering for each of those events.

The choice of ξ is a hyperparameter that should be tuned per league. Top-five European leagues with relatively stable rosters work well around ξ = 0.005 to 0.008. Leagues with more volatile rosters (MLS, Saudi Pro League) might benefit from faster decay (ξ = 0.010+) because team identity changes faster relative to match volume.

Given a training set of historical matches, we estimate the parameters (α, β, γ, ρ) by maximizing the time-weighted log-likelihood:

def negative_log_likelihood(params, matches):
    α, β, γ, ρ = unpack(params)
    nll = 0
    for match in matches:
        h, a = match.home_team, match.away_team
        x, y = match.home_goals, match.away_goals
        days_ago = (today - match.date).days
        w = exp(-ξ × days_ago)

        λ = exp(α[h] + β[a] + γ)
        μ = exp(α[a] + β[h])

        log_p = (
            poisson_logpmf(x, λ) +
            poisson_logpmf(y, μ) +
            log(τ(x, y, λ, μ, ρ))
        )
        nll -= w × log_p
    return nll

We minimize this NLL using L-BFGS-B (a quasi-Newton method that handles bound constraints well). For a standard EPL training set of ~3,000 historical matches with 20 teams, the optimization converges in roughly 100 iterations and 3-5 seconds on a modern CPU.

Critical implementation details that turn this from a 30-minute fit into a 5-second fit:

• Vectorize everything. Use NumPy arrays for the inner loop. Compute λ and μ as vectors across all matches simultaneously, not in a Python for-loop.
• Use scipy.special.gammaln for log-factorials. The Poisson PMF involves log(x!) which becomes a bottleneck for naive implementations. gammaln(x+1) is the vectorized equivalent of log(x!) and is dramatically faster.
• Vectorize the tau correction. Mask the four low-score cells with boolean arrays rather than checking each match in a loop.
• Identifiability constraint via penalty. Adding a soft penalty for sum(α) deviating from zero is simpler than removing a parameter, and lets you keep all teams symmetric in the parameter vector.

Once you have fitted (α, β, γ, ρ), prediction for any future match is a closed-form computation. Build the joint goal distribution by computing P(X=i, Y=j) for i, j in [0, 8] (nine possible scores per side covers essentially 100% of probability mass):

def predict_match(model, home_team, away_team):
    λ = exp(model.α[home_team] + model.β[away_team] + model.γ)
    μ = exp(model.α[away_team] + model.β[home_team])

    # Build 9×9 joint distribution
    joint = np.zeros((9, 9))
    for i in range(9):
        for j in range(9):
            joint[i, j] = (
                poisson.pmf(i, λ) ×
                poisson.pmf(j, μ) ×
                τ(i, j, λ, μ, model.ρ)
            )
    joint /= joint.sum()  # renormalize

    # Aggregate into market probabilities
    P_home = sum(joint[i, j] for i, j in pairs if i > j)
    P_draw = sum(joint[i, i] for i in range(9))
    P_away = sum(joint[i, j] for i, j in pairs if i < j)

    return P_home, P_draw, P_away, joint

The joint distribution gives you every market for free. From the same 9×9 grid you can derive:

• Moneyline (1X2): sum cells by margin sign
• Asian handicap at line h: sum cells where (home_margin + h) > 0
• Totals over/under N: sum cells where (i + j) > N
• Both teams to score: 1 - P(home shutout) - P(away shutout)
• Correct score: joint[i, j] directly
• Half-time / full-time: requires fitting separate model on first-half goals

This is one of the most underrated properties of Dixon-Coles. Train one model, generate every soccer betting market. Most ML approaches require training a separate classifier for each market type.

We tested Dixon-Coles on EPL and Serie A using a strict walk-forward methodology: for each test month from January 2024 through May 2026, we trained the model on all matches that occurred strictly before that month, then predicted every match within the test month. This produces a clean, unbiased estimate of how the model would perform in real-time deployment.

Aggregate results across 1,216 graded predictions:

The accuracy number alone (49.3%) sounds modest until you compare it to the structural ceiling. Industry literature places competitive soccer ML between 52% and 58%. Our pure Dixon-Coles is in the lower end of that range, and we beat the naive home-favorite baseline by +3.5pp without any feature engineering beyond goals and time-decay.

Dixon-Coles is not a universal solution. Three known failure modes: