source('../Chapter13/utils.R')

mu = 0.12 ## annual
sigma = 0.2 ## annual
u = 1.1
d = 1/u
deltaT = 1/12

if( T ){
    ## solve "exactly":
    probs = trinomial_lattice_probabilities(u, mu, sigma, deltaT)
    print(probs)
}else{
    ## Use the books rounded numbers: 
    p1 = 0.228
    p2 = 0.632
    p3 = 0.140
    probs = c(p1, p2, p3)
}

r = 0.1 ## per year rate of return 
R0 = 1 + r*deltaT ## per month total return

duplicate_results_from_Example_16_2 = FALSE ## duplicates the alpha claimed in Example 16.2 else solve this exercise

max_fn = function(x){
    alpha = x[1]

    U_arg = alpha * c(u, 1, d) + (1-alpha) * R0
    U_arg[U_arg <= 0] = 1.e-6

    if( duplicate_results_from_Example_16_2 ){
        EU = sum(log(U_arg) * probs)
    }else{
        EU = sum(sqrt(U_arg) * probs)
    }

    return(EU)
}
## Test our function once: 
## print(max_fn(c(0.5)))

## Find the maximum with an optimization routine:
##
res = optim(c(0.5), max_fn, method='BFGS', control=list(fnscale=-1))
print(sprintf('alpha= %f; E(U)= %f', res$par, res$value))

## Verify that maximum:
##
alphas = seq(0, 2, length.out=100)
EUs = sapply(alphas, max_fn)
plot(alphas, EUs, type='l', xlab='alpha', ylab='E(U)')
abline(v=res$par, col='red')
grid()

## Compute the risk neutral probabilities:
##
alpha = res$par
U_arg = alpha * c(u, 1, d) + (1-alpha) * R0
if( duplicate_results_from_Example_16_2 ){
    qs = ( 1/U_arg ) * probs
}else{
    qs = ( 1/(2*sqrt(U_arg)) ) * probs
}
qs = qs/sum(qs) ## the risk neutral probabilities
print('Risk neutral probabilities:')
print(qs)