R/monte_carlo.R

In shill1729/OptionPricer: An option pricer for European and American vanilla stock options

#' Comptue a conditional expectation via Monte-Carlo simulation/averaging
#' 
#' @param spot the spot price
#' @param strike the strike sprice
#' @param maturity the maturity of the option contract
#' @param rate the risk-neutral rate
#' @param div the dividend yield rate
#' @param volat the volatility
#' @param what the payoff to use
#' @param samples number of variates to sample
#' 
#' @description {Compute a conditional expectation via discounting and averaging}
#' @export monte_carlo_black_scholes
monte_carlo_black_scholes <- function(spot, strike, maturity, rate, div, volat, what = "call", samples = 20000)
{
  variates <- spot*exp((rate-div-0.5*volat^2)*maturity+volat*sqrt(maturity)*stats::rnorm(samples))
  if(what == "call")
  {
    variates <- pmax(variates-strike, 0.0)
  } else if(what == "put")
  {
    variates <- pmax(strike - variates, 0.0)
  } else if(what == "straddle")
  {
    variates <- pmax(strike-variates, 0.0)+pmax(variates-strike, 0.0)
  } else if(what == "cdf")
  {
    variates <- ifelse(variates <= strike, 1, 0)
  }
  if(what != "cdf")
  {
    variates <- exp(-(rate-div)*maturity)*variates  
  }
  estimate <- data.frame(point = mean(variates), stderr = stats::qnorm(1-0.05/2)*stats::sd(variates)/sqrt(samples))
  return(estimate)
}


#' Comptue a conditional expectation via Monte-Carlo simulation/averaging
#' 
#' @param spot the spot price
#' @param strike the strike sprice
#' @param maturity the maturity of the option contract
#' @param parameters vector of model parameters
#' @param what the payoff to use
#' @param samples number of variates to sample
#' 
#' @description {Compute a conditional expectation via discounting and averaging}
#' @export monte_carlo_merton
monte_carlo_merton <- function(spot, strike, maturity, parameters, what = "call", samples = 20000)
{
  rate <- parameters[1]
  div <- parameters[2]
  volat <- parameters[3]
  lambda <- parameters[4]
  jm <- parameters[5]
  jv <- parameters[6]
  numjumps <- stats::rpois(1, lambda = lambda*maturity)
  jumpsizes <- stats::rnorm(numjumps, mean = jm, sd = jv)
  jumps <- sum(jumpsizes)
  eta <- exp(jm+0.5*jv^2)-1
  variates <- spot*exp((rate-div-lambda*eta-0.5*volat^2)*maturity+volat*sqrt(maturity)*stats::rnorm(samples) + jumps)
  if(what == "call")
  {
    variates <- pmax(variates-strike, 0.0)
  } else if(what == "put")
  {
    variates <- pmax(strike - variates, 0.0)
  } else if(what == "straddle")
  {
    variates <- pmax(strike-variates, 0.0)+pmax(variates-strike, 0.0)
  } else if(what == "cdf")
  {
    variates <- ifelse(variates <= strike, 1, 0)
  }
  if(what != "cdf")
  {
    variates <- exp(-(rate-div)*maturity)*variates  
  }
  estimate <- data.frame(point = mean(variates), stderr = stats::qnorm(1-0.05/2)*stats::sd(variates)/sqrt(samples))
  return(estimate)
}



#' Comptue a conditional expectation via Monte-Carlo simulation/averaging
#' 
#' @param spot the spot price
#' @param strike the strike sprice
#' @param maturity the maturity of the option contract
#' @param parameters vector of model parameters
#' @param what the payoff to use
#' @param samples number of variates to sample
#' 
#' @description {Compute a conditional expectation via discounting and averaging}
#' @export monte_carlo_kou
#' @importFrom sde rdkou
monte_carlo_kou <- function(spot, strike, maturity, parameters, what = "call", samples = 20000)
{
  rate <- parameters[1]
  div <- parameters[2]
  volat <- parameters[3]
  lambda <- parameters[4]
  prob <- parameters[5]
  alpha <- parameters[6]
  beta <- parameters[7]
  ku <- parameters[8]
  kd <- parameters[9]
  
  numjumps <- stats::rpois(1, lambda = lambda*maturity)
  jumpsizes <- sde::rdkou(numjumps, prob, alpha, beta, ku, kd)
  jumps <- sum(jumpsizes)
  eta <- mgf_dkou(1, prob, alpha, beta, ku, kd)-1
  variates <- spot*exp((rate-div-lambda*eta-0.5*volat^2)*maturity+volat*sqrt(maturity)*stats::rnorm(samples) + jumps)
  if(what == "call")
  {
    variates <- pmax(variates-strike, 0.0)
  } else if(what == "put")
  {
    variates <- pmax(strike - variates, 0.0)
  } else if(what == "straddle")
  {
    variates <- pmax(strike-variates, 0.0)+pmax(variates-strike, 0.0)
  } else if(what == "cdf")
  {
    variates <- ifelse(variates <= strike, 1, 0)
  }
  if(what != "cdf")
  {
    variates <- exp(-(rate-div)*maturity)*variates  
  }
  estimate <- data.frame(point = mean(variates), stderr = stats::qnorm(1-0.05/2)*stats::sd(variates)/sqrt(samples))
  return(estimate)
}


#' Monte-Carlo European option pricing under exponential compensated Poisson process for prices
#'
#' @param spot current underlying share price
#' @param strike the agreed upon strike price of the option
#' @param maturity the years until expiration (in trading years)
#' @param rate the discounting rate (i.e. the risk-neutral rate)
#' @param div the dividend yield rate
#' @param a the jump coefficient
#' @param b the compensator coefficient
#' @param what the type of option to price, call or put
#' @param t the current time
#' @param n the number of simulations
#' @param withCI whether to return the confidence intervals
#' @param alpha the significance level
#' 
#' @description {Monte-Carlo simulated European option prices under exponential Levy processes. Specifically log returns are compensated Poisson processes.}
#' @return vector or data.frame
#' @export monte_carlo_pois
monte_carlo_pois <- function(spot, strike, maturity, rate, div, a, b, what = "put", t = 0, n = 1000, withCI = TRUE, alpha = 0.05)
{
  duration <- (maturity-t)
  lambda <- (rate-div+b)/(exp(a)-1)
  if(what == "put")
  {
    h <- exp(-(rate-div)*duration)*pmax(strike-spot*exp(a*stats::rpois(n, lambda = lambda*duration)-b*duration),0)
  } else if(what == "call")
  {
    h <- exp(-(rate-div)*duration)*pmax(spot*exp(a*stats::rpois(n, lambda = lambda*duration)-b*duration)-strike,0)
  }
  if(!withCI)
  {
    return(mean(h))
  } else if(withCI)
  {
    value <- mean(h)
    std_err <- stats::qnorm(1-alpha/2)*stats::sd(h)/sqrt(n)
    LB <- value-std_err
    UB <- value+std_err
    return(data.frame(LB, value, UB, std_err, n))
  }
}