R/check_convergence.R

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

#' Compare prices for European vanilla options produced by the Monte-Carlo routine vs PDE solver
#'
#' @param spot the underlying spot price
#' @param maturity the maturity of the option in trading years
#' @param rate the risk-neutral rate of return
#' @param div the dividend yield of the underlying
#' @param volat the annual volatility
#' @param what what type of option
#' @param N time resolution
#' @param M space resolution
#' @param threshold threshold for MSE error for acceptance of convergence
#' @param samples number of samples to use in MC
#' @param strike_range range of strikes to price
#' 
#' @description {Check for mutual convergence of prices generated by both PDE and MC methods to a common value under Black-Scholes
#' dynamics for European vanilla options.}
#' @return list containing a data.frame of the strike and genertaed prices, the mse, boolean for success and success message if successful.
#' @export check_gbm_convergence
check_gbm_convergence <- function(spot, maturity, rate, div, volat, what, N, M, threshold = 0.001, samples = 10^6, strike_range = 5)
{
  style <- "european"
  output <- "price"
  if(what == "call")
  {
    K <- seq(spot-strike_range/2, spot+strike_range, by = 0.5)  
  } else if(what == "put")
  {
    K <- seq(spot-strike_range, spot+strike_range/2, by = 0.5)  
  }
  analytic <- matrix(0, nrow = length(K))
  pde <- matrix(0, nrow = length(K))
  monte <- matrix(0, nrow = length(K))
  # TODO vectorize pde_gbm and monte_carlo...
  for(i in 1:length(K))
  {
    analytic[i] <- analytic_gbm(spot, K[i], maturity, rate, div, volat, what, 0, output = output)
    pde[i] <- pde_gbm(spot, K[i], maturity, rate, div, volat, what, style, output = output, N, M)  
    monte[i] <- monte_carlo_black_scholes(spot, K[i], maturity, rate, div, volat, what, samples)$point
  }
  graphics::par(mfrow = c(1, 1))
  graphics::plot(K, analytic, type = "l", xlab = "strike", ylab = "price", main = "Analytic vs PDE vs MC", ylim = c(min(analytic, pde, monte), max(analytic, pde, monte)))
  graphics::lines(K, pde, col = "blue", lty = "dashed")
  graphics::lines(K, monte, col = "green", lty = "dashed")
  graphics::legend(x = "topleft", legend = c("Analytic", "PDE", "Monte-Carlo"), col = c("black", "blue", "green"), lty = 1, cex = 0.6)
  mse_pde <- mean((pde-analytic)^2)
  mse_monte <- mean((monte-analytic)^2)
  if(mse_pde < threshold)
  {
    pde_success <- TRUE
    pde_msg <- "MSE between PDE and true value is within threshold"
  } else{
    pde_success <- FALSE
    warning("MSE between PDE and true value is greater than threshold")
  }
  if(mse_monte < threshold)
  {
    monte_success <- TRUE
    monte_msg <- "MSE between MC and true value is within threshold"
  } else{
    monte_success <- FALSE
    warning("MSE between MC and true value is greater than threshold")
  }
  return(list(model_prices = data.frame(strike = K, analytic = analytic, pde = pde, monte = monte), 
              what = what, 
              pde_error = mse_pde, 
              mc_error = mse_monte,
              threshold = threshold,
              pde_success = pde_success, 
              monte_success = monte_success))
  
}


#' Compare prices for European vanilla options produced by the Monte-Carlo routine vs PDE solver
#'
#' @param spot the underlying spot price
#' @param maturity the maturity of the option in trading years
#' @param parameters vector of parameters defining the Merton jump diffusion dynamics, see details
#' @param what what type of option
#' @param N time resolution
#' @param M space resolution
#' @param threshold threshold for MSE error for acceptance of convergence
#' @param samples number of samples to use in MC
#' @param strike_range range of strikes to price
#' 
#' @description {Check for mutual convergence of prices generated by both PDE and MC methods to a common value under Black-Scholes
#' dynamics for European vanilla options.}
#' @details {The argument \code{parameters} should be a vector containing the values (not necessarily named)
#' \itemize{
#' \item \code{rate} the risk-free rate,
#' \item \code{div} the dividend yield rate,
#' \item \code{volat} the annual log volatility,
#' \item \code{lambda} the mean rate of jumps per year,
#' \item \code{jm} the mean size of jumps,
#' \item \code{jv} the mean volatility of jumps}}
#' @return list containing a data.frame of the strike and genertaed prices, the mse, boolean for success and success message if successful.
#' @export check_merton_convergence
check_merton_convergence <- function(spot, maturity, parameters, what, N = 100, M = 100, threshold = 0.001, samples = 1000, strike_range = 5)
{
  if(length(parameters) != 6)
  {
    stop("argument 'parameters' must be a vector of 6 coordinates defining the jump dynamics, see documentation")
  }
  rate <- parameters[1]
  div <- parameters[2]
  volat <- parameters[3]
  lambda <- parameters[4]
  jm <- parameters[5]
  jv <- parameters[6]
  
  
  style <- "european"
  output <- "price"
  if(what == "call")
  {
    K <- seq(spot-strike_range/2, spot+strike_range, by = 0.5)  
  } else if(what == "put")
  {
    K <- seq(spot-strike_range, spot+strike_range/2, by = 0.5)  
  }
  analytic <- matrix(0, nrow = length(K))
  pide <- matrix(0, nrow = length(K))
  monte <- matrix(0, nrow = length(K))
  # TODO vectorize pide_gbm and monte_carlo...
  for(i in 1:length(K))
  {
    analytic[i] <- analytic_merton(spot, K[i], maturity, rate, div, volat, lambda, jm, jv, what, output)
    pide[i] <- pide_merton(spot, K[i], maturity, parameters, what, style, output, N, M, L = M/2+1)
    monte[i] <- monte_carlo_merton(spot, K[i], maturity, parameters, what, samples = samples)$point
  }
  graphics::par(mfrow = c(1, 1))
  graphics::plot(K, analytic, type = "l", xlab = "strike", ylab = "price", main = "Analytic vs PIDE vs MC", ylim = c(min(analytic, pide, monte), max(analytic, pide, monte)))
  graphics::lines(K, pide, col = "blue", lty = "dashed")
  graphics::lines(K, monte, col = "green", lty = "dashed")
  graphics::legend(x = "topleft", legend = c("Analytic", "PIDE", "Monte-Carlo"), col = c("black", "blue", "green"), lty = 1, cex = 0.6)
  mse_pide <- mean((pide-analytic)^2)
  mse_monte <- mean((monte-analytic)^2)
  if(mse_pide < threshold)
  {
    pide_success <- TRUE
    pide_msg <- "MSE between PIDE and true value is within threshold"
  } else{
    pide_success <- FALSE
    warning("MSE between PIDE and true value is greater than threshold")
  }
  if(mse_monte < threshold)
  {
    monte_success <- TRUE
    monte_msg <- "MSE between MC and true value is within threshold"
  } else{
    monte_success <- FALSE
    warning("MSE between MC and true value is greater than threshold")
  }
  return(list(model_prices = data.frame(strike = K, analytic = analytic, pide = pide, monte = monte), 
              what = what, 
              pide_error = mse_pide, 
              mc_error = mse_monte,
              threshold = threshold,
              pide_success = pide_success, 
              monte_success = monte_success))
  
}