
####################################################
#    Function definitions for Chapter 12


####################################################
# Exact delta for a vanilla call option
VanillaCallDelta <- function(S0,K,T,r,sigma){
  d1 <- (log(S0/K)+(r+sigma^2/2)*T)/(sqrt(T)*sigma)
  return(pnorm(d1))  
}

####################################################
# NAIVE (very) estimation of delta by MC
VanillaDeltaMCNaive <- function(S0,K,T,r,sigma,dS,numRepl,
                                c.level=0.95){
  nuT <- (r - 0.5*sigma^2)*T
  sigmaT <- sigma*sqrt(T)
  payoff1 <- pmax(0,S0*exp(rnorm(numRepl,nuT,sigmaT))-K)
  payoff2 <- pmax(0,(S0+dS)*
                    exp(rnorm(numRepl, nuT, sigmaT))-K)
  SampleDiff <- exp(-r*T)*(payoff2 - payoff1)/dS
  aux <- t.test(SampleDiff,conf.level=c.level)
  value <- as.numeric(aux$estimate)
  ci <- as.numeric(aux$conf.int)
  return(list(estimate=value,conf.int=ci))
}

####################################################
# Using common random numbers in the estimation of delta 
VanillaDeltaMC <- function(S0,K,T,r,sigma,dS,numRepl,
                           c.level=0.95){
  nuT <- (r - 0.5*sigma^2)*T
  sigmaT <- sigma*sqrt(T)
  S1 <- S0 - dS
  S2 <- S0 + dS
  Veps <- rnorm(numRepl)
  payoff1 <- pmax(0,S1*exp(nuT+sigmaT*Veps)-K)
  payoff2 <- pmax(0,S2*exp(nuT+sigmaT*Veps)-K)
  SampleDiff <- exp(-r*T)*(payoff2 - payoff1)/(2*dS)
  aux <- t.test(SampleDiff,conf.level=c.level)
  value <- as.numeric(aux$estimate)
  ci <- as.numeric(aux$conf.int)
  return(list(estimate=value,conf.int=ci))
}

####################################################
# Delta of an Asian option by central differences
AsianDeltaCentral <- function(S0,K,r,T,sigma,dS,numRepl,
                              numSamples,c.level=0.95){
  df <- exp(-r*T)
  dt <- T/numSamples
  nuT <- (r-sigma^2/2)*dt
  sigmaT <- sqrt(dt)*sigma
  diffs <- numeric(numRepl)
  for (k in 1:numRepl){
    # store normals for CRN
    normVet <- rnorm(numSamples,nuT,sigmaT)    
    logIncrements1 <- c(log(S0-ds),normVet)
    path1 <- exp(cumsum(logIncrements1))
    logIncrements2 <- c(log(S0+ds),normVet)
    path2 <- exp(cumsum(logIncrements2))
    payoff1 <- max(0,mean(path1[-1]) - K)
    payoff2 <- max(0,mean(path2[-1]) - K)
    diffs[k] <- df*(payoff2 - payoff1)/(2*ds)
  }
  aux <- t.test(diffs,conf.level=c.level)
  value <- as.numeric(aux$estimate)
  ci <- as.numeric(aux$conf.int)
  return(list(estimate=value,conf.int=ci))
}

####################################################
# AUXILIARY functions for geometric Asian control variates
# and path generation

# scalar Geometric Brownian Motions (vectorized version)
simvGBM <- function (S0,mu,sigma,T,numSteps,numRepl){
  dt <- T/numSteps
  nuT <- (mu-sigma^2/2)*dt
  sigmaT <- sqrt(dt)*sigma
  normMatrix <- matrix(rnorm(numSteps*numRepl,nuT,sigmaT),nrow=numRepl)
  logIncrements <- cbind(log(S0),normMatrix)
  logPath <- t(apply(logIncrements,1,cumsum))
  return(exp(logPath))
}

# Geometric average Asian option
EuAsianGeometric <- function(S0,K,T,r,sigma,numSamples,q=0){
  dT <- T/numSamples
  nu <- r - q - sigma^2/2
  a <- log(S0)+nu*dT+0.5*nu*(T-dT)
  b <- sigma^2*dT + sigma^2*(T-dT)*(2*numSamples-1)/6/numSamples
  x <- (a-log(K)+b)/sqrt(b)
  value <- exp(-r*T)*(exp(a+b/2)*pnorm(x) - K*pnorm(x-sqrt(b)))
  return(value)  
}

