Bradley-Terry model for soccer, incorporating team strength and draw factors

Extended Bradley-Terry model

data {
    int<lower=0> N; // N games
    int<lower=0> P; // P teams

    // Each team is referred to by an integer that acts as an index for the ratings vector.
    int team1[N]; // Indicator arrays for team 1
    int team2[N]; // Indicator arrays for team 1
    int results[N]; // Results. 1 if home win, 2 if draw, 3 if away win.

    real<lower=0> nu_sigma;

    vector[P] alpha; // Parameters for Dirichlet prior.
}

parameters {
    // Vector of ratings for each player
    // The simplex constrains the ratings to sum to 1
    simplex[P] ratings;

    // Parameter adjusting the probability of draw for each team
    // Parameter in logspace
    simplex[P] nu;
}


model {

    // Array of length 3 vectors for the three outcome probabilies for each game.
    vector[3] result_probabilities[N];
    real nu_game;
    real nu_rating_prod;

    ratings ~ dirichlet(alpha); // Dirichlet prior on the ratings.

    nu ~ normal(0, nu_sigma); // exponential prior on nu.

    for (i in 1:N) {

        # Find the draw factor for this game and shift out of logspace
        nu_game = exp(nu[team1[i]] + nu[team1[i]]);

        // nu multiplied by the harmonic mean of the ratings.
        nu_rating_prod = sqrt(ratings[team1[i]] * ratings[team2[i]]) * nu_game;

        result_probabilities[i][3] = nu_rating_prod / (ratings[team1[i]] + ratings[team2[i]] + nu_rating_prod);
        result_probabilities[i][1] = ratings[team1[i]] / (ratings[team1[i]] + ratings[team2[i]] + nu_rating_prod);
        result_probabilities[i][2] = 1 - (result_probabilities[i][1] + result_probabilities[i][3]);

        results[i] ~  categorical(result_probabilities[i]);
    }

}