R/sample_path2.R

In shill1729/FeynmanKacSolver: Feynman-Kac PDE solver for Ito diffusions and processes

#' Simulate Ito-Levy process by model name and parameters
#'
#' @param t terminal time
#' @param n number of nodes-1
#' @param continuous.model "gbm" or "lvm"
#' @param jump.model "kou", "norm", or "unif"
#' @param parameters parameter list depending on models, see details
#' @param measure either "P" or "Q"
#'
#' @description {TO BE DEPRECATED--------Simulate stock prices under either the real or the risk-neutral measure for a combination of jump-diffusion models, including geometric Brownian motion, local volatility Gaussian mixtures,
#' for the continuous-diffusion part, and the double exponential Kou model, the Merton normal jump-diffusion model, and uniformly distributed jumps for the
#' discontinuous jump part. See details for specifics on \code{parameters}.}
#' @details {Regardless of the models \code{parameters} should contain
#' \itemize{
#' \item \code{spot} the current stock price or interest rate level
#' \item \code{rate} the risk-neutral rate
#' \item \code{lambda} the mean rate of jumps function (see below)}
#' where \code{lambda} is a function of \code{(x,t)} describing the mean rate of arrivals of jumps,
#'
#' for \code{continuous.model == "gbm"} \code{parameters} should additionally contain
#' \itemize{
#' \item \code{mu} the annual mean return of the stock in the geometric Brownian motion model for stock prices
#' \item \code{volat} the annual volatility of the stock in the geometric Brownian motion model for stock prices}
#' for \code{continuous.model = "lvm"}, \code{parameters} should additionally contain
#' \itemize{
#' \item \code{probs} the vector of probabilities for each mixture component
#' \item \code{mus} the vector of annual mean-returns for each mixture component
#' \item \code{sigmas} the vector of annual volatilitys for each mixture component}
#' for \code{jump.model = "kou"}, \code{parameters} should additionally contain
#' \itemize{
#' \item \code{p} the probability of positive jump
#' \item \code{alpha} the inverse of the mean size of positive jumps
#' \item \code{beta} the inverse of the mean size of negative jumps}
#' for \code{jump.model = "runif"} \code{parameters} should additionally contain
#' \itemize{
#' \item \code{min} the smallest possible jump size
#' \item \code{max} the largest possible jump size}
#' and finally for \code{jump.model = "norm"} \code{parameters} should additionally contain
#' \itemize{
#' \item \code{mean} the mean jump-size
#' \item \code{sd} the volatility of the jump-size}
#' }
#' @return data.frame
#' @export sample_path2
sample_path2 <- function(t, n, continuous.model, jump.model, parameters, measure = "Q")
{

  # Missing input error handling
  if(!all(c("spot", "lambda") %in% names(parameters)))
  {
    stop("'spot', and 'lambda' must always be passed in 'parameters'")
  }

  # Missing input error handling for jump-parameters
  if(jump.model == "kou")
  {
    if(!all(c("p", "alpha", "beta") %in% names(parameters)))
    {
      stop("'parameters' must contain 'p', 'alpha', and 'beta' for 'kou' jump-model")
    }
  } else if(jump.model == "unif")
  {
    if(!all(c("min", "max") %in% names(parameters)))
    {
      stop("'parameters' must contain 'min' and 'max' for 'unif' jump-model")
    }
  } else if(jump.model == "norm")
  {
    if(!all(c("mean", "sd") %in% names(parameters)))
    {
      stop("'parameters' must contain 'mean' and 'sd' for 'norm' jump-model")
    }
  }
  # exponential for stock models, FALSE for "cir" for interest rates
  exponential <- ifelse(continuous.model == "cir", FALSE, TRUE)
  x0 <- ifelse(continuous.model == "cir", parameters$spot, 0)



  # Assign jump parameters
  switch (jump.model,
          kou = param <- list(p = parameters$p, alpha = parameters$alpha, beta = parameters$beta),
          unif = param <- list(min = parameters$min, max = parameters$max),
          norm = param <- list(mean = parameters$mean, sd = parameters$sd)
  )


  # Jump-size component definitions
  jumps <- list(distr = jump.model,
                param = param
  )

  if(continuous.model == "gbm")
  {
    if(measure == "Q")
    {
      if(!all(c("volat", "rate") %in% names(parameters)))
      {
        stop("'rate' and volat' must be passed in 'parameters' for 'gbm' model under Q")
      }
      spot <- parameters$spot
      rate <- parameters$rate
      volat <- parameters$volat
      lambda <- parameters$lambda
      # Coefficient functions
      coeff <- list(f.mu = function(x, t) rate-0.5*volat^2,
                    f.vo = function(x, t) volat,
                    f.lambda = lambda
      )
    } else if(measure == "P")
    {
      if(!all(c("volat", "mu") %in% names(parameters)))
      {
        stop("'mu' and 'volat' must be passed in 'parameters' for 'gbm' model under P")
      }
      spot <- parameters$spot
      mu <- parameters$mu
      volat <- parameters$volat
      lambda <- parameters$lambda
      # Coefficient functions
      coeff <- list(f.mu = function(x, t) mu-0.5*volat^2,
                    f.vo = function(x, t) volat,
                    f.lambda = lambda
      )
    }

  } else if(continuous.model == "lvm")
  {

    if(measure == "Q")
    {
      if(!all(c("rate","probs", "sigmas") %in% names(parameters)))
      {
        stop("'probs' and 'sigmas' must be passed in 'parameters' for 'lvm' model under Q")
      }
      spot <- parameters$spot
      rate <- parameters$rate
      probs <- parameters$probs
      sigmas <- parameters$sigmas
      lambda <- parameters$lambda
      # Coefficient functions
      coeff <- list(f.mu = function(x, t) rate-0.5*volat_lvm(spot*exp(x), t, probs, sigmas, rate, spot)^2,
                    f.vo = function(x, t) volat_lvm(spot*exp(x), t, probs, sigmas, rate, spot),
                    f.lambda = lambda
      )
    } else if(measure == "P")
    {
      if(!all(c("probs","mus", "sigmas") %in% names(parameters)))
      {
        stop("'probs', 'mus', and 'sigmas' must be passed in 'parameters' for 'lvm' model under P")
      }
      spot <- parameters$spot
      mus <- parameters$mus
      probs <- parameters$probs
      sigmas <- parameters$sigmas
      lambda <- parameters$lambda
      # Coefficient functions
      coeff <- list(f.mu = function(x, t) drift_lvm(spot*exp(x), t, probs, mus, sigmas, spot)-0.5*volat_lvm(spot*exp(x), t, probs, sigmas, mus, spot)^2,
                    f.vo = function(x, t) volat_lvm(spot*exp(x), t, probs, sigmas, mus, spot),
                    f.lambda = lambda
      )
    }
  } else if(continuous.model == "cir")
  {

    if(!all(c("theta", "mu", "volat") %in% names(parameters)))
    {
      stop("'theta', 'mu', and 'volat' must be passed in 'parameters' for 'cir'")
    }
    spot <- parameters$spot
    theta <- parameters$theta
    mu <- parameters$mu
    volat <- parameters$volat
    lambda <- parameters$lambda
    if(2*theta*mu < volat^2)
    {
      stop("Does not preclude zero interest rates: check parameter values so that 2*'theta'*'mu'>='volat'^2")
    }
    # Coefficient functions
    coeff <- list(f.mu = function(x, t) drift_vasicek(x, t, theta = theta, mu = mu),
                  f.vo = function(x, t) volat_cir(x, t, vcf = volat),
                  f.lambda = lambda
    )
  }
  x <- euler_maruyama(t = t, coeff = coeff, IC = list(x0 = x0, spot = spot), jumps = jumps, n = n, exponential = exponential)
  return(x)
}