################################################################
#### Examples 5.6 and 5.7: Online sports betting
################################################################

###
# Definitions
###

pvec <- c(1/3,1/3,1/3)
qvec <- c(1/2.50,1/3.25,1/2.70)
V0 <- 1
pvec/qvec


p1 <- pvec[1]
pX <- pvec[2]
p2 <- pvec[3]
q1 <- qvec[1]
qX <- qvec[2]
q2<- qvec[3]

###
# Example 5.6
###

gamvec <- seq(from = 0.1, to = 25, by= 0.1)
weights <- matrix(0,length(gamvec),3)

# Compute portfolio weights

for(k in 1:length(gamvec)){
gamma <- gamvec[k]	
w1 <- V0*q1^(1-gamma)/(q1^(1-gamma)+qX^(1-gamma)+q2^(1-gamma))
wX <- V0*qX^(1-gamma)/(q1^(1-gamma)+qX^(1-gamma)+q2^(1-gamma))
w2 <- V0*q2^(1-gamma)/(q1^(1-gamma)+qX^(1-gamma)+q2^(1-gamma))
weights[k,] <- c(w1,wX,w2)
}

# Plot portfolio weights

plot(gamvec,weights[,1],type="l",xlab="",ylab="", ylim = c(0,1))
lines(gamvec,weights[,1]+weights[,2])

###
# Synthetic risk-free asset
###

# Re-compute the weights when we use the synthetic risk free asset

weights.new <- matrix(0,length(gamvec),3)
for(k in 1:length(gamvec)){
	nofbonds <- weights[k,1]*sum(qvec)/q1
	weights.new[k,] <- c(nofbonds, weights[k,2]-nofbonds*qX/sum(qvec), weights[k,3]-nofbonds*q2/sum(qvec)) 
	}


plot(gamvec,weights.new[,1],type="l",xlab="",ylab="", ylim = c(0,1))
lines(gamvec,weights.new[,1]+weights.new[,3])

###
# Example 5.7: Allow saving money in "the pocket"
###

# Try first with both wX and w2 positive and some different gammas

gamma <- 0.1
#gamma <- 1
#gamma <- 25

lambda <- (gamma*V0)^(-1/gamma)*(pvec[1]^gamma*(1-qX-q2)^(1-gamma)+pX^gamma*qX^(1-gamma)+p2^gamma*q2^(1-gamma))^(1/gamma)

w0 <- (1/gamma)*(p1/(lambda*(1-qX-q2)))^gamma
wX <- qX*((1/gamma)*(lambda*qX/pX)^(-gamma)-w0)
w2 <- q2*((1/gamma)*(lambda*q2/p2)^(-gamma)-w0)
w0
wX
w2 # NOTE: w2 becomes negative
w0+wX+w2

# This leads us to wake w2=0

gamvec <- seq(from = 1, to = 25, by= 0.1)
weights <- matrix(0,length(gamvec),2)

for(k in 1:length(gamvec)){
gamma <- gamvec[k]	
nlambda <- (gamma*V0)^(-1/gamma)*(((1-qX)/(p1+p2))^(-gamma)*(1-qX)+(qX/pX)^(-gamma)*qX)^(1/gamma)
w0 <- (1/gamma)*(nlambda*(1-qX)/(p1+p2))^(-gamma)
wX <- qX*((1/gamma)*(nlambda*qX/pX)^(-gamma)-w0)
weights[k,] <- c(w0,wX)
}

plot(gamvec,weights[,1],type="l",xlab="",ylab="", ylim=c(0,1))