State-Space Models for Football Matches¶

Sam Duffield deployed a factorial state-space model for predicting the 2026 World Cup games, cuthberto-carlos. During this period, the model was consistently making predictions close to, if not better than market predictions on Polymarket. I decided to do a complete writeup to have a better understanding of the model and the workflow.

The model is built off the factorial state-space model from Duffield, Power and Rimella: A state-space perspective on modelling and inference for online skill rating. Details on how to use the exact model can be found in cuthbert.factorial (link) and other methods associated with factorial state-space models.

1 Factorial State-Space Models (fSSMs)¶

Duffield, Power and Remilla (2024) formally introduces the concept of factorial state-space models (fSSMs). In a general state space model, the joint distribution is defined as

where \(m_0\) is the intial distribution, \(\mathcal{M}_{t-1, t}\) is transition distribution and \(G_t\) is the likelihood function given the current state. We can simplify this model by making a few assumptions:

1. assuming that the latent state of indexes \(k\) are independent \(x_0^k \sim m_0^k\),
2. evolution of indexes are independent and Markovian \(x_{t-1}^{k} \mid x_t^{i} \sim \mathcal{M}_{t, t'}(x_{t}^k, \cdot)\) and
3. Given the state of index \(i\), the outcome \(y_t\) is conditionally independent of all other indexes \(-i\) and is drawn from \(G_t (y_k \mid x_t^i)\).

This gives us a simplified model known as the fSSM,

The fSSM only involves \(\mathcal{O}(N + K)\) terms as compared to the original \(\mathcal{O}(N \cdot K)\) in equation (1).

2 Football Matches with fSSMs¶

We define the number of matches \(N\) and number of teams as \(K\). We also denote that each football match consist of 2 teams, \(h(t)\) and \(a(t)\), representing the home and away teams respectively. We utilize the model defined in cuthberto-carlos,

where \(\phi_t = \exp(- \kappa \Delta t)\), \(Q_t = \Sigma_0 - \Phi_t \Sigma_0 \Phi_t\) and \(\Phi_t = \text{diag}(\phi_t)\).

\(y_t = (y_t^{\text{h}}, y_t^{\text{a}})\) is the observed goals for the home \(X_t^{\text{h}}\) and away \(X_t^{\text{a}}\) teams at time \(t\).

where \(\lambda_1 = \exp(\alpha + x_t^{\text{att}, \text{h}} - x_t^{\text{def}, \text{a}})\), \(\lambda_2 = \exp(\alpha + x_t^{\text{att}, \text{a}} - x_t^{\text{def}, \text{h}})\), \(\lambda_3 = \exp(\beta)\).

The static parameters are \(\Theta = \{\mu_0 \in \mathbb{R}^{2}, \Sigma_0 \in \mathbb{R}^{2 \times 2}, \alpha \in \mathbb{R}, \beta \in \mathbb{R}\}\). This model implicitly says that the teams skill level revert to a common level \(\mu_0\) and uncertainty \(\Sigma_0\) over time.

3 Filtering and Smoothing¶

Inference is performed with a linearized moments Kalman filter (from cuthbert.gaussian.moments), wrapped in the factorial structure of the fSSM (via cuthbert.factorial.gaussian). By assuming that the the non-linear likelihood \(G_t\) produces an approximate distribution, the filter can be computed analytically - predict and update steps are closed-form Kalman recursions.

where \(H\) is the Jacobian of the conditional mean with respect to the state, \(d\) is the intercept term and \(R\) is the Cholesky factor of the observation covariance, used for numerical stability during optimization. This method is extremely helpful with high dimensional problems and where we believe that the underlying state uncertainty is distributed normally.

3.1 Linearization of bivariate-poisson model¶

To linearize the bivariate Poisson likelihood function, we perform the following steps.

Step 1 — Conditional moments. The observation model is summarized by its conditional mean and covariance given the state. For the bivariate Poisson, the mean is

and the covariance is

Step 2 — First-order Taylor expansion. We expand the conditional mean around the predicted state \(\bar{x}_{t \mid t-1}\):

where \(H\) is the Jacobian of the conditional mean with respect to the state, evaluated at \(\bar{x}_{t \mid t-1}\).

Step 3 — Observation matrix \(H\). For a match between home team \(i\) and away team \(j\), the joint state is \(x_t = (x_t^{\text{att},i}, x_t^{\text{def},i}, x_t^{\text{att},j}, x_t^{\text{def},j})\). Using \(\frac{\partial}{\partial z} e^{a + z/s} = \frac{\lambda}{s}\), the Jacobian is

where \(s\) is the scale (1 for competitive matches, friendly_scale for friendlies).

Step 4 — Intercept \(d\). The intercept is chosen so the linearization matches the true mean at the linearization point:

Step 5 — Observation covariance \(R\). Putting these together gives the linear-Gaussian approximation

where \(R = \text{chol}\big(\text{Cov}[y_t \mid x_t]\big)\). In the implementation, \(H\) is obtained automatically via JAX automatic differentiation of the conditional mean function rather than by hand.

3.2 Prediction and Update¶

Each match applies the standard Kalman recursion for the specific teams. Given the filtered state \((\mu_{t-1 \mid t-1}, \Sigma_{t-1 \mid t-1})\) at the previous time step, the predict step propagates the state forward through the (exactly linear) OU dynamics:

Note that as \(\Delta t \to \infty\), \(\phi_t \to 0\) and the predicted state reverts to the prior \((\mu_0, \Sigma_0)\).

The update step then conditions on the observed goals \(y_t\) using the linearized observation model \(G_t(y_t \mid x_t) \approx \mathcal{N}(y_t \mid H x_t + d, RR^\top)\):

where \(K_t\) is the Kalman gain.

4 Parameter Estimation¶

Since the marginal likelihood \(p(y_{1:T} \mid \Theta)\) is tractable, the optimal parameters can be easily computed. Furthermore, the linearization of the likelihood function introduces a bias on the posterior, so "exact" EM isnt readily available.

where \(Z_t\) is the normalizing constant of the observation update.

An alternative to parameter estimation is the expectation-maximization (EM) algorithm, which iteratively maximizes the joint likelihood \(p(x_{0:T}, y_{1:T})\) given a set of parameters denoted by \(\theta \rightarrow Q(\theta \mid \hat{\theta})\).