# Using Geometric option as a control variate
EuAsianCVGeom <- function(S0,K,T,r,q=0,sigma,numSamples,numRepl,numPilot,c.level=0.95){
  # precompute discount factor and expected value of control variate
  df <- exp(-r*T)
  geomExact <- EuAsianGeometric(S0,K,T,r,sigma,numSamples,q)
  # pilot replications to set control parameter
  geomPrices <- numeric(numPilot)
  arithPrices <- numeric(numPilot)
  for (k in 1:numPilot){
    path <- simvGBM(S0,r-q,sigma,T,numSamples,1)
    geomPrices[k] <- df*max(0,(prod(path[-1]))^(1/numSamples) - K)
    arithPrices[k] <- df*max(0,mean(path[-1]) - K)
  }
  covar <- cov(geomPrices, arithPrices)
  cstar <- -covar/var(geomPrices)
  # MC run
  cVars = numeric(numRepl)
  for (k in 1:numRepl){
    path <- simvGBM(S0,r-q,sigma,T,numSamples,1)
    geomPrice <- df*max(0,(prod(path[-1]))^(1/numSamples) - K)
    arithPrice <- df*max(0,mean(path[-1]) - K)
    cVars[k] <- arithPrice + cstar*(geomPrice-geomExact)
  }
  aux <- t.test(cVars,conf.level=c.level)
  value=as.numeric(aux$estimate)
  ci=as.numeric(aux$conf.int)
  return(list(estimate=value,conf.int=ci))
}

########################################################
# Estimating the vanilla call delta by a pathwise estimator
VanillaDeltaMCPath <- function(S0,K,T,r,sigma,numRepl,
                               c.level=0.95){
  nuT <- (r - 0.5*sigma^2)*T
  sigmaT <- sigma*sqrt(T)
  VLogn <- exp(nuT+sigmaT*rnorm(numRepl))
  SampleDelta <- exp(-r*T) * VLogn * (S0*VLogn > K)
  aux <- t.test(SampleDelta,conf.level=c.level)
  value <- as.numeric(aux$estimate)
  ci <- as.numeric(aux$conf.int)
  return(list(estimate=value,conf.int=ci))
}


########################################################
# Estimating the Asian call delta by a pathwise estimator
AsianDeltaPath <- function(S0,K,r,T,sigma,numRepl,
                           numSamples,c.level=0.95){
  df <- exp(-r*T)
  dt <- T/numSamples
  nuT <- (r-sigma^2/2)*dt
  sigmaT <- sqrt(dt)*sigma
  sampleDiff <- numeric(numRepl)
  for (k in 1:numRepl){
    # store normals 
    normVet <- rnorm(numSamples,nuT,sigmaT)    
    logIncrements <- c(log(S0),normVet)
    path <- exp(cumsum(logIncrements))
    meanS <- mean(path[-1])
    sampleDiff[k] <- df*(meanS > K)*meanS/S0
  }
  aux <- t.test(sampleDiff,conf.level=c.level)
  value <- as.numeric(aux$estimate)
  ci <- as.numeric(aux$conf.int)
  return(list(estimate=value,conf.int=ci))
}

########################################################
# Estimating the vanilla call delta by likelihood ratios
VanillaDeltaMCLR <- function(S0,K,T,r,sigma,numRepl,
                             c.level=0.95){
  nuT <- (r - 0.5*sigma^2)*T
  sigmaT <- sigma*sqrt(T)
  Z <- rnorm(numRepl)
  ST <- S0*exp(nuT+sigmaT*Z)
  SampleDelta <- exp(-r*T)*pmax(0,ST-K)*Z/(S0*sigmaT)
  aux <- t.test(SampleDelta,conf.level=c.level)
  value <- as.numeric(aux$estimate)
  ci <- as.numeric(aux$conf.int)
  return(list(estimate=value,conf.int=ci))
}


########################################################
# Estimating the Asian call delta by likelihood ratios
AsianDeltaLR <- function(S0,K,r,T,sigma,numRepl,numSamples,
                         c.level=0.95){
  df <- exp(-r*T)
  dt <- T/numSamples
  nuT <- (r-sigma^2/2)*dt
  sigmaT <- sqrt(dt)*sigma
  sampleDiff <- numeric(numRepl)
  for (k in 1:numRepl){
    # store STANDARD normals 
    Veps <- rnorm(numSamples)
    logIncrements <- c(log(S0),nuT+Veps*sigmaT)
    path <- exp(cumsum(logIncrements))
    payoff <- max(0,mean(path[-1]) - K)
    sampleDiff[k] <- df*payoff*Veps[1]/(S0*sigmaT)
  }
  aux <- t.test(sampleDiff,conf.level=c.level)
  value <- as.numeric(aux$estimate)
  ci <- as.numeric(aux$conf.int)
  return(list(estimate=value,conf.int=ci))
}
