multifractal: Multifractal volatility models

Documented in fitMSM

# Check starting values for MSM without leverage
# INPUT
#   parms0: starting values
# OUTPUT
#   unconstrained starting value 
#' @importFrom stats var
f_init_msm <- function(parms0, y, kbar, p1, n) {
  init <- NULL
  if (is.null(parms0)) {
    nstates <- 2^kbar
    g       <- numeric(nstates)
    w       <- numeric(nstates)
    P       <- matrix(0, nrow = nstates, ncol = nstates)
    linit   <- expand.grid(list(seq(0.1, 0.9, 0.1),
                                var(y),
                                seq(0.01, 0.4, 0.025),
                                seq(1.5, 5.5, 0.5)))    
    
    llh.ini <- apply(linit, 1, function(x) floglikMSM(y, ftrans_c_nc(x), p1, kbar, n, g, w, P)) 
    init    <- as.numeric(linit[which.min(llh.ini),])
  }
  else{
    names.parms  <- toupper(names(parms0))
    ex.nam.parms <- c("M0", "SIGMA2", "GAMMA1", "B")
    if ((length((names.parms)) != 4) | (length(parms0) != 4)) {
      stop("unexpected names for parms0. the components should be named m0, sigma2, gamma1, b.")
    }
    
    if (!all(ex.nam.parms %in% names.parms)) {
      stop("unexpected names for parms0. the components should be named m0, sigma2, gamma1, b.")
    }
    init <- c("m0"     = parms0[names.parms == "M0"],
              "sigma2" = parms0[names.parms == "SIGMA2"],
              "gamma1" = parms0[names.parms == "GAMMA1"],
              "b"      = parms0[names.parms == "B"])
  }
  ncparms <- ftrans_c_nc(init)
  ncparms
}


# Check starting values for MSM without leverage
# INPUT
#   parms0: starting values
# OUTPUT
#   unconstrained starting value 

#' @importFrom stats var
f_init_msm_lev <- function(parms0, y, kbar, p1, n, NL) {
  init    <- NULL
  if (is.null(parms0)) {
    nstates <- 2^kbar
    g       <- numeric(nstates)
    w       <- numeric(nstates)
    P       <- matrix(0, nrow = nstates, ncol = nstates)
    Lt      <- numeric(n)
    linit   <- expand.grid(list(seq(0.1, 0.9, 0.1),
                                var(y),
                                seq(0.01, 0.4, 0.05),
                                seq(1.5, 5.5, 0.8),
                                seq(0.1, 1, 0.15),
                                seq(0.6, 0.95, 0.07)))    
    
    llh.ini <- apply(linit, 1, function(x) flevloglikMSM(y, ftrans_c_nc_lev(x), p1, kbar, n, NL, g, w, P, Lt)) 
    init    <- as.numeric(linit[which.min(llh.ini),])
  }
  else{
    names.parms  <- toupper(names(parms0))
    ex.nam.parms <- c("M0", "SIGMA2", "GAMMA1", "B", "L1", "THETAL")
    if ((length((names.parms)) != 6) | (length(parms0) != 6)) {
      stop("unexpected names for parms0. the components should be named m0, sigma2, gamma1, b, l1, thetaL.")
    }
    if (!all(ex.nam.parms %in% names.parms)) {
      stop("unexpected names for parms0. the components should be named m0, sigma2, gamma1, b, l1, thetaL.")
    }
    
    init <- c("m0"     = parms0[names.parms == "M0"],
              "sigma2" = parms0[names.parms == "SIGMA2"],
              "gamma1" = parms0[names.parms == "GAMMA1"],
              "b"      = parms0[names.parms == "B"],
              "l1"     = parms0[names.parms == "L1"],
              "thetaL" = parms0[names.parms == "THETAL"])
  }
  
  ncparms <- ftrans_c_nc_lev(init)
  ncparms
}

