R/betategarch-source.R

In betategarch: Simulation, Estimation and Forecasting of Beta-Skew-t-EGARCH Models

###########################################################
## This R file contains the functions of the betategarch
## package by Genaro Sucarrat.
##
## First created 15 September 2011, Cambridge.
## This version: 25 March 2025, Oslo.
##
## CONTENTS:
##
## rST         #generate random draws from uncentred skew t
## dST         #pdf of uncentred skew t
## STmean      #analytical mean of uncentred skew t
## STvar       #analytical variance of uncentred skew t
## STskewness  #analytical 3rd moment of standardised skew t
## STkurtosis  #analytical 4th moment of standardised skew t
##
## tegarchSim         #simulate from 1 or 2 comp tegarch
## tegarchRecursion   #recursion of 1 component
## tegarchRecursion2  #recursion of 2 component
## tegarchLogl        #log-likelihood of 1 component
## tegarchLogl2       #log-likelihood of 2 component
## tegarch            #estimate beta-skew-t-egarch
##
## coef.tegarch       #estimated coefficients
## fitted.tegarch     #fitted values (wrapper to tegarchRecursion/2)
## logLik.tegarch     #log-likelihood
## predict.tegarch    #forecast volatility n-steps ahead
## print.tegarch      #print method
## residuals.tegarch  #standardised residuals
## summary.tegarch    #summarise
## vcov.tegarch       #variance-covariance matrix
##
## tegarch-demo       #demo: estimation and forecasting of Nasdaq
###########################################################


###########################################################
## 1 DENSITY FUNCTIONS
###########################################################

##=========================================================
## Generate random draws from skew t
rST <- function(n, df=10, skew=1)
{
  zstar <- rt(n=n,df=df)
  weight <- skew/(skew + 1/skew)
  z <- runif(n, -weight, 1 - weight)
  signz <- sign(z)
  epsilon <- skew^signz
  result <- -abs(zstar)/epsilon * signz
  return(result)
} #close rST()

##=========================================================
## Skew t pdf
dST <- function(y, df=10, sd=1, skew=1, log=FALSE)
{
  if(log){
    lnumerator <- log(2) - log(skew + 1/skew)
    ldenom1 <- lbeta(0.5, df/2) + log(sd) + log(df)/2
    ldenom2 <- (df+1) * log( 1+y^2/(skew^(2*sign(y))*df*sd^2) )/2
    result <-  lnumerator - ldenom1 - ldenom2
  }else{
    numerator <- (2/(skew + 1/skew))
    denom1 <- beta(0.5, df/2) * sd * sqrt(df)
    denom2 <- (1 + y^2/(skew^(2*sign(y)) * df * sd^2 ) )^( (df+1)/2 )
    result <-  numerator/(denom1 * denom2)
  }
  return(result)
} #close dST()

##=========================================================
## Analytical mean of skew t (uncentred)
STmean <- function(df, skew=1)
{
  #(skew-1/skew) * sqrt(df) * gamma((df-1)/2)/(gamma(df/2) * sqrt(pi))
  as.numeric( (skew - 1/skew) * sqrt(df) * beta((df - 1)/2, 0.5)/pi )
} #close STmean()

##=========================================================
## Analytical mean of skew t (uncentred)
STvar <- function(df, skew=1)
{
  mueps <- STmean(df,skew=skew)
  #Eabseps2 <- df*gamma(3/2)*gamma((df-2)/2)/(gamma(1/2)*gamma(df/2))
  Eabseps2 <- df * gamma(3/2) * beta((df-2)/2,1)/sqrt(pi)
  Eeps2 <- (skew^3+1/skew^3)*Eabseps2/(skew+1/skew)
  vareps <- Eeps2 - mueps^2
  return(as.numeric(vareps))
} #close STvar()

##=========================================================
## 3rd moment of standardised skew t
STskewness <- function(df, skew=1)
{
  mueps <- STmean(df,skew=skew)
  vareps <- STvar(df,skew=skew)
  sdeps <- sqrt(vareps)
  #Eabseps3 <- df^(1.5)*gamma(2)*gamma((df-3)/2)/(gamma(1/2)*gamma(df/2))
  Eabseps3 <- df^(1.5) * beta((df-3)/2,1.5)*2/pi
  Eeps3 <- (skew^4-1/skew^4)*Eabseps3/(skew+1/skew)
  skewness <- (Eeps3-3*mueps*sdeps^2*mueps^3)/sdeps^3
  return(as.numeric(skewness))
} #close STskewness()

##=========================================================
## 4th moment of standardised skew t
STkurtosis <- function(df, skew=1)
{
  mueps <- STmean(df,skew=skew)
  vareps <- STvar(df,skew=skew)
  sdeps <- sqrt(vareps)
  Eabseps2 <- df * gamma(3/2) * beta((df-2)/2,1)/sqrt(pi)
  Eeps2 <- (skew^3+1/skew^3)*Eabseps2/(skew+1/skew)
  Eabseps3 <- df^(1.5) * beta((df-3)/2,1.5)*2/pi
  Eeps3 <- (skew^4-1/skew^4)*Eabseps3/(skew+1/skew)
  #Eabseps4 <- df^2*gamma(5/2)*gamma((df-4)/2)/(gamma(1/2)*gamma(df/2))
  Eabseps4 <- df^2 * 3/4 * beta((df-4)/2,2)
  Eeps4 <- (skew^5+1/skew^5)*Eabseps4/(skew+1/skew)
  kurtosis <- (Eeps4-4*mueps*Eeps3+6*mueps^2-3*mueps^4)/vareps^2
  return(as.numeric(kurtosis))
} #close STkurtosis()

