Code for simulating the outcome of the Elam Ending

elam_ending_simulation.R

# Code for creating an Elam Ending simulation model #
rm(list = ls(all = TRUE))

library(dplyr)
library(tidyr)

loc<-"" # enter file path for exporting model output #

# point probabilities for each possession : 0-1-2-3-4 points #
tp<-c(0.505,0.031,0.325,0.137,0.002) # team with possession #
op<-c(0.505,0.031,0.325,0.137,0.002) # team without possession #

tpc<-cumsum(tp)
opc<-cumsum(op)

# simulate each ame state 10,000 times to create win probabilities # 
# distance to target score: from 1 to 24 points #
# score differential: from -20 to +20 points #
gol<-vector('list')
mol<-vector('list')
for (str in 24:1) {
for (pos in 0:1) {
  for (mgn in -20:20) {
    print(paste(mgn,str,pos))
    flush.console()
    gs<-c(mgn,str,pos) # score margin, distance to target, possession #
    
    gml<-vector('list')
    edl<-vector('list')
    for (j in 1:10000) {
      m<-matrix(nrow=1001,ncol=7)
      colnames(m)<-c("i","m","dst","pos","ts","os","s")
      m[1,]<-c(0,gs,gs[1],0,NA)
      for (i in 1:1000) {
        g<-m[i,2:4]
        if (g[3]==1) {p<-tpc} else {p<-opc}
        r<-runif(n=1)
        s<-findInterval(r,p)
        m[i,7]<-s
        g1<-ifelse(g[3]==1,g[1]+s,g[1]-s)
        g2<-ifelse(g[3]==1,ifelse(g[1]<0,g[2]-max(0,s+g[1]),g[2]-s),ifelse(g[1]>0,g[2]-max(0,s-g[1]),g[2]-s))
        g3<-ifelse(g[3]==1,0,1)
        m[i+1,1]<-i
        m[i+1,2]<-g1
        m[i+1,3]<-g2
        m[i+1,4]<-g3
        m[i+1,5]<-ifelse(g[3]==1,m[i,5]+s,m[i,5])
        m[i+1,6]<-ifelse(g[3]==0,m[i,6]+s,m[i,6])
        if (g2<=0) {
          m<-m[!is.na(m[,1]),]
          break}
      }
      gml[[j]]<-m
      edl[[j]]<-m[nrow(m),1:(ncol(m)-1)]
    }
    edd<-data.frame(do.call('rbind',edl),stringsAsFactors = F)
    gmd<-data.frame(do.call('rbind',gml),stringsAsFactors = F)
    gmd$cnt<-1
    gma<-aggregate(gmd$cnt,by=list(gmd$m,gmd$dst),FUN=sum)
    colnames(gma)<-c("m","dst","cnt")
    gma<-cbind(mgn,str,pos,j,gma)
    wp<-mean(edd$m>0)
    mpos<-mean(edd$i)
    mmgn<-mean(edd$m)
    mmgnw<-mean(edd$m[edd$m>0])
    mmgnl<-mean(edd$m[edd$m<0])
    mgnz<-sum(edd$m==0)
    mo<-matrix(nrow=1,ncol=10,data=c(mgn,str,pos,wp,mpos,mmgn,mmgnw,mmgnl,mgnz,j))
    colnames(mo)<-c('mgn','str','pos','wp','mpos','mmgn','mmgnw','mmgnl','mgnz','j')
    mol[[paste(mgn,str,pos)]]<-mo
    gol[[paste(mgn,str,pos)]]<-gma
  }
}
}
mod<-data.frame(do.call('rbind',mol),stringsAsFactors = F)
god<-data.frame(do.call('rbind',gol),stringsAsFactors = F)

out<-mod %>%
  filter(mgn>=0, str>0) %>%
  select(str,mgn,pos,wp,mpos,mmgn) %>%
  arrange(-str,mgn,pos)

# raw data used for Elam Ending calculator #
write_csv(out,path=paste(loc,"elamcalc.csv"))