Calculates the Brier Score for the 2020 MLB season using Skiena and Bayes predictors

MLB_brier.r

# Skiena probability in 1-on-1 match
# Based on Skiena "Calculated Bets" 
#
# Parameters
# ----------
#   a,b: Probability of win against all competitors for Teams A, B
#       
# 
# Returns
# -------
#   p: Probability that the stronger team will win against weaker team

skiena_func <- function(a,b){
  (1 + abs(a-b)^0.4)/2
}

# Bayes probability in 1-on-1 match
#
# Parameters
# ----------
#   a,b: Probability of win against all competitors for Teams A, B
#       
# 
# Returns
# -------
#   p: Probability that the stronger team will win against weaker team

bayes_func <- function(a,b){
    c <- a
    a <- max(a,b)
    b <- min(b,c)
    a*(1-b)/(a*(1-b)+(1-a)*b)
}

# Lookup probability tables for Skiena and Bayes estimators 
# Save standings from MLB https://www.mlb.com/standings/2020 in Excel
# Add column for team nicknames: https://www.retrosheet.org/CurrentNames.csv
# Delete blank rows, extraneous columns then add headings "Team" and "Standings"
# Import data: File -> Import Dataset -> From Excel
#
# Parameters
# ----------
#   MLB: Team standings at end of 2020 season
#       
# 
# Returns
# -------
#   lookup: Probability table that the stronger team will win against weaker 
#           team using Skiena and Bayes estimators
#
# Example
# -------
# lookup$skiena["WAS","TEX"]
# 0.6685728

make_lookup_tables <- function(MLB){
  # Extract columns from MLB
  teams <- unlist(MLB[,c(1)])
  standings <- unlist(MLB[,c(2)])
  
  # Tables
  skiena <- outer(standings,standings,skiena_func)
  rownames(skiena) <- teams
  colnames(skiena) <- teams
  
  bayes <- outer(standings,standings,bayes_func)
  rownames(bayes) <- teams
  colnames(bayes) <- teams  
  
  # Return lookup tables
  lookup <- list("skiena" = skiena, "bayes" = bayes)
  return(lookup)
  
}

brier_score <- function(lookup,MLB,games){
  # Initialize Brier Score for each estimator
  BS_Skiena <- 0
  BS_Bayes <- 0
  
  # Number of games
  N <- lengths(games[,c(1)])
  
  # Evaluate Brier Score for each game
  for (i in 1:N){
    # Game data
    game <- games[i,c(4,7,10,11)]
    Team_A <- game[1]$X4
    Team_B <- game[2]$X7
    score_A <- game[3]$X10
    score_B <- game[4]$X11
    
    # Team standings at end of season
    pA <- MLB %>% select(Standing) %>% filter(MLB[,c(1)] == Team_A)
    pA <- pA$Standing
    pB <- MLB %>% select(Standing) %>% filter(MLB[,c(1)] == Team_B)
    pB <- pB$Standing
    
    # Determine outcome
    # Outcome is 1 if (pA < pB) and (score_A < score_B) or
    #                 (pA > pB) and (score_A > score_B) 
    #       otherwise 0
    # Multiply boolean by 1 to convert to numeric (0/1)
    o <- 1*( (pA < pB) & (score_A < score_B) | (pA > pB) & (score_A > score_B) )
    
    # Skiena forecast
    if (pA > pB){
      f_Skiena <- lookup$skiena[Team_A,Team_B]
    } else {
      f_Skiena <- 1 - lookup$skiena[Team_A,Team_B]
    }
    
    
    # Bayes forecast
    if (pA > pB){
      f_Bayes <- lookup$bayes[Team_A,Team_B]
    } else {
      f_Bayes <- lookup$bayes[Team_B,Team_A]
    }
    
    # Update Brier scores
    BS_Skiena <- BS_Skiena + (f_Skiena - o)^2
    BS_Bayes <- BS_Bayes + (f_Bayes - o)^2
    
  }
  
  # Normalize Brier scores by number of games played
  BS_Skiena <- BS_Skiena / N
  BS_Bayes <- BS_Bayes / N
  
  # Return Brier scores for each predictor
  BS <- list("skiena" = BS_Skiena, "bayes" = BS_Bayes)
  return(BS)
  
}