############################################################
### OLD: C FUNCTIONS
############################################################
#
###replace exp with std::exp in "mycode" as requested by the CRAN team?
#
##1 component:
#mysig <- signature(iN="integer", omega="numeric", phi1="numeric",
#  kappa1="numeric", kappastar="numeric", df="numeric", skew2="numeric",
#  dfpluss1="numeric", mueps="numeric", y="numeric", signnegy="numeric",
#  signarg="numeric", skewterm="numeric", lambda="numeric",
#  lambdadagg="numeric", u="numeric")
#mycode <- "
#  for (int i=0; i < *iN-1; i++) {
#    signarg[i] = y[i] + *mueps * exp(lambda[i]);
#    if(signarg[i]==0){ skewterm[i] = 1; }
#    if(signarg[i] > 0){ skewterm[i] = *skew2; }
#    if(signarg[i] < 0){ skewterm[i] = pow(*skew2,-1); }
#    u[i] = *dfpluss1 * signarg[i] * y[i]/(*df * exp(2 * lambda[i]) * skewterm[i] + pow(signarg[i], 2) ) - 1;
#    lambdadagg[i+1] = *phi1 * lambdadagg[i] + *kappa1 * u[i] + *kappastar * signnegy[i] * (u[i]+1);
#    lambda[i+1] = *omega + lambdadagg[i+1];
#  }
#"
#myfns <- cfunction( list(fn=mysig), list(mycode), convention=".C") #compile
#.tegarchRecursion <- myfns[["fn"]]
#rm(mysig, mycode, myfns) #remove files
#
##2 component:
#mysig <- signature(iN="integer", omega="numeric", phi1="numeric",
#  phi2="numeric", kappa1="numeric", kappa2="numeric", kappastar="numeric",
#  df="numeric", skew2="numeric", dfpluss1="numeric", mueps="numeric",
#  y="numeric", y2="numeric", signnegy="numeric", signarg="numeric",
#  skewterm="numeric", lambda="numeric", lambda1dagg="numeric",
#  lambda2dagg="numeric", u="numeric")
#mycode <- "
#  for (int i=0; i < *iN-1; i++) {
#    signarg[i] = y[i] + *mueps * exp(lambda[i]);
#    if(signarg[i]==0){ skewterm[i] = 1; }
#    if(signarg[i] > 0){ skewterm[i] = *skew2; }
#    if(signarg[i] < 0){ skewterm[i] = pow(*skew2,-1); }
#    u[i] = *dfpluss1 * signarg[i] * y[i]/(*df * exp(2 * lambda[i]) * skewterm[i] + pow(signarg[i], 2) ) - 1;
#    lambda1dagg[i+1] = *phi1 * lambda1dagg[i] + *kappa1 * u[i];
#    lambda2dagg[i+1] = *phi2 * lambda2dagg[i] + *kappa2 * u[i] + *kappastar * signnegy[i]*(u[i]+1);
#    lambda[i+1] = *omega + lambda1dagg[i+1] + lambda2dagg[i+1];
#  }
#"
#myfns <- cfunction( list(fn=mysig), list(mycode), convention=".C") #compile
#.tegarchRecursion2 <- myfns[["fn"]]
#rm(mysig, mycode, myfns) #remove files


############################################################
## 2 BETA-SKEW-T-EGARCH FUNCTIONS
############################################################

