####################################################
#    Scripts for Chapter 11 of
#    P. Brandimarte
#    Handbook in Monte Carlo Simulation
#    Applications in Financial Engineering,
#    Risk Management, and Economics
#    Wiley, 2014

#  Required packages
require(randtoolbox)

#  Source the function file
source("Chapter11_Functs.R")

####################################################
#   Chunk 1: European-style exchange option
####################################################

### Experiment 1: Compare naive MC and exact price
V0 <- 50; U0 <- 60; sigmaV <- 0.3
sigmaU <- 0.4; rho <- 0.7; T <- 5/12
r <- 0.05; K <- 0
EuExchange(V0,U0,sigmaV,sigmaU,rho,T,r)
numRepl <- 200000
set.seed(55555)
EuSpreadNaiveMC(V0,U0,0,sigmaV,sigmaU,rho,T,r,numRepl)

### Experiment 2: Compare naive MC and single/double stratification
V0 <- 50; U0 <- 45; sigmaV <- 0.2
sigmaU <- 0.3; rho <- 0.5; T <- 1
r <- 0.05; K <- 0
exactPrice <- EuExchange(V0,U0,sigmaV,sigmaU,rho,T,r)
numRepl <- 10000
set.seed(55555)
out <- EuSpreadNaiveMC(V0,U0,K,sigmaV,sigmaU,rho,T,r,numRepl)
numStrats <- 25
out1 <- EuSpreadStrat1(V0,U0,K,sigmaV,sigmaU,rho,T,r,
                 numRepl,numStrats)
numStrats <- 5
out2 <- EuSpreadStrat2(V0,U0,K,sigmaV,sigmaU,rho,T,r,
                 numRepl,numStrats)
buf <- sprintf('Exact price = %f\n', exactPrice)
cat(buf)
p <- out$estimate; ci <- out$conf.int
buf <- sprintf('Naive MC: Price = %f, CI = (%f, %f) [%f]\n', 
               p,ci[1],ci[2],(ci[2]-ci[1])/p);
cat(buf)
p <- out1$estimate; ci <- out1$conf.int
buf <- sprintf('Strat1:   Price = %f, CI = (%f, %f) [%f]\n', 
               p,ci[1],ci[2],(ci[2]-ci[1])/p);
cat(buf)
p <- out2$estimate; ci <- out2$conf.int
buf <- sprintf('Strat2:   Price = %f, CI = (%f, %f) [%f]\n', 
               p,ci[1],ci[2],(ci[2]-ci[1])/p);
cat(buf)


####################################################
#   Chunk 2: European-style barrier options
####################################################

### Experiment 1: Exact down-and-out put option price
###     with continuous barrier monitoring and approximation
###     for discrete monitoring
S0 <- 50; K <- 45; Sb <- 40
r <- 0.05; sigma <- 0.4; T <- 2/12
numSteps <- 60; numDays <- 360
EuVanillaPut(S0,K,T,r,sigma)
EuDOPut(S0,K,T,r,sigma,Sb)
Sbc <- Sb*exp(-0.5826*sigma*sqrt(1/numDays))
EuDOPut(S0,K,T,r,sigma,Sbc)

### Experiment 2: Down-and-out put option with naive MC
EuDOPut(S0,K,T,r,sigma,Sbc)
numRepl <- 50000
set.seed(55555)  
DownOutPutNaive(S0,K,T,r,sigma,Sb,numSteps,numRepl)

### Experiment 3: Down-and-in option with importance sampling
S0 <- 50; K <- 50; Sb <- 40; r <- 0.05; sigma <- 0.4; 
T <- 5/12; numSteps <- 5*30; numDays <- 360
Sbc <- Sb*exp(-0.5826*sigma*sqrt(1/numDays))
numRepl <- 50000
# first no correction to drift -> naive MC
set.seed(55555)
DownInPutIS(S0,K,T,r,sigma,Sb,numSteps,numRepl,0)
set.seed(55555)
DownInPutIS(S0,K,T,r,sigma,Sb,numSteps,numRepl,0.8)


####################################################
#   Chunk 3: European-style Asian options
####################################################

### Experiment 1: Naive MC (and comparison with vanilla)
S0 <- 50; K <- 50; r <- 0.05; sigma <- 0.4
T <- 1; numSamples <- 12; numRepl <- 50000
EuVanillaCall(S0,K,T,r,sigma)
set.seed(55555)
EuAsianNaive(S0,K,T,r,sigma,numSamples,numRepl)

### Experiment 2: Trying different control variates
S0 <- 50; K <- 50; r <- 0.05; sigma <- 0.4
T <- 1; numSamples <- 12; numRepl <- 50000
set.seed(55555)
EuAsianNaive(S0,K,T,r,sigma,numSamples,numRepl)
numPilot <- 5000; numRepl <- 45000
# CV is average price
alpha <- exp(r*T/numSamples)
expCV1 <- S0*alpha*(1-alpha^numSamples)/
  (numSamples*(1-alpha))
