
####################################################
#    Function definitions for Chapter 13
####################################################

####################################################
# AUXILIARY
####################################################

### Black-Scholes-Merton formulas for vanilla call
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)  
}

### Black-Scholes-Merton formula for vanilla put
EuVanillaPut <- function(S0,K,T,r,sigma){
  d1 <- (log(S0/K)+(r+sigma^2/2)*T)/(sqrt(T)*sigma)
  d2 <- d1-sqrt(T)*sigma
  value <- K*exp(-r*T)*pnorm(-d2) - S0*pnorm(-d1)
  return(value)  
}

### Black-Scholes-Merton formulas for vanilla call/put greeks
BSMdelta <- function(S0,K,T,r,sigma){
  d1 <- (log(S0/K)+(r+sigma^2/2)*T)/(sqrt(T)*sigma)
  return(pnorm(d1))  
}

BSMgamma <- function(S0,K,T,r,sigma){
  d1 <- (log(S0/K)+(r+sigma^2/2)*T)/(sqrt(T)*sigma)
  return(dnorm(d1)/(S0*sigma*sqrt(T)))  
}

BSMvega <- function(S0,K,T,r,sigma){
  d1 <- (log(S0/K)+(r+sigma^2/2)*T)/(sqrt(T)*sigma)
  return(dnorm(d1)*S0*sqrt(T))
}

BSMthetacall <- 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*dnorm(d1)*sigma/(2*sqrt(T))-
    r*K*exp(-r*T)*pnorm(d2)
  return(value)  
}

BSMthetaput <- 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*dnorm(d1)*sigma/(2*sqrt(T))+
    r*K*exp(-r*T)*pnorm(-d2)
  return(value)  
}

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


####################################################
# Estimating probability of shortfall for a portfolio
# of vanilla calls and puts
####################################################

### Crude Monte Carlo and full repricing
ProbShortfallCrudeMC<-function(S0,K,T,mu,r,sigma,numStocks,
            numCalls,numPuts,dt,numRepl,thresh,c.level=0.95){
  # evaluate initial portfolio value
  Call0 <- EuVanillaCall(S0,K,T,r,sigma)
  Put0 <- EuVanillaPut(S0,K,T,r,sigma)
  Value0 <- numStocks*(numCalls*Call0+numPuts*Put0)
  # check if loss is larger than threshold
  numAbove <- 0
  loss <- numeric(numRepl)
  for (k in 1:numRepl){
    # scenarios are generated according to normal
    # distribution for comparison purposes
    newS <- S0*(1+rnorm(numStocks,mu*dt,sigma*sqrt(dt)))
    newCall <- EuVanillaCall(newS,K,T-dt,r,sigma)
    newPut <- EuVanillaPut(newS,K,T-dt,r,sigma)
    newValue <- numCalls*sum(newCall)+numPuts*sum(newPut)
    loss[k] <- Value0-newValue
    if (loss[k]>thresh) numAbove <- numAbove+1
  }
  aux <- binom.test(numAbove,numRepl,conf.level=c.level)
  value <- as.numeric(aux$estimate)
  ci <- as.numeric(aux$conf.int)
  return(list(estimate=value,conf.int=ci,loss=loss))
}


### Crude Monte Carlo and delta-gamma approximation
ProbShortfallDG<-function(S0,K,T,mu,r,sigma,numStocks,
          numCalls,numPuts,dt,numRepl,thresh,c.level=0.95){
  # evaluate greeks
  DeltaCall <- BSMdelta(S0,K,T,r,sigma)
  DeltaPut <- DeltaCall-1
  Gamma <- BSMgamma(S0,K,T,r,sigma)
  ThetaCall <- BSMthetacall(S0,K,T,r,sigma)
  ThetaPut <- BSMthetaput(S0,K,T,r,sigma)
  DeltaVet <- 
    rep(numCalls*DeltaCall+numPuts*DeltaPut,numStocks)
  GammaMat <- diag((numCalls+numPuts)*Gamma,numStocks)
  ThetaTot<-numStocks*(numCalls*ThetaCall+numPuts*ThetaPut)
  # check if loss is larger than threshold
  numAbove <- 0
  loss <- numeric(numRepl)
  for (k in 1:numRepl){
    # scenarios are generated according to normal
    # distribution for comparison purposes
    dS <- S0*rnorm(numStocks,mu*dt,sigma*sqrt(dt))
    loss[k] <- -(ThetaTot*dt+crossprod(dS,DeltaVet)+
                   0.5*dS%*%GammaMat%*%dS)
    if ( loss[k] > thresh) numAbove<-numAbove+1
  }
  aux <- binom.test(numAbove,numRepl,conf.level=c.level)
  value <- as.numeric(aux$estimate)
  ci <- as.numeric(aux$conf.int)
  return(list(estimate=value,conf.int=ci,loss=loss))
}


