State-Space Models for Tennis Matches¶

I built a Gaussian factorial state-space model for predicting WTA tennis match outcomes, tennis_ssm. The model tracks a single latent skill per player and produces win probabilities for upcoming fixtures, which are compared against market prices on Polymarket. This writeup documents 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 Tennis Matches with fSSMs¶

We define the number of matches \(N\) and number of players as \(K\). We also denote that each tennis match consist of 2 players, \(i(t)\) and \(j(t)\), representing the two players respectively. We utilize the model defined in tennis_ssm,

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

\(y_t \in \{0, 1\}\) is the observed outcome of the match at time \(t\), where \(y_t = 1\) indicates that player \(i(t)\) won and \(y_t = 0\) indicates that player \(j(t)\) won. The outcome is observed through a sigmoid (logistic) link on the skill difference,

where \(\sigma(x) = (1 + e^{-x})^{-1}\) is the sigmoid function and \(s_d \in \mathbb{R}^+\) is the observation scale. A larger \(s_d\) flattens the win probability toward uniform (0.5), meaning the outcome is more random.

The static parameters are \(\Theta = \{\mu_0 \in \mathbb{R}, \Sigma_0 \in \mathbb{R}, \tau_d \in \mathbb{R}^+, s_d \in \mathbb{R}^+\}\). We fix \(\mu_0 = 0\) so skills are centered at zero. This model implicitly says that the players' skill levels revert to a common level \(\mu_0\) and uncertainty \(\Sigma_0\) over time. Unlike TrueSkill2 and Elo, which use Wiener processes (unbounded variance growth), this model uses an Ornstein-Uhlenbeck process where skills revert to \(\mathcal{N}(\mu_0, \Sigma_0)\) during inactive periods.

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 the sigmoid model¶

To linearize the sigmoid 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 Bernoulli outcome, the conditional mean is the win probability

and the conditional variance is the Bernoulli variance

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 player \(i\) and player \(j\), the joint state is \(x_t = (x_t^{i}, x_t^{j})\). Using \(\frac{\partial}{\partial z} \sigma(z/s) = \frac{1}{s} \sigma(z/s)(1 - \sigma(z/s))\), the Jacobian is

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 = \sqrt{h(x_t)(1 - h(x_t))}\) is the Cholesky factor of the Bernoulli variance. 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 players. 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 outcome \(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.

3.3 Prediction¶

To predict the outcome of a future match, we propagate the latest filtered skill distributions of the two players to the match time and apply the sigmoid link. Given the propagated skill difference \(d = \mu^{i} - \mu^{j}\) with variance \(\sigma_d^2 = \Sigma^{i} + \Sigma^{j}\), the win probability is computed via a probit-style approximation for the sigmoid marginal:

This accounts for both the observation scale \(s_d\) and the uncertainty in the skill difference \(\sigma_d^2\).

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.

tennis_ssm performs parameter estimation via the maximum likelihood estimation (MLE) method. This is done by performing gradient descent on the negative log-likelihood to find the best parameters that maximizes the log-likelihood. This is done using the optax.adam package. The trainable parameters are \(\tau_d\), \(s_d\) and \(\Sigma_0\), each constrained to be positive via a softplus transform, while \(\mu_0\) is fixed at \(0\) (skills are relative). The log normalizing constant of the filter is the log-likelihood of the observed match outcomes under the model, and maximizing it optimizes the parameters.