R/pide_solve.R

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

#' PIDE solver
#'
#' @param pide.setup output object from pide_setup
#' @param variational boolean for optional stopping time
#' @param output output choice
#'
#' @return vector, data.frame or matrix
#' @export pide_solve
pide_solve <- function(pide.setup, variational = FALSE, output = "price")
{
  # list2env(pide.setup, envir = environment())
  N <- pide.setup$N
  M <- pide.setup$M
  L <- pide.setup$L
  h <- pide.setup$grid.x[2]-pide.setup$grid.x[1]
  k <- pide.setup$grid.t[2]-pide.setup$grid.t[1]

  Q <- M+1+2*L
  g <- pide.setup$payoff
  u <- matrix(0, nrow = N+1, ncol = Q)
  # IC
  u[1, ] <- g
  # Extended BC
  for(i in 1:(N+1))
  {
    u[i, 1:(L+1)] <- g[1:(L+1)]
    u[i, (M+L+1):(M+2*L+1)] <- g[(M+L+1):(M+2*L+1)]
  }
  # Jump-size density
  jumpden <- pide.setup$jumpden

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

    # Taking diagonals as vecotrs and using thomas algorithm
    aa <- -k*pide.setup$alpha[i, 3:M+L]
    bb <- 1-k*pide.setup$beta[i, 2:M+L]
    cc <- -k*pide.setup$delta[i, 2:(M-1)+L]

    # Populate jump-matrix
    jump_matrix <- matrix(0, M-1, 2*L+1)
    for(j in 2:(M))
    {
      for(l in (-L):L)
      {
        jump_matrix[j-1, l+L+1] <- u[i-1, j+l+L]
      }
    }

    Ju <- 0.5*h*pide.setup$lambda[i, 2:M+L]*jump_matrix%*%jumpden
    if(i == 2)
    {
      ju1 <- Ju
    }

    if(variational){
      target <- k*(b.vec+Ju)+u[i-1, 2:M+L]
      u[i, 2:M+L] <- pmax(thomas_algorithm(aa, bb, cc, target), g[2:M+L])
    } else{
      target <- k*(b.vec+Ju)+u[i-1, 2:M+L]
      u[i, 2:M+L] <- thomas_algorithm(aa, bb, cc, target)
    }

  }
  if(output == "grid")
  {
    return(u)
  } else if(output == "price")
  {
    return(u[N+1, Q/2+1])
  } else if(output == "greeks")
  {
    i <- (Q/2+1)
    spot <- pide.setup$spot
    r <- pide.setup$r
    div <- pide.setup$div
    v <- pide.setup$v
    la <- pide.setup$la
    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))
  }

}