R/pde_solve.R

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

#' Solve Feynman-Kac PDE/variational inequality for Ito diffusions
#'
#' @param pde.setup list containing problem setup. See details
#' @param variational boolean for variational inequality vs PDE proper
#' @param output type of output to return: price, greeks, or grid
#' @description {Solves the PDE or variational inequality for an associated Ito diffusion via Feynman Kac formula.}
#' @details {An implicit scheme is implemented to solve the PDE and dynamic programming used for the variational inequality.
#' The pde.setup list should contain "grid.x", "grid.t", "alpha", "beta", "delta".
#' The last three are the coefficients of the linear system, i.e. the tridiagonal entries of the matrix.}
#' @return depending on the \code{output} input, this function either returns the price (numeric), a vector of the price and greeks, or
#' the entire solution grid.
#' @export pde_solve
pde_solve <- function(pde.setup, variational = FALSE, output = "price")
{

  # list2env(pde.setup, envir = environment())
  N <- length(pde.setup$grid.t)-1
  M <- length(pde.setup$grid.x)-1
  h <- pde.setup$grid.x[2]-pde.setup$grid.x[1]
  k <- pde.setup$grid.t[2]-pde.setup$grid.t[1]
  g <- pde.setup$payoff
  u <- matrix(0, nrow = N+1, ncol = M+1)
  # IC
  u[1, ] <- g
  # BC
  u[,1] <- g[1]
  u[,M+1] <- g[M+1]

  # Create tridiagonal system
  for(i in 2:(N+1))
  {
    # auxillary vector
    b.vec <- c(pde.setup$alpha[i, 2]*u[i, 1], rep(0, M-3), pde.setup$delta[i, M]*u[i, M+1])

    # Taking diagonals as vecotrs and using thomas algorithm
    aa <- -k*pde.setup$alpha[i, 3:M]
    bb <- 1-k*pde.setup$beta[i, 2:M]
    cc <- -k*pde.setup$delta[i, 2:(M-1)]
    if(variational){
      u[i, 2:M] <- pmax(thomas_algorithm(aa, bb, cc, k*b.vec+u[i-1, 2:M]), g[2:M])
    } else{
      u[i, 2:M] <- thomas_algorithm(aa, bb, cc, k*b.vec+u[i-1, 2:M])
    }

  }
  if(output == "grid")
  {
    return(u)
  } else if(output == "price")
  {
    return(u[N+1, M/2+1])
  } else if(output == "greeks")
  {
    i <- (M/2+1)
    spot <- pde.setup$spot
    r <- pde.setup$r
    div <- pde.setup$div
    v <- pde.setup$v
    fee <- u[N+1, i]
    delta <- (u[N+1, i+1]-u[N+1, i-1])/(2*h)
    gamma <- (u[N+1, i+1]-2*u[N+1, i]+u[N+1, i-1])/(h^2)
    gamma <- (gamma-delta)/(spot^2)
    delta <- delta/spot
    theta <- (u[N, i]-u[N+1, i])/k
    theta <- theta/360
    return(data.frame(fee = fee, delta = delta, gamma = gamma, theta = theta))
  }

}