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

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

####################################################
#   Chunk 1: Antithetic sampling
####################################################

### Example 1: integral of the exponential function
# Naive MC
set.seed(55555)
X <- exp(runif(100))
T <- t.test(X)
T$estimate
CI <- as.numeric(T$conf.int)
CI
CI[2]-CI[1]
(CI[2]-CI[1])/2/(exp(1)-1)
# Antithetic sampling
set.seed(55555)
U1 <- runif(50)
U2 <- 1-U1
X <- 0.5*(exp(U1)+exp(U2))
T <- t.test(X)
T$estimate
CI <- as.numeric(T$conf.int)
CI
CI[2]-CI[1]
(CI[2]-CI[1])/2/(exp(1)-1)

### Example 2: pricing a vanilla call
S0 <- 50
K <- 50
sigma <- 0.4
r <- 0.05
T <- 5/12
numRepl <- 50000
exact <- EuVanillaCall(S0,K,T,r,sigma); exact
set.seed(55555)
out1 <- BlsMCNaive(S0,K,T,r,sigma,numRepl); out1
set.seed(55555)
out2 <- BlsMCAV(S0,K,T,r,sigma,numRepl/2); out2
# Repeat for an in-the-money call
K <- 25
exact <- EuVanillaCall(S0,K,T,r,sigma); exact
set.seed(55555)
out1 <- BlsMCNaive(S0,K,T,r,sigma,numRepl); out1
set.seed(55555)
out2 <- BlsMCAV(S0,K,T,r,sigma,numRepl/2); out2


####################################################
#   Chunk 2: Control variates
####################################################

### Example 1: Integral of the exponential function
set.seed(55555)
numPilot <- 50
U1 <- runif(numPilot)
exp1 <- exp(U1)
c <- -cov(U1,exp1)/var(U1)
numReplication <- 150
U2 <- runif(numReplication)
exp2 <- exp(U2)
Xc <- exp2 + c*(U2 - 0.5)
outCV <- t.test(Xc)
out <- t.test(exp(runif(numPilot+numReplication)))
as.numeric(out$estimate)
as.numeric(out$conf.int)
as.numeric(outCV$estimate)
as.numeric(outCV$conf.int)
cor(U1,exp1)

### Example 2: Pricing a vanilla call
S0 <- 50
K <- 50
sigma <- 0.4
r <- 0.05
T <- 5/12
numRepl <- 50000
numPilot <- 10000
exact <- EuVanillaCall(S0,K,T,r,sigma); exact
set.seed(55555)
out1 <- BlsMCNaive(S0,K,T,r,sigma,numRepl+numPilot); out1
set.seed(55555)
out2 <- BlsMCCV(S0,K,T,r,sigma,numRepl,numPilot); out2


####################################################
#   Chunk 3: Conditional Monte Carlo
####################################################

### An example from Robert and Casella
### A function of a Student's t variable 
numRepl <- 100000
df <- 5
mu <- 3
sigma <- 0.5
set.seed(55555)
# crude MC
X <- mu+sigma*rt(numRepl,df)
expX <- exp(-X^2)
out <- t.test(expX)
# Conditional MC
Y <- sqrt(rchisq(numRepl,df)/df)
aux <- 1+2*(sigma/Y)^2
expXc <- exp(-mu^2/aux)/sqrt(aux)
outC <- t.test(expXc)
as.numeric(out$estimate)
as.numeric(out$conf.int)
as.numeric(outC$estimate)
as.numeric(outC$conf.int)

plot1 <- cumsum(expX)/(1:numRepl)
plot2 <- cumsum(expXc)/(1:numRepl)
plot(plot1, type='l', ylim=c(0.006,0.009))
lines(plot2, col='red')

####################################################
#   Chunk 4: Stratified sampling
####################################################

### Example 1: estimating pi
set.seed(55555)
estPiNaiveMC(5000)
estPiStratified(5000)
estPiStratifiedAnti(5000)

### Example 2: Pricing a vanilla call
S0 <- 50; K <- 50; sigma <- 0.4; 
r <- 0.05; T <- 5/12
numRepl <- 10000
EuVanillaCall(S0,K,T,r,sigma)
set.seed(55555)
BlsMCNaive(S0,K,T,r,sigma,numRepl)
numStrats <- 25
set.seed(55555)
BsmStrat(S0,K,T,r,sigma,numRepl,numStrats)

####################################################
#   Chunk 5: Importance sampling
####################################################

### Example 1: estimating pi
pi
set.seed(55555)
estpiIS(1000,10)
estpiIS(1000,10)
estpiIS(1000,10)
estpiIS(1000,10)
estpiIS(1000,10)

### Example 2: expected value of a triangular
###            function of a normal variable
# Auxiliary triangular  function 
triangle <- function(x){
  out <- numeric(length(x))
  idx1 <- which((x>3)&(x<=4))
  idx2 <- which((x>4)&(x<5))
  out[idx1] <- 1000*(x[idx1]-3)
  out[idx2] <- -1000*(x[idx2]-5)
  return(out)
}

set.seed(55555)
z <- rnorm(10000)
out <- t.test(triangle(z))
likRatio <- function(x,theta){exp(-theta*x+theta^2/2)}
theta <- 4
x <- z+theta
outC <- t.test(triangle(x)*likRatio(x,theta))
as.numeric(out$estimate)
as.numeric(out$conf.int)
as.numeric(outC$estimate)
as.numeric(outC$conf.int)

### Example 3: Pricing a deeply out-of-the-money call
S0 <- 50; K <- 80; r <- 0.05;
sigma <- 0.4; T <- 5/12
numRepl <- 100000
numTimes <- 100
MCError <- numeric(numTimes)
MCISError <- numeric(numTimes)
truePrice <- EuVanillaCall(S0,K,T,r,sigma)
set.seed(55555)
for (k in 1:numTimes){
  MCPrice <- BlsMCNaive(S0,K,T,r,sigma,numRepl)$estimate
  MCError[k] <- abs(MCPrice - truePrice)/truePrice
}
set.seed(55555)
for (k in 1:numTimes){
  MCISPrice <- BlsMCIS(S0,K,T,r,sigma,numRepl)$estimate
  MCISError[k] <- abs(MCISPrice - truePrice)/truePrice
}
cat('Average Percentage Error:\n')
cat(sprintf(' MC    = %6.3f%%', 100*mean(MCError)), '\n')
cat(sprintf(' MC+IS = %6.3f%%', 100*mean(MCISError)), '\n')
cat('Worst Case Percentage Error:\n')
cat(sprintf(' MC    = %6.3f%%', 100*max(MCError)), '\n')
cat(sprintf(' MC+IS = %6.3f%%', 100*max(MCISError)),'\n')

