
####################################################
#    Function definitions for Chapter 8


####################################################
# AUXILIARY functions for comparison in option pricing
# exact BSM formula
EuVanillaCall <- function(S0,K,T,r,sigma){
  d1 <- (log(S0/K)+(r+sigma^2/2)*T)/(sqrt(T)*sigma)
  d2 <- d1-sqrt(T)*sigma
  value <- S0*pnorm(d1)-K*exp(-r*T)*pnorm(d2)
  return(value)  
}
# Naive MC
BlsMCNaive <- function(S0,K,T,r,sigma,numRepl,c.level=0.95){
  nuT <- (r - 0.5*sigma^2)*T
  sigmaT <- sigma*sqrt(T)
  ST <- S0*exp(rnorm(numRepl, nuT, sigmaT))
  dpayoff <- exp(-r*T)*pmax(0,ST-K)
  aux <- t.test(dpayoff,conf.level=c.level)
  value <- as.numeric(aux$estimate)
  ci <- as.numeric(aux$conf.int)
  return(list(estimate=value,conf.int=ci))
}


################################################
# Pricing a vanilla call option with antithetic sampling
BlsMCAV <- function(S0,K,T,r,sigma,numPairs,c.level=0.95){
  nuT <- (r - 0.5*sigma^2)*T
  sigmaT <- sigma*sqrt(T)
  veteps <- rnorm(numPairs)
  payoff1 <- pmax(0,S0*exp(nuT+sigmaT*veteps)-K)
  payoff2 <- pmax(0,S0*exp(nuT-sigmaT*veteps)-K)
  dpayoff <- exp(-r*T)*(payoff1+payoff2)/2  
  aux <- t.test(dpayoff,conf.level=c.level)
  value <- as.numeric(aux$estimate)
  ci <- as.numeric(aux$conf.int)
  return(list(estimate=value,conf.int=ci))
}


################################################
# Pricing a vanilla call option with control variates
BlsMCCV <- function(S0,K,T,r,sigma,numRepl,numPilot,
                    c.level=0.95){
  # invariant quantities
  nuT <- (r - 0.5*sigma^2)*T
  sigmaT <- sigma*sqrt(T)
  # pilot replication first 
  ST <- S0*exp(rnorm(numPilot,nuT,sigmaT))
  optVals <- exp(-r*T) * pmax( 0 , ST-K)
  expOptVal <- S0 * exp(r*T)
  varOptVal <- S0^2 * exp(2*r*T) * (exp(T * sigma^2) - 1)
  cstar <- -cov(ST,optVals)/varOptVal
  # now run MC with controlled estimator
  newST <- S0*exp(rnorm(numRepl,nuT,sigmaT))
  newOptVals <- exp(-r*T) * pmax( 0 , newST-K)
  controlledEst <- newOptVals+cstar*(newST-expOptVal)
  aux <- t.test(controlledEst,conf.level=c.level)
  value <- as.numeric(aux$estimate)
  ci <- as.numeric(aux$conf.int)
  return(list(estimate=value,conf.int=ci))
}

################################################
# Pricing a vanilla call option with statification
BsmStrat <- function(S0,K,T,r,sigma,numRepl,numStrats,
                     c.level=0.95){
  nuT <- (r - 0.5*sigma^2)*T
  sigmaT <- sigma*sqrt(T)
  numPerStrat <- ceiling(numRepl/numStrats)
  meanStrat <- numeric(numStrats)
  varStrat <- numeric(numStrats)
  for (k in 1:numStrats) {
    U <- (k-1)/numStrats+runif(numPerStrat)/numStrats
    ST <- S0*exp(nuT+sigmaT*qnorm(U))
    dPayoff <- exp(-r*T)*pmax(0,ST-K)
    meanStrat[k] <- mean(dPayoff)
    varStrat[k] <- var(dPayoff)
  }
  meanTot <- mean(meanStrat)
  varTot <- mean(varStrat)
  alpha <- (100-c.level)/100
  half <- qnorm(1-alpha/2)*sqrt(varTot)
  ci <- c(meanTot-half,meanTot+half)
  return(list(estimate=meanTot,conf.int=ci))
}

################################################
# Pricing a vanilla call option with importance sampling
BlsMCIS <- function(S0,K,T,r,sigma,numRepl,c.level=0.95){
  nuT <- (r - 0.5*sigma^2)*T
  sigmaT <- sigma*sqrt(T)
  ISnuT <- log(K/S0) - 0.5*sigma^2*T
  norms <- rnorm(numRepl,ISnuT,sigmaT)
  ISRatios <- exp( (2*(nuT - ISnuT)*norms - nuT^2 + 
                      ISnuT^2)/2/sigmaT^2)
  dpayoff <- exp(-r*T)*pmax(0, (S0*exp(norms)-K))
  aux <- t.test(dpayoff*ISRatios,conf.level=c.level)
  value <- as.numeric(aux$estimate)
  ci <- as.numeric(aux$conf.int)
  return(list(estimate=value,conf.int=ci))
}


################################################
# Estimating pi by random sampling
# Naive Monte Carlo
estPiNaiveMC <- function(numRepl, c.level=0.95){
  U <- runif(numRepl)
  X <- 4*sqrt(1-U^2)
  aux <- t.test(X,conf.level=c.level)
  value <- as.numeric(aux$estimate)
  ci <- as.numeric(aux$conf.int)
  return(list(estimate=value,conf.int=ci))
}

# Stratified sampling
estPiStratified <- function(numRepl, c.level=0.95){
  U <- runif(numRepl)
  jVet <- (1:numRepl)
  X <- 4*sqrt(1-((U+jVet-1)/numRepl)^2)
  aux <- t.test(X,conf.level=c.level)
  value <- as.numeric(aux$estimate)
  ci <- as.numeric(aux$conf.int)
  return(list(estimate=value,conf.int=ci))
}

# Stratified and antithetic sampling 
estPiStratifiedAnti <- function(numRepl, c.level=0.95){
  U <- runif(numRepl)
  jVet <- (1:numRepl)
  X <- 2*(sqrt(1-((U+jVet-1)/numRepl)^2) + 
            sqrt(1-((jVet-U)/numRepl)^2))
  aux <- t.test(X,conf.level=c.level)
  value <- as.numeric(aux$estimate)
  ci <- as.numeric(aux$conf.int)
  return(list(estimate=value,conf.int=ci))
}

# Importance sampling 
estpiIS <- function(m,L){
  # define left end points of subintervals
  s <- seq(from=0, by=1/L, to=1-1/L) + 1/(2*L)
  hvals <- sqrt(1 - s^2)
  # get cumulative probabilities
  cs <- cumsum(hvals)
  est <- numeric(m)
  for (j in 1:m){
    # locate subinterval
    loc <- sum(runif(1)*cs[L] > cs) + 1
    # sample uniformly within subinterval
    x <- (loc-1)/L + runif(1)/L
    p <- hvals[loc]/cs[L]
    est[j] <- sqrt(1 - x^2)/(p*L)
  }
  z <- 4*sum(est)/m
  return(z)
}