##=========================================================
## simulate from 1 or 2 component spec
tegarchSim <- function(n, omega=0, phi1=0.95, phi2=0,
  kappa1=0.01, kappa2=0, kappastar=0, df=10, skew=1,
  lambda.initial=NULL, verbose=FALSE)
{
  if(phi2==0 && kappa2==0){
    lambda <- rep(0,n) #lambda
    lambdadagg <- rep(0,n) #lambda dagger
    if(!is.null(lambda.initial)) lambdadagg[1] <- lambda.initial[2]
  
    epsilon <- rST(n, df=df, skew=skew)
    epsilon2 <- epsilon^2
    mueps <- STmean(df, skew=skew)
    eps2muepseps <- epsilon2 - mueps*epsilon
    signeps <- sign(epsilon)
    u <- (df+1)*eps2muepseps/(df*skew^(2*signeps) + epsilon2) - 1
    signnegyupluss1 <- sign(mueps-epsilon)*(u+1)
  
    #recursion:
    fn <- function(i){
      lambdadagg[i+1] <<- phi1*lambdadagg[i] + kappa1*u[i] + kappastar*signnegyupluss1[i]
    }
    indx <- 1:I(n-1)
    lambda1long <- sapply(indx,fn)
    lambda <- omega + lambdadagg
  
    #output:
    if(verbose){
      u[n] <- NA
      sigma <- exp(lambda)
      stdev <- sigma*sqrt(STvar(df,skew=skew))
      epsilon <- epsilon - mueps
      y <- sigma*epsilon
      result <- cbind(y, sigma, stdev, lambda, lambdadagg, u,
        epsilon)
    }else{
      sigma <- exp(lambda)
      result <- sigma*(epsilon-mueps)
    } #end 1-component switch
  }else{
    lambda <- rep(0,n) #lambda
    lambda1dagg <- rep(0,n) #lambda1 dagger
    lambda2dagg <- lambda1dagg #lambda2 dagger
    if(!is.null(lambda.initial)){
      lambda1dagg[1] <- lambda.initial[2]
      lambda2dagg[1] <- lambda.initial[3]
    }
  
    epsilon <- rST(n, df=df, skew=skew)
    epsilon2 <- epsilon^2
    mueps <- STmean(df, skew=skew)
    eps2muepseps <- epsilon2 - mueps*epsilon
    signeps <- sign(epsilon)
    u <- (df+1)*eps2muepseps/(df*skew^(2*signeps) + epsilon2) - 1
    signnegyupluss1 <- sign(mueps-epsilon)*(u+1)
  
    #recursion:
    fn <- function(i){
      lambda1dagg[i+1] <<- phi1*lambda1dagg[i] + kappa1*u[i]
      lambda2dagg[i+1] <<- phi2*lambda2dagg[i] + kappa2*u[i] + kappastar*signnegyupluss1[i]
    }
    indx <- 1:I(n-1)
    lambda1long <- sapply(indx,fn)
    lambda <- omega + lambda1dagg + lambda2dagg
    if(!is.null(lambda.initial)) lambda[1] <- lambda.initial[1]
  
    #output:
    if(verbose){
      u[n] <- NA
      sigma <- exp(lambda)
      stdev <- sigma*sqrt(STvar(df,skew=skew))
      epsilon <- epsilon - mueps
      y <- sigma*epsilon
      result <- cbind(y, sigma, lambda, lambda1dagg, lambda2dagg,
        u, epsilon)
    }else{
      sigma <- exp(lambda)
      epsilon <- epsilon - mueps
      result <- sigma*epsilon
    }
  } #end 2-component switch
  return(as.zoo(result))
} #close tegarchSim()

##=========================================================
## recursion of 1 component spec
tegarchRecursion <- function(y, omega=0.1, phi1=0.4,
  kappa1=0.2, kappastar=0.1, df=10, skew=0.6,
  lambda.initial=NULL, c.code=TRUE, verbose=FALSE, aux=NULL)
{
  if( is.null(aux) ){
    aux <- NULL
    aux$iN <- length(y)
    aux$signnegy <- sign(-y)
    aux$u <- rep(0,aux$iN)
  }
  dfpluss1 <- (df+1)
  u <- aux$u
  mueps <- STmean(df, skew=skew)
  lambda <- u
  lambdadagg <- u
  if(is.null(lambda.initial)){
    lambda[1] <- omega
  }else{
    lambda[1] <- lambda.initial[1]
    lambdadagg[1] <- lambda.initial[2]
  }
  
  if( c.code ){
    signarg <- u
    skewterm <- u
    skew2 <- skew^2
    tmp <- .C("tegarchrecursion", as.integer(aux$iN), as.double(omega),
      as.double(phi1), as.double(kappa1), as.double(kappastar),
      as.double(df), as.double(skew2), as.double(dfpluss1),
      as.double(mueps), as.double(y), as.double(aux$signnegy),
      as.double(signarg), as.double(skewterm), as.double(lambda),
      as.double(lambdadagg), as.double(u), PACKAGE="betategarch")
    names(tmp) <- c("iN", "omega", "phi1", "kappa1", "kappastar", "df",
      "skew2", "dfpluss1", "mueps", "y", "signnegy", "signarg", "skewterm",
      "lambda", "lambdadagg", "u")
    u <- tmp$u
    lambda <- tmp$lambda
    lambdadagg <- tmp$lambdadagg
  }else{
    y2 <- y^2
    fn <- function(i){
      u[i] <<- dfpluss1*(y2[i]+y[i]*mueps*exp(lambda[i]))/(df*exp(2*lambda[i]) * skew^(2*sign( y[i]+mueps*exp(lambda[i]) )) + (y[i]+mueps*exp(lambda[i]))^2) - 1
      lambdadagg[i+1] <<- phi1*lambdadagg[i] + kappa1*u[i] + kappastar*aux$signnegy[i]*(u[i]+1)
      lambda[i+1] <<- omega + lambdadagg[i+1]
    }
    indx <- 1:I(aux$iN-1)
    tmp <- sapply(indx,fn)
  }
  
  #output:
  if(verbose){
    u[aux$iN] <- NA
    sigma <- exp(lambda)
    epsilon <- y/sigma
  #  sdepsilon <- sd(epsilon)
    sdepsilon <- sqrt(STvar(df, skew=skew))
    stdev <- sigma*sdepsilon
    residstd <- epsilon/sdepsilon
    result <- cbind(y,sigma,stdev,lambda,lambdadagg,u,epsilon,residstd)
  }else{ result <- lambda }
  
  return(result)

} #close tegarchRecursion()