####################################################
# Estimating probability of shortfall for a portfolio
# of vanilla calls and puts, using IMPORTANCE SAMPLING
####################################################

### Some auxiliary functions related to cumulant
### generation functions and optimization of tilting 

psi <- function(tilt,a,b,lambda){
  # auxiliary function to compute the psi function
  # see Eq. (9.5) in Glasserman
  aux <- 1-2*tilt*lambda
  out <- a*tilt + 0.5*sum((tilt*b)^2/aux-log(aux))
  return(out)
}

psiPrime <- function(tilt,a,b,lambda){
  # auxiliary function to compute the derivative of 
  # the psi function
  aux <- 1-2*tilt*lambda
  out <- a + sum(tilt*b^2*(1-tilt*lambda)/(aux^2) + 
                   lambda/aux)
  return(out)
}

findTilt <- function(thresh,a,b,lambda){
  if (any(lambda > 0))
    tiltMax <- 1/(2*max(lambda))
  else
    tiltMax <- 1000 # some BIG number
  fun <- function(x) psiPrime(x,a,b,lambda)-thresh
  tiltStar <- uniroot(fun,c(0,tiltMax))$root
  return(tiltStar)
}

### The main function implementing IS by exponential tilting
ProbShortfallIS <- function(S0,K,T,mu,r,sigma,numStocks,
                            numCalls,numPuts,dt,numRepl,thresh){
  # evaluate initial portfolio value
  Call0 <- EuVanillaCall(S0,K,T,r,sigma)
  Put0 <- EuVanillaPut(S0,K,T,r,sigma)
  Value0 <- numStocks*(numCalls*Call0+numPuts*Put0)
  # evaluate greeks
  DeltaCall <- BSMdelta(S0,K,T,r,sigma)
  DeltaPut <- DeltaCall-1
  Gamma <- BSMgamma(S0,K,T,r,sigma)
  ThetaCall <- BSMthetacall(S0,K,T,r,sigma)
  ThetaPut <- BSMthetaput(S0,K,T,r,sigma)
  DeltaVet<-rep(numCalls*DeltaCall+numPuts*DeltaPut,numStocks)
  GammaMat<-diag((numCalls+numPuts)*Gamma,numStocks)
  ThetaTot<-numStocks*(numCalls*ThetaCall+numPuts*ThetaPut)
  # covariance matrix is trivial here
  varDeltaS <- S0^2*sigma^2*dt
  sigmaMat <- diag(varDeltaS,numStocks)
  # find a, b, and lambdas in quadratic approximation
  G <- t(chol(sigmaMat)) # upper triangle
  hessian <- -0.5 * G %*% GammaMat %*% t(G)
  aux <- eigen(hessian);   lambda <- aux$values
  U <- aux$vectors; C <- G %*% U
  a <- -dt*ThetaTot; b <- -C %*% DeltaVet
  # find optimal tilting parameter
  tilt <- findTilt(thresh,a,b,lambda)
  # find parameters of tilted normals
  muTheta <- tilt*b/(1-2*lambda*tilt)
  sigmaTheta <- 1/sqrt((1-2*lambda*tilt))
  psiVal <- psi(tilt,a,b,lambda)
  # Now sample scenarios
  numAbove <- 0
  estimate <- numeric(numRepl)
  for (k in 1:numRepl){
    Z <- rnorm(numStocks,muTheta,sigmaTheta)
    Q <- a+sum(b*Z+lambda*Z^2)
    likRatio <- exp(-tilt*Q+psiVal)
    newS <- S0 + muTheta + C %*% Z
    newCall <- EuVanillaCall(newS,K,T-dt,r,sigma)
    newPut <- EuVanillaPut(newS,K,T-dt,r,sigma)
    newValue <- numCalls*sum(newCall)+numPuts*sum(newPut)
    loss <- Value0-newValue
    estimate[k] <- likRatio * (loss > thresh)
  }
  return(mean(estimate))
}


