# Modelling time series where model output feeds into data

I’m trying to model a rating system for a competitive game that comprises of 2 parts: a match prediction network and a rating update network, each having several parameters.

Each team starts off with the same rating.
Iterating through the matches in ascending chronological order, a logistic regression predicts the match outcome based on a single predictor: the difference in the 2 team’s ratings, so that there is an intercept (b0) and a coefficient (b1) to be identified.
The predicted outcome from the logistic regression and the actual outcome are then used with a gain parameter (alpha) to detail how much rating should be assigned to the winner (and consequently removed from the loser).
Each team has their rating accordingly updated.
After several rounds of competition the ratings have steadied out.

i.e. for game i in 1…N:
team1 = games[i][‘team1’]
team2 = games[i][‘team2’]
logit(Pr(team1win)) = b0 + b1 * (rating[team1] - rating[team2])
rating_update = alpha * (outcome[i] - Pr(team1win))
ratings[team1] += rating_update
ratings[team2] -= rating_update

I’ve framed it as a general optimisation problem by using the total accuracy over the dataset as the fitness value and fitting the 3 parameters using a genetic algorithm works well, but I’d be interested in treating the parameters as random variables using PyMC.
However, I can’t work out how to parameterise it since the data inputs (ratings and Pr(team1win)) aren’t fixed in advance and are dependent upon the output of the 2 models at each timestep.

Instead of looking at individual rounds I’d just write this down as one logistic regression.

``````with pm.Model():
raw = pm.Normal('rating_raw', sd=1, shape=n_teams)
rating = pm.Deterministic('rating', rating_raw * rating_sd)
pm.Deterministic('elo', 1400 + 400 * rating / np.log(10))
pm.Bernoulli('winner', logit_p=rating[team1] - rating[team2], observed=team1_won)
``````

where `team1` and `team2` are arrays with `n_games` entries, that each contain the index of the corresponding team in game `i`. `team1_won` is 1 if `team1[i]` won game `i` and 0 otherwise.

2 Likes

I’ve implemented this model and adapted it slightly, as in the third line both ‘rating_raw’ and ‘rating_sd’ were undefined. What is the purpose of this line actually, I don’t understand the difference in the standard deviation of `rating_raw` and multiplying `rating_raw` by a (different) sd?

``````with pm.Model() as mod_logistic:
rating_raw = pm.Normal('rating_raw', sd=1, shape=n_teams)
rating = pm.Deterministic('rating', rating_raw * 2)
pm.Deterministic('elo', 1400 + 400 * rating / np.log(10))
pm.Bernoulli('winner', logit_p=rating[team1] - rating[team2], observed=team1_won)
``````

On my dataset (I’ll try and upload a simulated test one soon) it has worked and produced sensible ratings. However, these ratings are averages over the timeframe covered by the training data. The reason I implemented this as a time series in my Genetic Algorithm method was so that I can track each how each team’s rating changes over time.

I forgot to define those vars…

``````with pm.Model() as mod_logistic:
rating_sd = pm.HalfNormal('rating_sd', sd=2)
rating_raw = pm.Normal('rating_raw', sd=1, shape=n_teams)
rating = pm.Deterministic('rating', rating_sd * rating_raw)
pm.Deterministic('elo', 1400 + 400 * rating / np.log(10))
pm.Bernoulli('winner', logit_p=rating[team1] - rating[team2], observed=team1_won)
``````

You could also write that as

``````with pm.Model() as mod_logistic:
rating_sd = pm.HalfNormal('rating_sd', sd=2)
rating = pm.Normal('rating', sd=rating_sd, shape=n_teams)
pm.Deterministic('elo', 1400 + 400 * rating / np.log(10))
pm.Bernoulli('winner', logit_p=rating[team1] - rating[team2], observed=team1_won)
``````

The first one is a centered, the second a non-centered parametrization. They give the same results, but the sampler usually likes one of them much more than the other.
A nice explanation is here.

About the time dependence. That is I think a bit more tricky, but with a bit of experimentation it should be doable.

First, I guess you have to pick a definition for “time”. Is it the number of games played? Actual time? The log of the game number (people usually improve faster at the beginning)? I’d start with a simple linear change of rating over time (whichever). While not particularly realistic, it should be pretty simple to implement and gives you something to compare other models to. Something like that:

``````with pm.Model() as mod_logistic:
rating_mu_sd = pm.HalfNormal('rating_sd', sd=2)
rating_mu = pm.Normal('rating', sd=rating_mu_sd, shape=n_teams)