##=========================================================
## recursion of 2 component spec
tegarchRecursion2 <- function(y, omega=0.1, phi1=0.4,
  phi2=0.2, kappa1=0.05, kappa2=0.1, kappastar=0.02, df=10,
  skew=0.6, lambda.initial=NULL, c.code=TRUE, verbose=FALSE,
  aux=NULL)
{
  if(is.null(aux)){
    aux <- NULL
    aux$iN <- length(y)
    aux$signnegy <- sign(-y)
    aux$u <- rep(0,aux$iN)
  }
  dfpluss1 <- (df+1)
  y2 <- y^2
  mueps <- STmean(df, skew=skew)
  u <- aux$u
  lambda <- u #lambda
  lambda1dagg <- rep(0,aux$iN) #lambda1dagg
  lambda2dagg <- lambda1dagg #lambda2dagg
  if(is.null(lambda.initial)){
    lambda[1] <- omega
  }else{
    lambda[1] <- lambda.initial[1]
    lambda1dagg[1] <- lambda.initial[2]
    lambda2dagg[1] <- lambda.initial[3]
  }
  
  if(c.code){
    signarg <- u
    skewterm <- u
    skew2 <- skew^2
    tmp <- .C("tegarchrecursion2", as.integer(aux$iN), as.double(omega),
      as.double(phi1), as.double(phi2), as.double(kappa1), as.double(kappa2),
      as.double(kappastar), as.double(df), as.double(skew2),
      as.double(dfpluss1), as.double(mueps), as.double(y), as.double(y2),
      as.double(aux$signnegy), as.double(signarg), as.double(skewterm),
      as.double(lambda), as.double(lambda1dagg), as.double(lambda2dagg),
      as.double(u), PACKAGE="betategarch")
    names(tmp) <- c("iN", "omega", "phi1", "phi2", "kappa1", "kappa2",
      "kappastar", "df", "skew2", "dfpluss1", "mueps", "y", "y2", "signnegy",
      "signarg", "skewterm", "lambda", "lambda1dagg", "lambda2dagg", "u")
    u <- tmp$u
    lambda1dagg <- tmp$lambda1dagg
    lambda2dagg <- tmp$lambda2dagg
    lambda <- tmp$lambda
  }else{
    fn <- function(i){
      u[i] <<- dfpluss1*(y2[i]+y[i]*mueps*exp(lambda[i]))/(df*exp(2*lambda[i]) * skew^(2*sign( y[i]+mueps*exp(lambda[i]) )) + (y[i]+mueps*exp(lambda[i]))^2) - 1
      lambda1dagg[i+1] <<- phi1*lambda1dagg[i] + kappa1*u[i]
      lambda2dagg[i+1] <<- phi2*lambda2dagg[i] + kappa2*u[i] + kappastar*aux$signnegy[i]*(u[i]+1)
      lambda[i+1] <<- omega + lambda1dagg[i+1] + lambda2dagg[i+1]
    }
    indx <- 1:I(aux$iN-1)
    tmp <- sapply(indx,fn)
  }
  
  #output:
  if(verbose){
    u[aux$iN] <- NA
    sigma <- exp(lambda)
    epsilon <- y/sigma
  #  sdepsilon <- sd(epsilon)
    sdepsilon <- sqrt(STvar(df, skew=skew))
    stdev <- sigma*sdepsilon
    residstd <- epsilon/sdepsilon
    result <- cbind(y,sigma,stdev,lambda,lambda1dagg,lambda2dagg,u,epsilon,residstd)
  }else{ result <- lambda }
  
  return(result)
} #close tegarchRecursion2()

##=========================================================
## loglikelihood of 1 component spec
tegarchLogl <- function(y, pars, lower=-Inf, upper=Inf,
  lambda.initial=NULL, logl.penalty=-1e+100, c.code=TRUE, aux=NULL)
{
  if( any(is.na(pars)) || any(pars<=aux$lower) || any(pars>=aux$upper) ){
    chk.conds <- FALSE
  }else{
    chk.conds <- TRUE
  }
  if(!aux$skew){ pars <- c(pars,1) }
  if(!aux$asym){ pars <- c(pars[1:3],0,pars[4:5]) }
  
  if( chk.conds ){
    lambda <- tegarchRecursion(y, omega=pars[1], phi1=pars[2],
      kappa1=pars[3], kappastar=pars[4], df=pars[5],
      skew=pars[6], lambda.initial=lambda.initial, c.code=c.code,
      verbose=FALSE, aux=aux)
  
    term1 <- aux$iN*( log(2)-log(pars[6]+1/pars[6])+lgamma((pars[5]+1)/2)-lgamma(pars[5]/2)-log(pi*pars[5])/2 )
  
    term2 <- sum(lambda)
  
    yterm <- y + STmean(pars[5], skew=pars[6])*exp(lambda)
    num.term <- yterm^2
    denom.term <- pars[6]^(2*sign(yterm))*pars[5]*exp(2*lambda)
    term3 <- sum( (pars[5]+1)*log(1 + num.term/denom.term)/2 )
  
    logl <- term1 - term2 - term3
  
    if(is.nan(logl) || is.na(logl) || abs(logl) == Inf) logl <- logl.penalty
  }else{ logl <- logl.penalty }
  
  return(logl)

} #close tegarchLogl()

