R/simul.R

In bpvgoncalves/DynHedge: Dynamic Hedge Simulation

#' Simulation of correlated series
#'
#' @param n      Number of steps to be computed
#' @param S0     Vector of initial prices. The number of series will be deduced from the number of prices.
#' @param mu     Vector of Gaussian drift factors
#' @param sd     Vector of Standard Deviations
#' @param corr   Correlation Matrix
#' @param jumps  Optional. Use Jump-Diffusion process. Default: False
#' @param delta  Optional. Jump distribution parameter. Default: 0.001
#' @param lambda Optional. Jump frequency parameter. Default: 2
#'
#' @return       A matrix of n rows by #series columns
#' @import       stats
#' @export
#'
#' @examples
#' sim <- simulate_series(10,
#'                        c(8.88, 5.55),
#'                        c(0, 0),
#'                        c(0.02, 0.05),
#'                        matrix(c(1, 0.9, 0.9, 1), nrow = 2))
simulate_series <- function(n, S0, mu, sd, corr, jumps = FALSE, delta = 0.001, lambda = 2){

  S0 <- as.matrix(S0)
  mu <- as.matrix(mu)
  sd <- as.matrix(sd)
  num_series <- nrow(mu)
  GBMs <- matrix(nrow = n, ncol = num_series)
  colnames(GBMs) <- colnames(S0)

  Wt <- matrix(stats::rnorm(n * num_series, 0, 1), ncol = num_series)
  Wt <- apply(Wt, 2, cumsum)
  # upper triangular Cholesky decomposition
  Wt <- Wt %*% chol(corr)   # key trick for creating correlated paths
  for (i in 1:num_series) {
    if (jumps) {
      GBMs[, i] <- GBMJ(n, sd[i], mu[i] , S0[i], Wt[, i], delta, lambda)
    } else {
      GBMs[, i] <- GBM(n, sd[i], mu[i] , S0[i], Wt[, i])
    }
  }
  invisible(GBMs)
}


GBM <- function(N, sigma, mu, S0, Wt = NULL) {
  if (is.null(Wt)) {
    Wt <- cumsum(rnorm(N, 0, 1))
  }
  t <- (1:N) / 252
  p1 <- (mu - 0.5 * (sigma ** 2)) * t
  p2 <- sigma * Wt
  St = S0 * exp(p1 + p2)
  return(St)
}

GBMJ <- function(N, sigma, mu, S0, Wt = NULL, delta, lambda){
  # parameters for Poisson Process
  Nt <- stats::rpois(N, lambda)
  jump <- stats::runif(N, -delta, delta)
  p3 <- cumprod((1 + jump) ** Nt)

  # And those for the Brownian Motion
  if (is.null(Wt)) {
    Wt <- cumsum(stats::rnorm(N, 0, 1))
  }
  t <- (1:N) / 252
  p1 <- (mu - 0.5 * (sigma ** 2)) * t
  p2 <- sigma * Wt

  St <- S0 * exp(p1 + p2) * p3
  return(St)
}