#!/usr/local/bin/Rscript

# This R code sets up eight six-prize latteries,
# with the prizes chosen randomly between +10 and -10,
# and determines which lattery a clairvoyant would choose,
# and which a PV decision maker
# and which a maxmax and a maxmin DM would choose
# with the chosen payoff.
#
# batch execute: Rscript riskmethods.r
# see
# http://www.cureffi.org/2014/01/15/running-r-batch-mode-linux/

# Written by Robert Marks, 2019. robert.marks@gmail.com

# Parameters to changed for different simulation experiments:

# CARA, gamma
gamma1 <- 0.11

# CRRA, rho
rho <- .0001

# DRP, delta and beta
dee <- 1.4 # >= 1
beta <- 0


# Next: to simulate multiple times and calculate the mean payoff.

set.seed(42)

count <- 0
c1 <- 0 # Clairvoyant
c2 <- 0 # EV
c3 <- 0 # maxmax
c4 <- 0 # maxmin
c5 <- 0 # Laplace
c6 <- 0 # random
c7 <- 0 # CARA
c8 <- 0 # CRRA
c9 <- 0 # DRP

while(count <10000) {

count <-  count + 1
# set up the eight lotteries

# first, 48 probabilities (we use prob*1000)

probl1a <- runif(6)
probl1 <- probl1a/sum(probl1a)
probl2a <- runif(6)
probl2 <- probl2a/sum(probl2a)
probl3a <- runif(6)
probl3 <- probl3a/sum(probl3a)
probl4a <- runif(6)
probl4 <- probl4a/sum(probl4a)
probl5a <- runif(6)
probl5 <- probl5a/sum(probl5a)
probl6a <- runif(6)
probl6 <- probl6a/sum(probl6a)
probl7a <- runif(6)
probl7 <- probl7a/sum(probl7a)
probl8a <- runif(6)
probl8 <- probl8a/sum(probl8a)
#cat("\n")
prob <- c(probl1, probl2, probl3, probl4, probl5, probl6, probl7, probl8)
#cat(length(prob)); cat("\n")

# second, the 48 prizes

priz1 <- runif(48, -10, 10)   # the 48 prizes of the 8 lotteries
mp <- sum(priz1)
#cat(mp/48); cat(" mean of the 48 prizes, close to 0\n")

# maxmax and maxmin
maxlot1 <- max(priz1[1:6])
maxlot2 <- max(priz1[7:12])
maxlot3 <- max(priz1[13:18])
maxlot4 <- max(priz1[19:24])
maxlot5 <- max(priz1[25:30])
maxlot6 <- max(priz1[31:36])
maxlot7 <- max(priz1[37:42])
maxlot8 <- max(priz1[43:48])
maxmaxlott <- which.max(c(maxlot1, maxlot2, maxlot3, maxlot4, maxlot5, maxlot6, maxlot7, maxlot8))
#cat(maxmaxlott); cat(" maxmaxlott\n")
minlot1 <- min(priz1[1:6])
minlot2 <- min(priz1[7:12])
minlot3 <- min(priz1[13:18])
minlot4 <- min(priz1[19:24])
minlot5 <- min(priz1[25:30])
minlot6 <- min(priz1[31:36])
minlot7 <- min(priz1[37:42])
minlot8 <- min(priz1[43:48])
maxminlott <- which.max(c(minlot1, minlot2, minlot3, minlot4, minlot5, minlot6, minlot7, minlot8))
#cat(maxminlott); cat(" maxminlott\n")

# the eight EVs

ev1 <- sum(priz1[1:6]*prob[1:6])
ev2 <- sum(priz1[7:12]*prob[7:12])
ev3 <- sum(priz1[13:18]*prob[13:18])
ev4 <- sum(priz1[19:24]*prob[19:24])
ev5 <- sum(priz1[25:30]*prob[25:30])
ev6 <- sum(priz1[31:36]*prob[31:36])
ev7 <- sum(priz1[37:42]*prob[37:42])
ev8 <- sum(priz1[43:48]*prob[43:48])

eight <- c(ev1, ev2, ev3, ev4, ev5, ev6, ev7, ev8)
#cat (eight); cat("\n")

# then the EV decision maker chooses the lottery with the hightes EV

evchoice <- which.max(eight)
#cat (evchoice); cat(" evchoice\n")

# the Laplace criterion of equal probabilities

lp1 <- sum(priz1[1:6])/6
lp2 <- sum(priz1[7:12])/6
lp3 <- sum(priz1[13:18])/6
lp4 <- sum(priz1[19:24])/6
lp5 <- sum(priz1[25:30])/6
lp6 <- sum(priz1[31:36])/6
lp7 <- sum(priz1[37:42])/6
lp8 <- sum(priz1[43:48])/6
eightlp <- c(lp1, lp2, lp3, lp4, lp5, lp6, lp7, lp8)
#cat(eightlp); cat(" eightlp\n")

lpchoice <- which.max(eightlp)
#cat (lpchoice); cat(" lpchoice\n")

#CARA rankings

# gamma1 <- 0.11
cara1 <- sum((1 - exp(- gamma1 * (priz1[1:6]+10))/gamma1)*prob[1:6])
cara2 <- sum((1 - exp(- gamma1 * (priz1[7:12]+10))/gamma1)*prob[7:12])
cara3 <- sum((1 - exp(- gamma1 * (priz1[13:18]+10))/gamma1)*prob[13:18])
cara4 <- sum((1 - exp(- gamma1 * (priz1[19:24]+10))/gamma1)*prob[19:24])
cara5 <- sum((1 - exp(- gamma1 * (priz1[25:30]+10))/gamma1)*prob[25:30])
cara6 <- sum((1 - exp(- gamma1 * (priz1[31:36]+10))/gamma1)*prob[31:36])
cara7 <- sum((1 - exp(- gamma1 * (priz1[37:42]+10))/gamma1)*prob[37:42])
cara8 <- sum((1 - exp(- gamma1 * (priz1[43:48]+10))/gamma1)*prob[43:48])

eightcara = c(cara1, cara2, cara3, cara4, cara5, cara6, cara7, cara8)
#cat(eightcara); cat("\n")
carachoice <- which.max(eightcara)
#cat(carachoice); cat(" carachoice\n")

# CRRA
# rho <- .0001
w <- (c8/count) + 10 #starting average wealth >0
if (rho == 1) {
crra1 <- sum(log(w + (priz1[1:6]+10))*prob[1:6])
crra2 <- sum(log(w + (priz1[7:12]+10))*prob[7:12])
crra3 <- sum(log(w + (priz1[13:18]+10))*prob[13:18])
crra4 <- sum(log(w + (priz1[19:24]+10))*prob[19:24])
crra5 <- sum(log(w + (priz1[25:30]+10))*prob[25:30])
crra6 <- sum(log(w + (priz1[31:36]+10))*prob[31:36])
crra7 <- sum(log(w + (priz1[37:42]+10))*prob[37:42])
crra8 <- sum(log(w + (priz1[43:48]+10))*prob[43:48])
} else {
crra1 <- sum((w + (priz1[1:6]+10))^(1 - rho)/(1 - rho)*prob[1:6])
crra2 <- sum((w + (priz1[7:12]+10))^(1 - rho)/(1 - rho)*prob[7:12])
crra3 <- sum((w + (priz1[13:18]+10))^(1 - rho)/(1 - rho)*prob[13:18])
crra4 <- sum((w + (priz1[19:24]+10))^(1 - rho)/(1 - rho)*prob[19:24])
crra5 <- sum((w + (priz1[25:30]+10))^(1 - rho)/(1 - rho)*prob[25:30])
crra6 <- sum((w + (priz1[31:36]+10))^(1 - rho)/(1 - rho)*prob[31:36])
crra7 <- sum((w + (priz1[37:42]+10))^(1 - rho)/(1 - rho)*prob[37:42])
crra8 <- sum((w + (priz1[43:48]+10))^(1 - rho)/(1 - rho)*prob[43:48])
}

eightcrra = c(crra1, crra2, crra3, crra4, crra5, crra6, crra7, crra8)
#cat(eightcrra); cat("\n")
crrachoice <- which.max(eightcrra)
#cat(crrachoice); cat(" crrachoice\n")


# DRP from Kahnemann & Tversky

# dee <- 1.4 # >= 1
# beta <- 0

#calculate the values of the DRP Vs,
#for both +ve and -ve prizes, non-zero betas
# then for zero beta

if (beta != 0) {
	drpvpos <- ((1 - exp(-beta * priz1))/(1 - exp(-100 * beta)))
	drpvneg <- -dee * ((1 - exp( beta * priz1))/(1 - exp(-100 * beta)))
	drpv <- ifelse (priz1 < 0, drpvneg, drpvpos)
} else {
	drpv <- ifelse (priz1 < 0, priz1*dee, priz1)
}

drp1 <- sum((drpv[1:6]+100)*prob[1:6])
drp2 <- sum((drpv[7:12]+100)*prob[7:12])
drp3 <- sum((drpv[13:18]+100)*prob[13:18])
drp4 <- sum((drpv[19:24]+100)*prob[19:24])
drp5 <- sum((drpv[25:30]+100)*prob[25:30])
drp6 <- sum((drpv[31:36]+100)*prob[31:36])
drp7 <- sum((drpv[37:42]+100)*prob[37:42])
drp8 <- sum((drpv[43:48]+100)*prob[43:48])

eightdrp = c(drp1, drp2, drp3, drp4, drp5, drp6, drp7, drp8)
#cat(eightdrp); cat("\n")
drpchoice <- which.max(eightdrp)
#cat(drpchoice); cat(" drpchoice\n")

#eight realisations

ra1 <- sample(6, 1, probl1, replace=T)
rl1 <- priz1[ra1]
ra2 <- sample(6, 1, probl2, replace=T)
rl2 <- priz1[ra2+6]
ra3 <- sample(6, 1, probl3, replace=T)
rl3 <- priz1[ra3+12]
ra4 <- sample(6, 1, probl4, replace=T)
rl4 <- priz1[ra4+18]
ra5 <- sample(6, 1, probl5, replace=T)
rl5 <- priz1[ra5+24]
ra6 <- sample(6, 1, probl6, replace=T)
rl6 <- priz1[ra6+30]
ra7 <- sample(6, 1, probl7, replace=T)
rl7 <- priz1[ra7+36]
ra8 <- sample(6, 1, probl8, replace=T)
rl8 <- priz1[ra8+42]

#actions to values needed
raa <- c(ra1, ra2, ra3, ra4, ra5, ra6, ra7, ra8)
#cat(raa); cat("\n")
#realised lottery payoffs
rla <- c(rl1, rl2, rl3, rl4, rl5, rl6, rl7, rl8)
#cat(rla); cat("\n")
clpayoff <- max(rla) # the clairvoyant's payoff
which.clch <- which.max(rla)
#cat ("Claivoyant payoff "); cat (clpayoff, which.clch); cat("\n")
c1 <- c1 + clpayoff

evpayoff <- rla[evchoice]
#cat ("EV payoff "); cat (evpayoff); cat("\n")
c2 <- c2 + evpayoff

maxmaxpayoff <- rla[maxmaxlott]
#cat ("maxmax payoff "); cat (maxmaxpayoff); cat("\n")
c3 <- c3 + maxmaxpayoff
maxminpayoff <- rla[maxminlott]
#cat ("maxmin payoff "); cat (maxminpayoff); cat("\n")
c4 <- c4 + maxminpayoff

lppayoff <- rla[lpchoice]
c5 <- c5 + lppayoff

chrand <- sample(1:8, 1)
c6 <- c6 + rla[chrand]

carapayoff <- rla[carachoice]
c7 <- c7 + carapayoff

crrapayoff <- rla[crrachoice]
#cat ("CRRApayoff "); cat (crrapayoff); cat("\n")
c8 <- c8 + crrapayoff
#cat ("CRRApayoff "); cat (c8); cat("\n")

drppayoff <- rla[drpchoice]
c9 <- c9 + drppayoff

}