sd = pm.HalfNormal('rating_time_sd', sd=2)
raw = pm.Normal('rating_time', sd=1, shape=n_teams)
rating_time = pm.Deterministic('rating_time', sd * raw)

rating_time1 = rating_time[team1] * team1_log_gametime
rating_time2 = rating_time[team2] * team2_log_gametime
theta = (rating[team1] + rating_time1) - (rating[team2] + rating_time2)

max_rating = rating + max_log_gametime
pm.Deterministic('max_elo', 1400 + 400 * max_rating / np.log(10))
pm.Bernoulli('winner', logit_p=theta, observed=team1_won)
``````

Where `team1_log_gametime` is an array/dataframe of shape `(n_games)` that contains in at index `i` the log of how many games team1 played up to that particular game. `max_log_gametime` contains the log of the number of games each team has played up to now. You might want to standardize all `log_gametime` values so that `mean(max_log_gametime) == 0` and `std(max_log_gametime) == 1`.

It’s probably also a good idea to think a bit about the mean of `rating_time`. Does it even make sense to have negative values here? Somehow we are measuring the performance of teams not absolute but in relation to other teams, so maybe?

If you are lucky the linear log_gametime works mostly. If you want to look at non-linear changes over time you could also try GPs, but with more than a few teams that probably takes some work to get right.

1 Like

That’s really useful, I hadn’t thought about including time in that manner.
The only issue now is that I don’t follow what `max_log_gametime` is and am getting an input mismatch where it should be added to `rating_mu`, where the latter has `nteams` on dimension 0 and `max_log_gametime` has `ngames` dimension 0. This is because I thought it should be an `ngames` vector with each entry `i` being `max(team1_log_gametime[i], team2_log_gametime[i])`, but I can’t work out from your description or from playing around what it should be.

``````with pm.Model() as mod_time:
rating_mu_sd = pm.HalfNormal('rating_mu_sd', sd=2)
rating_mu = pm.Normal('rating_mu', sd=rating_mu_sd, shape=n_teams)

time_sd = pm.HalfNormal('time_sd', sd=2)
time_raw = pm.Normal('time_raw', sd=1, shape=n_teams)
rating_time = pm.Deterministic('rating_time', time_sd * time_raw)

rating_time1 = rating_time[team1] * team1_log_gametime
rating_time2 = rating_time[team2] * team2_log_gametime
theta = (rating_mu[team1] + rating_time1) - (rating_mu[team2] + rating_time2)

max_rating = rating_mu + max_log_gametime
pm.Deterministic('max_elo', 1400 + 400 * max_rating / np.log(10))
pm.Bernoulli('winner', logit_p=theta, observed=team1_won)``````

Does that help?

``````df = pd.DataFrame({
'team1': ['blueberry', 'bilberry', 'strawberry', 'blueberry'],
'team2': ['bilberry', 'strawberry', 'blueberry', 'strawberry'],
'team1_won': [False, True, False, True],
'time': ['2019-07-01', '2019-07-02', '2019-07-03', '2019-07-04'],
'team1_gamenum': [1, 2, 2, 3], # How many games did team1 already play?
}).rename_axis(index='game')
df['time'] = pd.to_datetime(df['time'])
# How many games were played by

max_log_gametime = np.log(df.team1.value_counts() + df.team2.value_counts())

max_log_gametime.index = max_log_gametime.index.astype('category')
df = df.astype({
'team1': max_log_gametime.index.dtype,
'team2': max_log_gametime.index.dtype})
n_teams = len(max_log_gametime)

with pm.Model() as model:
rating_mu_sd = pm.HalfNormal('rating_mu_sd', sd=2)
rating_mu = pm.Normal('rating_mu', sd=rating_mu_sd, shape=n_teams)

sd = pm.HalfNormal('rating_time_sd', sd=2)
raw = pm.Normal('rating_time_raw', sd=1, shape=n_teams)
rating_time = pm.Deterministic('rating_time', sd * raw)

theta = (
(rating_mu[df.team1.cat.codes] + rating_time1)
- (rating_mu[df.team2.cat.codes] + rating_time2))

max_rating = rating_mu + max_log_gametime.values
pm.Deterministic('max_elo', 1400 + 400 * max_rating / np.log(10))
pm.Bernoulli('winner', logit_p=theta, observed=df.team1_won)

trace = pm.sample()
``````
2 Likes

Yes it does indeed, thank you so much for your help!

1 Like