
####################################################
#    Function definitions for Chapter 9
####################################################


####################################################
# 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))
}

####################################################
# Finding direction numbers to define Sobol sequences
####################################################
GetDirNumbers <- function(p,m0,n){
  degree <- length(p)-1
  p <- p[2:degree]
  m <- c(m0, numeric(n-degree))
  for (i in ((degree+1):n)){
    m[i] <- bitXor(m[i-degree], 2^degree * m[i-degree])
    for (j in 1:(degree-1))
      m[i] <- bitXor(m[i], 2^j * p[j]*m[i-j])
  }
  v <- m/(2^(1:length(m)))
  return(list(v=v,m=m))
}

####################################################
# Pricing a vanilla call with low-discrepancy seqs
####################################################
BlsLowD <- function(S0,K,T,r,sigma,LDsequence){
  nuT <- (r - 0.5*sigma^2)*T
  sigmaT <- sigma*sqrt(T)
  vetEps <- qnorm(LDsequence)
  ST <- S0*exp(nuT+sigmaT*vetEps)
  dpayoff <- exp(-r*T)*pmax(0,ST-K)
  return(mean(dpayoff))
}

####################################################
# Sampling GBMs with low-discrepancy seqs
####################################################
lowDiscrGBM <- function (S0,mu,sigma,T,numSteps,numRepl,
                         type="sobol",scrambling=0){
  dt <- T/numSteps
  nuT <- (mu-sigma^2/2)*dt
  sigmaT <- sqrt(dt)*sigma
  if (toupper(type) == "HALTON")
    uSeq <- halton(numRepl,dim=numSteps)
  else
    uSeq <- sobol(numRepl,dim=numSteps,scrambling=scrambling)
  normMatrix <- nuT + sigmaT*qnorm(uSeq)
  logIncrements <- cbind(log(S0),normMatrix)
  logPath <- t(apply(logIncrements,1,cumsum))
  return(exp(logPath))
}