SampleCV1 <- function(path) mean(path[-1])
set.seed(55555)
EuAsianCV(S0,K,T,r,sigma,numSamples,numRepl,
          SampleCV1,expCV1,numPilot)
# CV is vanilla call
expCV2 <- EuVanillaCall(S0,K,T,r,sigma)
SampleCV2 <- function(path) exp(-r*T)*
  max(path[length(path)]-K,0)
set.seed(55555)
EuAsianCV(S0,K,T,r,sigma,numSamples,numRepl,
          SampleCV2,expCV2,numPilot)
# CV is geometric
SampleCV3 <- function(path) exp(-r*T)*
  max(0,(prod(path[-1]))^(1/numSamples) - K)
expCV3 <- EuAsianGeometric(S0,K,T,r,sigma,numSamples)
set.seed(55555)
EuAsianCV(S0,K,T,r,sigma,numSamples,numRepl,
          SampleCV3,expCV3,numPilot)

### Experiment 3: Checking correlations with different CVs
set.seed(55555)
EuAsianCompareCorr(S0,K,T,r,sigma,numSamples,50000)

### Experiment 4: Using low-discrepancy sequences
S0 <- 50; K <- 50; r <- 0.05; sigma <- 0.4
T <- 1; numSamples <- 12; numRepl <- 10000
EuAsianLOWD(S0,K,T,r,sigma,numSamples,
                           numRepl,"halton")
EuAsianLOWD(S0,K,T,r,sigma,numSamples,
                           numRepl,"sobol",0)
EuAsianLOWD(S0,K,T,r,sigma,numSamples,
                           numRepl,"sobol",1)
EuAsianLOWD(S0,K,T,r,sigma,numSamples,
                           numRepl,"sobol",2)

####################################################
#   Chunk 4: Pricing American-style options by 
#            linear regression and ADP
####################################################

### Replicating the original Longstaff-Schwartz example
paths <- cbind(1, c(1.09,1.16,1.22,.93,1.11,.76,.92,.88),
               c(1.08,1.26,1.07,.97,1.56,.77,.84,1.22),
               c(1.34,1.54,1.03,.92,1.52,.9,1.01,1.34))         
putF <- function(x) max(0,1.1-x)
r <- 0.06
T <- 3
b1 <- function(x) rep(1,length(x))
b2 <- function(x) x
b3 <- function(x) x^2
basisFuns <- list(b1,b2,b3)
genericLS(paths,T,r,putF,basisFuns)

####################################################
#   Chunk 5: Options with stochastic volatility (Heston model)
####################################################

### Pricing a call option
S0 <- 50; K <- 50; T <- 1; r <- 0.03; V0 <- 0.07
alpha <- 1; Vbar <- 0.05; xi <- 0.1; rho <- 0.5
numRepl <- 10000; numSteps <- 50
set.seed(55555)
CallStochVol(S0,K,T,r,V0,alpha,Vbar,xi,rho,numRepl,numSteps)

####################################################
#   Chunk 6: Bond and bond options (Vasicek model)
####################################################

### Experiment 1: Compare exact ZCB price under Vasicek against MC
r0 <- 0.02; rbar <- 0.06; gamma <- 0.45; sigma <- 0.03
T <- 1; F <- 100; numRepl <- 1000; numSteps <- 200
ZCBVasicek(r0,rbar,gamma,sigma,T,F)
set.seed(55555)
ZCBVasicekMC(r0,rbar,gamma,sigma,T,mumRepl,numSteps,F)


### Experiment 2: Compare exact call option price 
###               on a ZCB under Vasicek against MC
F <- 100; K <- 80; TO <- 1; TB <- 5; rbar <- 0.06
gamma <- 0.46; sigma <- 0.025; r0 <- 0.016
numRepl <- 10000; numSteps <- 200
set.seed(55555)
ZCBCallVasicek(r0,rbar,gamma,sigma,TB,K,TO,F)
ZCBCallVasicekMC(r0,rbar,gamma,sigma,TB,K,TO,
                 numRepl,numSteps,F)


####################################################
#   Chunk 7: Bond and bond options (CIR model)
####################################################

### Compare exact ZCB price under CIR with MC
r0 <- 0.02
rbar <- 0.06
gamma <- 0.45
alpha <- 0.03^2
T <- 1
F <- 100
numRepl <- 1000
numSteps <- 200
ZCBCIR(r0,rbar,gamma,alpha,T,F)
ZCBCIRMC(r0,rbar,gamma,alpha,T,mumRepl,numSteps,F)
