R/pide_setup.R

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

#' Set up grids for PIDE
#'
#' @param f list of coefficient functions and mean-rate of jumps function
#' @param contract_specs list of variables defining contract specifications, see details
#' @param jump_specs list of variables defining the jump dynamics, see details
#' @param numeric_specs list of variables defining numerical specifications, see details
#'
#' @description {Implicit-Explicit PIDE solver for the extended Feynman-Kac connection for Ito-Levy processes with compound Poisson jump dynamics.
#' Pass a list of coefficient functions and the mean-rate of jumps, a list of contract specifications, a list of jump-dynamics specifications,
#' and finally a list of numerical specifications to solve a given PIDE problem.}
#' @details {The list \code{contract_specs} must contain
#' \itemize{
#' \item \code{spot} numeric
#' \item \code{div} numeric
#' \item \code{maturity} numeric
#' \item \code{payoff} the payoff name
#' \item \code{payoff_param} list of parameters for the payoff function}, the list \code{jump_specs} should contain
#' \itemize{
#' \item \code{distr} the jump-size distribution name: "kou", "norm" or "unif",
#' \item \code{param} the list of parameters for the above distribution}
#' \code{numeric_specs} should contain
#' \itemize{
#' \item \code{N} time-resolution (integer)
#' \item \code{M} space-resolution (integer)
#' \item \code{L} jump-resolution (integer)
#' \item \code{B} boundary}}
#' @return list
#' @export pide_setup
pide_setup <- function(f, contract_specs, jump_specs, numeric_specs)
{
  if(!contract_specs$payoff %in% c("put", "call", "put_debit", "call_debit", "iron_condor", "straddle", "strangle", "indicator", "cdf"))
  {
    stop("'payoff' in 'contract_specs' must be one of: put, call, put_debit, call_debit, iron_condor, straddle, strangle, indicator, cdf")
  }

  # Write function for calling ddistr(x,...)
  distr <- jump_specs$distr
  ddistr_name <- paste("d", distr, sep = "")
  ddistr <- function(y) do.call(what = ddistr_name, args = c(y, jump_specs$param))
  if(distr == "kou")
  {
    if(!all(c("p", "alpha", "beta") %in% names(jump_specs$param)))
    {
      stop("list 'param' in list 'jump_specs' must contain 'p', 'alpha', 'beta' for 'kou' model")
    }
  } else if(distr == "norm")
  {
    if(!all(c("mean", "sd") %in% names(jump_specs$param)))
    {
      stop("list 'param' in list 'jump_specs' must contain 'mean' and 'sd' for 'norm' model")
    }
  } else if(distr == "unif")
  {
    if(!all(c("min", "max") %in% names(jump_specs$param)))
    {
      stop("list 'param' in list 'jump_specs' must contain 'min' and 'max' for 'unif' model")
    }
  } else if(distr == "dkou")
  {
    if(!all(c("p", "alpha", "beta", "ku", "kd") %in% names(jump_specs$param)))
    {
      stop("list 'param' in list 'jump_specs' must contain 'p', 'alpha', 'beta', 'ku', 'kd', for 'dkou' model")
    }
  }


  spot <- contract_specs$spot
  div <- contract_specs$div
  maturity <- contract_specs$maturity
  N <- numeric_specs$N
  M <- numeric_specs$M
  L <- numeric_specs$L
  B <- numeric_specs$B
  if(is.null(B))
  {
    B <- 0.5*M * f$f.vo(spot,0)*sqrt(3 * maturity / N)
  }


  h <- (2*B)/M
  k <- maturity/N

  Bext <- B+L*h
  Q <- M+1+2*L
  grid.x <- seq(-Bext, Bext, length.out = Q)
  grid.t <- seq(0, maturity, length.out = N+1)

  # Grids, mu, volat.  1:N+1, 1:M+1; on paper 0:N, 0:M
  grids <- lapply(f, function(x) grid_function(x, grid.x, grid.t))
  names(grids) <- c("m", "v", "r", "la")

  eta <- mean_jump_size(jump_specs)
  g <- grids$m-div-grids$la*eta
  alpha <- grids$v^2/(2*h^2)-g/(2*h)
  beta <- -(grids$r-div)-grids$la-grids$v^2/(h^2)
  delta <- grids$v^2/(2*h^2)+g/(2*h)

  payoff_name_call <- paste("payoff_", contract_specs$payoff, sep = "")
  payoff <- function(y) do.call(what = payoff_name_call, args = append(contract_specs$payoff_param, spot*exp(y), after = 0))
  payoff <- grid_function(payoff, x = grid.x)

  jumpden <- matrix(0, 2*L+1)
  jumpden <- unlist(lapply(((-L):L)*h, ddistr))
  jumpden[-c(1, 2*L+1)] <- 2*jumpden[-c(1, 2*L+1)]
  return(list(grid.x = grid.x,
              grid.t = grid.t,
              alpha = alpha,
              beta = beta,
              delta = delta,
              spot = spot,
              div = div,
              eta = eta,
              m = grids$m[N+1, Q/2+1],
              v = grids$v[N+1, Q/2+1],
              r = grids$r[N+1, Q/2+1],
              la = grids$la[N+1, Q/2+1],
              lambda = grids$la,
              N = N,
              M = M,
              L = L,
              Bext = Bext,
              B = B,
              jumpden = jumpden,
              payoff = payoff,
              type = "pide"
              )
         )
}