cat ("Claivoyant payoff "); cat (c1/count); cat("\n")
cat ("EV payoff "); cat (c2/count); cat(" i.e.%Cl "); cat (c2/c1 * 100); cat("\n")
cat ("maxmax payoff "); cat (c3/count);cat(" i.e.%Cl "); cat (c3/c1 * 100);  cat("\n")
cat ("maxmin payoff "); cat (c4/count); cat(" i.e.%Cl "); cat (c4/c1 * 100); cat("\n")
cat ("Laplace payoff "); cat (c5/count); cat(" i.e.%Cl "); cat (c5/c1 * 100); cat("\n")
cat("gamma = "); cat(gamma1)
cat ("CARA payoff "); cat (c7/count); cat(" i.e.%Cl "); cat (c7/c1 * 100);
cat(" i.e.%EV "); cat (c7/c2 *100); cat("\n")
cat("rho = "); cat(rho); cat("\n")
cat ("CRRA payoff "); cat (c8/count); cat(" i.e.%Cl "); cat (c8/c1 * 100);
cat(" i.e.%EV "); cat (c8/c2 *100); cat("\n")
cat("beta = "); cat(beta); cat(" dee = "); cat(dee); cat("\n")
cat ("DRP payoff "); cat (c9/count); cat(" i.e.%Cl "); cat (c9/c1 * 100);
cat(" i.e.%EV "); cat (c9/c2 *100); cat("\n")
cat ("random payoff "); cat (c6/count); cat("\n")
cat ("iterations = "); cat (count); cat("\n")