#' @title Fitting Multifractal Models
#' @description `fitMSM` is used to fits multifractal model. Unlike \pkg{MSM} package fitting function, `fitMSM` allows leverage.
#' @param y the vector of returns.
#' @param kbar the number of components.
#' @param parms0 the vector of starting values for the optimization with `NULL` as default value. The vector imputs should be named (see details).
#' @param p1 the vector of initial probabilities.
#' @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 details).
#' @param optim.method the optimisation method. The expected values are \code{"nlm"}, \code{"Nelder-Mead"}, \code{"BFGS"}, \code{"CG"}, 
#' \code{"L-BFGS-B"}, \code{"nlm"} and \code{"solnp"}. It is possible to combine severals optimization methods.
#' @return list with the following elements:
#'         \item{estimates}{a list containing the parameters estimations and other settings of the model.}
#'         \item{Lt}{A vector of the components \eqn{L_t}{Lt} of the latent variance.}
#'         \item{iterations}{the number of iterations the solver computes the likelihood before convergence.}
#'         \item{convergence}{a boolean value indicating if the estimation converges or not}
#'         \item{NL}{The set value for the parameter \eqn{N_L}{NL}.}
#'         \item{y}{the vector of returns.}
#' @details Let \eqn{\{r_t\}}{{rt}} a returns process. The multifractal model is given by
#' \deqn{r_t = V_t \varepsilon_t}{rt = Vt*et}
#' where \eqn{V_t}{Vt} is the latent variance defined by \eqn{V_t = \sigma^2 M_t}{Vt = sigma2*Mt} if the model is without leverage and 
#' \eqn{\varepsilon_t}{et} is an independent and identically distributed (iid)
#' gaussian process with mean 0 and variance 1. The process \eqn{\{M_t\}}{{Mt}} is constructed as a product of \code{Kbar} independent two-state
#' Markov chains, denoted by \eqn{M^{(i)}_t, ~ i = 1,...,\bar{K}}{M(i)t, i = 1,...,Kbar} and defined on a common state space comprised of the values 
#' \eqn{\{m_0, 2 - m_0\}}{{m0, 2 - m0}}, where \eqn{m_0 \in (0, 1)}{m0 belongs to (0,1)}. The transition probability matrix of 
#' the Markov chain \eqn{M^{(i)}_t}{M(i)t} is such that the diagonal entries are \eqn{p_i}{pi} and the off-diagonal entries \eqn{1 - p_i}{1 - pi}. 
#' The probability \eqn{p_i}{pi} is defined by \eqn{\frac{1 + (1 - \gamma_1)^{b^i}}{2}}{0.5 + 0.5*(1 - gamma1)^(b^i)}, where 
#' \eqn{\gamma_1 \in (0, 1)}{gamma1 belongs to (0, 1)} and \eqn{b} is a positive parameter.\cr
#' 
#' When the model is with leverage effect, \eqn{V_t = \sigma^2 M_t L_t}{Vt = sigma2*Mt*Lt} where, 
#' \eqn{L_t}{Lt} definition involves two new unknown parameters, \eqn{l_1}{l1} and \eqn{\theta_L}{thetaL} (see Augustyniak et al., 2019).\cr
#' 
#' The starting values should be named as \code{"m0"}, \code{"sigma2"}, \code{"gamma1"} and \code{"b"} when \code{leverage = FALSE}.
#' Additional components such that \code{"l1"} and \code{"thetaL"} should be added for \code{leverage = TRUE}. 
#' @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")
#' 
#' # Estimation
#' fit1 <- fitMSM(y1, 7, leverage = FALSE)
#' fit1$estimates$parms
#' fit1$estimates$log.likelihood
#' 
#' fit2 <- fitMSM(y2, 7, leverage = TRUE)
#' fit2$estimates$parms
#' fit2$estimates$log.likelihood
#' }
#' @importFrom utils capture.output
#' @export

fitMSM <- function(y, kbar, parms0 = NULL, p1 = NULL, leverage = TRUE, NL = NULL, optim.method = "nlm") {
  n <- length(y)
  
  # initial probabilities
  if (is.null(p1)) {
    p1 <- rep(1, 2^kbar)
  }
  if (length(p1) != 2^kbar) {
    stop("the length of P1 should be 2^kbar")
  }
  p1 <- p1 / sum(p1)
  
  # initial solution
  ncparms <- NULL
  if (leverage) {
    if(is.null(NL)) {
      NL    <- 70
    }
    ncparms <- f_init_msm_lev(parms0, y, kbar, p1, n, NL)
  }
  else {
    if(!is.null(NL)) {
      NL    <- NULL
      warning("unused argument: NL, because leverage = FALSE")
    }
    ncparms <- f_init_msm(parms0, y, kbar, p1, n)
  }
  
  # Available Methods
  AVAI.ME <- c("NELDER-MEAD", "BFGS", "CG", "L-BFGS-B", "NLM", "SOLNP")
  
  if (!all(toupper(optim.method) %in% AVAI.ME)) {
    stop("Available optimization methods are Nelder-Mead, BFGS, CG, L-BFGS-B, nlm, solnp")
  }
  
  id.me   <- which(AVAI.ME %in% toupper(optim.method))
  met.opt <- c("Nelder-Mead", "BFGS", "CG", "L-BFGS-B", "nlm", "solnp")
  n.met   <- length(id.me)
  if (length(n.met) == 0) {
    stop("No optimization method declared")
  }
  
  # optimizer function
  fun_opt <- c(rep(list(f_opt_optim), 4), list(f_opt_nlm, f_opt_solnp))[id.me]
  fun.lik <- list(floglikMSM, flevloglikMSM)
  
  # optimization
  tmp <- list()
  j   <- 1
  invisible(capture.output(
    try({
      for (i in 1:n.met) {
        tmp[[j]] <- fun_opt[[i]](ncparms, fun.lik,  y, p1, kbar, n, NL, met.opt[id.me[i]], leverage)
        tmp[[j]]$method <- met.opt[id.me[i]]
        j <- j + 1
      }
    })
  ))
  
  n.met.tmp      <- length(tmp)
  log.likelihood <- unlist(lapply(1:n.met.tmp, function(x) tmp[[x]]$estimates$log.likelihood))
  id.best.met    <- which.max(log.likelihood)
  out            <- tmp[[id.best.met]]
  out$NL         <- NL
  out$y          <- y
  class(out)     <- "fitMSM"
  out
}