R/pde_setup.R

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

#' Setting up PDE problem
#'
#' @param f the list of coefficient functions (f.mu, f.vo, f.rate)
#' @param contract_specs list of variables defining contract specifications, see details
#' @param numeric_specs list of variables defining numerical specifications, see details
#'
#' @description {Sets up grid and function evaluations over grid to create tridiagonal system coefficients for the solver to take.}
#' @details {A basic discretization scheme is used, functions are evaluated via loops then coefficients are returned in a list. For risk neutral pricing
#' the drift and rate functions should match.
#' The list \code{contract_specs} must contain
#' \itemize{
#' \item \code{spot} numeric, the spot price
#' \item \code{div} numeric, the dividend yield rate
#' \item \code{maturity} numeric, the time until expiration in trading years (DTE/252)
#' \item \code{payoff} the payoff name: "put", "call", "strangle", "straddle", "iron_condor", "put_debit", "call_debit", "indicator" and "cdf"
#' \item \code{payoff_param} list of parameters for the payoff function},
#' while the list \code{numeric_specs} should contain
#' \itemize{
#' \item \code{N} time-resolution (integer)
#' \item \code{M} space-resolution (integer)
#' \item \code{B} boundary for the log-price space, which can be optional}}
#' @return list containing
#' \itemize{
#' \item grid.x, a discretized interval of the spatial/price variable
#' \item grid.t, a discretized interval of the time variable
#' \item alpha, the lower diagonal in the tridiagonal matrix in the linear system resulting from the discretized PDE
#' \item beta, the diagonal in the tridiagonal matrix in the linear system resulting from the discretized PDE
#' \item delta, the upper diagonal in the tridiagonal matrix in the linear system resulting from the discretized PDE
#' \item spot, the given spot price, needed to be passed to \code{pde_solve}
#' \item div, the dividend rate, also needed to be passed to \code{pde_solve}
#' \item m, the drift rate at the current time and price level
#' \item v, the volatility at the current time and price level
#' \item r, the risk-free rate, at the current time and price level
#' }
#' @export pde_setup
pde_setup <- function(f, contract_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")
  }

  spot <- contract_specs$spot
  div <- contract_specs$div
  maturity <- contract_specs$maturity
  N <- numeric_specs$N
  M <- numeric_specs$M
  B <- numeric_specs$B
  if(is.null(B))
  {
    B <- 0.5*M * f$f.vo(spot,0)*sqrt(3 * maturity / N)
  }
  grid.x <- seq(-B, B, length.out = M+1)
  grid.t <- seq(0, maturity, length.out = N+1)
  h <- (2*B)/M
  k <- maturity/N
  # 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")
  m <- grids$m
  v <- grids$v
  r <- grids$r
  alpha <- v^2/(2*h^2)-m/(2*h)
  beta <- -(r-div)-v^2/(h^2)
  delta <- v^2/(2*h^2)+m/(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)

  return(list(grid.x = grid.x,
              grid.t = grid.t,
              alpha = alpha,
              beta = beta,
              delta = delta,
              spot = spot,
              div = div,
              m = m[N+1, M/2+1],
              v = v[N+1, M/2+1],
              r = r[N+1, M/2+1],
              payoff = payoff,
              type = "pde"
              )
         )
}