R/predict.R

In ahoundetoungan/multifractal: Multifractal volatility models

#' @title Forecasting from multifractal fits
#' @description Forecasts from multifractal fits.
#' @param object a fitted object of class inheriting from `fitMSM`.
#' @param h the prediction horizon.
#' @param n.sim the number of simulation to generate for each period of forecasting.
#' @param ...	further arguments passed to or from other methods.
#' @return a matrix of dimension \code{h}\eqn{\times}{*}\code{nsim} of simulations. 
#' @examples 
#' \dontrun{
#' pred <- predict(fit, h = 5, n.sim = 3e3) }
#' @importFrom stats rnorm
#' @export
predict.fitMSM <- function(object, h, n.sim, ...) {
  parms     <- object$estimates$parms
  g         <- object$estimates$g
  w         <- object$estimates$w
  P         <- object$estimates$P
  n.stes    <- length(g)
  stes      <- 1:n.stes
  Lt        <- object$Lt
  n         <- length(Lt)
  y         <- object$y
  sigma2    <- parms["sigma2"]
  NL        <- object$NL
  
  sim.st     <- matrix(data = 0, nrow = h, ncol = n.sim)
  sim.st[1,] <- sample(x = stes, size = n.sim, replace = TRUE, prob = w)
  
  # simule the states
  if (h > 1) {
    for (i in 2:h) {
      for (j in 1:n.sim) {
        sim.st[i,j] <- sample(x = stes, size = 1L, replace = FALSE, prob = P[sim.st[i - 1,j],])
      }
    }
  }
  
  # simule the returns
  # without leverage
  sims <- NULL
  if (is.null(Lt)) {
    sims <- rnorm(n = h * n.sim, mean = 0, sd = sqrt(sigma2*g[sim.st]))
    sims <- matrix(data = sims, nrow = h, ncol = n.sim)
  }
  # with leverage
  else {
    sims          <- matrix(data = NA, nrow = n + h, ncol = n.sim)
    sim.lt        <- matrix(data = NA, nrow = n + h, ncol = n.sim)
    sims[1:n,]    <- y
    sim.lt[1:n,]  <- Lt
    theta         <- parms["thetaL"]
    l1            <- parms["l1"]
    for (j in 1:n.sim) {
      t                <- n + 1
      while (t <= min(NL, n + h)) {
        Lti            <- numeric(t - 1)
        for (i in 1:(t - 1)) {
          li           <- l1*theta^(i - 1)
          Lti[i]       <- ifelse(sims[t - i, j] >= 0, 1, 1 + li*abs(sims[t - i, j])/sqrt(sim.lt[t - i, j]))
        }
        sim.lt[t, j]   <- prod(Lti)
        sims[t, j] <- rnorm(1, 0, sqrt(sigma2*g[sim.st[t - n, j]]*sim.lt[t, j]))
        t <- t + 1
      }
      while (t <= (n + h)) {
        Lti <- numeric(NL) 
        for (i in 1:NL) {
          li      <- l1*theta^(i - 1)
          Lti[i]  <- ifelse(sims[t - i, j] >= 0, 1, 1 + li*abs(sims[t - i, j])/sqrt(sim.lt[t - i, j]))
        }
        sim.lt[t, j] <- prod(Lti)
        sims[t, j]   <- rnorm(1, 0, sqrt(sigma2*g[sim.st[t - n, j]]*sim.lt[t, j]))
        t            <- t + 1
      }
    }
    sims <- as.matrix(sims[-(1:n),,drop=FALSE])
  }
  sims
}