NL_eredivisie_2014_2019.Rmd

---
title: "Bayesian estimation of a Poisson model for Dutch football matches odds"
author: "Piet Stam"
date: "August 15th, 2020"
output: 
  github_document:
    toc: true
    toc_depth: 4
---

```{r knitr_setup, include=FALSE, purl=FALSE}
knitr::opts_chunk$set(message=FALSE, warning=FALSE, fig.width=9,
                      fig.height=5, dpi=90, fig.path='results/')

knitr::knit_hooks$set(purl = knitr::hook_purl)
```

> *Copyright 2019 [Piet Stam](http://www.pietstam.nl). The code and the documentation are licensed under the Creative Commons [Attribution 4.0 International license](http://creativecommons.org/licenses/by/4.0/).* 

## Introduction

### Quick summary
We applied the original work of [Rasmus Baath](http://www.sumsar.net/about.html) to the Dutch Eredivisie football competition. With `r-bayesian-football-odds` we estimated the odds of football matches in the last two weeks of the 2018/2019 Dutch Eredivisie football competiton. We provide the code and evaluate the results of our predictions. 

### Acknowledgements
This piece of work is based on the works of [Rasmus Baath](http://www.sumsar.net/blog/2013/07/modeling-match-results-in-la-liga-part-one/). Rasmus Baath submitted his code to the [UseR 2013 data analysis contest](https://www.r-project.org/conferences/useR-2013/) and licensed it under the Creative Commons [Attribution 3.0 Unported license](http://creativecommons.org/licenses/by/3.0/). 

He predicted the results of the 50 last matches of the 2012/2013 Spanish LaLiga season. He used data of the 2008/09-2012/13 seasons (5 seasons in total) to estimate his regression model in a [Bayesian](https://en.wikipedia.org/wiki/Bayes_estimator) way. See [this thread](https://stats.stackexchange.com/questions/252577/bayes-regression-how-is-it-done-in-comparison-to-standard-regression) for an intuitive explanation of the difference between the bayesian approach and the conventional approaches of linear regression and maximum likelihood.

I slightly adpated his code to predict the results of the last two competition rounds (that is, the last 18 matches) of the 2018/2019 Dutch Eredivisie season. These predictions are based on soccer match data of the 2014/15-2018/19 seasons (5 seasons in total). The source of these data is [here](http://www.football-data.co.uk/netherlandsm.php). Out of the three model specifications that Rasmus developed, I used the most sophisticated model that allowed for year-to-year variability in team skill (called "iteration 3" by Rasmus).

You can find my code at [GitHub](https://github.com/pjastam/r-bayesian-football-odds). Rasmus deserves all the credits, I deserve all the blame in case of any errors in my application to the Dutch football competition.

## Data and methods

### Theoretical description of the model

<!-- The source of the folowing function is  [https://github.com/STAT545-UBC/Discussion/issues/102](https://github.com/STAT545-UBC/Discussion/issues/102). -->

```{r, echo=FALSE, warning=FALSE, purl=FALSE}
latexImg = function(latex){

    link = paste0('http://latex.codecogs.com/gif.latex?',
           gsub('\\=','%3D',URLencode(latex)))

    link = gsub("(%..)","\\U\\1",link,perl=TRUE)
    return(paste0('![](',link,')'))
}
```

#### Basic model

The first thing to notice is that not all teams are equally good. Therefore, it will be assumed that all teams have a latent skill variable and the skill of a team *minus* the skill of the opposing team defines the predicted outcome of a game. As the number of goals are assumed to be Poisson distributed it is natural that the skills of the teams are on the log scale of the mean of the distribution.

In its simplest form, the distribution of the number of goals for team $i$ when facing team $j$ is then

`r latexImg('{goals}_{i,j} \\sim \\text{Poisson}(\\lambda_{i,j})')`

`r latexImg('\\text{log}(\\lambda_{i,j}) = \\text{baseline} + \\text{skill}_i - \\text{skill}_j')`

where <code>baseline</code> is the log average number of goals when both teams are equally good. Note that this model description does not capture the variation in the number of goals among football seasons and between home vs away teams.

#### General model

In order to allow for variation in the number of goals among football seasons and between home vs away teams, we refine the distribution of the goal outcome of a match between home team $i$ and away team $j$ in season $s$ as follows: 

`r latexImg('{goals}^\\text{home}_{s,i,j} \\sim \\text{Poisson}(\\lambda^\\text{home}_{s,i,j})')` 

with the <code>lambdas</code> defined as follows 

`r latexImg('\\lambda^\\text{home}_{s,i,j} = \\exp(\\text{baseline}^\\text{home}_s + \\text{skill}_{s,i} - \\text{skill}_{s,j})')`

`r latexImg('\\lambda^\\text{away}_{s,i,j} = \\exp(\\text{baseline}^\\text{away}_s + \\text{skill}_{s,j} - \\text{skill}_{s,i})')`

Note that the <code>baseline</code> is split into <code>home_baseline</code> and <code>away_baseline</code> in order to account for the home advantage. Furthermore, we introduced the index t for the baseline and skill parameters to allow for variation among seasons.

##### Defining the baseline distributions

I set the prior distributions of the baselines in season $s$ to:

`r latexImg('\\text{baseline}^\\text{home}_s \\sim \\text{Normal}(\\text{baseline}^\\text{home}_{s-1}, \\sigma_{\\text{seasons}}^2)')`

`r latexImg('\\text{baseline}^\\text{away}_s \\sim \\text{Normal}(\\text{baseline}^\\text{away}_{s-1}, \\sigma_{\\text{seasons}}^2)')`

and in the *first* season to:

`r latexImg('\\text{baseline}^\\text{home}_1 \\sim \\text{Normal}(0, 4^2)')`

`r latexImg('\\text{baseline}^\\text{away}_1 \\sim \\text{Normal}(0, 4^2)')`

with <code>sigma-seasons</code> defined as:

`r latexImg('\\sigma_\\text{seasons} \\sim \\text{Uniform}(0, 3)')`

##### Defining the team skill distributions

I set the prior distributions over the skills of team $i$ (or $j$, denoted as i|j) in season $s$ to:

`r latexImg('\\text{skill}_{s,i|j} \\sim \\text{Normal}(\\text{skill}_{s-1,i|j}, \\sigma_{\\text{seasons}}^2)')`

and in the *first* season to:

`r latexImg('\\text{skill}_{1,i|j} \\sim \\text{Normal}(\\mu_\\text{teams}, \\sigma_{\\text{teams}}^2)')`

with the <code>sigma-seasons</code> defined as above and <code>mu-teams</code> and <code>sigma-teams</code> defined as:

`r latexImg('\\mu_\\text{teams} \\sim \\text{Normal}(0, 4^2)')`

`r latexImg('\\sigma_\\text{teams} \\sim \\text{Uniform}(0, 3)')`

We apply a normalizing restriction with respect to the (arbitrarily chosen) *first* team in each season $s$ as follows

`r latexImg('\\text{skill}_{s,1} = 0')`

We choose very vague priors. For example, the prior on the baseline have a SD of 4 but since this is on the log scale of the mean number of goals it corresponds to one SD from the mean $0$ covering the range of $[0.02, 54.6]$ goals. A very wide prior indeed.

#### Probabilistic Graphical Model

We graphed the dependencies described above with the help of a probabilistic graphical model. To this end, we make use of the CRAN package [DiagrammeR](https://cran.r-project.org/web/packages/DiagrammeR/index.html) with the help of which you can use the [Graph Visualization Software](https://graphviz.gitlab.io/).

```{r pgm, echo=FALSE, warning=FALSE, purl=FALSE}
library(DiagrammeR)

grViz("
digraph G {
  rankdir=LR
  bgcolor=grey
  
  labelloc='b'
  labeljust='r'

  node [shape=rectangle]

  a [label='N(0,4^2)', color=none, fontcolor=white]
  b [label='N(0,4^2)', color=none, fontcolor=white]
  c [label='mu_teams', shape=ellipse, color=none, fontcolor=black]
  d [label='sigma_teams', shape=ellipse, color=none, fontcolor=black]
  e [label='sigma_seasons', shape=ellipse, color=none, fontcolor=black]
  
  a -> q [color=white]
  b -> r [color=white]

  c -> t [color=white]
  d -> t [color=white]

  e -> h [color=white]
  e -> s [color=white]
  
  e -> f [color=white]
  e -> g [color=white]
  e -> o [color=white]
  e -> p [color=white]
  
  l [label='U(0,3)', color=none, fontcolor=white]
  m [label='N(0,4^2)', color=none, fontcolor=white]
  n [label='U(0,3)', color=none, fontcolor=white]
  
  l -> e [color=white]
  m -> c [color=white]
  n -> d [color=white]
  
  subgraph cluster_seasons {
    label='Season s &isin; {1,...,S}'

    q [label='baseline_home 1', color=none, fontcolor=black, style=filled, fillcolor=none]
    r [label='baseline_away 1', color=none, fontcolor=black, style=filled, fillcolor=none]
    o [label='baseline_home s-1', color=none, fontcolor=black, style=filled, fillcolor=none]
    p [label='baseline_away s-1', color=none, fontcolor=black, style=filled, fillcolor=none]
    f [label='baseline_home s', color=none, fontcolor=black, style=filled, fillcolor=none]
    g [label='baseline_away s', color=none, fontcolor=black, style=filled, fillcolor=none]
  
    q -> o [color=white]
    r -> p [color=white]
    o -> f [color=white]
    p -> g [color=white]
  
    f -> j [color=white]
    g -> k [color=white]
  
    subgraph cluster_teams {
      label='Team i,j &isin; {1,...,T | i < > j}'

      h [label='skill s,i|j', color=none, fontcolor=black, style=filled, fillcolor=none]
      s [label='skill s-1,i|j', color=none, fontcolor=black, style=filled, fillcolor=none]
      t [label='skill 1,i|j', color=none, fontcolor=black, style=filled, fillcolor=none]
      j [label='goals_home s,i,j', color=none, fontcolor=black, style=filled, fillcolor=none]
      k [label='goals_away s,i,j', color=none, fontcolor=black, style=filled, fillcolor=none]
  
      t -> s [color=white]
      s -> h [color=white]
      h -> j [color=white]
      h -> k [color=white]
    }
  }
}
")
```

### Read data

```{r setup, include=FALSE}
library(runjags)
library(coda)
library(mcmcplots)
library(stringr)
library(plyr)
library(dplyr)
library(xtable)
library(ggplot2)

source("functions/plotPost.R")

set.seed(12345)

# Convenience function to generate the type of column names Jags outputs.
col_name <- function(name, ...) {
  paste0(name, "[", paste(..., sep=",") , "]")
}
```

We first read the Dutch soccer match data of the 2014/15-2018/19 Dutch Eredivisie seasons from the original csv-files and cache them. The result is a database called `eredivisie`.

```{r scan_data_dir, cache = TRUE, cache.extra = file.info(list.files(path = "data", pattern = ".csv", full.names = TRUE)), include=FALSE, purl=FALSE}
```

```{r raw_data, cache = 2, dependson = "scan_data_dir"}
from_year <- 2014
to_year <- 2019
source(paste0("functions/Import_Data_Eredivisie.R"))
```

Then the cached `eredivisie` data are cleaned and new variables are created.

```{r processed_data, cache = 2, dependson = "raw_data"}
eredivisie <- eredivisie %>%
  mutate(MatchResult = sign(HomeGoals - AwayGoals)) # -1 Away win, 0 Draw, 1 Home win

# Creating a data frame d with only the complete match results
d <- na.omit(eredivisie)

# Lists with the unique names of all teams and all seasons in the database
teams <- unique(c(d$HomeTeam, d$AwayTeam))
seasons <- unique(d$Season)

# A list for JAGS with the data from d where the strings are coded as integers
data_list <- list(HomeGoals = d$HomeGoals, AwayGoals = d$AwayGoals, 
                  HomeTeam = as.numeric(factor(d$HomeTeam, levels=teams)),
                  AwayTeam = as.numeric(factor(d$AwayTeam, levels=teams)),
                  Season = as.numeric(factor(d$Season, levels=seasons)),
                  n_teams = length(teams), n_games = nrow(d), 
                  n_seasons = length(seasons))
```

The data set <code>eredivisie</code> contains data from 5 different seasons. In this model we allow for variability in year-to-year team performance. This variablitity in team performance can be demonstrated by the following diagram, which shows that some teams do not even participate in all seasons in the <code>eredivisie</code> data set, as a result of dropping out of the first division:

```{r participation_by_season, fig.height=5, fig.width=7}
qplot(Season, HomeTeam, data=d, ylab="Team", xlab = "Season")
```

## Estimation, simulation and validation

### Model estimation

Turning this into a JAGS model results in the following string. Note that the model loops over all seasons and all match results. JAGS parameterizes the normal distribution with precision (the reciprocal of the variance) instead of variance so the hyper priors have to be converted. Finally, we \"anchor\" the skill of one team to a constant otherwise the mean skill can drift away freely. Doing these adjustments results in the following model description:

```{r jags_model_description, cache=2, tidy=FALSE}
m3_string <- "model {
for(i in 1:n_games) {
HomeGoals[i] ~ dpois(lambda_home[Season[i], HomeTeam[i],AwayTeam[i]])
AwayGoals[i] ~ dpois(lambda_away[Season[i], HomeTeam[i],AwayTeam[i]])
}

for(season_i in 1:n_seasons) {
for(home_i in 1:n_teams) {
for(away_i in 1:n_teams) {
lambda_home[season_i, home_i, away_i] <- exp( home_baseline[season_i] + skill[season_i, home_i] - skill[season_i, away_i])
lambda_away[season_i, home_i, away_i] <- exp( away_baseline[season_i] + skill[season_i, away_i] - skill[season_i, home_i])
}
}
}

skill[1, 1] <- 0 
for(j in 2:n_teams) {
skill[1, j] ~ dnorm(group_skill, group_tau)
}

group_skill ~ dnorm(0, 0.0625)
group_tau <- 1/pow(group_sigma, 2)
group_sigma ~ dunif(0, 3)

home_baseline[1] ~ dnorm(0, 0.0625)
away_baseline[1] ~ dnorm(0, 0.0625)

for(season_i in 2:n_seasons) {
skill[season_i, 1] <- 0 
for(j in 2:n_teams) {
skill[season_i, j] ~ dnorm(skill[season_i - 1, j], season_tau)
}
home_baseline[season_i] ~ dnorm(home_baseline[season_i - 1], season_tau)
away_baseline[season_i] ~ dnorm(away_baseline[season_i - 1], season_tau)
}

season_tau <- 1/pow(season_sigma, 2) 
season_sigma ~ dunif(0, 3) 
}"
```

We then run this model directly from R using RJAGS and the <code>textConnection</code> function. This takes about half an hour on my computer, but of course this depends on the configuration at hand.

```{r monte_carlo, cache=2, dependson=list("processed_data","jags_model_description"), results='hide'}
# Compiling the model
m3 <- run.jags(method="parallel",
               model=m3_string,
               monitor=c("home_baseline", "away_baseline","skill", "season_sigma", "group_sigma", "group_skill"),
               data=data_list,
               n.chains=3,
               adapt=10000,
               burnin=10000,
               sample=15000,
               thin=8,
               summarise=FALSE,
               plots=FALSE)
# Generating MCMC samples
s3 <- as.mcmc.list(m3)
# Merging the three MCMC chains into one matrix
ms3 <- as.matrix(s3)
```

The following graphs shows the trace plots and probability distributions of the team mean, team sigma and season sigma parameters, respectively.

```{r mu_sigma_params, fig.height=3}
plot(s3[, "group_skill"])
plot(s3[, "group_sigma"])
plot(s3[, "season_sigma"])
```

We can also calculate the default home advantage by looking at the difference between <code>exp(home_baseline) - exp(away_baseline)</code>. The next graph shows that there is a home advantage of more than 0.4 goals, on average, and it differs significantly from zero.

```{r overall_home_advantage}
plotPost(exp(ms3[,col_name("home_baseline",to_year-from_year)]) - exp(ms3[,col_name("away_baseline",to_year-from_year)]), compVal = 0, xlab = "Home advantage in number of goals")
```

### Generating predictions (in- and out-of-sample)

In the <code>eredivisie</code> data set included in this project, the results of the 18 last games of the 2018/2019 season are missing. Using our model we can now both predict and simulate the outcomes of these 18 games. The R code below calculates a number of measures for each game (both the games with known and unknown outcomes):

  - The mode of the simulated number of goals, that is, the *most likely* number of scored goals. If we were asked to bet on the number of goals in a game this is what we would use.
  - The mean of the simulated number of goals, this is our best guess of the average number of goals in a game.
  - The most likely match result for each game.
  - A random sample from the distributions of credible home scores, away scores and match results. This is how the Eredivisie actually could have played out in an alternative reality.

```{r model_predictions, cache=2, dependson="monte_carlo", tidy=FALSE}
n <- nrow(ms3)
m3_pred <- sapply(1:nrow(eredivisie), function(i) {
  home_team <- which(teams == eredivisie$HomeTeam[i])
  away_team <- which(teams == eredivisie$AwayTeam[i])
  season <- which(seasons == eredivisie$Season[i])
  home_skill <- ms3[, col_name("skill", season, home_team)] 
  away_skill <- ms3[, col_name("skill", season, away_team)]
  home_baseline <- ms3[, col_name("home_baseline", season)]
  away_baseline <- ms3[, col_name("away_baseline", season)]
  
  home_goals <- rpois(n, exp(home_baseline + home_skill - away_skill))
  away_goals <- rpois(n, exp(away_baseline + away_skill - home_skill))
  home_goals_table <- table(home_goals)
  away_goals_table <- table(away_goals)
  match_results <- sign(home_goals - away_goals)
  match_results_table <- table(match_results)
  
  mode_home_goal <- as.numeric(names(home_goals_table)[ which.max(home_goals_table)])
  mode_away_goal <- as.numeric(names(away_goals_table)[ which.max(away_goals_table)])
  match_result <-  as.numeric(names(match_results_table)[which.max(match_results_table)])
  rand_i <- sample(seq_along(home_goals), 1)
  
  c(mode_home_goal = mode_home_goal, mode_away_goal = mode_away_goal, match_result = match_result,
    mean_home_goal = mean(home_goals), mean_away_goal = mean(away_goals),
    rand_home_goal = home_goals[rand_i], rand_away_goal = away_goals[rand_i],
    rand_match_result = match_results[rand_i])
})
m3_pred <- t(m3_pred)
```

### Model validation

First let\'s compare the distribution of the actual number of goals in the data with the predicted mode, mean and randomized number of goals for all the games (focusing on the number of goals for the home team).  

First the actual distribution of the number of goals for the home teams.

```{r hist_home_goal, fig.height=3, fig.width=4}
hist(eredivisie$HomeGoals, breaks= (-1:max(eredivisie$HomeGoals, na.rm=TRUE)) + 0.5, xlim=c(-0.5, 10), main = "Distribution of the number of goals\nscored by a home team in a match",
    xlab = "")
```

This next plot shows the distribution of the modes from the predicted distribution of home goals from each game. That is, what is the most probable outcome, for the home team, in each game.

```{r mode_home_goal, fig.height=3, fig.width=4}
hist(m3_pred[ , "mode_home_goal"], breaks= (-1:max(m3_pred[ , "mode_home_goal"])) + 0.5, xlim=c(-0.5, 10),
    main = "Distribution of predicted most \nprobable score by a home team in\na match",
    xlab = "")
```

For almost all games the single most likely number of goals is one. Actually, if you know nothing about an Eredivisie game, betting on one goal for the home team is 78 % of the times the best bet.

Let\'s instead look at the distribution of the predicted mean number of home goals in each game.

```{r mean_home_goal, fig.height=3, fig.width=4}
hist(m3_pred[ , "mean_home_goal"], breaks= (-1:max(m3_pred[ , "mean_home_goal"])) + 0.5, xlim=c(-0.5, 10),
    main = "Distribution of predicted mean \n score by a home team in a match",
    xlab = "")
```

For most games the expected number of goals are 2. That is, even if your safest bet is one goal you would expect to see around two goals.

The distribution of the mode and the mean number of goals doesn't look remotely like the actual number of goals. This was not to be expected, we would however expect the distribution of randomized goals (where for each match the number of goals has been randomly drawn from that match's predicted home goal distribution) to look similar to the actual number of home goals. Looking at the histogram below, this seems to be the case.

```{r rand_home_goal, fig.height=3, fig.width=4}
hist(m3_pred[ , "rand_home_goal"], breaks= (-1:max(m3_pred[ , "rand_home_goal"])) + 0.5, xlim=c(-0.5, 10),
    main = "Distribution of randomly drawn \n score by a home team in a match",
    xlab = "")
```

We can also look at how well the model predicts the data. This should probably be done using cross validation, but as the number of effective parameters are much smaller than the number of data points a direct comparison should at least give an estimated prediction accuracy in the right ballpark. 

```{r validation_mode_home_goal}
mean(eredivisie$HomeGoals == m3_pred[ , "mode_home_goal"], na.rm=T)
```
```{r validation_mean_home_goal} 
mean((eredivisie$HomeGoals - m3_pred[ , "mean_home_goal"])^2, na.rm=T)
```

So on average the model predicts the correct number of home goals 31% of the time and guesses the average number of goals with a mean squared error of 1.51. Now we'll look at the actual and predicted match outcomes. The graph below shows the match outcomes in the data with 1 being a home win, 0 being a draw and -1 being a win for the away team.

```{r hist_actual_match_result, fig.height=3, fig.width=4}
hist(eredivisie$MatchResult, breaks= (-2:1) + 0.5, xlim=c(-1.5, 1.5), ylim=c(0, 1000), main = "Actual match results",
    xlab = "")
```

Now looking at the most probable outcomes of the matches according to the model.

```{r hist_pred_match_result, fig.height=3, fig.width=4}
hist(m3_pred[ , "match_result"], breaks= (-2:1) + 0.5, xlim=c(-1.5, 1.5), ylim=c(0, 1000), main = "Predicted match results",
    xlab = "")
```

For almost all matches the safest bet is to bet on the home team. While draws are not uncommon it is *never* the safest bet.

As in the case with the number of home goals, the randomized match outcomes have a distribution similar to the actual match outcomes: 

```{r hist_rand_match_result, fig.height=3, fig.width=4}
hist(m3_pred[ , "rand_match_result"], breaks= (-2:1) + 0.5, xlim=c(-1.5, 1.5), ylim=c(0, 1000), main = "Randomized match results",
    xlab = "")
```

```{r validation_match_result}
mean(eredivisie$MatchResult == m3_pred[ , "match_result"], na.rm=T)
```
The model predicts the correct match outcome (i.e. home team wins / a draw / away team wins) 57% of the time. Pretty good!

## Results

Disclaimer: my comments below may be out of sync with the empirical results and graphs, because these comments (as well as the description of my betting experience in the last section) are based on the results of running [VERSION 1.0](https://github.com/pjastam/r-bayesian-football-odds/releases/tag/v1.0) instead of the current version of the app.

### The ranking of the teams

We'll start by ranking the teams of the <code>Eredivisie</code> using the estimated skill parameters from the 2018/2019 season, which are based on the estimation sample of the five seasons 2014/2015-2018/2019. Note that for one of the teams the skill parameter is "anchored at zero". This "anchoring" is done for the very same \"identification\" reason that one of the parameters in a traditional logit analysis is always set to zero by default: the value of a parameter automatically follows if you already know all the other parameters in your model and given the fact that probabilities always sum up to 1 in total.

Consequently, as Rasmus noted before, the skill parameters are difficult to interpret as they are relative to the skill of the team that had its skill parameter "anchored" at zero. To put them on a more interpretable scale the skill paramters are first zero centered by subtracting the mean skill of all teams. Then he added the home baseline and exponentiated the resulting values. These rescaled skill parameters are now on the scale of expected number of goals when playing as a home team.

```{r team_skill, fig.height=8, dpi=90}
team_skill <- ms3[, str_detect(string=colnames(ms3), paste0("skill\\[",to_year-from_year,","))]
team_skill <- (team_skill - rowMeans(team_skill)) + ms3[, paste0("home_baseline[",to_year-from_year,"]")]
team_skill <- exp(team_skill)
colnames(team_skill) <- teams
team_skill <- team_skill[,order(colMeans(team_skill), decreasing=T)]
old_par <- par(mar=c(2,0.7,0.7,0.7), xaxs='i')
caterplot(team_skill, labels.loc="above", val.lim=c(0.7, 3.8))
par(old_par)
```

Two teams are clearly ahead of the rest, Ajax and PSV. Let\'s look at the credible difference between these two teams. Ajax is a better team than PSV with a probabilty of 74%, i.e. the odds in favor of Ajax are 74% / 26% = 3. So, on average, PSV only wins one out of four games that they play against Ajax.

```{r team_skill_PSV_Ajax, fig.height=3, fig.width=7}
plotPost(team_skill[, "Ajax"] - team_skill[, "PSV Eindhoven"], compVal = 0, xlab = "<- PSV     vs     Ajax ->")
```

### Predicting the future

Now that we've checked that the model reasonably predicts the Eredivisie history let\'s predict the Eredivisie endgame! 

At the time when I executed my version of this model applied to the Dutch Eredivisie competition (2019-05-10), most of the matches in the 2018/2019 season had already been played. Yet two out of 34 competition rounds had to be played (that is, competition rounds 33 and 34). With these two rounds still to go, Ajax and PSV both have 80 points, but Ajax leads the competition as their goal difference is larger (111-30 = 81) than that of PSV (95-24 = 71). The code below displays the predicted mean and mode number of goals for the endgame and the predicted winner of each game.

```{r predict_future_matches, results='asis'}
eredivisie_forecast <- eredivisie[is.na(eredivisie$HomeGoals), c("Season", "Week", "HomeTeam", "AwayTeam")]
m3_forecast <- m3_pred[is.na(eredivisie$HomeGoals),] 
eredivisie_forecast$mean_home_goals <- round(m3_forecast[,"mean_home_goal"], 1) 
eredivisie_forecast$mean_away_goals <- round(m3_forecast[,"mean_away_goal"], 1)
eredivisie_forecast$mode_home_goals <- m3_forecast[,"mode_home_goal"] 
eredivisie_forecast$mode_away_goals <- m3_forecast[,"mode_away_goal"]
eredivisie_forecast$predicted_winner <- ifelse(m3_forecast[ , "match_result"] == 1, eredivisie_forecast$HomeTeam, 
                                           ifelse(m3_forecast[ , "match_result"] == -1, eredivisie_forecast$AwayTeam, "Draw"))

rownames(eredivisie_forecast) <- NULL
print(xtable(eredivisie_forecast, align="cccccccccc"), type="html")
```

These predictions are perfectly useful if you want to bet on the likely winner of each game. However, they do not reflect how the actual endgame will play out, e.g., there is not a single draw in the <code>predicted_winner</code> column. So at last let\'s look at a *possible* version of the Eredivisie endgame by displaying the simulated match results calculated earlier. 

```{r predict_past_matches, results='asis'}
eredivisie_sim <- eredivisie[is.na(eredivisie$HomeGoals), c("Season", "Week", "HomeTeam", "AwayTeam")]
eredivisie_sim$home_goals <- m3_forecast[,"rand_home_goal"] 
eredivisie_sim$away_goals <- m3_forecast[,"rand_away_goal"]
eredivisie_sim$winner <- ifelse(m3_forecast[ , "rand_match_result"] == 1, eredivisie_forecast$HomeTeam, 
                            ifelse(m3_forecast[ , "rand_match_result"] == -1, eredivisie_forecast$AwayTeam, "Draw"))

rownames(eredivisie_sim) <- NULL
print(xtable(eredivisie_sim, align="cccccccc"), type="html")
```

Now we see a number of games resulting in a draw. We also see that Ajax and FC Utrecht tie in round 33, which puts PSV on top of the leaderboard! However, in round 34 the image is reversed when PSV and Heracles tie, against all odds. So, in the end, Ajax wins the competition in this *possible* version of the Eredivisie endgame by their better goal difference.

### Betting on the match outcome

One of the powers with using Bayesian modeling and MCMC sampling is that once you have the MCMC samples of the parameters it is straight forward to calculate any quantity resulting from these estimates while still retaining the uncertainty of the parameter estimates. So let\'s look at the predicted distribution of the number of goals for AZ Alkmaar vs PSV Eindhoven game and see if I can use my model to make some money. I'll start by using the MCMC samples to calculate the distribution of the number of goals for AZ Alkmaar and PSV Eindhoven.

```{r function_plot_goals}
plot_goals <- function(home_goals, away_goals) { 
  old_par <- par(mar = c(0, 0, 0, 0))
  par(mfrow = c(2, 2), mar=rep(2.2, 4))
  
	n_matches <- length(home_goals) 
	goal_diff <- home_goals - away_goals 
	match_result <- ifelse(goal_diff < 0, "away_win", ifelse(goal_diff > 0, "home_win", "equal")) 
	hist(home_goals, xlim = c(-0.5, 10), breaks = (0:100) - 0.5)
	hist(away_goals, xlim = c(-0.5, 10), breaks = (0:100) - 0.5) 
	hist(goal_diff, xlim = c(-6, 6), breaks = (-100:100) - 0.5)
	barplot(table(match_result)/n_matches, ylim = c(0, 1))
	par(old_par)
}
```

```{r setup_top2teams_comparison}
n <- nrow(ms3)
home_team <- which(teams == "AZ Alkmaar")
away_team <- which(teams == "PSV Eindhoven")
season <- which(seasons == paste0(to_year-1,"/",to_year))
home_skill <- ms3[, col_name("skill", season, home_team)] 
away_skill <- ms3[, col_name("skill", season, away_team)]
home_baseline <- ms3[, col_name("home_baseline", season)]
away_baseline <- ms3[, col_name("away_baseline", season)]

home_goals <- rpois(n, exp(home_baseline + home_skill - away_skill))
away_goals <- rpois(n, exp(away_baseline + away_skill - home_skill))
```

Looking at summary of these two distributions in the first two graphs below, it shows that AZ and PSV both have the biggest chance to score one goal (as the modus of both distributions equals 1). From the third graph it follows that the most likely goal difference is 0 or -1: either AZ and PSV draw (0), or PSV scores just one more goal than AZ (-1). In case of the latter, PSV turns out to be the match winner.

The fourth graph shows the probability distribution of a PSV win ('away_win'), a draw ('equal') and AZ win ('home_win'). This graph underlines that a PSV win is a likely scenario: it has a probability of more than 50%. The fact that the balance topples in favor of PSV should then be due to the one goal difference that is attributed a great chance according to the third graph. Note, however, that the probability that PSV will *not* turn out as the match winner (i.e. a draw or a loss) is still almost 50%.

```{r plot_goals, fig.height=5}
old_par <- par(mfrow = c(2, 2), mar=rep(2.2, 4))
plot_goals(home_goals, away_goals)
par(old_par)
```

At May 10th, that is just before the start of competition round 33, you got the following payouts (that is, how much would I get back if my bet was successful) for betting on the outcome of this game, after 288 bets being placed on the betting site [William Hill](http://www.williamhill.com/)

  |  AZ  | Draw | PSV  |
  |:----:|:----:|-----:|
  | 3.90 | 4.00 | 1.78 |

Using my simulated distribution of the number of goals I can calculate the predicted payouts of the model. It appears that the payouts of the model are very close to the payouts that William Hill offers.

```{r}
1 / c(AZ =  mean(home_goals > away_goals), Draw = mean(home_goals == away_goals), PSV = mean(home_goals < away_goals))
```

The most likely result is 1 - 1 with a predicted payout of 9.70, which can be compared to the William Hill payout of 7.50 for this bet. Thus, William Hill thinks that a 1 - 1 draw is even likier than our model predicts. If we want to earn some extra money, we should bet on a 1 - 0 win for AZ, as the William Hill payout is 19 and our model predicts 17.50.

### Betting on the correct score

It is also possible to bet on the final goal outcome so let\'s calculate what payouts my model predicts for different goal outcomes. The payouts that William Hill reports are

  |      | PSV 0 |  PSV 1 | PSV 2 | PSV 3 | PSV 4 |
  |------|:-----:|:------:|:-----:|:-----:|:-----:|
  | AZ 0 |  21   |   12   |  12   |  17   |  29   |
  | AZ 1 |  19   |   7.5  |  12   |  12   |  21   |
  | AZ 2 |  23   |   13   |  11   |  17   |  29   |
  | AZ 3 |  41   |   26   |  23   |  29   |  51   |
  | AZ 4 |  81   |   51   |  66   |  126  |  81   |

It follows that the 1 - 1 draw is also the most likely result at Wiliam Hill. Now, we are going to calculate the payouts that our model predicts.

```{r payouts, results='asis'}
goals_payout <- laply(0:6, function(home_goal) {
  laply(0:6, function(away_goal) {
    1 / mean(home_goals == home_goal & away_goals  == away_goal)
  })
})

colnames(goals_payout) <- paste("PSV Eindhoven", 0:6, sep=" - ")
rownames(goals_payout) <- paste("AZ Alkmaar", 0:6, sep=" - ")
goals_payout <- round(goals_payout, 1)
print(xtable(goals_payout, align="cccccccc"), type="html")
```

The most likely result is 1 - 1 with a predicted payout of 9.70, which can be compared to the William Hill payout of 7.50 for this bet. This, we can earn some extra money if we bet on this end score.

### Betting results

As the English say: "put your money where your mouth is". So, at the time, I placed a 1 euro bet on three potential outcomes of the AZ-PSV match:

* a 0-0 equalizer
* a 1-1 equalizer
* an AZ win over PSV

I uploaded my official betting ticket to GitHub as proof:

![](https://user-images.githubusercontent.com/8238853/61595945-b8fa5580-abfd-11e9-83ea-4933f6bcd9ee.PNG)

As the true outcome of this match was a 1-0 victory of AZ over PSV, the results of my betting adventure are a loss, a loss and a win, respectively. Here you can find the official ticket with my betting results (in reverse order):

![](https://user-images.githubusercontent.com/8238853/61595949-c283bd80-abfd-11e9-8f93-8efb7131be60.PNG)

So, with a 3 euro stake I earned 4.15 euro, which means a profit margin of 1.15 / 3 = 38%.

For those interested, my bets were placed on [the Dutch Lottery website](https://toto.nederlandseloterij.nl/) (disclosure: I do *not* own shares in that company) and I posted about my one-off betting experience [on LinkedIn](https://www.linkedin.com/pulse/de-bal-rond-piet-stam/) (in Dutch) at the time.

## Appendix

### Who do you call?

You can contact me at my GitHub email address if you would like to share your thoughts.