
####################################################
#   Function definitions for Chapter 6


################################################
# Sampling Geometric Brownian motion
# scalar GBM (loop version)
simGBM <- function (S0,mu,sigma,T,numSteps,numRepl){
  dt <- T/numSteps
  # precompute invariant quantities
  nuT <- (mu-sigma^2/2)*dt
  sigmaT <- sqrt(dt)*sigma
  # allocate matrix to store sample paths
  pathMatrix <- matrix(nrow=numRepl,ncol=numSteps+1)
  pathMatrix[,1] <- S0
  for (i in 1:numRepl)
    for (j in 2:(numSteps+1))
      pathMatrix[i,j] <- pathMatrix[i,j-1]*
    exp(rnorm(1,nuT,sigmaT))
  return(pathMatrix)
}

# scalar GBM (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))
}


################################################
# Pricing a vanilla call
# 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))
}


################################################
# Multidimensional GBM

# bidimensional GBM: receives drifts, volatilities, 
# instantaneous correlation and number of steps
# naive implementation
sim2GBM <- function(S0,mus,sigmas,rho,T,numSteps){
  dt <- T/numSteps
  nuTs <- (mus-sigmas^2/2)*dt
  sigmaTs <- sqrt(dt)*sigmas
  pathMatrix <- matrix(nrow=2,ncol=numSteps+1)
  pathMatrix[,1] <- S0
  sqrtRho <- sqrt(1-rho^2)
  for (j in 2:(numSteps+1)){
    eps1 <- rnorm(1)
    pathMatrix[1,j] <- pathMatrix[1,j-1]*
      exp(nuTs[1]+sigmaTs[1]*eps1)
    eps2 <- rho*eps1+sqrtRho*rnorm(1)
    pathMatrix[2,j] <- pathMatrix[2,j-1]*
      exp(nuTs[2]+sigmaTs[2]*eps2)
  }
  return(pathMatrix)
}

# multidimensional GBM: input vector of drifts, vector of
# volatilities, and matrix of instantaneous correlations
simMVGBM <- function(S0Vet,muVet,sigmaVet,corrMatrix,T,
                     numSteps,numRepl){
  # precompute invariant quantities
  dt <- T/numSteps
  nuT <- (muVet-sigmaVet^2/2)*dt
  sigmaT <- sqrt(dt)*sigmaVet
  # find Cholesky factor: transposition is necessary
  # to obtain the lower triangular factor
  lowerCholFactor <- t(chol(corrMatrix))
  numFactors <- length(muVet)
  # prepare three-dimensional array
  pathArray <- array(dim=c(numRepl,numSteps+1,numFactors))
  for (k in 1:numFactors) pathArray[,1,k] <- S0Vet[k]
  for (i in 1:numRepl){
    for (j in 2:(numSteps+1)){
      # generate independent standard normals 
      # and then correlate them
      zVet <- rnorm(numFactors)
      epsVet <- lowerCholFactor %*% zVet
      # new sample vector
      pathArray[i,j,] <- pathArray[i,j-1,]*
        exp(nuT+sigmaT*epsVet)
    }
  }
  return(pathArray)
}


################################################
# Path generation by the Brownian brdge

# Standard Wiener process 
WienerBridge <- function(T, numSteps){
  numBisections <- log2(numSteps)
  if (round(numBisections) != numBisections){
    cat('ERROR: numSteps must be a power of 2\n')
    return(NaN)
  }
  path <- numeric(numSteps+1)
  path[1] <- 0
  path[numSteps+1] <- sqrt(T)*rnorm(1)
  TJump <- T
  IJump <- numSteps
  for (k in 1:numBisections) {
    left <- 1
    i <- IJump/2 + 1
    right <- IJump + 1
    for (j in 1:2^(k-1)){
      a <- 0.5*(path[left] + path[right])
      b <- 0.5*sqrt(TJump)
      path[i] <- a + b*rnorm(1)
      right <- right + IJump
      left <- left + IJump
      i <- i + IJump
    }
    IJump <- IJump/2
    TJump <- TJump/2
  }
  return(path)
}

# The case of GBM
GBMBridge <- function(S0, mu, sigma, T, numSteps, numRepl){
  numBisections <- log2(numSteps)
  if (round(numBisections) != numBisections){
    cat('ERROR: numSteps must be a power of 2\n')
    return(NaN)
  }
  dt <- T/numSteps
  nudt <- (mu-0.5*sigma^2)*dt
  Paths <- matrix(0, nrow=numRepl, ncol=numSteps+1)
  for (k in 1:numRepl){
    W <- WienerBridge(T,numSteps)
    Increments <- nudt + sigma*diff(W)
    LogPath <- cumsum(c(log(S0) , Increments))
    Paths[k,] <- exp(LogPath)
  }
  Paths[,1] <- S0
  return(Paths)
}

################################################
# Path generation for short rate models

# The Vasicek model
simVasicek<-function(r0,rbar,gamma,sigma,T,numSteps,numRepl){
  dt <- T/numSteps
  # precompute time invariant quantities
  expg <- exp(-gamma*dt)
  auxc <- rbar*(1-expg)
  vol <- sigma*sqrt((1-expg^2)/(2*gamma))
  paths <- matrix(0,nrow=numRepl,ncol=numSteps+1)
  paths[,1] <- r0
  epsMat <- matrix(rnorm(numRepl*numSteps),nrow=numRepl)
  for (i in 1:numRepl){
    for (t in 1:numSteps)
      paths[i,t+1] <- paths[i,t]*expg+auxc+vol*epsMat[i,t]
  }
  return(paths)
}

# THE CIR model: Exact approach
simCIR <- function(r0,rbar,gamma,alpha,T,
                   numSteps,numRepl){
  dt <- T/numSteps
  # precompute time invariant quantities
  expg <- alpha*(1-exp(-gamma*dt))/(4*gamma)
  nu <- 4*rbar*gamma/alpha
  lambdaC <- exp(-gamma*dt)/expg
  paths <- matrix(0,nrow=numRepl,ncol=numSteps+1)
  paths[,1] <- r0
  for (i in 1:numRepl){
    for (t in 1:numSteps){
      chi <- rchisq(1,df=nu,ncp=lambdaC*paths[i,t])
      paths[i,t+1] <- expg*chi}}
  return(paths)
}

# The CIR model: Euler approximation
simEulerCIR <- function(r0,rbar,gamma,alpha,T,
                        numSteps,numRepl){
  dt <- T/numSteps
  # precompute time invariant quantities
  aux1 <- gamma*rbar*dt
  aux2 <- 1-gamma*dt
  aux3 <- sqrt(alpha*dt)
  paths <- matrix(0,nrow=numRepl,ncol=numSteps+1)
  paths[,1] <- r0
  epsMat <- matrix(rnorm(numRepl*numSteps),nrow=numRepl)
  for (i in 1:numRepl){
    for (t in 1:numSteps){
      paths[i,t+1] <- aux1+aux2*paths[i,t]+
        sqrt(paths[i,t])*aux3*epsMat[i,t]}}
  return(paths)
}
