R/merton_family.R

In shill1729/OptionPricer: An option pricer for European and American vanilla stock options

#' probability density function for Merton dynamics
#'
#' @param x log-return
#' @param t time in trading years
#' @param rate rate of return
#' @param volat annual volatility
#' @param lambda mean rate of jumps per year
#' @param jm mean size of jumps
#' @param jv mean volatility of jumps
#' @param N number of terms to use in infinite summation
#' 
#' @description {The PDF for the log-returns of an asset following a Merton jump diffusion.}
#' @return numeric
dmerton1 <- function(x, t, rate, volat, lambda, jm, jv, N = 100)
{
  n <- 0:N
  eta <- exp(jm+0.5*jv^2)-1
  drift <- (rate-0.5*volat^2-lambda*eta)
  p <- stats::dpois(n, lambda = lambda*t)
  phi <- stats::dnorm(x, drift*t+n*jm, sqrt(t*volat^2+n*jv^2))
  return(sum(p*phi))
}


#' probability density function for Merton dynamics
#'
#' @param x log-return
#' @param t time in trading years
#' @param rate rate of return
#' @param volat annual volatility
#' @param lambda mean rate of jumps per year
#' @param jm mean size of jumps
#' @param jv mean volatility of jumps
#' @param N number of terms to use in infinite summation
#' 
#' @description {The semi-analytic PDF for the log-returns of an asset following a Merton jump diffusion.}
#' @return numeric
#' @export dmerton
dmerton <- function(x, t, rate, volat, lambda, jm, jv, N = 100)
{
  if(length(x) == 1)
  {
    return(dmerton1(x, t, rate, volat, lambda, jm, jv, N = 100))
  } else if(length(x) > 1)
  {
    m <- mapply(function(X) {
      dmerton1(X, t, rate, volat, lambda, jm, jv, N)
    }, X = x)
    return(m)
  } else{
    stop("empty argument passed in x variable")
  }
}



#' cumulative distribution function for Merton dynamics
#'
#' @param x log-return
#' @param t time in trading years
#' @param rate rate of return
#' @param volat annual volatility
#' @param lambda mean rate of jumps per year
#' @param jm mean size of jumps
#' @param jv mean volatility of jumps
#' @param N number of terms to use in infinite summation
#' 
#' @description {The semi-analytic CDF for the log-returns of an asset following a Merton jump diffusion.}
#' @return numeric
pmerton1 <- function(x, t, rate, volat, lambda, jm, jv, N = 100)
{
  n <- 0:N
  eta <- exp(jm+0.5*jv^2)-1
  drift <- (rate-0.5*volat^2-lambda*eta)
  p <- stats::dpois(n, lambda = lambda*t)
  phi <- stats::pnorm(x, drift*t+n*jm, sqrt(t*volat^2+n*jv^2))
  return(sum(p*phi))
}

#' cumulative distribution functio for Merton dynamics
#'
#' @param x log-return
#' @param t time in trading years
#' @param rate rate of return
#' @param volat annual volatility
#' @param lambda mean rate of jumps per year
#' @param jm mean size of jumps
#' @param jv mean volatility of jumps
#' @param N number of terms to use in infinite summation
#' 
#' @description {The semi-analytic CDF for the log-returns of an asset following a Merton jump diffusion.}
#' @return numeric
#' @export pmerton
pmerton <- function(x, t, rate, volat, lambda, jm, jv, N = 100)
{
  if(length(x) == 1)
  {
    return(pmerton1(x, t, rate, volat, lambda, jm, jv, N = 100))
  } else if(length(x) > 1)
  {
    m <- mapply(function(X) {
      pmerton1(X, t, rate, volat, lambda, jm, jv, N)
    }, X = x)
    return(m)
  } else{
    stop("empty argument passed in x variable")
  }
}