# Becoming Warren Buffett simulator
# Foster et al TAS
# Nicholas Horton, nhorton@smith.edu Thu Oct 23 18:05:57 EDT 2008

green <- c(0.8,  0.9, 1.1, 1.1, 1.2, 1.4)
red   <- c(0.05, 0.2, 1,   3,   3,   3)
white <- c(0.95, 1,   1,   1,   1,   1.1)

years <- 20
numsims <- 5000
n <- years*numsims

library(Hmisc)
process <- function(color) {
   xmat <- matrix(rMultinom(matrix(rep(1/6,6), 1 , 6), n), nrow=numsims)
   res <- rep(1000,numsims)
   for (i in 1:years) {
	res <- res*color[xmat[,i]]
   }
   return(res)
}

pink <- function(col1,col2) {
   xmat1 <- matrix(rMultinom(matrix(rep(1/6,6), 1 , 6), n), nrow=numsims)
   xmat2 <- matrix(rMultinom(matrix(rep(1/6,6), 1 , 6), n), nrow=numsims)
   res <- rep(500,numsims)
   for (i in 1:years) {
	redtmp <- res*col1[xmat1[,i]]
	whitetmp <- res*col2[xmat2[,i]]
	res <- (redtmp+whitetmp)/2
   }
   return(res*2)
}
   
      
greenres <- process(green)
redres <- process(red)
whiteres <- process(white)
pinkres <- pink(red,white)
boxplot(greenres,redres,whiteres,pinkres)

results <- data.frame(greenres,redres,whiteres,pinkres)
write.csv(results,"res.csv")