present_value_of_CD = function(r, u, d, n_MC, S0, T){
    ##
    ## Computes the present value of the Great Western CD
    ##
    T_months = 12*T ## time to expiration in months

    ## What is the risk neutra1 probabiilty of each Monte Carlo path:
    ##
    R = 1 + r/12
    q = (R - d)/(u - d)

    MC = sample(c(+1, -1), T_months*n_MC, prob=c(q, 1-q), replace=TRUE) ## 1 -> an up; -1 -> a down
    MC = matrix(MC, nrow=n_MC, ncol=T_months)

    ## Compute the cumulative sum along the rows of MC (this will compute the value of n_u - n_d):
    ##
    MC_cumsum = t(apply(MC, 1, cumsum))
    ##MC_n_up = t(apply(MC>0, 1, cumsum)) # the number of up moves
    ##MC_n_down = t(apply(MC<0, 1, cumsum)) # the number of down moves

    ## What are the SNP prices at each time:
    ##
    St = S0 * u^MC_cumsum

    Account_Balance = rep(S0, n_MC) ## how much money is in the account

    ## What are the three yearly averages of the SNP:
    ##
    Px_0 = rep(S0, n_MC) ## the index price at the start of the first year
    year_1_A = apply(St[, 1:12], 1, mean) ## year 1 values of A
    R_1 = (year_1_A - Px_0)/Px_0
    I_1 = pmax(0, R_1)
    Account_Balance = Account_Balance + Account_Balance*I_1
    
    Px_0 = St[, 12] ## the index price at the start of the second year
    year_2_A = apply(St[, 12 + (1:12)], 1, mean)
    R_2 = (year_2_A - Px_0)/Px_0
    I_2 = pmax(0, R_2)
    Account_Balance = Account_Balance + Account_Balance*I_2
    
    Px_0 = St[, 24] ## the index price at the start of the third year
    year_3_A = apply(St[, 24 + (1:12)], 1, mean)
    R_3 = (year_3_A - Px_0)/Px_0
    I_3 = pmax(0, R_3)
    Account_Balance = Account_Balance + Account_Balance*I_3 ## this is the final account Balance

    E_PT = mean(Account_Balance) ## this is \hat{E}[P(T)]

    return(exp(-r*T) * E_PT) ## discount
}


##
## Using the above function lets try to solve this exercise:
##
set.seed(12345)

sigma = 0.2
S0 = 1.0

T = 3

## The risk neutra1 parameters of the stock:
##
delta_t = 1/12 ## monthly timesteps
u = exp(sigma*sqrt(delta_t))
d = 1/u
print(sprintf('u= %f; d= %f', u, d))

## present_value_of_CD(0.05, u, d, n_MC, S0, T) ## How to call the above function

n_MC = 1000 ## the number of Monte Carlo
n_rs = 100
rs = seq(0.05, 0.09, length.out=n_rs)
PV = rep(NA, n_rs)
for( ri in 1:n_rs ){
    PV[ri] = present_value_of_CD(rs[ri], u, d, n_MC, S0, T)
}

##postscript('../../WriteUp/Graphics/Chapter13/chap_13_ex_8_PVs.eps', onefile=FALSE, horizontal=FALSE)
plot(rs, PV, xlab='r (interest rate)', ylab='Great Western CD Present Value')
abline(h=1.0, col='green')
grid()
##dev.off()