####################################################
# Estimating the cost of hedging strategies for 
# a vanilla call option
####################################################

### Simulating a stop-loss strategy (scalar version)
stopLoss <- function(paths,K,r,T){
  # paths is a matrix
  numRepl <- dim(paths)[1]
  numSteps <- dim(paths)[2]
  numSteps <-  numSteps - 1 # true number of steps
  cost <- numeric(numRepl)
  dt <- T/numSteps
  discountFactors <- exp(-r*dt*(0:numSteps))
  for (k in 1:numRepl){
    # clear cash flow vector
    cashFlows <- numeric(numSteps+1)
    if (paths[k,1] >= K){ # should I cover at beginning?
      covered <- TRUE
      cashFlows[1] <- -paths[k,1]
    }
    else
      covered <- FALSE
    # now go on until end of horizon
    for (t in 2:(numSteps+1)){
      if (covered  && (paths[k,t] < K)){
        # Sell
        covered <- FALSE
        cashFlows[t] <- paths[k,t]
      }
      else if (!covered && (paths[k,t] > K)){
        # Buy
        covered <- TRUE
        cashFlows[t] <- -paths[k,t]
      }
    }
    # what happens at maturity?
    if (paths[k,numSteps + 1] >= K) # Option is exercised
      cashFlows[numSteps + 1] <- cashFlows[numSteps + 1] + K
    # at the end of each replication get discounted cost
    cost[k] <- -crossprod(discountFactors,cashFlows)
  } # end of main loop 
  return(list(mean=mean(cost),sd=sd(cost)))
}

### Simulating a stop-loss strategy (vectorized version)
stopLossV <- function(paths,K,r,T){
  # paths is a matrix
  numRepl <- dim(paths)[1]
  numSteps <- dim(paths)[2]
  numSteps <-  numSteps - 1 # true number of steps
  cost <- numeric(numRepl)
  dt <- T/numSteps
  discountFactors <- exp(-r*dt*(0:numSteps))
  cashFlows <- matrix(0, nrow=numRepl, ncol=numSteps+1);
  oldPrice <- c(numeric(numRepl), paths[,1:numSteps])
  upTimes <- which((oldPrice < K) & (paths >= K))
  downTimes <- which((oldPrice >= K) & (paths < K))
  cashFlows[upTimes] <- -paths[upTimes]
  cashFlows[downTimes] <- paths[downTimes]
  exPaths <- which(paths[,numSteps+1] >= K)
  cashFlows[exPaths,numSteps+1] <- 
    cashFlows[exPaths,numSteps+1] + K
  cost = -cashFlows %*% discountFactors
  return(list(mean=mean(cost),sd=sd(as.numeric(cost))))
}

### Simulating a delta-hedging strategy
deltaHedging <- function(paths,K,sigma,r,T){
  # paths is a matrix
  numRepl <- dim(paths)[1]
  numSteps <- dim(paths)[2]
  numSteps <-  numSteps - 1 # true number of steps
  cost <- numeric(numRepl)
  dt <- T/numSteps
  cashFlows <- numeric(numSteps+1)
  discountFactors <- exp(-r*dt*(0:numSteps))
  for (k in 1:numRepl){
    path <- paths[k,]
    position <- 0 # initial stock holding is zero
    # get vector of deltas at each point on the sample path
    deltas <- BSMdelta(path[1:numSteps],K,
                       T-(0:(numSteps-1))*dt,r,sigma);
    for (j in 1:numSteps){
      # delta hedge along the sample path
      cashFlows[j] <- (position - deltas[j])*path[j]
      position <- deltas[j]
    }
    # check if option is exercised at maturity 
    # and liquidate position
    if (path[numSteps+1] > K)
      cashFlows[numSteps+1] <- 
      K-(1-position)*path[numSteps+1]
    else
      cashFlows[numSteps+1] <- position*path[numSteps+1]
    # at the end of each replication get discounted cost
    cost[k] <- -crossprod(discountFactors,cashFlows)
  } # end of main loop 
  return(list(mean=mean(cost),sd=sd(cost)))
}

