R/Igarch.R

In KevinKotze/tsm: Time Series Modelling

#' Estimation of a Gaussian IGARCH(1,1) model.
#'
#' @param rtn Time series variable.
#' @param include.mean flag for the constant in the mean equation.
#' @param volcnt flag for the constant term of the volatility equation.
#' @return Gaussian IGARCH(1,1) results.
#' @examples

"Igarch" <- function(rtn,
                     include.mean = F,
                     volcnt = F) {
  # Estimation of a Gaussian IGARCH(1,1) model.
  # rtn: return series
  # include.mean: flag for the constant in the mean equation.
  # volcnt: flag for the constant term of the volatility equation.
  #### default is the RiskMetrics model
  #
  Idata <<- rtn
  Flag <<- c(include.mean, volcnt)
  #
  Mean <- mean(Idata)
  Var <- var(Idata)
  S <- 1e-6
  if ((volcnt) && (include.mean)) {
    params <- c(mu = Mean,
               omega = 0.1 * Var,
               beta = 0.85)
    lowerBounds <- c(mu = -10 * abs(Mean),
                    omega = S ^ 2,
                    beta = S)
    upperBounds <- c(mu = 10 * abs(Mean),
                    omega = 100 * Var,
                    beta = 1 - S)
  }
  if ((volcnt) && (!include.mean)) {
    params <- c(omega = 0.1 * Var, beta = 0.85)
    lowerBounds <- c(omega = S ^ 2, beta = S)
    upperBounds <- c(omega = 100 * Var, beta = 1 - S)
  }
  #
  if ((!volcnt) && (include.mean)) {
    params <- c(mu = Mean, beta = 0.8)
    lowerBounds <- c(mu = -10 * abs(Mean), beta = S)
    upperBounds <- c(mu = 10 * abs(Mean), beta = 1 - S)
  }
  if ((!volcnt) && (!include.mean)) {
    params <- c(beta = 0.85)
    lowerBounds <- c(beta = S)
    upperBounds <- c(beta = 1 - S)
  }
  # Step 3: set conditional distribution function:
  igarchDist <- function(z, hh) {
    dnorm(x = z / hh) / hh
  }
  # Step 4: Compose log-likelihood function:
  igarchLLH <- function(parm) {
    include.mean = Flag[1]
    volcnt = Flag[2]
    mu = 0
    omega = 0
    if ((include.mean) && (volcnt)) {
      my = parm[1]
      omega = parm[2]
      beta = parm[3]
    }
    if ((!include.mean) && (volcnt)) {
      omega = parm[1]
      beta = parm[2]
    }
    if ((!include.mean) && (!volcnt))
      beta = parm[1]
    if ((include.mean) && (!volcnt)) {
      mu = parm[1]
      beta = parm[2]
    }
    #
    z <- (Idata - mu)
    Meanz <- mean(z ^ 2)
    e  <-  omega + (1 - beta) * c(Meanz, z[-length(Idata)] ^ 2)
    h <- stats::filter(e, beta, "r", init = Meanz)
    hh <- sqrt(abs(h))
    llh <- -sum(log(igarchDist(z, hh)))
    llh
  }
  # Step 5: Estimate Parameters and Compute Numerically Hessian:
  fit <- nlminb(
    start = params,
    objective = igarchLLH,
    lower = lowerBounds,
    upper = upperBounds
  )
  ##lower = lowerBounds, upper = upperBounds, control = list(trace=3))
  epsilon <- 0.0001 * fit$par
  cat("Estimates: ", fit$par, "\n")
  npar <- length(params)
  Hessian <- matrix(0, ncol = npar, nrow = npar)
  for (i in 1:npar) {
    for (j in 1:npar) {
      x1 = x2 = x3 = x4  = fit$par
      x1[i] = x1[i] + epsilon[i]
      x1[j] = x1[j] + epsilon[j]
      x2[i] = x2[i] + epsilon[i]
      x2[j] = x2[j] - epsilon[j]
      x3[i] = x3[i] - epsilon[i]
      x3[j] = x3[j] + epsilon[j]
      x4[i] = x4[i] - epsilon[i]
      x4[j] = x4[j] - epsilon[j]
      Hessian[i, j] = (igarchLLH(x1) - igarchLLH(x2) - igarchLLH(x3) + igarchLLH(x4)) /
        (4 * epsilon[i] * epsilon[j])
    }
  }
  cat("Maximized log-likehood: ", igarchLLH(fit$par), "\n")
  # Step 6: Create and Print Summary Report:
  se.coef <- sqrt(diag(solve(Hessian)))
  tval <- fit$par / se.coef
  matcoef <- cbind(fit$par, se.coef, tval, 2 * (1 - pnorm(abs(tval))))
  dimnames(matcoef) <- list(names(tval),
                           c(" Estimate",
                             " Std. Error", " t value", "Pr(>|t|)"))
  cat("\nCoefficient(s):\n")
  printCoefmat(matcoef, digits = 6, signif.stars = TRUE)

  if ((include.mean) && (volcnt)) {
    mu = fit$par[1]
    omega = fit$par[2]
    beta = fit$par[3]
  }
  if ((include.mean) && (!volcnt)) {
    mu = fit$par[1]
    beta = fit$par[2]
    omega = 0
  }
  if ((!include.mean) && (volcnt)) {
    mu = 0
    omega = fit$par[1]
    beta = fit$par[2]
  }
  if ((!include.mean) && (!volcnt)) {
    mu = 0
    omega = 0
    beta = fit$par[1]
  }
  z <- Idata - mu
  Mz <- mean(z ^ 2)
  e <- omega + (1 - beta) * c(Mz, z[-length(z)] ^ 2)
  h <- stats::filter(e, beta, "r", init = Mz)
  vol <- sqrt(abs(h))

  Igarch <- list(par = fit$par, volatility = vol)
}