##=========================================================
## log-likelihood of 2 component spec
tegarchLogl2 <- function(y, pars, lower=-Inf, upper=Inf,
  lambda.initial=NULL, logl.penalty=-10e+100, c.code=TRUE,
  aux=NULL)
{
  if( any(is.na(pars)) || any(pars<=aux$lower) || any(pars>=aux$upper) ){
    chk.conds <- FALSE
  }else{
    chk.conds <- TRUE
  }
  if(!aux$skew){ pars <- c(pars,1) }
  
  if(chk.conds){
    lambda <- tegarchRecursion2(y, omega=pars[1], phi1=pars[2], phi2=pars[3],
      kappa1=pars[4], kappa2=pars[5], kappastar=pars[6], df=pars[7], skew=pars[8],
      lambda.initial=lambda.initial, c.code=c.code, aux=aux)
  
    #iN <- length(y)
    term1 <- aux$iN*( log(2)-log(pars[8]+1/pars[8])+lgamma((pars[7]+1)/2)-lgamma(pars[7]/2)-log(pi*pars[7])/2 )
  
    term2 <- sum(lambda)
  
    yterm <- y + STmean(pars[7], skew=pars[8])*exp(lambda)
    num.term <- yterm^2
    denom.term <- pars[8]^(2*sign(yterm))*pars[7]*exp(2*lambda)
    term3 <- sum( (pars[7]+1)*log(1 + num.term/denom.term)/2 )
  
    logl <- term1 - term2 - term3
    if(is.nan(logl) || is.na(logl) || abs(logl) == Inf) logl <- logl.penalty
  }else{ logl <- logl.penalty }
  
  return(logl)
} #close tegarchLogl2()

##=========================================================
## estimation of tegarch: 1 and 2 comp
tegarch <- function(y, asym=TRUE, skew=TRUE, components=1,
  initial.values=NULL, lower=NULL, upper=NULL, hessian=TRUE,
  lambda.initial=NULL, c.code=TRUE, logl.penalty=NULL,
  aux=NULL, ...)
{

  ##y argument:
  y <- as.zoo(y)
  y <- na.trim(y)
  y.index <- index(y)
  y <- coredata(y)
  
  ##xts (package) related:
  if(is.matrix(y)){
    if (NCOL(y) > 1)
      stop("Dependent variable not a 1-dimensional matrix")
    y <- y[, 1]
  }
  y <- as.numeric(y)
  
  aux$asym <- asym
  aux$skew <- skew
  aux$iN <- length(y)
  aux$signnegy <- sign(-y)
  aux$u <- rep(0,aux$iN)
  
  ##1-component spec:
  if( components==1 ){

    if(is.null(initial.values)){
      initial.values <- c(0.02,0.95,0.05,0.01,10,0.98)
    }
    if(is.null(lower)){
      lower <- c(-Inf,-1+.Machine$double.eps,-Inf,-Inf,
        2+.Machine$double.eps,.Machine$double.eps)
    }
    if(is.null(upper)){
      upper <- c(Inf,1-.Machine$double.eps,Inf,Inf,Inf,Inf)
    }
    if(!aux$skew){
      initial.values <- initial.values[-6]
      lower <- lower[-6]
      upper <- upper[-6]
    }
    if(!aux$asym){
      initial.values <- initial.values[-4]
      lower <- lower[-4]
      upper <- upper[-4]
    }
    if( is.null(logl.penalty) ){
      logl.penalty <- tegarchLogl(y, initial.values,
        lower=lower, upper=upper, lambda.initial=lambda.initial,
        logl.penalty=-1e+100, c.code=TRUE, aux=aux)
    }
    objective.f <- function(pars, x=y){f <- -tegarchLogl(x,
      pars, lower=lower, upper=upper,
      lambda.initial=lambda.initial, logl.penalty=logl.penalty,
      c.code=c.code, aux=aux); f}

  } #close if( components==1 )
  
  ##2-component spec:
  if( components==2 ){

    if(is.null(initial.values)){
      initial.values <- c(0.02,0.95,0.9,0.001,0.01,0.005,10,0.98)
    }
    if(is.null(lower)){
      lower <- c(-Inf,-1+.Machine$double.eps,
        -1+.Machine$double.eps,-Inf,-Inf,-Inf,
        2+.Machine$double.eps,.Machine$double.eps)
    }
    if(is.null(upper)){
    upper <- c(Inf,1-.Machine$double.eps,
      1-.Machine$double.eps,Inf,Inf,Inf,Inf,Inf)
    }
    if(!aux$skew){
      initial.values <- initial.values[-8]
      lower <- lower[-8]
      upper <- upper[-8]
    }
    if(!aux$asym)
      stop("For identification, asym=TRUE is required when components=2")
  #  if(!aux$asym){
  #    asym=TRUE
  #  }
    if(is.null(logl.penalty)){
      logl.penalty <- tegarchLogl2(y, initial.values,
        lower=lower, upper=upper, lambda.initial=lambda.initial,
        logl.penalty=-1e+100, c.code=c.code, aux=aux)
    }
    objective.f <- function(pars, x=y){f <- -tegarchLogl2(x,
      pars, lower=lower, upper=upper, lambda.initial=lambda.initial,
      logl.penalty=logl.penalty, c.code=c.code, aux=aux); f}

  } #close if( components==2 )
  
  ##estimate:
  est <- nlminb(initial.values, objective.f, lower=lower,
    upper=upper, x=y, ...)
  est$objective <- -est$objective
  
  ##compute Hessian:
  if(hessian){
    hessian.numeric <- -optimHess(est$par, objective.f)
    est <- c(list(hessian=hessian.numeric), est)
  }
  
  ##type of model:
  model <- c(components, asym, skew)
  names(model) <- c("components", "asym", "skew")
  est <- c(list(y=zoo(y, order.by=y.index), date=date(),
    initial.values=initial.values, lower=lower, upper=upper,
    lambda.initial=lambda.initial, model=model), est)
  
  #add names:
  if(components==1){
    parnames <- c("omega", "phi1", "kappa1", "kappastar", "df", "skew")
    if(!aux$skew){ parnames <- parnames[-6] }
    if(!aux$asym){ parnames <- parnames[-4] }
  }else{
    parnames <- c("omega", "phi1", "phi2", "kappa1", "kappa2",
      "kappastar", "df", "skew")
    if(!aux$skew){ parnames <- parnames[-8] }
  #  if(!aux$asym){
  #    est$NOTE <- "2 comp spec without leverage/asymmetry not available"
  #  }
  }
  names(est$par) <- parnames
  if(hessian){
    colnames(est$hessian) <- parnames
    rownames(est$hessian) <- parnames
  }
  names(est$initial.values) <- parnames
  names(est$lower) <- parnames
  names(est$upper) <- parnames
  
  #out:
  class(est) <- "tegarch"
  return(est)

} #close tegarch()


