dl_p_min_fn = function(x, pA1, pB1, UA, DA, UB, DB, sigmaAB){
    ##
    ## While not needed (since the system is linear), we can also solve the system in the section
    ## "The Double Lattice" using minimization.  Towards that end this is a function we could minimize
    ## to finding solutions to p_11, p_12, p_21, and p_22.
    ##
    p11 = x[1]
    p12 = x[2]
    p21 = x[3]
    p22 = x[4]

    pA2 = 1-pA1
    pB2 = 1-pB1

    r1 = (p11 + p12) - pA1 ## the residuals 
    r2 = (p21 + p22) - pA2
    r3 = (p11 + p21) - pB1
    r4 = ((p11 - pA1*pB1) * UA * UB + (p12 - pA1*pB2) * UA * DB + (p21 - pA2*pB1) * DA * UB + (p22 - pA2*pB2) * DA * DB) - sigmaAB

    err = mean(c(r1, r2, r3, r4)^2)
}


the_double_lattice = function(pA1, pB1, uA, dA, uB, dB, rhoAB){
    ##
    ## Solves the equations for p_11, p_12, p_21, and p_22 given in the section "The Double Lattice"
    ##
    UA = log(uA)
    DA = log(dA)
    UB = log(uB)
    DB = log(dB)

    E_Log_SA = pA1*UA + (1-pA1)*DA ## Compute E(ln(S_A))
    E_Log_SB = pB1*UB + (1-pB1)*DB

    sigmaA2 = pA1*UA^2 + (1-pA1)*DA^2 - E_Log_SA^2 ## Compute Var(ln(S_A))
    sigmaB2 = pB1*UB^2 + (1-pB1)*DB^2 - E_Log_SB^2

    sigmaA = sqrt(sigmaA2)
    sigmaB = sqrt(sigmaB2)
    
    sigmaAB = sigmaA * sigmaB * rhoAB ## Compute sigma_AB
    
    ##print(c(E_Log_SA, E_Log_SB, sigmaA2, sigmaB2, sigmaAB))

    ## Set up and solve the linear system:
    ##
    A = matrix(data=c(1, 1, 0, 0, 0, 0, 1, 1, 1, 0, 1, 0, UA*UB, UA*DB, DA*UB, DA*DB), nrow=4, ncol=4, byrow=TRUE)
    b = c(pA1, 1-pA1, pB1, sigmaAB + (pA1*UA + (1-pA1)*DA)*(pB1*UB + (1-pB1)*DB))
    x0 = solve(A, b)

    if( FALSE ){
        ##
        ## If we seek the solution via minimization (not needed).
        ## 

        ## Test our optimization function: 
        ##print(dl_p_min_fn(x0, pA1, pB1, UA, DA, UB, DB, sigmaAB))

        x0 = c(0.25, 0.25, 0.25, 0.25) ## the inital guess at p11, p12, p21, and p22
        res = optim(x0,
                    function(x) dl_p_min_fn(x, pA1, pB1, UA, DA, UB, DB, sigmaAB),
                    method='L-BFGS-B',
                    lower=c(0, 0, 0, 0),
                    upper=c(1, 1, 1, 1))
        x0 = res$par
    }

    return(x0)
}


dl_q_min_fn = function(x, qA, qB, real_prob_ratio){
    ##
    ## Uses the ratio theorem to compute q_{ij} given the risk-neutral probabilities of the individual
    ## assets qA and qB along with the ratio (p_11*p_22)/(p_12*p_21).
    ##
    q11 = x[1]
    q12 = x[2]
    q21 = x[3]
    q22 = x[4]

    r1 = (q11 + q12) - qA ## the residuals 
    r2 = (q11 + q21) - qB
    r3 = (q11 + q12 + q21 + q22) - 1.0
    r4 = (q11*q22)/(q12*q21) - real_prob_ratio

    err = mean(c(r1, r2, r3, r4)^2)
}


dup_Example_16_3 = function(){
    ##
    ## This code will verify the codes above by duplicating the numbers found in Example 16.3.
    ##
    pA1 = 0.6
    pB1 = 0.6
    uA = 1.3
    dA = 0.9
    uB = uA
    dB = dA
    rhoAB = 0.3
    
    x = the_double_lattice(pA1, pB1, uA, dA, uB, dB, rhoAB)
    print('Probabilities:')
    print(x)

    p11 = x[1] ## lets see how well this solution satisifies the equations
    p12 = x[2]
    p21 = x[3]
    p22 = x[4]

    pA2 = 1-pA1
    pB2 = 1-pB1

    UA = log(uA)
    DA = log(dA)
    UB = log(uB)
    DB = log(dB)

    E_Log_SA = pA1*UA + (1-pA1)*DA ## Compute E(ln(S_A))
    E_Log_SB = pB1*UB + (1-pB1)*DB

    sigmaA2 = pA1*UA^2 + (1-pA1)*DA^2 - E_Log_SA^2 ## Compute Var(ln(S_A))
    sigmaB2 = pB1*UB^2 + (1-pB1)*DB^2 - E_Log_SB^2

    sigmaA = sqrt(sigmaA2)
    sigmaB = sqrt(sigmaB2)
    
    sigmaAB = sigmaA * sigmaB * rhoAB ## Compute sigma_AB
    
    r1 = (p11 + p12) - pA1 ## what are the residuals 
    r2 = (p21 + p22) - pA2
    r3 = (p11 + p21) - pB1
    r4 = ((p11 - pA1*pB1) * UA * UB + (p12 - pA1*pB2) * UA * DB + (p21 - pA2*pB1) * DA * UB + (p22 - pA2*pB2) * DA * DB) - sigmaAB
    
    print('Residuals:')
    print(c(r1, r2, r3, r4))
}
## dup_Example_16_3()


exponential_certainty_equivalent = function(a, p, X1, X2){
    ##
    ## Computes the certainty equivalent for an exponetial utility function: 
    ## okular ../EBook/Investment_Science.djvu -p 480 & 
    ##
    return( - log( p * exp(-a*X1) + (1-p) * exp(-a*X2) ) / a )
}