R/analytic_merton.R

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

#' Merton's jump-diffusion model analytic price formula for European options
#'
#' @param spot the spot price
#' @param strike The strike price
#' @param maturity the year until expirations
#' @param rate the drift rate
#' @param div The dividend yield rate
#' @param volat the diffusion coefficients of each component (a vector)
#' @param lambda the yearly jump-rate
#' @param jm the jump mean
#' @param jv the jump volatility
#' @param what What type of option, put or call
#' @param num_terms the number of terms to sum in the Poisson expectation
#' 
#' @description {Semi-nalytic series formula for European options under Merton's jump diffusion dynamics}
#' @return numeric
analytic_merton_price <- function(spot, strike, maturity, rate, div, volat, lambda, jm, jv, what = "call", num_terms = 100)
{
  eta <- exp(jm+jv^2/2)-1
  lam <- lambda*(1+eta)
  n  <- 0:num_terms
  sigma_n <- sqrt(volat^2+(n*jv^2)/maturity)
  rtilde <- rate-div-lambda*eta-0.5*volat^2+n*jm/maturity
  r_n <- rtilde+sigma_n^2/2
  bs_n <- mapply(function(X, Y){
    analytic_gbm_price(spot, strike, maturity, Y, 0, X, what, current_time = 0)
  }, X = sigma_n, Y = r_n)
  pp <- stats::dpois(n, lambda = lam*maturity)
  price <- round(sum(bs_n*pp),6)
  return(price)
}

#' Merton's jump-diffusion model delta for European options
#'
#' @param spot the spot price
#' @param strike The strike price
#' @param maturity the year until expirations
#' @param rate the drift rate
#' @param div The dividend yield rate
#' @param volat the diffusion coefficients of each component (a vector)
#' @param lambda the yearly jump-rate
#' @param jm the jump mean
#' @param jv the jump volatility
#' @param what What type of option, put or call
#' @param num_terms the number of terms to sum in the Poisson expectation
#' 
#' @description {Semi-analytic series formula for Delta greek under jump diffusion dynamics, i.e. the price sensitivity.}
#' @return Returns numeric
analytic_merton_delta <- function(spot, strike, maturity, rate, div, volat, lambda, jm, jv, what = "call", num_terms = 100)
{
  eta <- exp(jm+jv^2/2)-1
  lam <- lambda*(1+eta)
  n  <- 0:num_terms
  sigma_n <- sqrt(volat^2+(n*jv^2)/maturity)
  rtilde <- rate-div-lambda*eta-0.5*volat^2+n*jm/maturity
  r_n <- rtilde+sigma_n^2/2
  bs_n <- mapply(function(X, Y){
    analytic_gbm_delta(spot, strike, maturity, Y, 0, X, what, current_time = 0)
  }, X = sigma_n, Y = r_n)
  pp <- stats::dpois(n, lambda = lam*maturity)
  delta <- round(sum(bs_n*pp),6)
  return(delta)
}

#' Merton's jump-diffusion model gamma for European options
#'
#' @param spot the spot price
#' @param strike The strike price
#' @param maturity the year until expirations
#' @param rate the drift rate
#' @param div The dividend yield rate
#' @param volat the diffusion coefficients of each component (a vector)
#' @param lambda the yearly jump-rate
#' @param jm the jump mean
#' @param jv the jump volatility
#' @param num_terms the number of terms to sum in the Poisson expectation
#' 
#' @description {Semi-analytic series formula for Gamma greek under jump diffusion dynamics, i.e. the delta sensitivity.}
#' @return Returns numeric
analytic_merton_gamma <- function(spot, strike, maturity, rate, div, volat, lambda, jm, jv, num_terms = 100)
{
  eta <- exp(jm+jv^2/2)-1
  lam <- lambda*(1+eta)
  n  <- 0:num_terms
  sigma_n <- sqrt(volat^2+(n*jv^2)/maturity)
  rtilde <- rate-div-lambda*eta-0.5*volat^2+n*jm/maturity
  r_n <- rtilde+sigma_n^2/2
  bs_n <- mapply(function(X, Y){
    analytic_gbm_gamma(spot, strike, maturity, Y, 0, X, current_time = 0)
  }, X = sigma_n, Y = r_n)
  pp <- stats::dpois(n, lambda = lam*maturity)
  gamma <- round(sum(bs_n*pp),6)
  return(gamma)
}