###########################################################
## S3 METHODS
###########################################################

##=========================================================
coef.tegarch <- function(object, ...)
{
  return(object$par)
} #close coef.tegarch()

##=========================================================
## fitted values (wrapper to tegarchRecursion)
fitted.tegarch <- function(object, verbose=FALSE, ...)
{
  y.index <- index(object$y)
  y <- coredata(object$y)
  
  if(object$model["components"]==1){
    if(object$model["skew"]==0){
      pars <- c(object$par, 1)
    }else{ pars <- object$par }
    if(length(pars) < 6){
      pars <- c(pars[1:3], 0, pars[4:5])
    }
    out <- tegarchRecursion(y, omega=pars[1], phi1=pars[2],
      kappa1=pars[3], kappastar=pars[4], df=pars[5],
      skew=pars[6], lambda.initial=object$lambda.initial,
      c.code=TRUE, verbose=TRUE, aux=NULL)
    out <- zoo(out, order.by=y.index)
  }else{
    if(object$model["skew"]==0){
      pars <- c(object$par, 1)
    }else{ pars <- object$par }
    out <- tegarchRecursion2(y, omega=pars[1], phi1=pars[2],
      phi2=pars[3], kappa1=pars[4], kappa2=pars[5],
      kappastar=pars[6], df=pars[7], skew=pars[8],
      lambda.initial=object$lambda.initial, c.code=TRUE,
      verbose=TRUE, aux=NULL)
    out <- zoo(out, order.by=y.index)
  }
  if(!verbose){
    out <- out[,"stdev"]
  }
  return(out)
} #close fitted.tegarch()

##=========================================================
## return log-likelihood
logLik.tegarch <- function(object, ...)
{
  out <- object$objective
  attr(out, "nobs") <- length(object$y)
  attr(out, "df") <- length(object$par)
  class(out) <- "logLik"
  return(out)
} #close logLik.tegarch()

