R/sample_path_cir.R

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

#' Simulate a sample path of an Cox-Ingersoll-Ross short-rate model starting from a given point
#'
#' @param t length of time to simulate over
#' @param IC initial value of the Cox-Ingersoll-Ross short-rate model, no default
#' @param parameters a named list containing constants \code{theta}, \code{mu}, and \code{volat} describing the mean-reversion speed, drift and volatility coefficients
#' @param jumps a list of objects defining the jump-size distribution, can be NULL for purely continuous models, see details.
#' @param n number of time sub-intervals in the discretization
#'
#' @description {Non-exported Wrapper to \code{sample_path_em} to directly simulate a sample path of Cox-Ingersoll-Ross short-rate model without having to specify the constant
#' coefficient functions. This function should not be called directly but rather through \code{sample_path} with \code{continuous.model} set to \code{"cir"}.}
#' @details {The drift of the Cox-Ingersoll-Ross short-rate model is in the form of \eqn{\theta(mu-r_t)} where \eqn{\theta} is the mean-reversion speed and \eqn{\mu} is the long-term mean level.
#' The volatility coefficient is given by \eqn{\sigma \sqrt{r_t}}
#' \itemize{
#' \item \code{distr} the name of the distribution of the jump-sizes e.g. "norm", "unif", "kou"
#' \item \code{param} named list of parameters for the distribution matching the input in \code{rdistr} for a given "distr".
#' }}
#' @return data.frame containing \code{t} (time) and \code{X} (state).
sample_path_cir <- function(t, IC, parameters, jumps = NULL, n = 10000)
{
  if(is.null(IC))
  {
    stop("Need to pass initial value list 'IC' list for 'cir'")
  } else{
    if(is.null(IC$x0))
    {
      stop("Need to pass initial value 'x0' in 'IC' list for 'cir'")
    }
  }

  theta <- parameters$theta
  mu <- parameters$mu
  volat <- parameters$volat
  if(is.null(theta))
  {
    stop("'theta' must be passed in the 'parameters' list")
  }
  if(theta <= 0)
  {
    stop("'theta' must be positive")
  }
  if(is.null(mu))
  {
    stop("'mu' must be passed in the 'parameters' list")
  }
  if(is.null(volat))
  {
    stop("'volat' must be passed in the 'parameters' list")
  }
  if(volat <= 0)
  {
    stop("'volat' must be positive")
  }
  if(2*theta*mu < volat^2)
  {
    stop("Does not preclude zero interest rates: check parameter values so that 2*'theta'*'mu'>='volat'^2")
  }
  if(is.null(parameters$lambda))
  {
    lambda <- function(x, t) 0
  } else
  {
    lambda <- parameters$lambda
  }

  coeff <- list(f.mu = function(x, t) drift_vasicek(x, t, theta, mu),
                f.vo = function(x, t) volat_cir(x, t, vcf = volat),
                f.lambda = function(x, t) lambda(x, t))
  B <- euler_maruyama(t = t, coeff = coeff, IC = IC, jumps = jumps, exponential = FALSE, n = n)
  return(B)
}