S0 = 62
sigma = 0.2
r = 0.1
K = 60

delta_t = 1/12 ## monthly timesteps

T_months = 5 ## time to expiration

## Parameters of the stock:
##
u = exp(sigma*sqrt(delta_t))
d= 1/u
R = 1 + delta_t*r ## total return over the next month
q = (R - d)/(u - d)
print(sprintf('u= %f; d= %f; R= %f; q= %f', u, d, R, q))

## Price the standard call option using a couple of methods:
##
source('../Chapter12/utils.R') ## price using a binomial tree
SSpot = make_spot_table(S0, u, d, T_months)
CGrid = evaluate_call_option(SSpot, K, q, R)
C_BT = CGrid[1, 1]

source('../Chapter13/utils.R') ## price using the Black-Scholes model
C_BS = BS_Call(S0, 0, r, K, T_months/12, sigma)
print(sprintf('C_BT= %f; C_BS= %f', C_BT, C_BS))


## Lets simulate random wa1ks according to the risk neutral probabilities:
##
run_mc_trials = function(n_MC=10000, delta_t, S0, r, K, T_months, sigma){
    ## Parameters of the stock:
    ##
    u = exp(sigma*sqrt(delta_t))
    d = 1/u
    R = 1 + delta_t*r ## total return over the next month
    q = (R - d)/(u - d)
    ##print(c(u, d, R, q))

    ##
    ## Make a function for these ca1culations ..
    ##
    set.seed(12345)
    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))

    ## What are the stock prices on each month after t=0:
    ##
    St = S0 * u^MC_cumsum

    ## The value of a European call:
    ##
    C_fv = pmax(St[, T_months] - K, 0) ## future values
    C_pv = C_fv / R^T_months ## present values

    ## What is the average spot price across time:
    ##
    St_w_S0 = cbind(rep(S0, n_MC), St) ## include the original price
    S_avg = apply(St_w_S0, 1, mean)
    
    ## The value of the Asian option:
    ##
    A_fv = pmax(S_avg - K, 0) ## future values
    A_pv = A_fv / R^T_months

    return(list(C_pv=C_pv, S_avg=S_avg, A_pv=A_pv))
}

##res = run_mc_trials(100, delta_t, S0, r, K, T_months, sigma) ## how to call the above function 


## The Monte Carlo estimates of both options:
##
for( n_MC in c(10, 100, 1000, 10000, 10^5, 10^6, 10^7) ){
    res = run_mc_trials(n_MC, delta_t, S0, r, K, T_months, sigma)

    C_pv = res$C_pv ## the Monte Carlo estimates of the premium for the European call
    A_pv = res$A_pv ## the Monte Carlo estimates of the premium for the Asian call
    
    print(sprintf('n_MC= %8d; C_MC= %f (%f); A_MC= %f (%f)', n_MC, mean(C_pv), sd(C_pv)/sqrt(n_MC), mean(A_pv), sd(A_pv) /sqrt (n_MC)))
}


## Premiums estimated using control variants:
##
for( n_MC in c(10, 100, 1000, 10000, 10^5, 10^6, 10^7) ){
    res = run_mc_trials(n_MC, delta_t, S0, r, K, T_months, sigma)

    C_pv = res$C_pv ## the Monte Carlo estimates of the premium for the European call
    A_pv = res$A_pv ## the Monte Carlo estimates of the premium for the Asian call
    S_avg = res$S_avg

    ## Part (a): We'll now use the European calls as the control variant:
    ##
    X = A_pv ## MC samples of Asian options
    Y = C_pv ## MC samples of European call
    Y_bar = C_BT ## the expected value for the European call
    a = -cov(X, Y)/var(Y)
    x_hat = mean(X) + a * (mean(Y) - Y_bar)
    print(sprintf('n_MC= %8d; Control Variate= %15s; a= %f; x_hat= %f', n_MC, 'European Call', a, x_hat))

    ## Part (b): Use the average spot price as the control variate:
    ##
    X = A_pv ## MC samples of Asian options
    Y = (S_avg - K) / R^T_months ## MC samples of the control varient= ( S_sample_average - K )/R^5
    Y_bar = ( mean(R^seq(0, T_months) * S0) - K )/R^T_months
    a = -cov(X, Y)/var(Y)
    x_hat = mean(X) + a * (mean(Y) - Y_bar)
    print(sprintf('n_MC= %8d; Control Variate= %15s; a= %f; x_hat= %f', n_MC, 'S_average', a, x_hat))
    
}