##=========================================================
## forecast volatility
predict.tegarch <- function(object, n.ahead=1,
  initial.values=NULL, n.sim=10000, verbose=FALSE, ...)
{
  mFit <- coredata(fitted(object, verbose=TRUE))
  n.mFit <- dim(mFit)[1]
  
  if(object$model[1]==1){
    #1-comp model parameters:
    omega <- as.numeric(object$par["omega"])
    phi1 <- as.numeric(object$par["phi1"])
    kappa1 <- as.numeric(object$par["kappa1"])
    if(object$model[2]==0){kappastar<-0}else{kappastar<-as.numeric(object$par["kappastar"])}
    df <- as.numeric(object$par["df"])
    if(object$model[3]==0){skew<-1}else{skew<-as.numeric(object$par["skew"])}
  
    #initial values:
    if(is.null(initial.values)){
      y.initial <- mFit[n.mFit,"y"]
      lambda.initial <- mFit[n.mFit,"lambda"]
      lambdadagg.initial <- mFit[n.mFit,"lambdadagg"]
    }else{
      y.initial <- initial.values$y
      lambda.initial <- initial.values$lambda
      lambdadagg.initial <- initial.values$lambdadagg
    }
    y2.initial <- y.initial^2
    mueps <- STmean(df, skew=skew)
    dfpluss1 <- df+1
    u.initial <- dfpluss1*(y2.initial+y.initial*mueps*exp(lambda.initial))/(df*exp(2*lambda.initial) * skew^(2*sign( y.initial+mueps*exp(lambda.initial) )) + (y.initial+mueps*exp(lambda.initial))^2) - 1
    gterm.initial <- kappa1*u.initial + kappastar*sign(-y.initial)*(u.initial+1)
  
    #1-step ahead:
    sigma <- exp(omega + phi1*lambdadagg.initial + gterm.initial)
    stdev <- sigma*sqrt(STvar(df, skew=skew))
  
    if(n.ahead>1){
      #simulate g-term:
      epsilon <- rST(n.sim, df=df, skew=skew)
      epsilon2 <- epsilon^2
      mueps <- STmean(df, skew=skew)
      eps2muepseps <- epsilon2 - mueps*epsilon
      signeps <- sign(epsilon)
      u <- (df+1)*eps2muepseps/(df*skew^(2*signeps) + epsilon2) - 1
      signnegyupluss1 <- sign(mueps-epsilon)*(u+1)
      gterm.sim <- kappa1*u + kappastar*signnegyupluss1
  
      for(i in 2:n.ahead){
        Eexpgterm <- rep(NA,i)
        Eexpgterm[1] <- exp(phi1^I(n.ahead-1)*gterm.initial)
        for(j in 2:n.ahead){
          Eexpgterm[j] <- mean(exp(phi1^I(n.ahead-j)*gterm.sim))
        }
        sigma[i] <- exp(omega + phi1^i * lambdadagg.initial) * prod(Eexpgterm)
        stdev[i] <- sigma[i]*sqrt(STvar(df, skew=skew))
      } #end for(i in..)
    } #end if(n.ahead>1)
  }else{
    #2-comp model parameters:
    omega <- as.numeric(object$par["omega"])
    phi1 <- as.numeric(object$par["phi1"])
    phi2 <- as.numeric(object$par["phi2"])
    kappa1 <- as.numeric(object$par["kappa1"])
    kappa2 <- as.numeric(object$par["kappa2"])
    kappastar<-as.numeric(object$par["kappastar"])
    df <- as.numeric(object$par["df"])
    if(object$model[3]==0){skew<-1}else{skew<-as.numeric(object$par["skew"])}
  
    #initial values:
    if(is.null(initial.values)){
      y.initial <- mFit[n.mFit,"y"]
      lambda.initial <- mFit[n.mFit,"lambda"]
      lambda1dagg.initial <- mFit[n.mFit,"lambda1dagg"]
      lambda2dagg.initial <- mFit[n.mFit,"lambda2dagg"]
    }else{
      y.initial <- initial.values$y
      lambda.initial <- initial.values$lambda
      lambda1dagg.initial <- initial.values$lambda1dagg
      lambda2dagg.initial <- initial.values$lambda2dagg
    }
    y2.initial <- y.initial^2
    mueps <- STmean(df, skew=skew)
    dfpluss1 <- df+1
    u.initial <- dfpluss1*(y2.initial+y.initial*mueps*exp(lambda.initial))/(df*exp(2*lambda.initial) * skew^(2*sign( y.initial+mueps*exp(lambda.initial) )) + (y.initial+mueps*exp(lambda.initial))^2) - 1
    gterm1.initial <- kappa1*u.initial
    gterm2.initial <- kappa2*u.initial + kappastar*sign(-y.initial)*(u.initial+1)
  
    #1-step ahead:
    sigma <- exp(omega+phi1*lambda1dagg.initial+phi2*lambda2dagg.initial+gterm1.initial+gterm2.initial)
    stdev <- sigma*sqrt(STvar(df, skew=skew))
  
    if(n.ahead>1){
      #simulate g-term:
      epsilon <- rST(n.sim, df=df, skew=skew)
      epsilon2 <- epsilon^2
      mueps <- STmean(df, skew=skew)
      eps2muepseps <- epsilon2 - mueps*epsilon
      signeps <- sign(epsilon)
      u <- (df+1)*eps2muepseps/(df*skew^(2*signeps) + epsilon2) - 1
      signnegyupluss1 <- sign(mueps-epsilon)*(u+1)
      gterm1.sim <- kappa1*u
      gterm2.sim <- kappa2*u + kappastar*signnegyupluss1
  
      for(i in 2:n.ahead){
        Eexpgterm <- rep(NA,i)
        Eexpgterm[1] <- exp(phi1^I(n.ahead-1)*gterm1.initial + phi2^I(n.ahead-1)*gterm2.initial)
        for(j in 2:n.ahead){
          Eexpgterm[j] <- mean(exp(phi1^I(n.ahead-j)*gterm1.sim + phi2^I(n.ahead-j)*gterm2.sim))
        }
        sigma[i] <- exp(omega + phi1^i * lambda1dagg.initial + phi2^i * lambda2dagg.initial) * prod(Eexpgterm)
        stdev[i] <- sigma[i]*sqrt(STvar(df, skew=skew))
      } #end for(i in..)
    } #end if(n.ahead>1)
  } #end if(object$model[1]==...)
  
  #out:
  if(verbose){
    out <- as.zoo(cbind(sigma, stdev))
  }else{
    out <- as.zoo(stdev)
  }
  return(out)
} #close predict.tegarch()