#' Merton's jump-diffusion model theta for European options
#'
#' @param spot the spot price
#' @param strike The strike price
#' @param maturity the year until expirations
#' @param rate the drift rate
#' @param div The dividend yield rate
#' @param volat the diffusion coefficients of each component (a vector)
#' @param lambda the yearly jump-rate
#' @param jm the jump mean
#' @param jv the jump volatility
#' @param what What type of option, put or call
#' @param num_terms the number of terms to sum in the Poisson expectation
#' 
#' @description {Semi-analytic series formula for Theta greek under jump diffusion dynamics, i.e. the time until expiration sensitivity.}
#' @return Returns numeric
analytic_merton_theta <- function(spot, strike, maturity, rate, div, volat, lambda, jm, jv, what = "call", num_terms = 100)
{
  eta <- exp(jm+jv^2/2)-1
  lam <- lambda*(1+eta)
  n  <- 0:num_terms
  sigma_n <- sqrt(volat^2+(n*jv^2)/maturity)
  rtilde <- rate-div-lambda*eta-0.5*volat^2+n*jm/maturity
  r_n <- rtilde+sigma_n^2/2
  bs_n <- mapply(function(X, Y){
    analytic_gbm_price(spot, strike, maturity, Y, 0, X, what, current_time = 0)
  }, X = sigma_n, Y = r_n)
  theta_n <- mapply(function(X, Y){
    analytic_gbm_theta(spot, strike, maturity, Y, 0, X, what, current_time = 0)*360
  }, X = sigma_n, Y = r_n)
  pp <- stats::dpois(n, lambda = lam*maturity)
  drv <- pp*(lam-n/maturity)
  theta <- round(sum(drv*bs_n+pp*theta_n),6)/360
  return(theta)
}

#' Analytic option pricing and greeks under Merton's jump diffusion model
#'
#' @param spot the spot price
#' @param strike The strike price
#' @param maturity the year until expirations
#' @param rate the drift rate
#' @param div The dividend yield rate
#' @param volat the diffusion coefficients of each component (a vector)
#' @param lambda the yearly jump-rate
#' @param jm the jump mean
#' @param jv the jump volatility
#' @param what What type of option, put or call
#' @param output either "greeks" or "price"
#' @param num_terms the number of terms to sum in the Poisson expectation
#' 
#' @description {Calls all analytic formulas for option pricing under Merton's log-normal jump-diffusion}
#' @return data frame
#' @export analytic_merton
#'
#' @examples
#' spot <- 100
#' maturity <- 30/252
#' strike <- 100
#' rate <- 0.05
#' div <- 0
#' volat <- 0.19
#' lambda <- 5
#' jm <- 0.01
#' jv <- 0.03
#' what <- "call"
#' print(analytic_merton(spot, strike, maturity, rate, div, volat, lambda, jm, jv, what))
analytic_merton <- function(spot, strike, maturity, rate, div, volat, lambda, jm, jv, what = "call", output = "greeks", num_terms = 100)
{
  price <- analytic_merton_price(spot, strike, maturity, rate, div, volat, lambda, jm, jv, what, num_terms)
  if(output == "greeks")
  {
    delta <- analytic_merton_delta(spot, strike, maturity, rate, div, volat, lambda, jm, jv, what, num_terms)
    gamma <- analytic_merton_gamma(spot, strike, maturity, rate, div, volat, lambda, jm, jv, num_terms)
    theta <- analytic_merton_theta(spot, strike, maturity, rate, div, volat, lambda, jm, jv, what, num_terms)
    model_output <- data.frame(price = price, delta = delta, gamma = gamma, theta = theta)
    return(model_output)  
  } else if(output == "price")
  {
    return(price)
  }
  
}