gistfile1.r

RoR.sim <- function(prizes, probs, BI, BR, games=NA, CI=.95, sims=1000, plots=200, warns=TRUE) {
  if(sum(prizes < 0) > 0) {stop("prizes must not be below 0.")}
  if(length(prizes) != length(probs)) {stop("prizes and probs must be the same length (have the same number of elements).")}
  if(warns==T) {
    if(sum(probs) > 1.00001 | sum(probs) < .99999) {warning("Probs might not total 1.0. Check probabilities total or set warns==FALSE")}
  }
  if(CI <= 0 | CI >= 1.0) {stop("CI must be between .00 and 1.00")}
  LL <- (1 - CI) / 2
  UL <- 1 - LL

  RoR.iter <- function(prizes, probs, BI, BR) {
    prizes.BI <- prizes/BI
    BR.BI <- BR/BI
    ROI.BI <- sum((prizes.BI-1)*probs)
    ifelse(ROI.BI < 0, R <- 1, R <- uniroot(function(x) sum(probs*x^prizes.BI) - x, c(0,.9999999), tol=1e-70)$root)
    RoR <- R^BR.BI
    return(RoR) 
  }

  ITM <- sum(probs[prizes > 0])
  profits <- prizes - BI
  profits.BI <- profits / BI
  ROI <- sum(profits*probs)
  ROI.BI <- sum(profits.BI*probs)
  VAR <- sum(profits^2 * probs) - ROI^2
  VAR.BI <- sum(profits.BI^2 * probs) - ROI.BI^2
  SD <- sqrt(VAR)
  SD.BI <- sqrt(VAR.BI)
  
  if(is.na(games)) {
    RoR <- exp(-2*ROI*BR/VAR)
    if(RoR > 1.0) {RoR <- 1}
    highVarRoR <- RoR.iter(prizes, probs, BI, BR)
    out.reg <- rbind(BI, ITM*100, ROI, SD, RoR*100, highVarRoR*100) 
    out.BI <- rbind(1, ITM*100, ROI.BI, SD.BI, RoR*100, highVarRoR*100)
    out <- round(cbind(out.reg, out.BI),5)
    colnames(out) <- c("Results $", "Results BI")
    rownames(out) <- c("BuyIn", "ITM %", "ROI/game", "SD", "Risk of Ruin %", "High Variance RoR %") 
    return(out)
  }
  if(!is.na(games)) {
    norm.val1 <- pnorm((-BR-ROI*games)/(SD*sqrt(games)))
    norm.val2 <- pnorm((-BR+ROI*games)/(SD*sqrt(games)))
    RoR <- exp(-2*ROI*games*BR/(VAR*games))*norm.val2 + norm.val1
  }

  if(sims!=FALSE) {
    dswings <- rep(0,sims)
    uswings <- rep(0,sims)
    final <- rep(0,sims)
    ruin.cnt <- 0
    for (i in 1:sims) {
      runvec <- sample(profits, games, prob=probs, replace=T)
      dswings[i] <- max(rle(runvec)$lengths[rle(runvec)$values<0])
      uswings[i] <- max(rle(runvec)$lengths[rle(runvec)$values>0])
      final[i] <- sum(runvec)
      ruin.cnt <- ruin.cnt + any(cumsum(runvec) + BR < BI) 
    }
  }

  if (plots > 0) {

###### Simulated Profits Plots Only

      inc <- ifelse(games < 2000,1,floor(games/1000))
      pts <- seq(1,games,inc)
      num.pts <- length(pts)
      cummat <- matrix(NA, nrow=plots, ncol=num.pts)
      for(i in 1:plots)
        cummat[i,] <- cumsum(sample(profits, games, prob=probs, replace=T))[pts]

      op <- par(mfcol=c(2,2), font.main=1)

     
      plot(pts, cumsum(rep(ROI*inc,num.pts)), type="l", lty=5, lwd=4, xlim=c(1,games),
        ylim=c(min(cummat) - 3*BI, max(cummat) + 3*BI), main=paste(plots, " Simulated Profits"),
        xlab="Tournaments", ylab="Profits ($)")
      for(i in 1:plots) {
        lines(pts, cummat[i,], col=sample(sims,sims,T))
      }
      lines(pts, cumsum(rep(ROI*inc,num.pts)), type="l", lty=5, lwd=4)
      legend("topleft", legend="Expected Line", col="black", lty=5, bty="n", lwd=4)

      plot(pts, cumsum(rep(ROI.BI*inc,num.pts)), type="l", lty=5, lwd=4, xlim=c(1,games),
        ylim=c(min(cummat/BI) - 3, max(cummat/BI) + 3), main=paste(plots, " Simulated Profits"),
        xlab="Tournaments", ylab="Profits (BuyIns)")
      for(i in 1:plots) {
        lines(pts, cummat[i,]/BI, col=sample(sims,sims,T))
      }
      lines(pts, cumsum(rep(ROI.BI*inc,num.pts)), type="l", lty=5, lwd=4)
      legend("topleft", legend="Expected Line", col="black", lty=5, bty="n", lwd=4)

  }

######

  if (plots == 0) op <- par(mfcol=c(2,1), font.main=1)

  plot(density(final), main="Distribution of Profits", ylab="Density", xlab="Profits ($)", col="green")
  polygon(density(final), col="green")

  plot(density(final/BI), main="Distribution of Profits", ylab="Density", xlab="Profits (BuyIns)", col="green")
  polygon(density(final/BI), col="green")

  out.reg <- rbind(sims, games, BI, ITM*100, ROI, SD, RoR*100, ruin.cnt/sims*100,
  min(final), quantile(final,LL), mean(final),
  median(final), quantile(final,UL), max(final),
  max(dswings), max(uswings), ((sum(final < 0))/sims)*100)
  
  out.BI <- rbind(sims, games, 1, ITM*100, ROI.BI, SD.BI, RoR*100, ruin.cnt/sims*100,
  min(final/BI), quantile(final/BI,LL), mean(final/BI),
  median(final/BI), quantile(final/BI,UL), max(final/BI),
  max(dswings), max(uswings), ((sum(final/BI < 0))/sims)*100)
  
  out <- round(cbind(out.reg, out.BI),5)
  colnames(out) <- c("Results $", "Results BI")
  rownames(out) <- c("Simulations", "Tournies per Sim", "BuyIn", "ITM %", "ROI/game", 
  "SD", "Risk of Ruin %", "Simulated Ruin %", "Worst Profit", "CI Lowerbound", "Avg. Profit", 
  "Median Profit", "CI Upperbound", "Best Profit", "Longest OOTM Streak", "Longest ITM Streak", "% Finishes with Loss")
 
  return(out)
}