##=========================================================
## summarise output
print.tegarch <- function(x, ...)
{
  vcovmat <- vcov(x)
  out1 <- rbind(x$par, sqrt(diag(vcovmat)))
  rownames(out1) <- c("Estimate:", "Std. Error:")
  y.n <- length(x$y)
  bic <- (-2*x$objective + length(x$par)*log(y.n))/y.n
  out2 <- rbind(x$objective, bic)
  rownames(out2) <- c("Log-likelihood:", "BIC:")
  colnames(out2) <- ""
  
  cat("Date:", x$date, "\n")
  cat("Message (nlminb):", x$message, "\n")
  if(!is.null(x$NOTE)){
    cat("NOTE:", x$NOTE, "\n")
  }
  cat("\n")
  cat("Coefficients:\n")
  print(out1)
  print(out2)
} #close print.tegarch()

##=========================================================
## fitted values (wrapper to tegarchRecursion)
residuals.tegarch <- function(object, standardised=TRUE, ...)
{
  y.index <- index(object$y)
  y <- coredata(object$y)
  
  if(object$model["components"]==1){
    if(object$model["skew"]==0){
      pars <- c(object$par, 1)
    }else{ pars <- object$par }
    if(length(pars) < 6){
      pars <- c(pars[1:3], 0, pars[4:5])
    }
    out <- tegarchRecursion(y, omega=pars[1], phi1=pars[2],
      kappa1=pars[3], kappastar=pars[4], df=pars[5],
      skew=pars[6], lambda.initial=object$lambda.initial,
      c.code=TRUE, verbose=TRUE, aux=NULL)
    out <- zoo(out, order.by=y.index)
  }else{
    if(object$model["skew"]==0){
      pars <- c(object$par, 1)
    }else{ pars <- object$par }
    out <- tegarchRecursion2(y, omega=pars[1], phi1=pars[2],
      phi2=pars[3], kappa1=pars[4], kappa2=pars[5],
      kappastar=pars[6], df=pars[7], skew=pars[8],
      lambda.initial=object$lambda.initial, c.code=TRUE,
      verbose=TRUE, aux=NULL)
    out <- zoo(out, order.by=y.index)
  }
  if(standardised){
    out <- out[,"residstd"]
  }else{
    out <- out[,"epsilon"]
  }
  return(out)
} #close residuals.tegarch()

##=========================================================
## summarise output
summary.tegarch <- function(object, verbose=FALSE, ...)
{
  if(verbose){
    out <- object[-1]
  }else{
    if(is.null(object$hessian)){
      out <- object[-c(1,3:7)]
    }else{
      out <- object[-c(1,3:8)]
    }
  }
  return(out)
} #close summary.tegarch()

##=========================================================
## return variance-covariance matrix
vcov.tegarch <- function(object, ...)
{
  if(is.null(object$hessian)){
    aux <- list()
    aux$asym <- object$model[2]
    aux$skew <- object$model[3]
    aux$iN <- length(object$y)
    aux$signnegy <- sign(-object$y)
    aux$u <- rep(0,aux$iN)
  
    #construct objective.f:
    if(object$model[1]==1){
      logl.penalty <- tegarchLogl(object$y, object$initial.values,
        lower=object$lower, upper=object$upper,
        lambda.initial=object$lambda.initial,
        logl.penalty=-1e+100, c.code=TRUE, aux=aux)
      objective.f <- function(pars, x=object$y){f <- -tegarchLogl(x,
        pars, lower=object$lower, upper=object$upper,
        lambda.initial=object$lambda.initial,
        logl.penalty=logl.penalty, c.code=TRUE, aux=aux); f}
    }else{
      logl.penalty <- tegarchLogl2(object$y, object$initial.values,
        lower=object$lower, upper=object$upper,
        lambda.initial=object$lambda.initial,
        logl.penalty=-1e+100, c.code=TRUE, aux=aux)
      objective.f <- function(pars, x=object$y){f <- -tegarchLogl2(x,
        pars, lower=object$lower, upper=object$upper,
        lambda.initial=object$lambda.initial,
        logl.penalty=logl.penalty, c.code=TRUE, aux=aux); f}
    }
    #compute Hessian vcovmat:
    hessian.numeric <- -optimHess(object$par, objective.f)
    vcovmat <- solve(-hessian.numeric)
  
  }else{
     vcovmat <- solve(-object$hessian)
  }
  return(vcovmat)
} #close vcov.tegarch()