R/simul.R

In ahoundetoungan/multifractal: Multifractal volatility models

#' @title Simulation From Multifractal Model
#' @description Simulates returns following multifractal model with leverage.
#' @param parms the vector of model parameters.
#' @param kbar the number of components.
#' @param N the sample size.
#' @param leverage a boolean value indicating if the model has leverage or not.
#' @param NL A parameter \eqn{N_L}{NL} as integer required for the leverage effect. The default value \code{NULL} set \eqn{N_L}{NL} to 70, which is the suggested
#' value by Augustyniak et al. (2019) for daily return (see \code{\link{fitMSM}}). 
#' @return a vector of returns
#' @examples 
#' \dontrun{
#' # Model without leverage
#' parms1 <- c("m0" = 0.97, "sigma2" = 1e-2, "b" = 2.7, "gamma1" = 0.95)
#' y1     <- simMSM(parms = parms1, kbar = 5, leverage = FALSE)
#' 
#' plot(y1, type = "l")
#' 
#' # Model with leverage
#' parms2 <- c("m0" = 0.37, "sigma2" = 1e-2, "b" = 2.7, "gamma1" = 0.95, "l1" = 1.35, "thetal" = 0.99)
#' y2     <- simMSM(parms = parms2, kbar = 5, leverage = TRUE)
#' 
#' plot(y2, type = "l")
#' }
#' @export

simMSM <- function(parms, kbar, N = 1000L, leverage = TRUE, NL = NULL) {
  stopifnot(N > 1)
  nstates        <- 2^kbar
  
  # parameters
  msg            <- paste0("Unexpected names for parms. The components should be named m0, sigma2, gamma1, b",
                           ifelse(leverage, ", l1, thetaL.", "."))
  if (is.null(names(parms))) {
    stop(msg)
  } 
  
  ex.nam.parms   <- c("M0", "SIGMA2", "GAMMA1", "B")
  names.parms    <- toupper(names(parms))
  
  if (leverage) {
    if ((length((names.parms)) != 6) | (length(parms) != 6)) {
      stop(msg)
    }
    ex.nam.parms <- c(ex.nam.parms, "L1", "THETAL")
  } else {
    if ((length((names.parms)) != 4) | (length(parms) != 4)) {
      stop(msg)
    }
  }
  
  if (!all(ex.nam.parms %in% names.parms)) {
    stop(msg)
  }
  
  m0             <- parms[names.parms == "M0"]
  sigma2         <- parms[names.parms == "SIGMA2"]
  gamma1         <- parms[names.parms == "GAMMA1"]
  b              <- parms[names.parms == "B"]
  l1             <- parms[names.parms == "L1"]
  thetaL         <- parms[names.parms == "THETAL"]

  # initial probabilities
  p1             <- rep(1, nstates)
  p1             <- p1/sum(p1)
  # markov chains states
  g              <- fg(m0, kbar)
  # transition matrix
  P              <- fP(gamma1, b, kbar)
  
  M              <- numeric(1)
  y              <- numeric(0)
  Lt             <- numeric(0)
  Lt[1]          <- 1
  M[1]           <- 1
  y[1]           <- rnorm(1, 0, sqrt(sigma2*g[1]))

  # Marvol chain
  if(leverage) {
    if (is.null(NL)) {
      NL         <- 70
    }
    
    NLst         <- min(NL, N)
    for (t in 2:NLst) {
      Lti        <- numeric(t - 1)
      for (i in 1:(t - 1)) {
        li       <- l1*thetaL^(i - 1)
        Lti[i]   <- ifelse(y[t - i] >= 0, 1, 1 + li*abs(y[t - i])/sqrt(Lt[t - i]))
      }
      Lt[t]      <- prod(Lti)
      M[t]       <- sample(1:(2^kbar), 1, prob = P[M[t - 1],])
      y[t]       <- rnorm(1, 0, sqrt(sigma2*g[M[t]]*Lt[t]))
    }
    
    if (N > NLst) {
      for (t in (NLst + 1):N) {
        Lti      <- numeric(NLst)
        for (i in 1:NLst) {
          li     <- l1*thetaL^(i - 1)
          Lti[i] <- ifelse(y[t - i] >= 0, 1, 1 + li*abs(y[t - i])/sqrt(Lt[t - i]))
        }
        Lt[t]    <- prod(Lti)
        M[t]     <- sample(1:(2^kbar), 1, prob = P[M[t - 1],])
        y[t]     <- rnorm(1, 0, sqrt(sigma2*g[M[t]]*Lt[t]))
      }
    }
  } else {
    if(!is.null(NL)) {
      warning("unused argument: NL, because leverage = FALSE")
    }
    for (t in 2:N) {
      M[t]       <- sample(1:(2^kbar), 1, prob = P[M[t - 1],])
      y[t]       <- rnorm(1, 0, sqrt(sigma2*g[M[t]]))
    }
  }
  y
}