R/euler_maruyama.R

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

#' Simulate a sample path of a Ito-Levy process with the Euler-Maruyama scheme
#'
#' @param t the terminal time to simulate until
#' @param coeff a list of named coefficient functions for the diffusion part and named function for mean-rate of jumps, see details
#' @param IC a list of initial conditions
#' @param jumps a list of objects defining the jump-size distribution, can be NULL for purely continuous models
#' @param exponential boolean for exponential processes
#' @param n number of sub-intervals in time-grid
#' @param M number of samples to use in Monte-Carlo estimate of \code{eta} for the jump-diffusion drift correction.
#'
#' @details {\code{coefficients} must contain
#' \itemize{
#' \item \code{f.mu} the drift-coefficient function
#' \item \code{f.vo} the volatility-coefficient function
#' \item \code{f.lambda} the mean-rate function of the number of jumps
#' } while \code{jumps} should contain
#' \itemize{
#' \item \code{distr} the name of the distribution of the jump-sizes: "norm", "unif", "kou"
#' \item \code{param} named list of parameters for the distribution matching the input in \code{rdistr} for a given "distr".
#' }
#' finally, \code{IC} should contain initial conditions \code{x0} and \code{spot}}
#' @return data.frame
#' @export euler_maruyama
 euler_maruyama <- function(t, coeff, IC, jumps = NULL, exponential = TRUE, n = 1000, M = 20000)
 {
   # Extract coefficient-functions
   mu <- coeff$f.mu
   volat <- coeff$f.vo
   # Extract mean number of jumps per year
   lambda <- coeff$f.lambda
   if(!is.null(jumps))
   {
      # Extract jump-size distribution
      distr <- jumps$distr
      # Write function for calling rdistr(n)
      rdistr_name <- paste("r", distr, sep = "")
      rdistr <- function(y) do.call(what = rdistr_name, args = c(y, jumps$param))
      # Check if eta is passed and compute it if it is not
      if(is.null(jumps$eta))
      {

         # Error checking parameters for known distributions
         if(distr == "kou")
         {
            if(jumps$param$alpha >= 1)
            {
               stop("mean positive jump must be < 1")
            }
            p <- jumps$param$p
            q <- 1-p
            alpha <- 1/jumps$param$alpha
            beta <- 1/jumps$param$beta
            eta <- q*beta/(beta+1)+p*alpha/(alpha-1)-1
         }
         if(distr == "unif")
         {
            if(jumps$param$min > jumps$param$max)
            {
               stop("min must be less than max for uniform distribution")
            }
            eta <- (exp(jumps$param$max)-exp(jumps$param$min))/(jumps$param$max-jumps$param$min)-1
         }
         if(distr == "norm")
         {
            if(jumps$param$sd <= 0)
            {
               stop("sd of jump-size must be positive")
            }
            eta <- exp(jumps$param$mean+0.5*jumps$param$sd^2)-1
         }
         if(distr == "dkou")
         {
            if(jumps$param$alpha >= 1 || jumps$param$beta >= 1)
            {
               stop("mean jump sizes must be < 1")
            }
            if(jumps$param$p > 1 || jumps$param$p < 0)
            {
               stop("Probability of positive jumps must be in unit interval")
            }
            if(jumps$param$ku < 0 || jumps$param$kd > 0)
            {
               stop("Upward displacement must be positive and downward displacement must be negative")
            }
            p <- jumps$param$p
            q <- 1-p
            alpha <- 1/jumps$param$alpha
            beta <- 1/jumps$param$beta
            eta <- q*exp(jumps$param$kd)*beta/(beta+1)+p*exp(jumps$param$ku)*alpha/(alpha-1)-1
         }
         # If a known distr is not passed, use Monte-Carlo estimates
         if(!distr %in% c("kou", "unif", "norm", "dkou"))
         {
            eta <- mean(exp(rdistr(M)))-1
         }

      } else
      {
         eta <- jumps$eta
      }
   } else
   {
      eta <- 0
   }

   # Solution vector to populate
   x <- matrix(0, nrow = n+1)
   x[1] <- IC$x0
   spot <- IC$spot
   # Fixed time-step size
   k <- t/n
   # Time grid
   tt <- seq(0, t, k)
   # Loop through time grid
   for(i in 2:(n+1))
   {
      if(!is.null(jumps))
      {
         # Simulate number of jumps in single time-step
         m <- stats::rpois(1, lambda = lambda(x[i-1], tt[i-1])*k)
         # Compound Poisson sum of the random jump-amplitudes
         cps_pois <- sum(rdistr(m))
      } else{
         cps_pois <- 0
      }
      x[i] <- x[i-1] + (mu(x[i-1], tt[i-1])-eta*lambda(x[i-1], tt[i-1]))*k+volat(x[i-1], tt[i-1])*stats::rnorm(1, 0, sd = sqrt(k))+cps_pois

   }
   if(exponential)
   {
     if(!is.null(spot))
     {
       X <- data.frame(t = tt, X = spot*exp(x))
     } else
     {
       stop("'spot' is null but 'exponential' = TRUE; please enter valid spot price")
     }

   } else
   {
     X <- data.frame(t = tt, X = x)
   }
   return(X)
 }