Nothing
R/rgarch-gjrgarch.R
In rgarch: Flexible GARCH modelling in R
#################################################################################
##
## R package rgarch by Alexios Ghalanos Copyright (C) 2008, 2009, 2010, 2011
## This file is part of the R package rgarch.
##
## The R package rgarch is free software: you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation, either version 3 of the License, or
## (at your option) any later version.
##
## The R package rgarch is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
## GNU General Public License for more details.
##
#################################################################################
#---------------------------------------------------------------------------------
# SECTION gjrGARCH spec
#---------------------------------------------------------------------------------
.gjrgarchspec = function(vmodel, mmodel, dmodel, start.pars, fixed.pars)
{
# extract values:
# mean model
include.mean = mmodel$include.mean
armap = mmodel$armaOrder[1]
armaq = mmodel$armaOrder[2]
garchInMean = mmodel$garchInMean
inmeanType = mmodel$inMeanType
# inmeanType(1) means sd
# inmeanType(2) means var
#if(garchInMean && (inmeanType!=1 | inmeanType!=2)) stop("valid inmeanType values are 1 and 2")
im = inmeanType/2
arfima = mmodel$arfima
if(!is.null(mmodel$external.regressors)){
include.mex = TRUE
mexdata = as.matrix(mmodel$external.regressors)
mxn = dim(mexdata)[2]
} else{
include.mex = FALSE
mexdata = NULL
mxn = 0
}
# vmodel
garchp = vmodel$garchOrder[1]
garchq = vmodel$garchOrder[2]
if(!is.null(vmodel$external.regressors)){
include.vex = TRUE
vexdata = as.matrix(vmodel$external.regressors)
vxn = dim(vexdata)[2]
} else{
include.vex = FALSE
vexdata = NULL
vxn = 0
}
include.omega = ifelse(vmodel$targeting, FALSE, TRUE)
if(!include.omega) fixed.pars$omega = 0
# dmodel/skew/shape
include.dlambda = dmodel$include.dlambda
include.skew = dmodel$include.skew
include.shape = dmodel$include.shape
model = vmodel$model
submodel = vmodel$submodel
maxOrder = max(c(garchp,garchq))
maxOrder2 = max(c(garchp, garchq, armap, armaq))
vmodel = list(model = model, submodel = submodel, garchOrder=vmodel$garchOrder,
include.vex = include.vex, vexdata = vexdata, vxn = vxn, pow=2, include.omega = include.omega)
mmodel = list(include.mean = include.mean, armaOrder = mmodel$armaOrder,
garchInMean = garchInMean, inMeanType = inmeanType, im = im, arfima = arfima,
include.mex = include.mex, mexdata = mexdata, mxn = mxn)
#specnames = c("mu","ar","ma","darfima","mxreg","omega","alpha","beta","vxreg")
# match
modelnames = c(
if(include.mean) "mu",
if(armap > 0) paste("ar", 1:armap, sep = ""),
if(armaq > 0) paste("ma", 1:armaq, sep = ""),
if(garchInMean) paste("inmean"),
if(arfima) paste("darfima"),
if(include.mex) paste("mxreg", 1:mxn, sep = ""),
paste("omega"),
if(garchp > 0) paste("alpha", 1:garchp, sep = ""),
if(garchp > 0) paste("gamma", 1:garchp, sep = ""),
if(garchq > 0) paste("beta", 1:garchq, sep = ""),
if(include.vex) paste("vxreg", 1:vxn, sep = ""),
if(include.dlambda) "dlambda",
if(include.skew) "skew",
if(include.shape) "shape")
pos=1
pos.matrix = matrix(0, ncol = 3, nrow = 14)
colnames(pos.matrix) = c("start","stop","include")
rownames(pos.matrix) = c("mu", "ar", "ma", "inmean", "darfima", "mxreg", "omega",
"alpha", "gamma", "beta", "vxreg", "dlambda", "skew", "shape")
if(include.mean){ pos.matrix[1,1:3] = 1
pos=pos.matrix[1,2]+1}
if(armap>0){pos.matrix[2,1:3] = c(pos,pos+armap-1,1)
pos=max(pos.matrix[1:2,2])+1}
if(armaq>0){pos.matrix[3,1:3] = c(pos,pos+armaq-1,1)
pos=max(pos.matrix[1:3,2])+1}
if(garchInMean){pos.matrix[4,1:3] = c(pos,pos,1)
pos=max(pos.matrix[1:4,2])+1}
if(arfima){pos.matrix[5,1:3] = c(pos,pos,1)
pos=max(pos.matrix[1:5,2])+1}
if(include.mex){pos.matrix[6,1:3] = c(pos,pos+mxn-1,1)
pos=max(pos.matrix[1:6,2])+1}
pos.matrix[7,1:3] = c(pos,pos,1)
pos=max(pos.matrix[1:7,2])+1
if(garchp > 0){pos.matrix[8,1:3] = c(pos,pos+garchp-1,1)
pos=max(pos.matrix[1:8,2])+1}
if(garchp > 0){pos.matrix[9,1:3] = c(pos,pos+garchp-1,1)
pos=max(pos.matrix[1:9,2])+1}
if(garchq > 0){pos.matrix[10,1:3] = c(pos,pos+garchq-1,1)
pos=max(pos.matrix[1:10,2])+1}
if(include.vex){pos.matrix[11,1:3] = c(pos,pos+vxn-1,1)
pos=max(pos.matrix[1:11,2])+1}
if(include.dlambda){pos.matrix[12,1:3] = c(pos,pos,1)
pos=max(pos.matrix[1:12,2])+1}
if(include.skew){pos.matrix[13,1:3] = c(pos,pos,1)
pos=max(pos.matrix[1:13,2])+1}
if(include.shape){pos.matrix[14,1:3] = c(pos,pos,1)}
nn = length(modelnames)
modelmatrix = matrix(0, ncol = 3, nrow = nn)
rownames(modelmatrix) = modelnames
colnames(modelmatrix) = c("opt","fixed","start")
fixed.names = names(fixed.pars)
fp=charmatch(fixed.names,modelnames)
if(!is.null(fixed.names) && any(!is.na(fp))){
fixed = fp[!is.na(fp)]
modelmatrix[fixed,2]=1
fz=charmatch(modelnames,fixed.names)
fz=fz[!is.na(fz)]
fixed.pars=fixed.pars[fz]
names(fixed.pars) = fixed.names[fz]
} else{
fixed.pars=NULL
}
modelmatrix[,1] = 1
start.names = names(start.pars)
sp=charmatch(start.names,modelnames)
if(!is.null(start.names) && any(!is.na(sp))){
start = sp[!is.na(sp)]
modelmatrix[start,3]=1
sz=charmatch(modelnames,start.names)
sz=sz[!is.na(sz)]
start.pars=start.pars[sz]
} else{
start.pars=NULL
}
# the way we handle fixed pars will be through upper and lower bound constraints
opt.pars = vector(mode="numeric",length=length(modelnames))
names(opt.pars) = modelnames
# optimization model
omodel = list(modelnames = modelnames, opt.pars = opt.pars, start.pars = start.pars,
fixed.pars = fixed.pars, modelmatrix = modelmatrix, maxOrder = maxOrder,
maxOrder2 = maxOrder2, pos.matrix = pos.matrix)
ans = new("gjrGARCHspec",
mean.model = mmodel,
variance.model = vmodel,
distribution.model = dmodel,
optimization.model = omodel)
return(ans)
}
#---------------------------------------------------------------------------------
# SECTION gjrGARCH fit
#---------------------------------------------------------------------------------
.gjrgarchfit = function(spec, data, out.sample = 0, solver = "solnp", solver.control = list(),
fit.control = list(stationarity = 1, fixed.se = 0, scale = 0))
{
tic = Sys.time()
if(is.null(solver.control$trace)) trace = 0 else trace = solver.control$trace
if(is.null(fit.control$stationarity)) fit.control$stationarity = 1
if(is.null(fit.control$fixed.se)) fit.control$fixed.se = 0
if(is.null(fit.control$scale)) fit.control$scale = 0
# if we have external regressors in variance turn off scaling
if(spec@variance.model$include.vex) fit.control$scale = 0
# at the moment we do not allow scaling with fixed parameters
if(!is.null(spec@optimization.model$fixed.pars)) fit.control$scale = 0
xdata = .extractdata(data)
if(!is.numeric(out.sample)) stop("\nugarchfit-->error: out.sample must be numeric\n")
if(as.numeric(out.sample)<0) stop("\nugarchfit-->error: out.sample must be positive\n")
n.start = round(out.sample,0)
n = length(xdata$data)
if((n-n.start)<100) stop("\nugarchfit-->error: function requires at least 100 data\n points to run\n")
data = xdata$data[1:(n-n.start)]
dates = xdata$pos[1:(n-n.start)]
origdata = xdata$data
origdates = xdata$pos
dformat = xdata$dformat
# get current environment
tempenvir = new.env(hash = TRUE, parent = .GlobalEnv)
# expand the spec object and assign spec lists
vmodel = spec@variance.model
mmodel = spec@mean.model
dmodel = spec@distribution.model
omodel = spec@optimization.model
if(mmodel$include.mex) {
mmodel$origmexdata = mmodel$mexdata
mmodel$mexdata = mmodel$mexdata[1:(n-n.start),]
} else{
mmodel$origmexdata = NULL
}
if(vmodel$include.vex) {
vmodel$origvexdata = vmodel$vexdata
vmodel$vexdata = vmodel$vexdata[1:(n-n.start),]
} else{
vmodel$origvexdata = NULL
}
assign("vmodel", vmodel, envir = tempenvir)
assign("mmodel", mmodel, envir = tempenvir)
assign("dmodel", dmodel, envir = tempenvir)
assign("omodel", omodel, envir = tempenvir)
assign("dates", dates, envir = tempenvir)
assign("fit.control", fit.control, envir = tempenvir)
assign("solver", solver, envir = tempenvir)
assign("trace", trace, envir = tempenvir)
include.mean = mmodel$include.mean
include.omega = vmodel$include.omega
m = omodel$maxOrder
T = length(as.numeric(data))
dist = dmodel$distribution
model = vmodel$model
submodel = vmodel$submodel
if(fit.control$scale) dscale = sd(data) else dscale = 1
assign("dscale", dscale, envir = tempenvir)
zdata = data/dscale
# Optimization Starting Parameters Vector & Bounds
ipars = .garchstart(data = zdata, model = model, submodel = submodel, garchenv = tempenvir)
pars = ipars$pars
LB = ipars$modelLB
UB = ipars$modelUB
assign("npar", length(pars), envir = tempenvir)
# we now split out any fixed parameters
if(any(LB==UB)){
fixid = which(LB==UB)
if(length(fixid)==length(pars)){
if(fit.control$fixed.se==0) {
# if all parameters are fixed an no standard erros are to
# be calculated then we return a ugarchfilter object
cat("\nugarchfit-->warning: all parameters fixed...returning ugarchfilter
object instead\n")
return(ugarchfilter(data = data, spec = spec))
} else{
# if all parameters are fixed but we require standard errors, we
# skip the solver
use.solver = 0
fixed = pars[fixid]
pars = pars[-fixid]
}
} else{
# with some parameters fixed we extract them (to be rejoined at end)
# so that they do not enter the solver
use.solver = 1
fixed = pars[fixid]
pars = pars[-fixid]
LB = LB[-fixid]
UB = UB[-fixid]
}
} else{
# no fixed parameters, all normal
fixed = NULL
fixid = NULL
use.solver=1
}
assign("fixed", fixed, envir = tempenvir)
assign("fixid", fixid, envir = tempenvir)
# start counter
assign(".llh", 1, envir = tempenvir)
# scale parameters (for non-solnp solvers)
parscale = rep(1,length=length(pars))
names(parscale) = names(pars)
if(include.mean) parscale["mu"] = abs(mean(zdata))
if(include.omega) parscale["omega"] = var(zdata)
assign("omega.calculate", TRUE, envir = tempenvir)
# assisgn solver constraints (solnp directly else exterior type penalty
# for other solvers)
if(fit.control$stationarity == 1 && vmodel$include.vex == 0){
cb = .garchconbounds()
Ifn = .gjrgarchcon
ILB = cb$LB
IUB = cb$UB
if(solver == "solnp" | solver == "gosolnp") fit.control$stationarity = 0
} else{
Ifn = NULL
ILB = NULL
IUB = NULL
}
# conditions controls the non-solnp solver penalty
assign("fit.control", fit.control, , envir = tempenvir)
if(use.solver){
solution = .garchsolver(solver, pars, fun = .gjrgarchLLH, Ifn, ILB, IUB,
gr = NULL, hessian = NULL, parscale = parscale,
control = solver.control, LB = LB, UB = UB, ux = NULL, ci = NULL,
mu = NULL, data = zdata, returnType = "llh", garchenv = tempenvir)
sol = solution$sol
hess = solution$hess
timer = Sys.time()-tic
opt.pars = sol$par
if(!include.omega){
fixed = get("fixed", tempenvir)
fixed["omega"] = get("omega", tempenvir)
fixed["omega"] = (fixed["omega"]) * dscale^2
assign("fixed", fixed, envir = tempenvir)
}
assign("omega.calculate", FALSE, envir = tempenvir)
if(!is.null(opt.pars)) names(opt.pars) = names(pars)
if(is.null(fixed)){
if(include.mean) opt.pars["mu"] = opt.pars["mu"]*dscale
if(mmodel$include.mex){
opt.pars[omodel$pos[6,1]:omodel$pos[6,2]] = opt.pars[omodel$pos[6,1]:omodel$pos[6,2]] * dscale
}
opt.pars["omega"] = (opt.pars["omega"]) * dscale^2
}
convergence = sol$convergence
} else{
hess = NULL
timer = Sys.time()-tic
opt.pars = NULL
convergence = 0
sol = list()
sol$message = "all parameters fixed"
}
fit = list()
model = list()
assign("pow", 2, envir = tempenvir)
# check convergence else write message/return
if(convergence==0){
if(!is.null(fixed) && fit.control$fixed.se == 1){
fixed = get("fixed", tempenvir)
if(!include.omega){
vmodel$include.omega = TRUE
assign("vmodel", vmodel, envir = tempenvir)
}
npar = get("npar", tempenvir)
pall = vector(mode = "numeric", length = npar)
pall[fixid] = fixed
pall[-fixid] = opt.pars
opt.pars = pall
Names = omodel$modelnames
names(opt.pars) = Names
assign("fixed", NULL, envir = tempenvir)
assign("fixid", NULL, envir = tempenvir)
}
fit = .makefitmodel(garchmodel = "gjrGARCH", f = .gjrgarchLLH, opt = opt.pars,
data = data, T = T, m = m, timer = timer, convergence = convergence,
message = sol$message, hess, garchenv = tempenvir)
fit$dates = dates
fit$date.format = dformat
fit$origdata = origdata
fit$origdates = origdates
} else{
fit$message = sol$message
fit$convergence = 1
}
# make model list to return some usefule information which
# will be called by other functions (show, plot, sim etc)
model = .makemodel(tempenvir)
model$n.start = n.start
# return object gjrGARCHfit
ans = new("gjrGARCHfit",
fit = fit,
model = model)
rm(tempenvir)
return(ans)
}
#---------------------------------------------------------------------------------
# SECTION gjrGARCH LLH
#---------------------------------------------------------------------------------
.gjrgarchLLH = function(pars, data, returnType = "llh", garchenv)
{
# prepare inputs
# rejoin fixed and pars
fixed = get("fixed", garchenv)
fixid = get("fixid", garchenv)
npar = get("npar", garchenv)
trace = get("trace", garchenv)
if(!is.null(fixed)){
pall = vector(mode = "numeric", length = npar)
pall[fixid] = fixed
pall[-fixid] = pars
pars = pall
}
T = length(data)
vmodel = get("vmodel", garchenv)
mmodel = get("mmodel", garchenv)
dmodel = get("dmodel", garchenv)
omodel = get("omodel", garchenv)
fit.control = get("fit.control", garchenv)
armaOrder = c(mmodel$armaOrder,
mmodel$im,
mmodel$arfima,
mmodel$include.mex)
m = omodel$maxOrder
N = c(m,T)
garchOrder = c(vmodel$garchOrder, vmodel$include.vex)
include.omega = vmodel$include.omega
omega.calculate = get("omega.calculate", garchenv)
include.inmean = mmodel$garchInMean
mexdata = mmodel$mexdata
vexdata = vmodel$vexdata
distribution = dmodel$distribution
# distribution number
# 1 = norm, 2=snorm, 3=std, 4=sstd, 5=ged, 6=sged, 7=nig
dist = dmodel$distno
pos = omodel$pos.matrix
omega = 0.00001
mu = ar = ma = inmean = darfima = mxreg = alpha = beta = gm = vxreg = skew = shape = dlambda = 0
if(pos[1,3]) mu = pars[pos[1,1]:pos[1,2]]
if(pos[2,3]) ar = pars[pos[2,1]:pos[2,2]]
if(pos[3,3]) ma = pars[pos[3,1]:pos[3,2]]
if(pos[4,3]) inmean = pars[pos[4,1]:pos[4,2]]
if(pos[5,3]) darfima = pars[pos[5,1]:pos[5,2]]
if(pos[6,3]) mxreg = pars[pos[6,1]:pos[6,2]]
if(pos[8,3]) alpha = pars[pos[8,1]:pos[8,2]]
if(pos[9,3]) gm = pars[pos[9,1]:pos[9,2]]
if(pos[10,3]) beta = pars[pos[10,1]:pos[10,2]]
if(pos[11,3]) vxreg = pars[pos[11,1]:pos[11,2]]
if(pos[12,3]) dlambda = pars[pos[12,1]:pos[12,2]]
if(pos[13,3]) skew = pars[pos[13,1]:pos[13,2]]
if(pos[14,3]) shape = pars[pos[14,1]:pos[14,2]]
#rx = .getresiduals2(data = data, include.mean = mmodel$include.mean,
# mu = mu, armaOrder=armaOrder, ar = ar, ma = ma, include.inmean,
# inmean = inmean, darfima = darfima, mexdata = mexdata, mxreg = mxreg,
# h = get("meqspars", garchenv)$h)
hm = get("meqspars", garchenv)$h
rx = .arfimaxfilter(armaOrder = armaOrder, mu = mu, ar = ar, ma = ma, inmean = inmean,
darfima = darfima, mxreg = mxreg, mexdata = mexdata, h = hm,
data = data, N = N, garchenv)
res = rx$res
zrf = rx$zrf
res = res[!is.na(res) && is.finite(res) && !is.nan(res)]
#names(pars) = omodel$modelnames
names(pars) = omodel$modelnames
# gjrgarch persistence value
persist = sum(beta)+ sum(alpha)+sum(apply(as.data.frame(gm),1,FUN=function(x)
x*pneg(dlambda, shape, skew, distribution)))
# unconditional sigma value
mvar = mean(res*res)
if(pos[11,3] == 1) {
mv = sum(apply(matrix(vexdata, ncol = vmodel$vxn), 2, "mean")*vxreg)
} else{
mv = 0
}
if(include.omega){
omega = max(eps, pars[pos[7,1]:pos[7,2]])
hEst = mvar
} else{
if(omega.calculate){
persist = .persistgjrgarch(pars = pars, distribution = distribution)
hEst = mvar
omega = hEst * (1 - persist) - mv
assign("omega", omega, envir = garchenv)
} else{
omega = get("omega", envir = garchenv)
hEst = mvar
}
}
# if we have external regressors in variance equation we cannot have
# stationarity checks in likelihood
if(fit.control$stationarity == 1 && vmodel$include.vex == 0){
if(!is.na(persist) && persist >= 1) return(llh = get(".llh", garchenv) +
0.1*(abs(get(".llh", garchenv))))
}
# C Code [gjrgarch.c] via wrapper
ans = .gjrgarchfitC(armaOrder, garchOrder, mu, ar, ma, inmean, mxreg, mexdata,
omega, hEst, alpha, beta, gm, vxreg, vexdata, dlambda, skew, shape, dist, data,
zrf, N, res, garchenv)
if(get(".csol", garchenv) == 1){
if( trace > 0 ) cat(paste("\nugarchfit-->warning: ", get(".filtermessage",garchenv),"\n", sep=""))
return(llh = (get(".llh", garchenv) + 0.1*(abs(get(".llh", garchenv)))))
}
z = ans$z
h = ans$h
epsx = ans$res
llh = ans$llh
if(is.finite(llh) && !is.na(llh) && !is.nan(llh)){
assign(".llh", llh, envir = garchenv)
} else {
llh = (get(".llh", garchenv) + 0.1*(abs(get(".llh", garchenv))))
}
# LHT = raw scores
LHT = -ans$LHT
ans = switch(returnType,
llh = llh,
LHT = LHT,
all = list(llh = llh, h = h, epsx = epsx, z = z, kappa = kappa, LHT = LHT, persistence = persist))
return(ans)
}
#---------------------------------------------------------------------------------
# SECTION gjrGARCH filter
#---------------------------------------------------------------------------------
.gjrgarchfilter = function(spec, data, out.sample = 0, n.old = NULL)
{
.garchenv = environment()
if(!is.null(n.old)) Nx = n.old else Nx = length(data)
xdata = .extractdata(data)
data = xdata$data
dates = xdata$pos
dformat = xdata$dformat
origdata = data
origdates = dates
T = length(origdata) - out.sample
data = origdata[1:T]
dates = origdates[1:T]
vmodel = spec@variance.model
mmodel = spec@mean.model
dmodel = spec@distribution.model
omodel = spec@optimization.model
include.mean = mmodel$include.mean
modelnames = omodel$modelnames
pars = unlist(spec@optimization.model$fixed.pars)
parnames = names(pars)
if(length(modelnames) != length(pars))
stop("parameters do not match expected specification length")
if(is.na(all(match(modelnames, parnames), 1:length(modelnames)))) {
stop("parameters names do not match specification")
} else{
pars = pars[modelnames]
}
m = omodel$maxOrder
T = length(as.numeric(data))
dist = dmodel$distribution
model = vmodel$model
submodel = vmodel$submodel
armaOrder = c(mmodel$armaOrder,
mmodel$im,
mmodel$arfima,
mmodel$include.mex)
garchOrder = c(vmodel$garchOrder,
vmodel$include.vex)
include.omega = vmodel$include.omega
include.inmean = mmodel$garchInMean
mexdata = mmodel$mexdata
vexdata = vmodel$vexdata
distribution = dmodel$distribution
garchOrder = c(vmodel$garchOrder, vmodel$include.vex)
include.inmean = mmodel$garchInMean
mexdata = mmodel$mexdata
vexdata = vmodel$vexdata
distribution = dmodel$distribution
# distribution number
# 1 = norm, 2=snorm, 3=std, 4=sstd, 5=ged, 6=sged, 7=nig
dist = dmodel$distno
N = c(m,T)
pos = omodel$pos.matrix
omega = 0.00001
mu = ar = ma = inmean = darfima = mxreg = alpha = beta = gm = vxreg = skew = shape = dlambda = 0
if(pos[1,3]) mu = pars[pos[1,1]:pos[1,2]]
if(pos[2,3]) ar = pars[pos[2,1]:pos[2,2]]
if(pos[3,3]) ma = pars[pos[3,1]:pos[3,2]]
if(pos[4,3]) inmean = pars[pos[4,1]:pos[4,2]]
if(pos[5,3]) darfima = pars[pos[5,1]:pos[5,2]]
if(pos[6,3]) mxreg = pars[pos[6,1]:pos[6,2]]
if(pos[8,3]) alpha = pars[pos[8,1]:pos[8,2]]
if(pos[9,3]) gm = pars[pos[9,1]:pos[9,2]]
if(pos[10,3]) beta = pars[pos[10,1]:pos[10,2]]
if(pos[11,3]) vxreg = pars[pos[11,1]:pos[11,2]]
if(pos[12,3]) dlambda = pars[pos[12,1]:pos[12,2]]
if(pos[13,3]) skew = pars[pos[13,1]:pos[13,2]]
if(pos[14,3]) shape = pars[pos[14,1]:pos[14,2]]
rx = .arfimaxfilter(armaOrder = armaOrder, mu = mu, ar = ar, ma = ma, inmean = inmean,
darfima = darfima, mxreg = mxreg, mexdata = mexdata, h = 0,
data = data, N = N, garchenv = .garchenv)
res = rx$res
zrf = rx$zrf
if(!is.null(n.old)){
rx2 = .arfimaxfilter(armaOrder = armaOrder, mu = mu, ar = ar, ma = ma, inmean = inmean,
darfima = darfima, mxreg = mxreg, mexdata = mexdata[1:Nx, , drop = FALSE], h = 0,
data = origdata[1:Nx], N = c(m, Nx), garchenv = .garchenv)
res2 = rx2$res
# unconditional sigma value
mvar = mean(res2*res2)
} else{
mvar = mean(res*res)
}
# gjrgarch persistence value
persist = .persistgjrgarch(pars, distribution)
#if(persist>1) stop("\nPersistence of model is > 1...try again!\n")
#hEst = omega + persistence*mvar
if(pos[11,3] == 1) {
mv = sum(apply(matrix(vexdata, ncol = vmodel$vxn), 2, "mean")*vxreg)
} else{
mv = 0
}
if(include.omega){
omega = pars[pos[7,1]:pos[7,2]]
hEst = mvar
} else{
omega = hEst * (1 - persist) - mv
hEst = mvar
}
ans = .gjrgarchfitC(armaOrder, garchOrder, mu, ar, ma, inmean, mxreg, mexdata,
omega, hEst, alpha, beta, gm, vxreg, vexdata, dlambda, skew, shape, dist, data,
zrf, N, res, garchenv = .garchenv)
filter = list()
filter$z = ans$z
filter$sigma = sqrt(abs(ans$h))
filter$residuals = ans$res
filter$LLH = -ans$llh
filter$log.likelihoods = ans$LHT
filter$data = data
filter$dates = dates
filter$origdata = origdata
filter$origdates = origdates
filter$date.format = dformat
filter$persistence = persist
filter$distribution = distribution
filter$pars = pars
model = c(mmodel, vmodel, dmodel)
sol = new("gjrGARCHfilter",
filter = filter,
model = model)
return(sol)
}
#---------------------------------------------------------------------------------
# SECTION gjrGARCH forecast
#---------------------------------------------------------------------------------
.gjrgarchforecast = function(fitORspec, data = NULL, n.ahead = 10, n.roll = 0, out.sample = 0,
external.forecasts = list(mregfor = NULL, vregfor = NULL), ...)
{
fit = fitORspec
data = fit@fit$data
N = length(as.numeric(data))
origData = fit@fit$origdata
origDates = fit@fit$origdates
dformat = fit@fit$date.format
Nor = length(as.numeric(origData))
ns = fit@model$n.start
if( n.roll > ns )
stop("n.roll must not be greater than out.sample!")
include.mean = fit@model$include.mean
armaOrder = fit@model$armaOrder
garchInMean = fit@model$garchInMean
im = inMeanType = fit@model$inMeanType
arfima = fit@model$arfima
include.mex = fit@model$include.mex
mxn = fit@model$mxn
if(mxn>0) mexdata = matrix(fit@model$mexdata, ncol=mxn) else mexdata = NULL
pow = fit@model$pow
garchOrder = fit@model$garchOrder
include.vex = fit@model$include.vex
vxn = fit@model$vxn
if(vxn>0) vexdata = matrix(fit@model$vexdata, ncol=vxn) else vexdata = NULL
modelnames = fit@model$modelnames
pars = fit@fit$coef
distribution = fit@model$distribution
include.dlambda = fit@model$include.dlambda
include.skew = fit@model$include.skew
include.shape = fit@model$include.shape
# check if necessary the external regressor forecasts provided first
xreg = .forcregressors(include.mex, mxn, fit@model$mexdata,
mregfor = external.forecasts$mregfor, include.vex, vxn,
fit@model$vexdata, external.forecasts$vregfor, pars, n.ahead, Nor,
out.sample = ns, n.roll)
mxreg = xreg$mxreg
mxf = xreg$mxf
vxreg = xreg$vxreg
vxf = xreg$vxf
pos = fit@model$pos
ar = armap = ma = armaq = inmean = darfima = 0
alpha = garchp = beta = garchq = gm = 0
skew = shape = dlambda = 0
if(pos[2,3] == 1){
ar = pars[pos[2,1]:pos[2,2]]
armap = length(ar)
}
if(pos[3,3] == 1){
ma = pars[pos[3,1]:pos[3,2]]
armaq = length(ma)
}
if(pos[4,3] == 1) inmean = pars[pos[4,1]:pos[4,2]]
if(pos[5,3] == 1) darfima = pars[pos[5,1]:pos[5,2]]
if(pos[8,3] == 1){
alpha = pars[pos[8,1]:pos[8,2]]
garchp=garchOrder[1]
}
if(pos[9,3] == 1) gm = pars[pos[9,1]:pos[9,2]]
if(pos[10,3] == 1){
beta = pars[pos[10,1]:pos[10,2]]
garchq=garchOrder[2]
}
if(pos[12,3]) dlambda = pars[pos[12,1]:pos[12,2]]
if(pos[13,3]) skew = pars[pos[13,1]:pos[13,2]]
if(pos[14,3]) shape = pars[pos[14,1]:pos[14,2]]
kappa = 0
if(garchp>0) kappa = pneg(dlambda, shape, skew, distribution)
fcreq = ifelse(ns >= (n.ahead+n.roll), n.ahead+n.roll, ns)
fspec = ugarchspec(variance.model = list(model = "gjrGARCH",
garchOrder = c(garchp, garchq), submodel = NULL,
external.regressors = vxf[1:(N + fcreq), , drop = FALSE]),
mean.model = list(armaOrder = c(armap, armaq),
include.mean = include.mean,
garchInMean = garchInMean, inMeanType = inMeanType, arfima = arfima,
external.regressors = mxf[1:(N + fcreq), , drop = FALSE]),
distribution.model = distribution, fixed.pars = as.list(pars))
tmp = data.frame(origData[1:(N + fcreq)])
rownames(tmp) = as.character(origDates[1:(N + fcreq)])
colnames(tmp) = "data"
flt = .gjrgarchfilter(data = tmp[,1,drop=FALSE], spec = fspec,
n.old = length(data))
sigmafilter = flt@filter$sigma
resfilter = flt@filter$residuals
zfilter = flt@filter$z
# forecast GARCH process
forecasts = vector(mode="list", length = n.roll+1)
fwdd = vector(mode="list", length = n.roll+1)
for(i in 1:(n.roll+1)){
np = N + i - 1
if(pos[1,3] == 1){
mu = rep(pars[pos[1,1]:pos[1,2]], N+i+n.ahead-1)
} else{
mu = rep(0, N+i+n.ahead-1)
}
if(pos[7,3] == 1){
omega = rep(pars[pos[7,1]:pos[7,2]], N+i+n.ahead-1)
} else {
omega = rep(0, N+i+n.ahead-1)
}
# no look-ahead
h = c(sigmafilter[1:(N+i-1)], rep(0, n.ahead))
epsx = c(resfilter[1:(N+i-1)], rep(0, n.ahead))
x = c(origData[1:(N+i-1)], rep(0, n.ahead))
z = c(zfilter[1:(N+i-1)], rep(0, n.ahead))
# forecast of externals is provided outside the system
mxfi = mxf[1:(N+i-1+n.ahead), , drop=FALSE]
vxfi = vxf[1:(N+i-1+n.ahead), , drop=FALSE]
ans = .ngjrsgarchforecast(N = np, n.ahead, include.vex, vxfi, vxreg, vxn, h,
omega, alpha, gm, beta, garchp, garchq, epsx, z, data = x,
armaOrder, mu, ar, ma, arfima, inmean, im, include.mex, mxreg, mxfi,
darfima = darfima, kappa)
fdf = data.frame(sigma = ans$h, series = ans$x)
fwdd[[i]] = .forcdates( origDates, n.ahead, N, i, ns, dformat )
rownames(fdf) = as.character(fwdd[[i]])
forecasts[[i]] = fdf
}
fcst = list()
fcst$n.ahead = n.ahead
fcst$N = N+ns
fcst$n.start = ns
fcst$n.roll = n.roll
fcst$sigma = fit@fit$sigma
fcst$series = fit@fit$data
fcst$forecasts = forecasts
fcst$fdates = fwdd
model = fit@model
model$dates = fit@fit$dates
model$origdata = fit@fit$origdata
model$origdates = fit@fit$origdates
ans = new("gjrGARCHforecast",
forecast = fcst,
model = model,
flt)
return(ans)
}
.ngjrsgarchforecast = function(N, n.ahead, include.vex, vxf, vxreg, vxn, h, omega,
alpha, gm, beta, garchp, garchq, epsx, z, data, armaOrder, mu, ar, ma,
arfima, inmean, im, include.mex, mxreg, mxf, darfima, kappa)
{
if(include.vex){
omega = omega + vxf%*%t(matrix(vxreg, ncol = vxn))
}
for(i in 1:n.ahead){
h[N+i] = omega[N+i] + sum(beta*h[N+i-(1:garchq)]^2)
for (j in 1:garchp){
if (i-j > 0){
s1 = (h[N + i - j]^2)
s2 = kappa*gm[j]* (h[N + i - j]^2)
} else{
s1 = epsx[N + i - j]^2
s2 = gm[j]*(epsx[N + i - j]^2)*as.integer(epsx[N + i - j]<0)
}
h[N+i] = h[N+i] + alpha[j] * s1 + s2
}
h[N+i] = sqrt(h[N+i])
}
if(arfima){
res = arfimaf(data[1:N], h, armaOrder, mu, ar, ma, darfima, inmean, im, epsx,
include.mex, mxreg, mxf, n.ahead)
} else{
res = armaf(data[1:N], h, armaOrder, mu, ar, ma, inmean, im,
epsx, include.mex, mxreg, mxf, n.ahead)
}
return(list(h = h[(N+1):(N+n.ahead)], x = res[(N+1):(N+n.ahead)]))
}
#---------------------------------------------------------------------------------
# 2nd dispatch method for forecast
.gjrgarchforecast2 = function(fitORspec, data = NULL, n.ahead = 10, n.roll = 0, out.sample = 0,
external.forecasts = list(mregfor = NULL, vregfor = NULL), ...)
{
spec = fitORspec
if(is.null(data))
stop("\nugarchforecast-->error: data must not be NULL when using a specification!")
# we then extract the data/coefs etc
xdata = .extractdata(data)
Nor = length(as.numeric(xdata$data))
data = xdata$data[1:(Nor - out.sample)]
N = length(as.numeric(data))
dates = xdata$pos[1:(Nor - out.sample)]
origData = xdata$data
origDates = xdata$pos
dformat = xdata$dformat
vmodel = spec@variance.model
mmodel = spec@mean.model
dmodel = spec@distribution.model
omodel = spec@optimization.model
ns = out.sample
if( n.roll > ns )
stop("\nugarchforecast-->error: n.roll must not be greater than out.sample!")
include.mean = mmodel$include.mean
armaOrder = mmodel$armaOrder
garchInMean = mmodel$garchInMean
im = inMeanType = mmodel$inMeanType
arfima = mmodel$arfima
include.mex = mmodel$include.mex
mxn = mmodel$mxn
if(mxn>0) mexdata = matrix(mmodel$mexdata, ncol = mxn) else mexdata = NULL
pow = 2
garchOrder = vmodel$garchOrder
include.vex = vmodel$include.vex
vxn = vmodel$vxn
if(vxn>0) vexdata = matrix(vmodel$vexdata, ncol = vxn) else vexdata = NULL
modelnames = omodel$modelnames
pars = unlist(omodel$fixed.pars)
distribution = dmodel$distribution
include.dlambda = dmodel$include.dlambda
include.skew = dmodel$include.skew
include.shape = dmodel$include.shape
# check if necessary the external regressor forecasts provided first
xreg = .forcregressors(include.mex, mxn, mexdata,
mregfor = external.forecasts$mregfor, include.vex, vxn,
vexdata, external.forecasts$vregfor, pars, n.ahead, Nor,
out.sample = ns, n.roll)
mxreg = xreg$mxreg
mxf = xreg$mxf
vxreg = xreg$vxreg
vxf = xreg$vxf
pos = omodel$pos
ar = armap = ma = armaq = inmean = darfima = 0
alpha = garchp = beta = garchq = gm = 0
skew = shape = dlambda = 0
if(pos[2,3] == 1){
ar = pars[pos[2,1]:pos[2,2]]
armap = length(ar)
}
if(pos[3,3] == 1){
ma = pars[pos[3,1]:pos[3,2]]
armaq = length(ma)
}
if(pos[4,3] == 1) inmean = pars[pos[4,1]:pos[4,2]]
if(pos[5,3] == 1) darfima = pars[pos[5,1]:pos[5,2]]
if(pos[8,3] == 1){
alpha = pars[pos[8,1]:pos[8,2]]
garchp=garchOrder[1]
}
if(pos[9,3] == 1) gm = pars[pos[9,1]:pos[9,2]]
if(pos[10,3] == 1){
beta = pars[pos[10,1]:pos[10,2]]
garchq=garchOrder[2]
}
if(pos[12,3]) dlambda = pars[pos[12,1]:pos[12,2]]
if(pos[13,3]) skew = pars[pos[13,1]:pos[13,2]]
if(pos[14,3]) shape = pars[pos[14,1]:pos[14,2]]
kappa = 0
if(garchp>0) kappa = pneg(dlambda, shape, skew, distribution)
fcreq = ifelse(ns >= (n.ahead+n.roll), n.ahead+n.roll, ns)
fspec = ugarchspec(variance.model = list(model = "gjrGARCH",
garchOrder = c(garchp, garchq), submodel = NULL,
external.regressors = vxf[1:(N + fcreq), , drop = FALSE]),
mean.model = list(armaOrder = c(armap, armaq),
include.mean = include.mean,
garchInMean = garchInMean, inMeanType = inMeanType, arfima = arfima,
external.regressors = mxf[1:(N + fcreq), , drop = FALSE]),
distribution.model = distribution, fixed.pars = as.list(pars))
tmp = data.frame(origData[1:(N + fcreq)])
rownames(tmp) = as.character(origDates[1:(N + fcreq)])
flt = .gjrgarchfilter(data = tmp[,1,drop=FALSE], spec = fspec,
n.old = N)
sigmafilter = flt@filter$sigma
resfilter = flt@filter$residuals
zfilter = flt@filter$z
# forecast GARCH process
forecasts = vector(mode="list", length = n.roll+1)
fwdd = vector(mode="list", length = n.roll+1)
for(i in 1:(n.roll+1)){
np = N + i - 1
if(pos[1,3] == 1){
mu = rep(pars[pos[1,1]:pos[1,2]], N+i+n.ahead-1)
} else{
mu = rep(0, N+i+n.ahead-1)
}
if(pos[7,3] == 1){
omega = rep(pars[pos[7,1]:pos[7,2]], N+i+n.ahead-1)
} else {
omega = rep(0, N+i+n.ahead-1)
}
# no look-ahead
h = c(sigmafilter[1:(N+i-1)], rep(0, n.ahead))
epsx = c(resfilter[1:(N+i-1)], rep(0, n.ahead))
x = c(origData[1:(N+i-1)], rep(0, n.ahead))
z = c(zfilter[1:(N+i-1)], rep(0, n.ahead))
# forecast of externals is provided outside the system
mxfi = mxf[1:(N+i-1+n.ahead), , drop=FALSE]
vxfi = vxf[1:(N+i-1+n.ahead), , drop=FALSE]
ans = .ngjrsgarchforecast(N = np, n.ahead, include.vex, vxfi, vxreg, vxn, h,
omega, alpha, gm, beta, garchp, garchq, epsx, z, data = x,
armaOrder, mu, ar, ma, arfima, inmean, im, include.mex, mxreg, mxfi,
darfima = darfima, kappa)
fdf = data.frame(sigma = ans$h, series = ans$x)
fwdd[[i]] = .forcdates( origDates, n.ahead, N, i, ns, dformat )
rownames(fdf) = as.character(fwdd[[i]])
forecasts[[i]] = fdf
}
fcst = list()
fcst$n.ahead = n.ahead
fcst$N = N+ns
fcst$n.start = ns
fcst$n.roll = n.roll
fcst$sigma = flt@filter$sigma[1:N]
fcst$series = data
fcst$forecasts = forecasts
fcst$fdates = fwdd
model = c(mmodel, vmodel, omodel)
model$dates = dates
model$origdata = origData
model$origdates = origDates
ans = new("gjrGARCHforecast",
forecast = fcst,
model = model,
flt)
return(ans)
}
#---------------------------------------------------------------------------------
# SECTION gjrGARCH simulate
#---------------------------------------------------------------------------------
# choose the fastest method to use:
.gjrgarchsim = function(fit, n.sim = 1000, n.start = 0, m.sim = 1, startMethod =
c("unconditional","sample"), presigma = NA, prereturns = NA,
preresiduals = NA, rseed = NA, custom.dist = list(name = NA, distfit = NA),
mexsimdata = NULL, vexsimdata = NULL)
{
if( (n.sim+n.start) < 100 && m.sim > 100 ){
ans = .gjrgarchsim2(fit = fit, n.sim = n.sim, n.start = n.start, m.sim = m.sim,
startMethod = startMethod, presigma = presigma, prereturns = prereturns,
preresiduals = preresiduals, rseed = rseed, custom.dist = custom.dist,
mexsimdata = mexsimdata, vexsimdata = vexsimdata)
} else{
ans = .gjrgarchsim1(fit = fit, n.sim = n.sim, n.start = n.start, m.sim = m.sim,
startMethod = startMethod, presigma = presigma, prereturns = prereturns,
preresiduals = preresiduals, rseed = rseed, custom.dist = custom.dist,
mexsimdata = mexsimdata, vexsimdata = vexsimdata)
}
return( ans )
}
.gjrgarchsim1 = function(fit, n.sim = 1000, n.start = 0, m.sim = 1, startMethod =
c("unconditional","sample"), presigma = NA, prereturns = NA,
preresiduals = NA, rseed = NA, custom.dist = list(name = NA, distfit = NA),
mexsimdata = NULL, vexsimdata = NULL)
{
# some checks
if(is.na(rseed[1])){
sseed = as.integer(runif(1,0,as.integer(Sys.time())))
} else{
if(length(rseed) != m.sim) sseed = as.integer(rseed[1]) else sseed = rseed[1:m.sim]
}
# Enlarge Series:
n = n.sim + n.start
startMethod = startMethod[1]
data = fit@fit$data
N = length(as.numeric(data))
m = fit@model$maxOrder2
resids = fit@fit$residuals
sigma = fit@fit$sigma
include.mean = fit@model$include.mean
armaOrder = fit@model$armaOrder
garchInMean = fit@model$garchInMean
inMeanType = fit@model$inMeanType
im = inMeanType/2
arfima = fit@model$arfima
include.mex = fit@model$include.mex
mxn = fit@model$mxn
if(mxn>0) mexdata = matrix(fit@model$mexdata, ncol=mxn) else mexdata = NULL
pow = fit@model$pow
garchOrder = fit@model$garchOrder
include.vex = fit@model$include.vex
vxn = fit@model$vxn
if(vxn>0) vexdata = matrix(fit@model$vexdata, ncol=vxn) else vexdata = NULL
modelnames = fit@model$modelnames
pars = fit@fit$coef
distribution = fit@model$distribution
include.dlambda = fit@model$include.dlambda
include.skew = fit@model$include.skew
include.shape = fit@model$include.shape
# check if necessary the external regressor forecasts provided first
xreg = .simregressors(include.mex, mxn, mexsimdata, mexdata, include.vex,
vxn, vexsimdata, vexdata, N, pars, n, m.sim, m)
if(mxn>0) mxreg = matrix(xreg$mxreg, ncol = mxn)
mexsim = xreg$mexsimlist
vxreg = xreg$vxreg
vexsim = xreg$vexsimlist
pos = fit@model$pos
ar = armap = ma = armaq = inmean = darfima = mu = 0
alpha = garchp = beta = garchq = gm = 0
skew = shape = dlambda = 0
omega = 0.00001
if(pos[1,3] == 1) mu = pars[pos[1,1]:pos[1,2]]
if(pos[2,3] == 1){
ar = pars[pos[2,1]:pos[2,2]]
armap = length(ar)
}
if(pos[3,3] == 1){
ma = pars[pos[3,1]:pos[3,2]]
armaq = length(ma)
}
if(pos[4,3] == 1) inmean = pars[pos[4,1]:pos[4,2]]
if(pos[5,3] == 1) darfima = pars[pos[5,1]:pos[5,2]]
if(pos[7,3] == 1) omega = pars[pos[7,1]:pos[7,2]]
if(pos[8,3] == 1){
alpha = pars[pos[8,1]:pos[8,2]]
garchp=garchOrder[1]
}
if(pos[9,3] == 1) gm = pars[pos[9,1]:pos[9,2]]
if(pos[10,3] == 1){
beta = pars[pos[10,1]:pos[10,2]]
garchq=garchOrder[2]
}
if(pos[12,3]) dlambda = pars[pos[12,1]:pos[12,2]]
if(pos[13,3]) skew = pars[pos[13,1]:pos[13,2]]
if(pos[14,3]) shape = pars[pos[14,1]:pos[14,2]]
if(N < n.start){
startmethod[1] = "unconditional"
cat("\nugarchsim-->warning: n.start greater than length of data...using unconditional
start method...\n")
}
# Random Samples from the Distribution
#dst = match(distribution, c("norm","snorm","std","sstd","ged","sged","nig","ghyp","jsu"))
if(length(sseed) == 1){
zmatrix = data.frame(dist = distribution, dlambda = dlambda,
skew = skew, shape = shape, n = n * m.sim, seed = sseed[1])
z = .custzdist(custom.dist, zmatrix, m.sim, n)
} else{
zmatrix = data.frame(dist = rep(distribution, m.sim), dlambda = rep(dlambda, m.sim),
skew = rep(skew, m.sim), shape = rep(shape, m.sim), n = rep(n, m.sim), seed = sseed)
z = .custzdist(custom.dist, zmatrix, m.sim, n)
}
if(startMethod == "unconditional"){
z = rbind(matrix(0, nrow = m, ncol = m.sim), z)
} else{
z = rbind(matrix(tail(fit@fit$z, m), nrow = m, ncol = m.sim), z)
}
# create the presample information
if(!is.na(presigma[1])){
presigma = as.vector(presigma)
if(length(presigma)<m) stop(paste("\nugarchsim-->error: presigma must be of length ", m, sep=""))
}
if(!is.na(prereturns[1])){
prereturns = as.vector(prereturns)
if(length(prereturns)<m) stop(paste("\nugarchsim-->error: prereturns must be of length ", m, sep=""))
}
if(!is.na(preresiduals[1])){
preresiduals = as.vector(preresiduals)
if(length(preresiduals)<m) stop(paste("\nugarchsim-->error: preresiduals must be of length ", m, sep=""))
preres = matrix(preresiduals, nrow = m)
}
if(is.na(presigma[1])){
if(startMethod[1] == "unconditional"){
persist = persistence(fit)
hEst = sqrt(omega/abs(1-persist))
presigma = as.numeric(rep(hEst, times = m))}
else{
presigma = tail(sigma, m)
}
}
if(is.na(prereturns[1])){
if(startMethod[1] == "unconditional"){
prereturns = as.numeric(rep(uncmean(fit), times = m))
}
else{
prereturns = tail(data, m)
}
}
# input vectors/matrices
h = c(presigma * presigma, rep(0, n))
x = c(prereturns, rep(0, n))
constm = matrix(mu, ncol = m.sim, nrow = n + m)
# outpus matrices
sigmaSim = matrix(0, ncol = m.sim, nrow = n.sim)
seriesSim = matrix(0, ncol = m.sim, nrow = n.sim)
residSim = matrix(0, ncol = m.sim, nrow = n.sim)
for(i in 1:m.sim){
if(is.na(preresiduals[1])){
if(startMethod[1] == "unconditional"){
preres = as.numeric(z[1:m, i])*presigma
} else{
preres = tail(resids, m)
}
}
ngrd = which(preres<0)
tmpr = rep(0, length(preres))
tmpr[ngrd] = preres[ngrd] * preres[ngrd]
nres = c(tmpr, rep(0, n))
res = c(preres, rep(0, n))
ans1 = .gjrgarchsimC(garchp, garchq, incvex = vxn, omega, alpha,
beta, gm, vxreg, h, as.numeric(z[,i]), res, nres, vexsim[[i]], T = n + m, m)
sigmaSim[,i] = sqrt( ans1$h[(n.start + m + 1):(n+m)] )
residSim[,i] = ans1$res[(n.start + m + 1):(n+m)]
if(include.mex){
constm[,i] = constm[,i] + mxreg%*%t(matrix(mexsim[[i]], ncol = mxn))
}
if(garchInMean) constm[,i] = constm[,i] + inmean*(ans1$h^im)
if(arfima){
## if(constant) constm[,i] = constm[,i]*(1-sum(ar))
ans2 = .arfimaxsim(armap, armaq, constm[1:n, i], ar, ma, darfima, x, ans1$res[1:(n+armaq)], T = n)
seriesSim[,i] = ans2$s[(n.start + armaq + 1):(n + armaq)]
} else{
# if(constant) constm[,i] = constm[,i]*(1-sum(ar))
ans2 = .armaxsim(armap, armaq, constm[,i], ar, ma, x, ans1$res, T = n + m, m)
seriesSim[,i] = ans2$x[(n.start + m + 1):(n+m)]
}
}
smodel = fit@model
smodel$dates = fit@fit$dates
sim = list(sigmaSim = sigmaSim, seriesSim = seriesSim, residSim = residSim)
sim$series = data
sim$sigma = sigma
sol = new("gjrGARCHsim",
simulation = sim,
model = smodel,
seed = as.integer(sseed))
return(sol)
}
#---------------------------------------------------------------------------------
# SECTION gjrGARCH path
#---------------------------------------------------------------------------------
.gjrgarchpath = function(spec, n.sim = 1000, n.start = 0, m.sim = 1,
presigma = NA, prereturns = NA, preresiduals = NA, rseed = NA,
custom.dist = list(name = NA, distfit = NA), mexsimdata = NULL,
vexsimdata = NULL)
{
if( (n.sim+n.start) < 20 && m.sim > 100 ){
ans = .gjrgarchpath2(spec = spec, n.sim = n.sim, n.start = n.start, m.sim = m.sim,
presigma = presigma, prereturns = prereturns, preresiduals = preresiduals,
rseed = rseed, custom.dist = custom.dist, mexsimdata = mexsimdata,
vexsimdata = vexsimdata)
} else{
ans = .gjrgarchpath1(spec = spec, n.sim = n.sim, n.start = n.start, m.sim = m.sim,
presigma = presigma, prereturns = prereturns, preresiduals = preresiduals,
rseed = rseed, custom.dist = custom.dist, mexsimdata = mexsimdata,
vexsimdata = vexsimdata)
}
return( ans )
}
.gjrgarchpath1 = function(spec, n.sim = 1000, n.start = 0, m.sim = 1,
presigma = NA, prereturns = NA, preresiduals = NA, rseed = NA,
custom.dist = list(name = NA, distfit = NA), mexsimdata = NULL,
vexsimdata = NULL)
{
# some checks
if(is.na(rseed[1])){
sseed = as.integer(runif(1,0,as.integer(Sys.time())))
} else{
if(length(rseed) != m.sim) sseed = as.integer(rseed[1]) else sseed = rseed[1:m.sim]
}
# Enlarge Series:
n = n.sim + n.start
armaOrder = spec@mean.model$armaOrder
armap = armaOrder[1]
armaq = armaOrder[2]
include.mean=spec@mean.model$include.mean
garchOrder = spec@variance.model$garchOrder
garchp = garchOrder[1]
garchq = garchOrder[2]
include.mex = spec@mean.model$include.mex
mxn = spec@mean.model$mxn
N = 0
if(include.mex) {
mexdata = matrix(spec@mean.model$mexdata, ncol=mxn)
N = dim(mexdata)[1]
} else { mexdata = NULL }
include.vex = spec@variance.model$include.vex
vxn = spec@variance.model$vxn
if(include.vex) {
vexdata = matrix(spec@variance.model$vexdata, ncol=vxn)
N = dim(vexdata)[1]
} else { vexdata = NULL }
garchInMean = spec@mean.model$garchInMean
inMeanType = spec@mean.model$inMeanType
im = inMeanType/2
distribution = spec@distribution.model$distribution
pos = spec@optimization.model$pos.matrix
pars = unlist(spec@optimization.model$fixed.pars)
# perform checks on supplied parameters
exp.pars = rownames(spec@optimization.model$modelmatrix)
spec.pars = names(pars)
tmatch = match(sort(spec.pars), sort(exp.pars))
if(any(is.na(tmatch))){
cat("\nugarchpath-->error: supplied parameters do not match expected parameters\n")
cat(paste("epxected parameters:", sort(exp.pars), "\n", sep=""))
cat(paste("supplied parameters:", sort(spec.pars), "\n", sep=""))
stop("\n")
}
omega = 0.00001
mu = ar = ma = inmean = darfima = mxreg = alpha = beta = gm = vxreg = skew = shape = dlambda = 0
if(pos[1,3]) mu = pars[pos[1,1]:pos[1,2]]
if(pos[2,3]) ar = pars[pos[2,1]:pos[2,2]]
if(pos[3,3]) ma = pars[pos[3,1]:pos[3,2]]
if(pos[4,3]) inmean = pars[pos[4,1]:pos[4,2]]
if(pos[5,3]) darfima = pars[pos[5,1]:pos[5,2]]
if(pos[6,3]) mxreg = pars[pos[6,1]:pos[6,2]]
if(pos[7,3]) omega = pars[pos[7,1]:pos[7,2]]
if(pos[8,3]) alpha = pars[pos[8,1]:pos[8,2]]
if(pos[9,3]) gm = pars[pos[9,1]:pos[9,2]]
if(pos[10,3]) beta = pars[pos[10,1]:pos[10,2]]
if(pos[11,3]) vxreg = pars[pos[11,1]:pos[11,2]]
if(pos[12,3]) dlambda = pars[pos[12,1]:pos[12,2]]
if(pos[13,3]) skew = pars[pos[13,1]:pos[13,2]]
if(pos[14,3]) shape = pars[pos[14,1]:pos[14,2]]
arfima = spec@mean.model$arfima
m = spec@optimization.model$maxOrder2
modelnames = spec@optimization.model$modelnames
include.skew = spec@distribution.model$include.skew
include.shape = spec@distribution.model$include.shape
# check if necessary the external regressor forecasts provided first
xreg = .simregressors(include.mex, mxn, mexsimdata, mexdata, include.vex,
vxn, vexsimdata, vexdata, N, pars, n, m.sim, m)
if(mxn>0) mxreg = matrix(xreg$mxreg, ncol = mxn)
mexsim = xreg$mexsimlist
vxreg = xreg$vxreg
vexsim = xreg$vexsimlist
persist = persistence(pars = pars, distribution = distribution, model = "gjrGARCH")
if(persist >= 1)
cat(paste("\nugarchpath->warning: persitence :", round(persist,5), sep=""))
# Random Samples from the Distribution
#dst = match(distribution, c("norm","snorm","std","sstd","ged","sged","nig","ghyp","jsu"))
if(length(sseed) == 1){
zmatrix = data.frame(dist = distribution, dlambda = dlambda,
skew = skew, shape = shape, n = n * m.sim, seed = sseed[1])
z = .custzdist(custom.dist, zmatrix, m.sim, n)
} else{
zmatrix = data.frame(dist = rep(distribution, m.sim), dlambda = rep(dlambda, m.sim),
skew = rep(skew, m.sim), shape = rep(shape, m.sim), n = rep(n, m.sim), seed = sseed)
z = .custzdist(custom.dist, zmatrix, m.sim, n)
}
z = rbind(matrix(0, nrow = m, ncol = m.sim), z)
# create the presample information
if(!is.na(presigma[1])){
presigma = as.vector(presigma)
if(length(presigma)<m) stop(paste("\nugarchsim-->error: presigma must be of length ", m, sep=""))
}
if(!is.na(prereturns[1])){
prereturns = as.vector(prereturns)
if(length(prereturns)<m) stop(paste("\nugarchsim-->error: prereturns must be of length ", m, sep=""))
}
if(!is.na(preresiduals[1])){
preresiduals = as.vector(preresiduals)
if(length(preresiduals)<m) stop(paste("\nugarchsim-->error: preresiduals must be of length ", m, sep=""))
preres = matrix(preresiduals, nrow = m)
}
if(is.na(presigma[1])){
hEst = (omega/abs(1-persist))^(1/2)
presigma = as.numeric(rep(hEst, times = m))
}
if(is.na(prereturns[1])){
prereturns = as.numeric(rep(uncmean(spec), times = m))
}
# input vectors/matrices
h = c(presigma*presigma, rep(0, n))
x = c(prereturns, rep(0, n))
constm = matrix(mu, ncol = m.sim, nrow = n + m)
# outpus matrices
sigmaSim = matrix(0, ncol = m.sim, nrow = n.sim)
seriesSim = matrix(0, ncol = m.sim, nrow = n.sim)
residSim = matrix(0, ncol = m.sim, nrow = n.sim)
for(i in 1:m.sim){
if(is.na(preresiduals[1])){
preres = as.numeric(z[1:m, i])*presigma
}
res = c(preres, rep(0, n))
ngrd = which(preres<0)
tmpr = rep(0, length(preres))
tmpr[ngrd] = preres[ngrd]^2
nres = c(tmpr,rep(0, n))
ans1 = .gjrgarchsimC(garchp, garchq, incvex = vxn, omega, alpha,
beta, gm, vxreg, h, as.numeric(z[,i]), res, nres, vexsim[[i]], T = n + m, m)
sigmaSim[,i] = ans1$h[(n.start + m + 1):(n+m)]^(1/2)
residSim[,i] = ans1$res[(n.start + m + 1):(n+m)]
if(include.mex){
constm[,i] = constm[,i] + mxreg%*%t(matrix(mexsim[[i]], ncol = mxn))
}
if(garchInMean) constm[,i] = constm[,i] + inmean*(ans1$h^im)
if(arfima){
## if(constant) constm[,i] = constm[,i]*(1-sum(ar))
ans2 = .arfimaxsim(armap, armaq, constm[1:n, i], ar, ma, darfima, x, ans1$res[1:(n+armaq)], T = n)
seriesSim[,i] = ans2$s[(n.start + armaq + 1):(n + armaq)]
} else{
# if(constant) constm[,i] = constm[,i]*(1-sum(ar))
ans2 = .armaxsim(armap, armaq, constm[,i], ar, ma, x, ans1$res, T = n + m, m)
seriesSim[,i] = ans2$x[(n.start + m + 1):(n+m)]
}
}
path = list(sigmaSim = sigmaSim, seriesSim = seriesSim, residSim = residSim)
sol = new("gjrGARCHpath",
path = path,
seed = as.integer(sseed))
return(sol)
}
### Timing Tests suggest this is faster when n.sim is small and m.sim is large ###
.gjrgarchsim2 = function(fit, n.sim = 1000, n.start = 0, m.sim = 1, startMethod =
c("unconditional","sample"), presigma = NA, prereturns = NA,
preresiduals = NA, rseed = NA, custom.dist = list(name = NA, distfit = NA),
mexsimdata = NULL, vexsimdata = NULL)
{
# some checks
if(is.na(rseed[1])){
sseed = as.integer(runif(1,0,as.integer(Sys.time())))
} else{
if(length(rseed) != m.sim) sseed = as.integer(rseed[1]) else sseed = as.integer(rseed[1:m.sim])
}
# Enlarge Series:
n = n.sim + n.start
startMethod = startMethod[1]
data = fit@fit$data
N = length(as.numeric(data))
m = fit@model$maxOrder2
resids = fit@fit$residuals
sigma = fit@fit$sigma
include.mean = fit@model$include.mean
armaOrder = fit@model$armaOrder
garchInMean = fit@model$garchInMean
inMeanType = fit@model$inMeanType
im = inMeanType/2
arfima = fit@model$arfima
include.mex = fit@model$include.mex
mxn = fit@model$mxn
if(mxn>0) mexdata = matrix(fit@model$mexdata, ncol=mxn) else mexdata = NULL
pow = fit@model$pow
garchOrder = fit@model$garchOrder
include.vex = fit@model$include.vex
vxn = fit@model$vxn
if(vxn>0) vexdata = matrix(fit@model$vexdata, ncol=vxn) else vexdata = NULL
modelnames = fit@model$modelnames
pars = fit@fit$coef
distribution = fit@model$distribution
include.dlambda = fit@model$include.dlambda
include.skew = fit@model$include.skew
include.shape = fit@model$include.shape
# check if necessary the external regressor forecasts provided first
xreg = .simregressors(include.mex, mxn, mexsimdata, mexdata, include.vex,
vxn, vexsimdata, vexdata, N, pars, n, m.sim, m)
if(mxn>0) mxreg = matrix(xreg$mxreg, ncol = mxn)
mexsim = xreg$mexsimlist
vxreg = xreg$vxreg
vexsim = xreg$vexsimlist
pos = fit@model$pos
ar = armap = ma = armaq = inmean = darfima = mu = 0
alpha = garchp = beta = garchq = gm = 0
skew = shape = dlambda = 0
omega = 0.00001
if(pos[1,3] == 1) mu = pars[pos[1,1]:pos[1,2]]
if(pos[2,3] == 1){
ar = pars[pos[2,1]:pos[2,2]]
armap = length(ar)
}
if(pos[3,3] == 1){
ma = pars[pos[3,1]:pos[3,2]]
armaq = length(ma)
}
if(pos[4,3] == 1) inmean = pars[pos[4,1]:pos[4,2]]
if(pos[5,3] == 1) darfima = pars[pos[5,1]:pos[5,2]]
if(pos[7,3] == 1) omega = pars[pos[7,1]:pos[7,2]]
if(pos[8,3] == 1){
alpha = pars[pos[8,1]:pos[8,2]]
garchp=garchOrder[1]
}
if(pos[9,3] == 1) gm = pars[pos[9,1]:pos[9,2]]
if(pos[10,3] == 1){
beta = pars[pos[10,1]:pos[10,2]]
garchq=garchOrder[2]
}
if(pos[12,3]) dlambda = pars[pos[12,1]:pos[12,2]]
if(pos[13,3]) skew = pars[pos[13,1]:pos[13,2]]
if(pos[14,3]) shape = pars[pos[14,1]:pos[14,2]]
if(N < n.start){
startmethod[1] = "unconditional"
cat("\nugarchsim-->warning: n.start greater than length of data...using unconditional start method...\n")
}
# Random Samples from the Distribution # CHANGE for speed up
#dst = match(distribution, c("norm","snorm","std","sstd","ged","sged","nig","ghyp","jsu"))
if(length(sseed) == 1){
zmatrix = data.frame(dist = distribution, dlambda = dlambda,
skew = skew, shape = shape, n = n * m.sim, seed = sseed[1])
z = .custzdist(custom.dist, zmatrix, m.sim, n)
} else{
zmatrix = data.frame(dist = rep(distribution, m.sim), dlambda = rep(dlambda, m.sim),
skew = rep(skew, m.sim), shape = rep(shape, m.sim), n = rep(n, m.sim), seed = sseed)
z = .custzdist(custom.dist, zmatrix, m.sim, n)
}
if(startMethod == "unconditional"){
z = rbind(matrix(0, nrow = m, ncol = m.sim), z)
} else{
z = rbind(matrix(tail(fit@fit$z, m), nrow = m, ncol = m.sim), z)
}
# create the presample information
if(!is.na(presigma[1])){
presigma = as.vector(presigma)
if(length(presigma)<m) stop(paste("\nugarchsim-->error: presigma must be of length ", m, sep=""))
}
if(!is.na(prereturns[1])){
prereturns = as.vector(prereturns)
if(length(prereturns)<m) stop(paste("\nugarchsim-->error: prereturns must be of length ", m, sep=""))
}
if(!is.na(preresiduals[1])){
preresiduals = as.vector(preresiduals)
if(length(preresiduals)<m) stop(paste("\nugarchsim-->error: preresiduals must be of length ", m, sep=""))
preres = matrix(preresiduals, nrow = m)
}
if(is.na(presigma[1])){
if(startMethod[1] == "unconditional"){
persist = persistence(fit)
hEst = sqrt(omega/abs(1-persist))
presigma = as.numeric(rep(hEst, times = m))}
else{
presigma = tail(sigma, m)
}
}
if(is.na(prereturns[1])){
if(startMethod[1] == "unconditional"){
prereturns = as.numeric(rep(uncmean(fit), times = m))
}
else{
prereturns = tail(data, m)
}
}
# input vectors/matrices
h = matrix(c(presigma * presigma, rep(0, n)), nrow = n + m, ncol = m.sim)
x = matrix(c(prereturns, rep(0, n)), nrow = n + m, ncol = m.sim)
constm = matrix(mu, nrow = n + m, ncol = m.sim)
# outpus matrices
sigmaSim = matrix(0, ncol = m.sim, nrow = n.sim)
seriesSim = matrix(0, ncol = m.sim, nrow = n.sim)
residSim = matrix(0, ncol = m.sim, nrow = n.sim)
if(is.na(preresiduals[1])){
if(startMethod[1] == "unconditional"){
preres = matrix( z[1:m, 1:m.sim] * presigma, nrow = m, ncol = m.sim )
} else{
preres = matrix(tail(resids, m), nrow = m, ncol = m.sim)
}
}
ngrd = which(preres<0, arr.ind = TRUE)
tmpr = matrix(0, nrow = dim(preres)[1], ncol = m.sim)
if(dim(ngrd)[1]>0) tmpr[ngrd] = preres[ngrd]*preres[ngrd]
nres = rbind(tmpr, matrix(0, nrow = n, ncol = m.sim))
res = rbind(preres, matrix(0, nrow = n, ncol = m.sim))
nindx = matrix(1, nrow = n + m, ncol = m.sim)
nindx[which(z>0, arr.ind = TRUE)] = 0.0
garchindx = as.numeric( c(garchp, garchq, m) )
# we'll do the external regressors first for speed.
if(include.vex){
vxs = sapply(vexsim, FUN = function(x) vxreg%*%t(matrix(x, ncol = vxn)))
} else{
vxs = matrix(0, nrow = m + n, ncol = m.sim)
}
e = res * res
ans = .Call("mgjrgarchsim", garchindx = garchindx, omega = as.numeric(omega), alpha = as.numeric(alpha),
beta = as.numeric(beta), gm = as.numeric(gm), h = h, z = z, res = res, e = e, nres = nres,
nindx = nindx, vxs = vxs, PACKAGE = "rgarch")
# we have to re-check the matrix argument since when n.sim = 1 (i.e. 1 row), it returns a vector!
# (THIS IS SO ANNOYING!!!!)
sigmaSim = matrix(sqrt( ans$h[(n.start + m + 1):(n+m), ] ), ncol = m.sim)
residSim = matrix(ans$res[(n.start + m + 1):(n+m), ], ncol = m.sim)
if(include.mex){
mxs = sapply(mexsim, FUN = function(x) mxreg%*%t(matrix(x, ncol = mxn)))
} else{
mxs = 0
}
if(garchInMean){
imh = inmean*(sigmaSim^im)
} else{
imh = 0
}
constm = constm + mxs + imh
armaindx = c(armap, armaq, m, n)
if(arfima){
## if(constant) constm[,i] = constm[,i]*(1-sum(ar))
for(i in 1:m.sim){
tmp = .arfimaxsim(armap, armaq, constm[1:n, i], ar, ma, darfima, x, ans$res[1:(n+armaq), i], T = n)
seriesSim[,i] = tmp$s[(n.start + armaq + 1):(n + armaq)]
}
} else{
# if(constant) constm = constm * ( 1 - sum(ar) )
tmp = .Call("marmaxsim", armaindx = armaindx, ar = ar, ma = ma, mu = constm, x = x, res = ans$res,
PACKAGE = "rgarch")
seriesSim = matrix(tmp$x[(n.start + m + 1):(n+m), ], ncol = m.sim)
}
smodel = fit@model
smodel$dates = fit@fit$dates
sim = list(sigmaSim = sigmaSim, seriesSim = seriesSim, residSim = residSim)
sim$series = data
sim$sigma = sigma
sol = new("gjrGARCHsim",
simulation = sim,
model = smodel,
seed = as.integer(sseed))
return(sol)
}
.gjrgarchpath2 = function(spec, n.sim = 1000, n.start = 0, m.sim = 1,
presigma = NA, prereturns = NA, preresiduals = NA, rseed = NA,
custom.dist = list(name = NA, distfit = NA), mexsimdata = NULL,
vexsimdata = NULL)
{
# some checks
if(is.na(rseed[1])){
sseed = as.integer(runif(1,0,as.integer(Sys.time())))
} else{
if(length(rseed) != m.sim) sseed = as.integer(rseed[1]) else sseed = rseed[1:m.sim]
}
# Enlarge Series:
n = n.sim + n.start
armaOrder = spec@mean.model$armaOrder
armap = armaOrder[1]
armaq = armaOrder[2]
include.mean=spec@mean.model$include.mean
garchOrder = spec@variance.model$garchOrder
garchp = garchOrder[1]
garchq = garchOrder[2]
include.mex = spec@mean.model$include.mex
mxn = spec@mean.model$mxn
N = 0
if(include.mex) {
mexdata = matrix(spec@mean.model$mexdata, ncol=mxn)
N = dim(mexdata)[1]
} else { mexdata = NULL }
include.vex = spec@variance.model$include.vex
vxn = spec@variance.model$vxn
if(include.vex) {
vexdata = matrix(spec@variance.model$vexdata, ncol=vxn)
N = dim(vexdata)[1]
} else { vexdata = NULL }
garchInMean = spec@mean.model$garchInMean
inMeanType = spec@mean.model$inMeanType
im = inMeanType/2
distribution = spec@distribution.model$distribution
pos = spec@optimization.model$pos.matrix
pars = unlist(spec@optimization.model$fixed.pars)
# perform checks on supplied parameters
exp.pars = rownames(spec@optimization.model$modelmatrix)
spec.pars = names(pars)
tmatch = match(sort(spec.pars), sort(exp.pars))
if(any(is.na(tmatch))){
cat("\nugarchpath-->error: supplied parameters do not match expected parameters\n")
cat(paste("epxected parameters:", sort(exp.pars), "\n", sep=""))
cat(paste("supplied parameters:", sort(spec.pars), "\n", sep=""))
stop("\n")
}
omega = 0.00001
mu = ar = ma = inmean = darfima = mxreg = alpha = beta = gm = vxreg = skew = shape = dlambda = 0
if(pos[1,3]) mu = pars[pos[1,1]:pos[1,2]]
if(pos[2,3]) ar = pars[pos[2,1]:pos[2,2]]
if(pos[3,3]) ma = pars[pos[3,1]:pos[3,2]]
if(pos[4,3]) inmean = pars[pos[4,1]:pos[4,2]]
if(pos[5,3]) darfima = pars[pos[5,1]:pos[5,2]]
if(pos[6,3]) mxreg = pars[pos[6,1]:pos[6,2]]
if(pos[7,3]) omega = pars[pos[7,1]:pos[7,2]]
if(pos[8,3]) alpha = pars[pos[8,1]:pos[8,2]]
if(pos[9,3]) gm = pars[pos[9,1]:pos[9,2]]
if(pos[10,3]) beta = pars[pos[10,1]:pos[10,2]]
if(pos[11,3]) vxreg = pars[pos[11,1]:pos[11,2]]
if(pos[12,3]) dlambda = pars[pos[12,1]:pos[12,2]]
if(pos[13,3]) skew = pars[pos[13,1]:pos[13,2]]
if(pos[14,3]) shape = pars[pos[14,1]:pos[14,2]]
arfima = spec@mean.model$arfima
m = spec@optimization.model$maxOrder2
modelnames = spec@optimization.model$modelnames
include.skew = spec@distribution.model$include.skew
include.shape = spec@distribution.model$include.shape
# check if necessary the external regressor forecasts provided first
xreg = .simregressors(include.mex, mxn, mexsimdata, mexdata, include.vex,
vxn, vexsimdata, vexdata, N, pars, n, m.sim, m)
if(mxn>0) mxreg = matrix(xreg$mxreg, ncol = mxn)
mexsim = xreg$mexsimlist
vxreg = xreg$vxreg
vexsim = xreg$vexsimlist
persist = persistence(pars = pars, distribution = distribution, model = "gjrGARCH")
if(persist >= 1)
cat(paste("\nugarchpath->warning: persitence :", round(persist,5), sep=""))
# Random Samples from the Distribution
#dst = match(distribution, c("norm","snorm","std","sstd","ged","sged","nig","ghyp","jsu"))
if(length(sseed) == 1){
zmatrix = data.frame(dist = distribution, dlambda = dlambda,
skew = skew, shape = shape, n = n * m.sim, seed = sseed[1])
z = .custzdist(custom.dist, zmatrix, m.sim, n)
} else{
zmatrix = data.frame(dist = rep(distribution, m.sim), dlambda = rep(dlambda, m.sim),
skew = rep(skew, m.sim), shape = rep(shape, m.sim), n = rep(n, m.sim), seed = sseed)
z = .custzdist(custom.dist, zmatrix, m.sim, n)
}
z = rbind(matrix(0, nrow = m, ncol = m.sim), z)
# create the presample information
if(!is.na(presigma[1])){
presigma = as.vector(presigma)
if(length(presigma)<m) stop(paste("\nugarchsim-->error: presigma must be of length ", m, sep=""))
}
if(!is.na(prereturns[1])){
prereturns = as.vector(prereturns)
if(length(prereturns)<m) stop(paste("\nugarchsim-->error: prereturns must be of length ", m, sep=""))
}
if(!is.na(preresiduals[1])){
preresiduals = as.vector(preresiduals)
if(length(preresiduals)<m) stop(paste("\nugarchsim-->error: preresiduals must be of length ", m, sep=""))
preres = matrix(preresiduals, nrow = m)
}
if(is.na(presigma[1])){
hEst = (omega/abs(1-persist))^(1/2)
presigma = as.numeric(rep(hEst, times = m))
}
if(is.na(prereturns[1])){
prereturns = as.numeric(rep(uncmean(spec), times = m))
}
# input vectors/matrices
h = matrix(c(presigma * presigma, rep(0, n)), nrow = n + m, ncol = m.sim)
x = matrix(c(prereturns, rep(0, n)), nrow = n + m, ncol = m.sim)
constm = matrix(mu, nrow = n + m, ncol = m.sim)
# outpus matrices
sigmaSim = matrix(0, ncol = m.sim, nrow = n.sim)
seriesSim = matrix(0, ncol = m.sim, nrow = n.sim)
residSim = matrix(0, ncol = m.sim, nrow = n.sim)
if(is.na(preresiduals[1])){
preres = matrix( z[1:m, 1:m.sim] * presigma, nrow = m, ncol = m.sim )
} else{
preres = matrix( preresiduals, nrow = m, ncol = m.sim )
}
ngrd = which(preres<0, arr.ind = TRUE)
tmpr = matrix(0, nrow = dim(preres)[1], ncol = m.sim)
if(dim(ngrd)[1]>0) tmpr[ngrd] = preres[ngrd]*preres[ngrd]
nres = rbind(tmpr, matrix(0, nrow = n, ncol = m.sim))
res = rbind(preres, matrix(0, nrow = n, ncol = m.sim))
nindx = matrix(1, nrow = n + m, ncol = m.sim)
nindx[which(z>0, arr.ind = TRUE)] = 0.0
garchindx = as.numeric( c(garchp, garchq, m, n) )
# we'll do the external regressors first for speed.
if(include.vex){
vxs = sapply(vexsim, FUN = function(x) vxreg%*%t(matrix(x, ncol = vxn)))
} else{
vxs = matrix(0, nrow = m + n, ncol = m.sim)
}
e = res * res
ans = .Call("mgjrgarchsim", garchindx = garchindx, omega = as.numeric(omega), alpha = as.numeric(alpha),
beta = as.numeric(beta), gm = as.numeric(gm), h = h, z = z, res = res, e = e, nres = nres,
nindx = nindx, vxs = vxs, PACKAGE = "rgarch", DUP = FALSE)
sigmaSim = matrix(sqrt( ans$h[(n.start + m + 1):(n+m), ] ), ncol = m.sim)
residSim = matrix(ans$res[(n.start + m + 1):(n+m), ], ncol = m.sim)
if(include.mex){
mxs = sapply(mexsim, FUN = function(x) mxreg%*%t(matrix(x, ncol = mxn)))
} else{
mxs = 0
}
if(garchInMean){
imh = inmean*(sigmaSim^im)
} else{
imh = 0
}
constm = constm + mxs + imh
armaindx = c(armap, armaq, m, n)
if(arfima){
## if(constant) constm[,i] = constm[,i]*(1-sum(ar))
for(i in 1:m.sim){
tmp = .arfimaxsim(armap, armaq, constm[1:n, i], ar, ma, darfima, x, ans$res[1:(n+armaq), i], T = n)
seriesSim[,i] = tmp$s[(n.start + armaq + 1):(n + armaq)]
}
} else{
# if(constant) constm = constm * ( 1 - sum(ar) )
tmp = .Call("marmaxsim", armaindx = armaindx, ar = ar, ma = ma, mu = constm, x = x, res = ans$res)
seriesSim = matrix(tmp$x[(n.start + m + 1):(n+m), ], ncol = m.sim)
}
path = list(sigmaSim = sigmaSim, seriesSim = seriesSim, residSim = residSim)
sol = new("gjrGARCHpath",
path = path,
seed = as.integer(sseed))
return(sol)
}
Try the rgarch package in your browser
Any scripts or data that you put into this service are public.
rgarch documentation built on May 2, 2019, 5:22 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
#################################################################################
##
## R package rgarch by Alexios Ghalanos Copyright (C) 2008, 2009, 2010, 2011
## This file is part of the R package rgarch.
##
## The R package rgarch is free software: you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation, either version 3 of the License, or
## (at your option) any later version.
##
## The R package rgarch is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
## GNU General Public License for more details.
##
#################################################################################
#---------------------------------------------------------------------------------
# SECTION gjrGARCH spec
#---------------------------------------------------------------------------------
.gjrgarchspec = function(vmodel, mmodel, dmodel, start.pars, fixed.pars)
{
# extract values:
# mean model
include.mean = mmodel$include.mean
armap = mmodel$armaOrder[1]
armaq = mmodel$armaOrder[2]
garchInMean = mmodel$garchInMean
inmeanType = mmodel$inMeanType
# inmeanType(1) means sd
# inmeanType(2) means var
#if(garchInMean && (inmeanType!=1 | inmeanType!=2)) stop("valid inmeanType values are 1 and 2")
im = inmeanType/2
arfima = mmodel$arfima
if(!is.null(mmodel$external.regressors)){
include.mex = TRUE
mexdata = as.matrix(mmodel$external.regressors)
mxn = dim(mexdata)[2]
} else{
include.mex = FALSE
mexdata = NULL
mxn = 0
}
# vmodel
garchp = vmodel$garchOrder[1]
garchq = vmodel$garchOrder[2]
if(!is.null(vmodel$external.regressors)){
include.vex = TRUE
vexdata = as.matrix(vmodel$external.regressors)
vxn = dim(vexdata)[2]
} else{
include.vex = FALSE
vexdata = NULL
vxn = 0
}
include.omega = ifelse(vmodel$targeting, FALSE, TRUE)
if(!include.omega) fixed.pars$omega = 0
# dmodel/skew/shape
include.dlambda = dmodel$include.dlambda
include.skew = dmodel$include.skew
include.shape = dmodel$include.shape
model = vmodel$model
submodel = vmodel$submodel
maxOrder = max(c(garchp,garchq))
maxOrder2 = max(c(garchp, garchq, armap, armaq))
vmodel = list(model = model, submodel = submodel, garchOrder=vmodel$garchOrder,
include.vex = include.vex, vexdata = vexdata, vxn = vxn, pow=2, include.omega = include.omega)
mmodel = list(include.mean = include.mean, armaOrder = mmodel$armaOrder,
garchInMean = garchInMean, inMeanType = inmeanType, im = im, arfima = arfima,
include.mex = include.mex, mexdata = mexdata, mxn = mxn)
#specnames = c("mu","ar","ma","darfima","mxreg","omega","alpha","beta","vxreg")
# match
modelnames = c(
if(include.mean) "mu",
if(armap > 0) paste("ar", 1:armap, sep = ""),
if(armaq > 0) paste("ma", 1:armaq, sep = ""),
if(garchInMean) paste("inmean"),
if(arfima) paste("darfima"),
if(include.mex) paste("mxreg", 1:mxn, sep = ""),
paste("omega"),
if(garchp > 0) paste("alpha", 1:garchp, sep = ""),
if(garchp > 0) paste("gamma", 1:garchp, sep = ""),
if(garchq > 0) paste("beta", 1:garchq, sep = ""),
if(include.vex) paste("vxreg", 1:vxn, sep = ""),
if(include.dlambda) "dlambda",
if(include.skew) "skew",
if(include.shape) "shape")
pos=1
pos.matrix = matrix(0, ncol = 3, nrow = 14)
colnames(pos.matrix) = c("start","stop","include")
rownames(pos.matrix) = c("mu", "ar", "ma", "inmean", "darfima", "mxreg", "omega",
"alpha", "gamma", "beta", "vxreg", "dlambda", "skew", "shape")
if(include.mean){ pos.matrix[1,1:3] = 1
pos=pos.matrix[1,2]+1}
if(armap>0){pos.matrix[2,1:3] = c(pos,pos+armap-1,1)
pos=max(pos.matrix[1:2,2])+1}
if(armaq>0){pos.matrix[3,1:3] = c(pos,pos+armaq-1,1)
pos=max(pos.matrix[1:3,2])+1}
if(garchInMean){pos.matrix[4,1:3] = c(pos,pos,1)
pos=max(pos.matrix[1:4,2])+1}
if(arfima){pos.matrix[5,1:3] = c(pos,pos,1)
pos=max(pos.matrix[1:5,2])+1}
if(include.mex){pos.matrix[6,1:3] = c(pos,pos+mxn-1,1)
pos=max(pos.matrix[1:6,2])+1}
pos.matrix[7,1:3] = c(pos,pos,1)
pos=max(pos.matrix[1:7,2])+1
if(garchp > 0){pos.matrix[8,1:3] = c(pos,pos+garchp-1,1)
pos=max(pos.matrix[1:8,2])+1}
if(garchp > 0){pos.matrix[9,1:3] = c(pos,pos+garchp-1,1)
pos=max(pos.matrix[1:9,2])+1}
if(garchq > 0){pos.matrix[10,1:3] = c(pos,pos+garchq-1,1)
pos=max(pos.matrix[1:10,2])+1}
if(include.vex){pos.matrix[11,1:3] = c(pos,pos+vxn-1,1)
pos=max(pos.matrix[1:11,2])+1}
if(include.dlambda){pos.matrix[12,1:3] = c(pos,pos,1)
pos=max(pos.matrix[1:12,2])+1}
if(include.skew){pos.matrix[13,1:3] = c(pos,pos,1)
pos=max(pos.matrix[1:13,2])+1}
if(include.shape){pos.matrix[14,1:3] = c(pos,pos,1)}
nn = length(modelnames)
modelmatrix = matrix(0, ncol = 3, nrow = nn)
rownames(modelmatrix) = modelnames
colnames(modelmatrix) = c("opt","fixed","start")
fixed.names = names(fixed.pars)
fp=charmatch(fixed.names,modelnames)
if(!is.null(fixed.names) && any(!is.na(fp))){
fixed = fp[!is.na(fp)]
modelmatrix[fixed,2]=1
fz=charmatch(modelnames,fixed.names)
fz=fz[!is.na(fz)]
fixed.pars=fixed.pars[fz]
names(fixed.pars) = fixed.names[fz]
} else{
fixed.pars=NULL
}
modelmatrix[,1] = 1
start.names = names(start.pars)
sp=charmatch(start.names,modelnames)
if(!is.null(start.names) && any(!is.na(sp))){
start = sp[!is.na(sp)]
modelmatrix[start,3]=1
sz=charmatch(modelnames,start.names)
sz=sz[!is.na(sz)]
start.pars=start.pars[sz]
} else{
start.pars=NULL
}
# the way we handle fixed pars will be through upper and lower bound constraints
opt.pars = vector(mode="numeric",length=length(modelnames))
names(opt.pars) = modelnames
# optimization model
omodel = list(modelnames = modelnames, opt.pars = opt.pars, start.pars = start.pars,
fixed.pars = fixed.pars, modelmatrix = modelmatrix, maxOrder = maxOrder,
maxOrder2 = maxOrder2, pos.matrix = pos.matrix)
ans = new("gjrGARCHspec",
mean.model = mmodel,
variance.model = vmodel,
distribution.model = dmodel,
optimization.model = omodel)
return(ans)
}
#---------------------------------------------------------------------------------
# SECTION gjrGARCH fit
#---------------------------------------------------------------------------------
.gjrgarchfit = function(spec, data, out.sample = 0, solver = "solnp", solver.control = list(),
fit.control = list(stationarity = 1, fixed.se = 0, scale = 0))
{
tic = Sys.time()
if(is.null(solver.control$trace)) trace = 0 else trace = solver.control$trace
if(is.null(fit.control$stationarity)) fit.control$stationarity = 1
if(is.null(fit.control$fixed.se)) fit.control$fixed.se = 0
if(is.null(fit.control$scale)) fit.control$scale = 0
# if we have external regressors in variance turn off scaling
if(spec@variance.model$include.vex) fit.control$scale = 0
# at the moment we do not allow scaling with fixed parameters
if(!is.null(spec@optimization.model$fixed.pars)) fit.control$scale = 0
xdata = .extractdata(data)
if(!is.numeric(out.sample)) stop("\nugarchfit-->error: out.sample must be numeric\n")
if(as.numeric(out.sample)<0) stop("\nugarchfit-->error: out.sample must be positive\n")
n.start = round(out.sample,0)
n = length(xdata$data)
if((n-n.start)<100) stop("\nugarchfit-->error: function requires at least 100 data\n points to run\n")
data = xdata$data[1:(n-n.start)]
dates = xdata$pos[1:(n-n.start)]
origdata = xdata$data
origdates = xdata$pos
dformat = xdata$dformat
# get current environment
tempenvir = new.env(hash = TRUE, parent = .GlobalEnv)
# expand the spec object and assign spec lists
vmodel = spec@variance.model
mmodel = spec@mean.model
dmodel = spec@distribution.model
omodel = spec@optimization.model
if(mmodel$include.mex) {
mmodel$origmexdata = mmodel$mexdata
mmodel$mexdata = mmodel$mexdata[1:(n-n.start),]
} else{
mmodel$origmexdata = NULL
}
if(vmodel$include.vex) {
vmodel$origvexdata = vmodel$vexdata
vmodel$vexdata = vmodel$vexdata[1:(n-n.start),]
} else{
vmodel$origvexdata = NULL
}
assign("vmodel", vmodel, envir = tempenvir)
assign("mmodel", mmodel, envir = tempenvir)
assign("dmodel", dmodel, envir = tempenvir)
assign("omodel", omodel, envir = tempenvir)
assign("dates", dates, envir = tempenvir)
assign("fit.control", fit.control, envir = tempenvir)
assign("solver", solver, envir = tempenvir)
assign("trace", trace, envir = tempenvir)
include.mean = mmodel$include.mean
include.omega = vmodel$include.omega
m = omodel$maxOrder
T = length(as.numeric(data))
dist = dmodel$distribution
model = vmodel$model
submodel = vmodel$submodel
if(fit.control$scale) dscale = sd(data) else dscale = 1
assign("dscale", dscale, envir = tempenvir)
zdata = data/dscale
# Optimization Starting Parameters Vector & Bounds
ipars = .garchstart(data = zdata, model = model, submodel = submodel, garchenv = tempenvir)
pars = ipars$pars
LB = ipars$modelLB
UB = ipars$modelUB
assign("npar", length(pars), envir = tempenvir)
# we now split out any fixed parameters
if(any(LB==UB)){
fixid = which(LB==UB)
if(length(fixid)==length(pars)){
if(fit.control$fixed.se==0) {
# if all parameters are fixed an no standard erros are to
# be calculated then we return a ugarchfilter object
cat("\nugarchfit-->warning: all parameters fixed...returning ugarchfilter
object instead\n")
return(ugarchfilter(data = data, spec = spec))
} else{
# if all parameters are fixed but we require standard errors, we
# skip the solver
use.solver = 0
fixed = pars[fixid]
pars = pars[-fixid]
}
} else{
# with some parameters fixed we extract them (to be rejoined at end)
# so that they do not enter the solver
use.solver = 1
fixed = pars[fixid]
pars = pars[-fixid]
LB = LB[-fixid]
UB = UB[-fixid]
}
} else{
# no fixed parameters, all normal
fixed = NULL
fixid = NULL
use.solver=1
}
assign("fixed", fixed, envir = tempenvir)
assign("fixid", fixid, envir = tempenvir)
# start counter
assign(".llh", 1, envir = tempenvir)
# scale parameters (for non-solnp solvers)
parscale = rep(1,length=length(pars))
names(parscale) = names(pars)
if(include.mean) parscale["mu"] = abs(mean(zdata))
if(include.omega) parscale["omega"] = var(zdata)
assign("omega.calculate", TRUE, envir = tempenvir)
# assisgn solver constraints (solnp directly else exterior type penalty
# for other solvers)
if(fit.control$stationarity == 1 && vmodel$include.vex == 0){
cb = .garchconbounds()
Ifn = .gjrgarchcon
ILB = cb$LB
IUB = cb$UB
if(solver == "solnp" | solver == "gosolnp") fit.control$stationarity = 0
} else{
Ifn = NULL
ILB = NULL
IUB = NULL
}
# conditions controls the non-solnp solver penalty
assign("fit.control", fit.control, , envir = tempenvir)
if(use.solver){
solution = .garchsolver(solver, pars, fun = .gjrgarchLLH, Ifn, ILB, IUB,
gr = NULL, hessian = NULL, parscale = parscale,
control = solver.control, LB = LB, UB = UB, ux = NULL, ci = NULL,
mu = NULL, data = zdata, returnType = "llh", garchenv = tempenvir)
sol = solution$sol
hess = solution$hess
timer = Sys.time()-tic
opt.pars = sol$par
if(!include.omega){
fixed = get("fixed", tempenvir)
fixed["omega"] = get("omega", tempenvir)
fixed["omega"] = (fixed["omega"]) * dscale^2
assign("fixed", fixed, envir = tempenvir)
}
assign("omega.calculate", FALSE, envir = tempenvir)
if(!is.null(opt.pars)) names(opt.pars) = names(pars)
if(is.null(fixed)){
if(include.mean) opt.pars["mu"] = opt.pars["mu"]*dscale
if(mmodel$include.mex){
opt.pars[omodel$pos[6,1]:omodel$pos[6,2]] = opt.pars[omodel$pos[6,1]:omodel$pos[6,2]] * dscale
}
opt.pars["omega"] = (opt.pars["omega"]) * dscale^2
}
convergence = sol$convergence
} else{
hess = NULL
timer = Sys.time()-tic
opt.pars = NULL
convergence = 0
sol = list()
sol$message = "all parameters fixed"
}
fit = list()
model = list()
assign("pow", 2, envir = tempenvir)
# check convergence else write message/return
if(convergence==0){
if(!is.null(fixed) && fit.control$fixed.se == 1){
fixed = get("fixed", tempenvir)
if(!include.omega){
vmodel$include.omega = TRUE
assign("vmodel", vmodel, envir = tempenvir)
}
npar = get("npar", tempenvir)
pall = vector(mode = "numeric", length = npar)
pall[fixid] = fixed
pall[-fixid] = opt.pars
opt.pars = pall
Names = omodel$modelnames
names(opt.pars) = Names
assign("fixed", NULL, envir = tempenvir)
assign("fixid", NULL, envir = tempenvir)
}
fit = .makefitmodel(garchmodel = "gjrGARCH", f = .gjrgarchLLH, opt = opt.pars,
data = data, T = T, m = m, timer = timer, convergence = convergence,
message = sol$message, hess, garchenv = tempenvir)
fit$dates = dates
fit$date.format = dformat
fit$origdata = origdata
fit$origdates = origdates
} else{
fit$message = sol$message
fit$convergence = 1
}
# make model list to return some usefule information which
# will be called by other functions (show, plot, sim etc)
model = .makemodel(tempenvir)
model$n.start = n.start
# return object gjrGARCHfit
ans = new("gjrGARCHfit",
fit = fit,
model = model)
rm(tempenvir)
return(ans)
}
#---------------------------------------------------------------------------------
# SECTION gjrGARCH LLH
#---------------------------------------------------------------------------------
.gjrgarchLLH = function(pars, data, returnType = "llh", garchenv)
{
# prepare inputs
# rejoin fixed and pars
fixed = get("fixed", garchenv)
fixid = get("fixid", garchenv)
npar = get("npar", garchenv)
trace = get("trace", garchenv)
if(!is.null(fixed)){
pall = vector(mode = "numeric", length = npar)
pall[fixid] = fixed
pall[-fixid] = pars
pars = pall
}
T = length(data)
vmodel = get("vmodel", garchenv)
mmodel = get("mmodel", garchenv)
dmodel = get("dmodel", garchenv)
omodel = get("omodel", garchenv)
fit.control = get("fit.control", garchenv)
armaOrder = c(mmodel$armaOrder,
mmodel$im,
mmodel$arfima,
mmodel$include.mex)
m = omodel$maxOrder
N = c(m,T)
garchOrder = c(vmodel$garchOrder, vmodel$include.vex)
include.omega = vmodel$include.omega
omega.calculate = get("omega.calculate", garchenv)
include.inmean = mmodel$garchInMean
mexdata = mmodel$mexdata
vexdata = vmodel$vexdata
distribution = dmodel$distribution
# distribution number
# 1 = norm, 2=snorm, 3=std, 4=sstd, 5=ged, 6=sged, 7=nig
dist = dmodel$distno
pos = omodel$pos.matrix
omega = 0.00001
mu = ar = ma = inmean = darfima = mxreg = alpha = beta = gm = vxreg = skew = shape = dlambda = 0
if(pos[1,3]) mu = pars[pos[1,1]:pos[1,2]]
if(pos[2,3]) ar = pars[pos[2,1]:pos[2,2]]
if(pos[3,3]) ma = pars[pos[3,1]:pos[3,2]]
if(pos[4,3]) inmean = pars[pos[4,1]:pos[4,2]]
if(pos[5,3]) darfima = pars[pos[5,1]:pos[5,2]]
if(pos[6,3]) mxreg = pars[pos[6,1]:pos[6,2]]
if(pos[8,3]) alpha = pars[pos[8,1]:pos[8,2]]
if(pos[9,3]) gm = pars[pos[9,1]:pos[9,2]]
if(pos[10,3]) beta = pars[pos[10,1]:pos[10,2]]
if(pos[11,3]) vxreg = pars[pos[11,1]:pos[11,2]]
if(pos[12,3]) dlambda = pars[pos[12,1]:pos[12,2]]
if(pos[13,3]) skew = pars[pos[13,1]:pos[13,2]]
if(pos[14,3]) shape = pars[pos[14,1]:pos[14,2]]
#rx = .getresiduals2(data = data, include.mean = mmodel$include.mean,
# mu = mu, armaOrder=armaOrder, ar = ar, ma = ma, include.inmean,
# inmean = inmean, darfima = darfima, mexdata = mexdata, mxreg = mxreg,
# h = get("meqspars", garchenv)$h)
hm = get("meqspars", garchenv)$h
rx = .arfimaxfilter(armaOrder = armaOrder, mu = mu, ar = ar, ma = ma, inmean = inmean,
darfima = darfima, mxreg = mxreg, mexdata = mexdata, h = hm,
data = data, N = N, garchenv)
res = rx$res
zrf = rx$zrf
res = res[!is.na(res) && is.finite(res) && !is.nan(res)]
#names(pars) = omodel$modelnames
names(pars) = omodel$modelnames
# gjrgarch persistence value
persist = sum(beta)+ sum(alpha)+sum(apply(as.data.frame(gm),1,FUN=function(x)
x*pneg(dlambda, shape, skew, distribution)))
# unconditional sigma value
mvar = mean(res*res)
if(pos[11,3] == 1) {
mv = sum(apply(matrix(vexdata, ncol = vmodel$vxn), 2, "mean")*vxreg)
} else{
mv = 0
}
if(include.omega){
omega = max(eps, pars[pos[7,1]:pos[7,2]])
hEst = mvar
} else{
if(omega.calculate){
persist = .persistgjrgarch(pars = pars, distribution = distribution)
hEst = mvar
omega = hEst * (1 - persist) - mv
assign("omega", omega, envir = garchenv)
} else{
omega = get("omega", envir = garchenv)
hEst = mvar
}
}
# if we have external regressors in variance equation we cannot have
# stationarity checks in likelihood
if(fit.control$stationarity == 1 && vmodel$include.vex == 0){
if(!is.na(persist) && persist >= 1) return(llh = get(".llh", garchenv) +
0.1*(abs(get(".llh", garchenv))))
}
# C Code [gjrgarch.c] via wrapper
ans = .gjrgarchfitC(armaOrder, garchOrder, mu, ar, ma, inmean, mxreg, mexdata,
omega, hEst, alpha, beta, gm, vxreg, vexdata, dlambda, skew, shape, dist, data,
zrf, N, res, garchenv)
if(get(".csol", garchenv) == 1){
if( trace > 0 ) cat(paste("\nugarchfit-->warning: ", get(".filtermessage",garchenv),"\n", sep=""))
return(llh = (get(".llh", garchenv) + 0.1*(abs(get(".llh", garchenv)))))
}
z = ans$z
h = ans$h
epsx = ans$res
llh = ans$llh
if(is.finite(llh) && !is.na(llh) && !is.nan(llh)){
assign(".llh", llh, envir = garchenv)
} else {
llh = (get(".llh", garchenv) + 0.1*(abs(get(".llh", garchenv))))
}
# LHT = raw scores
LHT = -ans$LHT
ans = switch(returnType,
llh = llh,
LHT = LHT,
all = list(llh = llh, h = h, epsx = epsx, z = z, kappa = kappa, LHT = LHT, persistence = persist))
return(ans)
}
#---------------------------------------------------------------------------------
# SECTION gjrGARCH filter
#---------------------------------------------------------------------------------
.gjrgarchfilter = function(spec, data, out.sample = 0, n.old = NULL)
{
.garchenv = environment()
if(!is.null(n.old)) Nx = n.old else Nx = length(data)
xdata = .extractdata(data)
data = xdata$data
dates = xdata$pos
dformat = xdata$dformat
origdata = data
origdates = dates
T = length(origdata) - out.sample
data = origdata[1:T]
dates = origdates[1:T]
vmodel = spec@variance.model
mmodel = spec@mean.model
dmodel = spec@distribution.model
omodel = spec@optimization.model
include.mean = mmodel$include.mean
modelnames = omodel$modelnames
pars = unlist(spec@optimization.model$fixed.pars)
parnames = names(pars)
if(length(modelnames) != length(pars))
stop("parameters do not match expected specification length")
if(is.na(all(match(modelnames, parnames), 1:length(modelnames)))) {
stop("parameters names do not match specification")
} else{
pars = pars[modelnames]
}
m = omodel$maxOrder
T = length(as.numeric(data))
dist = dmodel$distribution
model = vmodel$model
submodel = vmodel$submodel
armaOrder = c(mmodel$armaOrder,
mmodel$im,
mmodel$arfima,
mmodel$include.mex)
garchOrder = c(vmodel$garchOrder,
vmodel$include.vex)
include.omega = vmodel$include.omega
include.inmean = mmodel$garchInMean
mexdata = mmodel$mexdata
vexdata = vmodel$vexdata
distribution = dmodel$distribution
garchOrder = c(vmodel$garchOrder, vmodel$include.vex)
include.inmean = mmodel$garchInMean
mexdata = mmodel$mexdata
vexdata = vmodel$vexdata
distribution = dmodel$distribution
# distribution number
# 1 = norm, 2=snorm, 3=std, 4=sstd, 5=ged, 6=sged, 7=nig
dist = dmodel$distno
N = c(m,T)
pos = omodel$pos.matrix
omega = 0.00001
mu = ar = ma = inmean = darfima = mxreg = alpha = beta = gm = vxreg = skew = shape = dlambda = 0
if(pos[1,3]) mu = pars[pos[1,1]:pos[1,2]]
if(pos[2,3]) ar = pars[pos[2,1]:pos[2,2]]
if(pos[3,3]) ma = pars[pos[3,1]:pos[3,2]]
if(pos[4,3]) inmean = pars[pos[4,1]:pos[4,2]]
if(pos[5,3]) darfima = pars[pos[5,1]:pos[5,2]]
if(pos[6,3]) mxreg = pars[pos[6,1]:pos[6,2]]
if(pos[8,3]) alpha = pars[pos[8,1]:pos[8,2]]
if(pos[9,3]) gm = pars[pos[9,1]:pos[9,2]]
if(pos[10,3]) beta = pars[pos[10,1]:pos[10,2]]
if(pos[11,3]) vxreg = pars[pos[11,1]:pos[11,2]]
if(pos[12,3]) dlambda = pars[pos[12,1]:pos[12,2]]
if(pos[13,3]) skew = pars[pos[13,1]:pos[13,2]]
if(pos[14,3]) shape = pars[pos[14,1]:pos[14,2]]
rx = .arfimaxfilter(armaOrder = armaOrder, mu = mu, ar = ar, ma = ma, inmean = inmean,
darfima = darfima, mxreg = mxreg, mexdata = mexdata, h = 0,
data = data, N = N, garchenv = .garchenv)
res = rx$res
zrf = rx$zrf
if(!is.null(n.old)){
rx2 = .arfimaxfilter(armaOrder = armaOrder, mu = mu, ar = ar, ma = ma, inmean = inmean,
darfima = darfima, mxreg = mxreg, mexdata = mexdata[1:Nx, , drop = FALSE], h = 0,
data = origdata[1:Nx], N = c(m, Nx), garchenv = .garchenv)
res2 = rx2$res
# unconditional sigma value
mvar = mean(res2*res2)
} else{
mvar = mean(res*res)
}
# gjrgarch persistence value
persist = .persistgjrgarch(pars, distribution)
#if(persist>1) stop("\nPersistence of model is > 1...try again!\n")
#hEst = omega + persistence*mvar
if(pos[11,3] == 1) {
mv = sum(apply(matrix(vexdata, ncol = vmodel$vxn), 2, "mean")*vxreg)
} else{
mv = 0
}
if(include.omega){
omega = pars[pos[7,1]:pos[7,2]]
hEst = mvar
} else{
omega = hEst * (1 - persist) - mv
hEst = mvar
}
ans = .gjrgarchfitC(armaOrder, garchOrder, mu, ar, ma, inmean, mxreg, mexdata,
omega, hEst, alpha, beta, gm, vxreg, vexdata, dlambda, skew, shape, dist, data,
zrf, N, res, garchenv = .garchenv)
filter = list()
filter$z = ans$z
filter$sigma = sqrt(abs(ans$h))
filter$residuals = ans$res
filter$LLH = -ans$llh
filter$log.likelihoods = ans$LHT
filter$data = data
filter$dates = dates
filter$origdata = origdata
filter$origdates = origdates
filter$date.format = dformat
filter$persistence = persist
filter$distribution = distribution
filter$pars = pars
model = c(mmodel, vmodel, dmodel)
sol = new("gjrGARCHfilter",
filter = filter,
model = model)
return(sol)
}
#---------------------------------------------------------------------------------
# SECTION gjrGARCH forecast
#---------------------------------------------------------------------------------
.gjrgarchforecast = function(fitORspec, data = NULL, n.ahead = 10, n.roll = 0, out.sample = 0,
external.forecasts = list(mregfor = NULL, vregfor = NULL), ...)
{
fit = fitORspec
data = fit@fit$data
N = length(as.numeric(data))
origData = fit@fit$origdata
origDates = fit@fit$origdates
dformat = fit@fit$date.format
Nor = length(as.numeric(origData))
ns = fit@model$n.start
if( n.roll > ns )
stop("n.roll must not be greater than out.sample!")
include.mean = fit@model$include.mean
armaOrder = fit@model$armaOrder
garchInMean = fit@model$garchInMean
im = inMeanType = fit@model$inMeanType
arfima = fit@model$arfima
include.mex = fit@model$include.mex
mxn = fit@model$mxn
if(mxn>0) mexdata = matrix(fit@model$mexdata, ncol=mxn) else mexdata = NULL
pow = fit@model$pow
garchOrder = fit@model$garchOrder
include.vex = fit@model$include.vex
vxn = fit@model$vxn
if(vxn>0) vexdata = matrix(fit@model$vexdata, ncol=vxn) else vexdata = NULL
modelnames = fit@model$modelnames
pars = fit@fit$coef
distribution = fit@model$distribution
include.dlambda = fit@model$include.dlambda
include.skew = fit@model$include.skew
include.shape = fit@model$include.shape
# check if necessary the external regressor forecasts provided first
xreg = .forcregressors(include.mex, mxn, fit@model$mexdata,
mregfor = external.forecasts$mregfor, include.vex, vxn,
fit@model$vexdata, external.forecasts$vregfor, pars, n.ahead, Nor,
out.sample = ns, n.roll)
mxreg = xreg$mxreg
mxf = xreg$mxf
vxreg = xreg$vxreg
vxf = xreg$vxf
pos = fit@model$pos
ar = armap = ma = armaq = inmean = darfima = 0
alpha = garchp = beta = garchq = gm = 0
skew = shape = dlambda = 0
if(pos[2,3] == 1){
ar = pars[pos[2,1]:pos[2,2]]
armap = length(ar)
}
if(pos[3,3] == 1){
ma = pars[pos[3,1]:pos[3,2]]
armaq = length(ma)
}
if(pos[4,3] == 1) inmean = pars[pos[4,1]:pos[4,2]]
if(pos[5,3] == 1) darfima = pars[pos[5,1]:pos[5,2]]
if(pos[8,3] == 1){
alpha = pars[pos[8,1]:pos[8,2]]
garchp=garchOrder[1]
}
if(pos[9,3] == 1) gm = pars[pos[9,1]:pos[9,2]]
if(pos[10,3] == 1){
beta = pars[pos[10,1]:pos[10,2]]
garchq=garchOrder[2]
}
if(pos[12,3]) dlambda = pars[pos[12,1]:pos[12,2]]
if(pos[13,3]) skew = pars[pos[13,1]:pos[13,2]]
if(pos[14,3]) shape = pars[pos[14,1]:pos[14,2]]
kappa = 0
if(garchp>0) kappa = pneg(dlambda, shape, skew, distribution)
fcreq = ifelse(ns >= (n.ahead+n.roll), n.ahead+n.roll, ns)
fspec = ugarchspec(variance.model = list(model = "gjrGARCH",
garchOrder = c(garchp, garchq), submodel = NULL,
external.regressors = vxf[1:(N + fcreq), , drop = FALSE]),
mean.model = list(armaOrder = c(armap, armaq),
include.mean = include.mean,
garchInMean = garchInMean, inMeanType = inMeanType, arfima = arfima,
external.regressors = mxf[1:(N + fcreq), , drop = FALSE]),
distribution.model = distribution, fixed.pars = as.list(pars))
tmp = data.frame(origData[1:(N + fcreq)])
rownames(tmp) = as.character(origDates[1:(N + fcreq)])
colnames(tmp) = "data"
flt = .gjrgarchfilter(data = tmp[,1,drop=FALSE], spec = fspec,
n.old = length(data))
sigmafilter = flt@filter$sigma
resfilter = flt@filter$residuals
zfilter = flt@filter$z
# forecast GARCH process
forecasts = vector(mode="list", length = n.roll+1)
fwdd = vector(mode="list", length = n.roll+1)
for(i in 1:(n.roll+1)){
np = N + i - 1
if(pos[1,3] == 1){
mu = rep(pars[pos[1,1]:pos[1,2]], N+i+n.ahead-1)
} else{
mu = rep(0, N+i+n.ahead-1)
}
if(pos[7,3] == 1){
omega = rep(pars[pos[7,1]:pos[7,2]], N+i+n.ahead-1)
} else {
omega = rep(0, N+i+n.ahead-1)
}
# no look-ahead
h = c(sigmafilter[1:(N+i-1)], rep(0, n.ahead))
epsx = c(resfilter[1:(N+i-1)], rep(0, n.ahead))
x = c(origData[1:(N+i-1)], rep(0, n.ahead))
z = c(zfilter[1:(N+i-1)], rep(0, n.ahead))
# forecast of externals is provided outside the system
mxfi = mxf[1:(N+i-1+n.ahead), , drop=FALSE]
vxfi = vxf[1:(N+i-1+n.ahead), , drop=FALSE]
ans = .ngjrsgarchforecast(N = np, n.ahead, include.vex, vxfi, vxreg, vxn, h,
omega, alpha, gm, beta, garchp, garchq, epsx, z, data = x,
armaOrder, mu, ar, ma, arfima, inmean, im, include.mex, mxreg, mxfi,
darfima = darfima, kappa)
fdf = data.frame(sigma = ans$h, series = ans$x)
fwdd[[i]] = .forcdates( origDates, n.ahead, N, i, ns, dformat )
rownames(fdf) = as.character(fwdd[[i]])
forecasts[[i]] = fdf
}
fcst = list()
fcst$n.ahead = n.ahead
fcst$N = N+ns
fcst$n.start = ns
fcst$n.roll = n.roll
fcst$sigma = fit@fit$sigma
fcst$series = fit@fit$data
fcst$forecasts = forecasts
fcst$fdates = fwdd
model = fit@model
model$dates = fit@fit$dates
model$origdata = fit@fit$origdata
model$origdates = fit@fit$origdates
ans = new("gjrGARCHforecast",
forecast = fcst,
model = model,
flt)
return(ans)
}
.ngjrsgarchforecast = function(N, n.ahead, include.vex, vxf, vxreg, vxn, h, omega,
alpha, gm, beta, garchp, garchq, epsx, z, data, armaOrder, mu, ar, ma,
arfima, inmean, im, include.mex, mxreg, mxf, darfima, kappa)
{
if(include.vex){
omega = omega + vxf%*%t(matrix(vxreg, ncol = vxn))
}
for(i in 1:n.ahead){
h[N+i] = omega[N+i] + sum(beta*h[N+i-(1:garchq)]^2)
for (j in 1:garchp){
if (i-j > 0){
s1 = (h[N + i - j]^2)
s2 = kappa*gm[j]* (h[N + i - j]^2)
} else{
s1 = epsx[N + i - j]^2
s2 = gm[j]*(epsx[N + i - j]^2)*as.integer(epsx[N + i - j]<0)
}
h[N+i] = h[N+i] + alpha[j] * s1 + s2
}
h[N+i] = sqrt(h[N+i])
}
if(arfima){
res = arfimaf(data[1:N], h, armaOrder, mu, ar, ma, darfima, inmean, im, epsx,
include.mex, mxreg, mxf, n.ahead)
} else{
res = armaf(data[1:N], h, armaOrder, mu, ar, ma, inmean, im,
epsx, include.mex, mxreg, mxf, n.ahead)
}
return(list(h = h[(N+1):(N+n.ahead)], x = res[(N+1):(N+n.ahead)]))
}
#---------------------------------------------------------------------------------
# 2nd dispatch method for forecast
.gjrgarchforecast2 = function(fitORspec, data = NULL, n.ahead = 10, n.roll = 0, out.sample = 0,
external.forecasts = list(mregfor = NULL, vregfor = NULL), ...)
{
spec = fitORspec
if(is.null(data))
stop("\nugarchforecast-->error: data must not be NULL when using a specification!")
# we then extract the data/coefs etc
xdata = .extractdata(data)
Nor = length(as.numeric(xdata$data))
data = xdata$data[1:(Nor - out.sample)]
N = length(as.numeric(data))
dates = xdata$pos[1:(Nor - out.sample)]
origData = xdata$data
origDates = xdata$pos
dformat = xdata$dformat
vmodel = spec@variance.model
mmodel = spec@mean.model
dmodel = spec@distribution.model
omodel = spec@optimization.model
ns = out.sample
if( n.roll > ns )
stop("\nugarchforecast-->error: n.roll must not be greater than out.sample!")
include.mean = mmodel$include.mean
armaOrder = mmodel$armaOrder
garchInMean = mmodel$garchInMean
im = inMeanType = mmodel$inMeanType
arfima = mmodel$arfima
include.mex = mmodel$include.mex
mxn = mmodel$mxn
if(mxn>0) mexdata = matrix(mmodel$mexdata, ncol = mxn) else mexdata = NULL
pow = 2
garchOrder = vmodel$garchOrder
include.vex = vmodel$include.vex
vxn = vmodel$vxn
if(vxn>0) vexdata = matrix(vmodel$vexdata, ncol = vxn) else vexdata = NULL
modelnames = omodel$modelnames
pars = unlist(omodel$fixed.pars)
distribution = dmodel$distribution
include.dlambda = dmodel$include.dlambda
include.skew = dmodel$include.skew
include.shape = dmodel$include.shape
# check if necessary the external regressor forecasts provided first
xreg = .forcregressors(include.mex, mxn, mexdata,
mregfor = external.forecasts$mregfor, include.vex, vxn,
vexdata, external.forecasts$vregfor, pars, n.ahead, Nor,
out.sample = ns, n.roll)
mxreg = xreg$mxreg
mxf = xreg$mxf
vxreg = xreg$vxreg
vxf = xreg$vxf
pos = omodel$pos
ar = armap = ma = armaq = inmean = darfima = 0
alpha = garchp = beta = garchq = gm = 0
skew = shape = dlambda = 0
if(pos[2,3] == 1){
ar = pars[pos[2,1]:pos[2,2]]
armap = length(ar)
}
if(pos[3,3] == 1){
ma = pars[pos[3,1]:pos[3,2]]
armaq = length(ma)
}
if(pos[4,3] == 1) inmean = pars[pos[4,1]:pos[4,2]]
if(pos[5,3] == 1) darfima = pars[pos[5,1]:pos[5,2]]
if(pos[8,3] == 1){
alpha = pars[pos[8,1]:pos[8,2]]
garchp=garchOrder[1]
}
if(pos[9,3] == 1) gm = pars[pos[9,1]:pos[9,2]]
if(pos[10,3] == 1){
beta = pars[pos[10,1]:pos[10,2]]
garchq=garchOrder[2]
}
if(pos[12,3]) dlambda = pars[pos[12,1]:pos[12,2]]
if(pos[13,3]) skew = pars[pos[13,1]:pos[13,2]]
if(pos[14,3]) shape = pars[pos[14,1]:pos[14,2]]
kappa = 0
if(garchp>0) kappa = pneg(dlambda, shape, skew, distribution)
fcreq = ifelse(ns >= (n.ahead+n.roll), n.ahead+n.roll, ns)
fspec = ugarchspec(variance.model = list(model = "gjrGARCH",
garchOrder = c(garchp, garchq), submodel = NULL,
external.regressors = vxf[1:(N + fcreq), , drop = FALSE]),
mean.model = list(armaOrder = c(armap, armaq),
include.mean = include.mean,
garchInMean = garchInMean, inMeanType = inMeanType, arfima = arfima,
external.regressors = mxf[1:(N + fcreq), , drop = FALSE]),
distribution.model = distribution, fixed.pars = as.list(pars))
tmp = data.frame(origData[1:(N + fcreq)])
rownames(tmp) = as.character(origDates[1:(N + fcreq)])
flt = .gjrgarchfilter(data = tmp[,1,drop=FALSE], spec = fspec,
n.old = N)
sigmafilter = flt@filter$sigma
resfilter = flt@filter$residuals
zfilter = flt@filter$z
# forecast GARCH process
forecasts = vector(mode="list", length = n.roll+1)
fwdd = vector(mode="list", length = n.roll+1)
for(i in 1:(n.roll+1)){
np = N + i - 1
if(pos[1,3] == 1){
mu = rep(pars[pos[1,1]:pos[1,2]], N+i+n.ahead-1)
} else{
mu = rep(0, N+i+n.ahead-1)
}
if(pos[7,3] == 1){
omega = rep(pars[pos[7,1]:pos[7,2]], N+i+n.ahead-1)
} else {
omega = rep(0, N+i+n.ahead-1)
}
# no look-ahead
h = c(sigmafilter[1:(N+i-1)], rep(0, n.ahead))
epsx = c(resfilter[1:(N+i-1)], rep(0, n.ahead))
x = c(origData[1:(N+i-1)], rep(0, n.ahead))
z = c(zfilter[1:(N+i-1)], rep(0, n.ahead))
# forecast of externals is provided outside the system
mxfi = mxf[1:(N+i-1+n.ahead), , drop=FALSE]
vxfi = vxf[1:(N+i-1+n.ahead), , drop=FALSE]
ans = .ngjrsgarchforecast(N = np, n.ahead, include.vex, vxfi, vxreg, vxn, h,
omega, alpha, gm, beta, garchp, garchq, epsx, z, data = x,
armaOrder, mu, ar, ma, arfima, inmean, im, include.mex, mxreg, mxfi,
darfima = darfima, kappa)
fdf = data.frame(sigma = ans$h, series = ans$x)
fwdd[[i]] = .forcdates( origDates, n.ahead, N, i, ns, dformat )
rownames(fdf) = as.character(fwdd[[i]])
forecasts[[i]] = fdf
}
fcst = list()
fcst$n.ahead = n.ahead
fcst$N = N+ns
fcst$n.start = ns
fcst$n.roll = n.roll
fcst$sigma = flt@filter$sigma[1:N]
fcst$series = data
fcst$forecasts = forecasts
fcst$fdates = fwdd
model = c(mmodel, vmodel, omodel)
model$dates = dates
model$origdata = origData
model$origdates = origDates
ans = new("gjrGARCHforecast",
forecast = fcst,
model = model,
flt)
return(ans)
}
#---------------------------------------------------------------------------------
# SECTION gjrGARCH simulate
#---------------------------------------------------------------------------------
# choose the fastest method to use:
.gjrgarchsim = function(fit, n.sim = 1000, n.start = 0, m.sim = 1, startMethod =
c("unconditional","sample"), presigma = NA, prereturns = NA,
preresiduals = NA, rseed = NA, custom.dist = list(name = NA, distfit = NA),
mexsimdata = NULL, vexsimdata = NULL)
{
if( (n.sim+n.start) < 100 && m.sim > 100 ){
ans = .gjrgarchsim2(fit = fit, n.sim = n.sim, n.start = n.start, m.sim = m.sim,
startMethod = startMethod, presigma = presigma, prereturns = prereturns,
preresiduals = preresiduals, rseed = rseed, custom.dist = custom.dist,
mexsimdata = mexsimdata, vexsimdata = vexsimdata)
} else{
ans = .gjrgarchsim1(fit = fit, n.sim = n.sim, n.start = n.start, m.sim = m.sim,
startMethod = startMethod, presigma = presigma, prereturns = prereturns,
preresiduals = preresiduals, rseed = rseed, custom.dist = custom.dist,
mexsimdata = mexsimdata, vexsimdata = vexsimdata)
}
return( ans )
}
.gjrgarchsim1 = function(fit, n.sim = 1000, n.start = 0, m.sim = 1, startMethod =
c("unconditional","sample"), presigma = NA, prereturns = NA,
preresiduals = NA, rseed = NA, custom.dist = list(name = NA, distfit = NA),
mexsimdata = NULL, vexsimdata = NULL)
{
# some checks
if(is.na(rseed[1])){
sseed = as.integer(runif(1,0,as.integer(Sys.time())))
} else{
if(length(rseed) != m.sim) sseed = as.integer(rseed[1]) else sseed = rseed[1:m.sim]
}
# Enlarge Series:
n = n.sim + n.start
startMethod = startMethod[1]
data = fit@fit$data
N = length(as.numeric(data))
m = fit@model$maxOrder2
resids = fit@fit$residuals
sigma = fit@fit$sigma
include.mean = fit@model$include.mean
armaOrder = fit@model$armaOrder
garchInMean = fit@model$garchInMean
inMeanType = fit@model$inMeanType
im = inMeanType/2
arfima = fit@model$arfima
include.mex = fit@model$include.mex
mxn = fit@model$mxn
if(mxn>0) mexdata = matrix(fit@model$mexdata, ncol=mxn) else mexdata = NULL
pow = fit@model$pow
garchOrder = fit@model$garchOrder
include.vex = fit@model$include.vex
vxn = fit@model$vxn
if(vxn>0) vexdata = matrix(fit@model$vexdata, ncol=vxn) else vexdata = NULL
modelnames = fit@model$modelnames
pars = fit@fit$coef
distribution = fit@model$distribution
include.dlambda = fit@model$include.dlambda
include.skew = fit@model$include.skew
include.shape = fit@model$include.shape
# check if necessary the external regressor forecasts provided first
xreg = .simregressors(include.mex, mxn, mexsimdata, mexdata, include.vex,
vxn, vexsimdata, vexdata, N, pars, n, m.sim, m)
if(mxn>0) mxreg = matrix(xreg$mxreg, ncol = mxn)
mexsim = xreg$mexsimlist
vxreg = xreg$vxreg
vexsim = xreg$vexsimlist
pos = fit@model$pos
ar = armap = ma = armaq = inmean = darfima = mu = 0
alpha = garchp = beta = garchq = gm = 0
skew = shape = dlambda = 0
omega = 0.00001
if(pos[1,3] == 1) mu = pars[pos[1,1]:pos[1,2]]
if(pos[2,3] == 1){
ar = pars[pos[2,1]:pos[2,2]]
armap = length(ar)
}
if(pos[3,3] == 1){
ma = pars[pos[3,1]:pos[3,2]]
armaq = length(ma)
}
if(pos[4,3] == 1) inmean = pars[pos[4,1]:pos[4,2]]
if(pos[5,3] == 1) darfima = pars[pos[5,1]:pos[5,2]]
if(pos[7,3] == 1) omega = pars[pos[7,1]:pos[7,2]]
if(pos[8,3] == 1){
alpha = pars[pos[8,1]:pos[8,2]]
garchp=garchOrder[1]
}
if(pos[9,3] == 1) gm = pars[pos[9,1]:pos[9,2]]
if(pos[10,3] == 1){
beta = pars[pos[10,1]:pos[10,2]]
garchq=garchOrder[2]
}
if(pos[12,3]) dlambda = pars[pos[12,1]:pos[12,2]]
if(pos[13,3]) skew = pars[pos[13,1]:pos[13,2]]
if(pos[14,3]) shape = pars[pos[14,1]:pos[14,2]]
if(N < n.start){
startmethod[1] = "unconditional"
cat("\nugarchsim-->warning: n.start greater than length of data...using unconditional
start method...\n")
}
# Random Samples from the Distribution
#dst = match(distribution, c("norm","snorm","std","sstd","ged","sged","nig","ghyp","jsu"))
if(length(sseed) == 1){
zmatrix = data.frame(dist = distribution, dlambda = dlambda,
skew = skew, shape = shape, n = n * m.sim, seed = sseed[1])
z = .custzdist(custom.dist, zmatrix, m.sim, n)
} else{
zmatrix = data.frame(dist = rep(distribution, m.sim), dlambda = rep(dlambda, m.sim),
skew = rep(skew, m.sim), shape = rep(shape, m.sim), n = rep(n, m.sim), seed = sseed)
z = .custzdist(custom.dist, zmatrix, m.sim, n)
}
if(startMethod == "unconditional"){
z = rbind(matrix(0, nrow = m, ncol = m.sim), z)
} else{
z = rbind(matrix(tail(fit@fit$z, m), nrow = m, ncol = m.sim), z)
}
# create the presample information
if(!is.na(presigma[1])){
presigma = as.vector(presigma)
if(length(presigma)<m) stop(paste("\nugarchsim-->error: presigma must be of length ", m, sep=""))
}
if(!is.na(prereturns[1])){
prereturns = as.vector(prereturns)
if(length(prereturns)<m) stop(paste("\nugarchsim-->error: prereturns must be of length ", m, sep=""))
}
if(!is.na(preresiduals[1])){
preresiduals = as.vector(preresiduals)
if(length(preresiduals)<m) stop(paste("\nugarchsim-->error: preresiduals must be of length ", m, sep=""))
preres = matrix(preresiduals, nrow = m)
}
if(is.na(presigma[1])){
if(startMethod[1] == "unconditional"){
persist = persistence(fit)
hEst = sqrt(omega/abs(1-persist))
presigma = as.numeric(rep(hEst, times = m))}
else{
presigma = tail(sigma, m)
}
}
if(is.na(prereturns[1])){
if(startMethod[1] == "unconditional"){
prereturns = as.numeric(rep(uncmean(fit), times = m))
}
else{
prereturns = tail(data, m)
}
}
# input vectors/matrices
h = c(presigma * presigma, rep(0, n))
x = c(prereturns, rep(0, n))
constm = matrix(mu, ncol = m.sim, nrow = n + m)
# outpus matrices
sigmaSim = matrix(0, ncol = m.sim, nrow = n.sim)
seriesSim = matrix(0, ncol = m.sim, nrow = n.sim)
residSim = matrix(0, ncol = m.sim, nrow = n.sim)
for(i in 1:m.sim){
if(is.na(preresiduals[1])){
if(startMethod[1] == "unconditional"){
preres = as.numeric(z[1:m, i])*presigma
} else{
preres = tail(resids, m)
}
}
ngrd = which(preres<0)
tmpr = rep(0, length(preres))
tmpr[ngrd] = preres[ngrd] * preres[ngrd]
nres = c(tmpr, rep(0, n))
res = c(preres, rep(0, n))
ans1 = .gjrgarchsimC(garchp, garchq, incvex = vxn, omega, alpha,
beta, gm, vxreg, h, as.numeric(z[,i]), res, nres, vexsim[[i]], T = n + m, m)
sigmaSim[,i] = sqrt( ans1$h[(n.start + m + 1):(n+m)] )
residSim[,i] = ans1$res[(n.start + m + 1):(n+m)]
if(include.mex){
constm[,i] = constm[,i] + mxreg%*%t(matrix(mexsim[[i]], ncol = mxn))
}
if(garchInMean) constm[,i] = constm[,i] + inmean*(ans1$h^im)
if(arfima){
## if(constant) constm[,i] = constm[,i]*(1-sum(ar))
ans2 = .arfimaxsim(armap, armaq, constm[1:n, i], ar, ma, darfima, x, ans1$res[1:(n+armaq)], T = n)
seriesSim[,i] = ans2$s[(n.start + armaq + 1):(n + armaq)]
} else{
# if(constant) constm[,i] = constm[,i]*(1-sum(ar))
ans2 = .armaxsim(armap, armaq, constm[,i], ar, ma, x, ans1$res, T = n + m, m)
seriesSim[,i] = ans2$x[(n.start + m + 1):(n+m)]
}
}
smodel = fit@model
smodel$dates = fit@fit$dates
sim = list(sigmaSim = sigmaSim, seriesSim = seriesSim, residSim = residSim)
sim$series = data
sim$sigma = sigma
sol = new("gjrGARCHsim",
simulation = sim,
model = smodel,
seed = as.integer(sseed))
return(sol)
}
#---------------------------------------------------------------------------------
# SECTION gjrGARCH path
#---------------------------------------------------------------------------------
.gjrgarchpath = function(spec, n.sim = 1000, n.start = 0, m.sim = 1,
presigma = NA, prereturns = NA, preresiduals = NA, rseed = NA,
custom.dist = list(name = NA, distfit = NA), mexsimdata = NULL,
vexsimdata = NULL)
{
if( (n.sim+n.start) < 20 && m.sim > 100 ){
ans = .gjrgarchpath2(spec = spec, n.sim = n.sim, n.start = n.start, m.sim = m.sim,
presigma = presigma, prereturns = prereturns, preresiduals = preresiduals,
rseed = rseed, custom.dist = custom.dist, mexsimdata = mexsimdata,
vexsimdata = vexsimdata)
} else{
ans = .gjrgarchpath1(spec = spec, n.sim = n.sim, n.start = n.start, m.sim = m.sim,
presigma = presigma, prereturns = prereturns, preresiduals = preresiduals,
rseed = rseed, custom.dist = custom.dist, mexsimdata = mexsimdata,
vexsimdata = vexsimdata)
}
return( ans )
}
.gjrgarchpath1 = function(spec, n.sim = 1000, n.start = 0, m.sim = 1,
presigma = NA, prereturns = NA, preresiduals = NA, rseed = NA,
custom.dist = list(name = NA, distfit = NA), mexsimdata = NULL,
vexsimdata = NULL)
{
# some checks
if(is.na(rseed[1])){
sseed = as.integer(runif(1,0,as.integer(Sys.time())))
} else{
if(length(rseed) != m.sim) sseed = as.integer(rseed[1]) else sseed = rseed[1:m.sim]
}
# Enlarge Series:
n = n.sim + n.start
armaOrder = spec@mean.model$armaOrder
armap = armaOrder[1]
armaq = armaOrder[2]
include.mean=spec@mean.model$include.mean
garchOrder = spec@variance.model$garchOrder
garchp = garchOrder[1]
garchq = garchOrder[2]
include.mex = spec@mean.model$include.mex
mxn = spec@mean.model$mxn
N = 0
if(include.mex) {
mexdata = matrix(spec@mean.model$mexdata, ncol=mxn)
N = dim(mexdata)[1]
} else { mexdata = NULL }
include.vex = spec@variance.model$include.vex
vxn = spec@variance.model$vxn
if(include.vex) {
vexdata = matrix(spec@variance.model$vexdata, ncol=vxn)
N = dim(vexdata)[1]
} else { vexdata = NULL }
garchInMean = spec@mean.model$garchInMean
inMeanType = spec@mean.model$inMeanType
im = inMeanType/2
distribution = spec@distribution.model$distribution
pos = spec@optimization.model$pos.matrix
pars = unlist(spec@optimization.model$fixed.pars)
# perform checks on supplied parameters
exp.pars = rownames(spec@optimization.model$modelmatrix)
spec.pars = names(pars)
tmatch = match(sort(spec.pars), sort(exp.pars))
if(any(is.na(tmatch))){
cat("\nugarchpath-->error: supplied parameters do not match expected parameters\n")
cat(paste("epxected parameters:", sort(exp.pars), "\n", sep=""))
cat(paste("supplied parameters:", sort(spec.pars), "\n", sep=""))
stop("\n")
}
omega = 0.00001
mu = ar = ma = inmean = darfima = mxreg = alpha = beta = gm = vxreg = skew = shape = dlambda = 0
if(pos[1,3]) mu = pars[pos[1,1]:pos[1,2]]
if(pos[2,3]) ar = pars[pos[2,1]:pos[2,2]]
if(pos[3,3]) ma = pars[pos[3,1]:pos[3,2]]
if(pos[4,3]) inmean = pars[pos[4,1]:pos[4,2]]
if(pos[5,3]) darfima = pars[pos[5,1]:pos[5,2]]
if(pos[6,3]) mxreg = pars[pos[6,1]:pos[6,2]]
if(pos[7,3]) omega = pars[pos[7,1]:pos[7,2]]
if(pos[8,3]) alpha = pars[pos[8,1]:pos[8,2]]
if(pos[9,3]) gm = pars[pos[9,1]:pos[9,2]]
if(pos[10,3]) beta = pars[pos[10,1]:pos[10,2]]
if(pos[11,3]) vxreg = pars[pos[11,1]:pos[11,2]]
if(pos[12,3]) dlambda = pars[pos[12,1]:pos[12,2]]
if(pos[13,3]) skew = pars[pos[13,1]:pos[13,2]]
if(pos[14,3]) shape = pars[pos[14,1]:pos[14,2]]
arfima = spec@mean.model$arfima
m = spec@optimization.model$maxOrder2
modelnames = spec@optimization.model$modelnames
include.skew = spec@distribution.model$include.skew
include.shape = spec@distribution.model$include.shape
# check if necessary the external regressor forecasts provided first
xreg = .simregressors(include.mex, mxn, mexsimdata, mexdata, include.vex,
vxn, vexsimdata, vexdata, N, pars, n, m.sim, m)
if(mxn>0) mxreg = matrix(xreg$mxreg, ncol = mxn)
mexsim = xreg$mexsimlist
vxreg = xreg$vxreg
vexsim = xreg$vexsimlist
persist = persistence(pars = pars, distribution = distribution, model = "gjrGARCH")
if(persist >= 1)
cat(paste("\nugarchpath->warning: persitence :", round(persist,5), sep=""))
# Random Samples from the Distribution
#dst = match(distribution, c("norm","snorm","std","sstd","ged","sged","nig","ghyp","jsu"))
if(length(sseed) == 1){
zmatrix = data.frame(dist = distribution, dlambda = dlambda,
skew = skew, shape = shape, n = n * m.sim, seed = sseed[1])
z = .custzdist(custom.dist, zmatrix, m.sim, n)
} else{
zmatrix = data.frame(dist = rep(distribution, m.sim), dlambda = rep(dlambda, m.sim),
skew = rep(skew, m.sim), shape = rep(shape, m.sim), n = rep(n, m.sim), seed = sseed)
z = .custzdist(custom.dist, zmatrix, m.sim, n)
}
z = rbind(matrix(0, nrow = m, ncol = m.sim), z)
# create the presample information
if(!is.na(presigma[1])){
presigma = as.vector(presigma)
if(length(presigma)<m) stop(paste("\nugarchsim-->error: presigma must be of length ", m, sep=""))
}
if(!is.na(prereturns[1])){
prereturns = as.vector(prereturns)
if(length(prereturns)<m) stop(paste("\nugarchsim-->error: prereturns must be of length ", m, sep=""))
}
if(!is.na(preresiduals[1])){
preresiduals = as.vector(preresiduals)
if(length(preresiduals)<m) stop(paste("\nugarchsim-->error: preresiduals must be of length ", m, sep=""))
preres = matrix(preresiduals, nrow = m)
}
if(is.na(presigma[1])){
hEst = (omega/abs(1-persist))^(1/2)
presigma = as.numeric(rep(hEst, times = m))
}
if(is.na(prereturns[1])){
prereturns = as.numeric(rep(uncmean(spec), times = m))
}
# input vectors/matrices
h = c(presigma*presigma, rep(0, n))
x = c(prereturns, rep(0, n))
constm = matrix(mu, ncol = m.sim, nrow = n + m)
# outpus matrices
sigmaSim = matrix(0, ncol = m.sim, nrow = n.sim)
seriesSim = matrix(0, ncol = m.sim, nrow = n.sim)
residSim = matrix(0, ncol = m.sim, nrow = n.sim)
for(i in 1:m.sim){
if(is.na(preresiduals[1])){
preres = as.numeric(z[1:m, i])*presigma
}
res = c(preres, rep(0, n))
ngrd = which(preres<0)
tmpr = rep(0, length(preres))
tmpr[ngrd] = preres[ngrd]^2
nres = c(tmpr,rep(0, n))
ans1 = .gjrgarchsimC(garchp, garchq, incvex = vxn, omega, alpha,
beta, gm, vxreg, h, as.numeric(z[,i]), res, nres, vexsim[[i]], T = n + m, m)
sigmaSim[,i] = ans1$h[(n.start + m + 1):(n+m)]^(1/2)
residSim[,i] = ans1$res[(n.start + m + 1):(n+m)]
if(include.mex){
constm[,i] = constm[,i] + mxreg%*%t(matrix(mexsim[[i]], ncol = mxn))
}
if(garchInMean) constm[,i] = constm[,i] + inmean*(ans1$h^im)
if(arfima){
## if(constant) constm[,i] = constm[,i]*(1-sum(ar))
ans2 = .arfimaxsim(armap, armaq, constm[1:n, i], ar, ma, darfima, x, ans1$res[1:(n+armaq)], T = n)
seriesSim[,i] = ans2$s[(n.start + armaq + 1):(n + armaq)]
} else{
# if(constant) constm[,i] = constm[,i]*(1-sum(ar))
ans2 = .armaxsim(armap, armaq, constm[,i], ar, ma, x, ans1$res, T = n + m, m)
seriesSim[,i] = ans2$x[(n.start + m + 1):(n+m)]
}
}
path = list(sigmaSim = sigmaSim, seriesSim = seriesSim, residSim = residSim)
sol = new("gjrGARCHpath",
path = path,
seed = as.integer(sseed))
return(sol)
}
### Timing Tests suggest this is faster when n.sim is small and m.sim is large ###
.gjrgarchsim2 = function(fit, n.sim = 1000, n.start = 0, m.sim = 1, startMethod =
c("unconditional","sample"), presigma = NA, prereturns = NA,
preresiduals = NA, rseed = NA, custom.dist = list(name = NA, distfit = NA),
mexsimdata = NULL, vexsimdata = NULL)
{
# some checks
if(is.na(rseed[1])){
sseed = as.integer(runif(1,0,as.integer(Sys.time())))
} else{
if(length(rseed) != m.sim) sseed = as.integer(rseed[1]) else sseed = as.integer(rseed[1:m.sim])
}
# Enlarge Series:
n = n.sim + n.start
startMethod = startMethod[1]
data = fit@fit$data
N = length(as.numeric(data))
m = fit@model$maxOrder2
resids = fit@fit$residuals
sigma = fit@fit$sigma
include.mean = fit@model$include.mean
armaOrder = fit@model$armaOrder
garchInMean = fit@model$garchInMean
inMeanType = fit@model$inMeanType
im = inMeanType/2
arfima = fit@model$arfima
include.mex = fit@model$include.mex
mxn = fit@model$mxn
if(mxn>0) mexdata = matrix(fit@model$mexdata, ncol=mxn) else mexdata = NULL
pow = fit@model$pow
garchOrder = fit@model$garchOrder
include.vex = fit@model$include.vex
vxn = fit@model$vxn
if(vxn>0) vexdata = matrix(fit@model$vexdata, ncol=vxn) else vexdata = NULL
modelnames = fit@model$modelnames
pars = fit@fit$coef
distribution = fit@model$distribution
include.dlambda = fit@model$include.dlambda
include.skew = fit@model$include.skew
include.shape = fit@model$include.shape
# check if necessary the external regressor forecasts provided first
xreg = .simregressors(include.mex, mxn, mexsimdata, mexdata, include.vex,
vxn, vexsimdata, vexdata, N, pars, n, m.sim, m)
if(mxn>0) mxreg = matrix(xreg$mxreg, ncol = mxn)
mexsim = xreg$mexsimlist
vxreg = xreg$vxreg
vexsim = xreg$vexsimlist
pos = fit@model$pos
ar = armap = ma = armaq = inmean = darfima = mu = 0
alpha = garchp = beta = garchq = gm = 0
skew = shape = dlambda = 0
omega = 0.00001
if(pos[1,3] == 1) mu = pars[pos[1,1]:pos[1,2]]
if(pos[2,3] == 1){
ar = pars[pos[2,1]:pos[2,2]]
armap = length(ar)
}
if(pos[3,3] == 1){
ma = pars[pos[3,1]:pos[3,2]]
armaq = length(ma)
}
if(pos[4,3] == 1) inmean = pars[pos[4,1]:pos[4,2]]
if(pos[5,3] == 1) darfima = pars[pos[5,1]:pos[5,2]]
if(pos[7,3] == 1) omega = pars[pos[7,1]:pos[7,2]]
if(pos[8,3] == 1){
alpha = pars[pos[8,1]:pos[8,2]]
garchp=garchOrder[1]
}
if(pos[9,3] == 1) gm = pars[pos[9,1]:pos[9,2]]
if(pos[10,3] == 1){
beta = pars[pos[10,1]:pos[10,2]]
garchq=garchOrder[2]
}
if(pos[12,3]) dlambda = pars[pos[12,1]:pos[12,2]]
if(pos[13,3]) skew = pars[pos[13,1]:pos[13,2]]
if(pos[14,3]) shape = pars[pos[14,1]:pos[14,2]]
if(N < n.start){
startmethod[1] = "unconditional"
cat("\nugarchsim-->warning: n.start greater than length of data...using unconditional start method...\n")
}
# Random Samples from the Distribution # CHANGE for speed up
#dst = match(distribution, c("norm","snorm","std","sstd","ged","sged","nig","ghyp","jsu"))
if(length(sseed) == 1){
zmatrix = data.frame(dist = distribution, dlambda = dlambda,
skew = skew, shape = shape, n = n * m.sim, seed = sseed[1])
z = .custzdist(custom.dist, zmatrix, m.sim, n)
} else{
zmatrix = data.frame(dist = rep(distribution, m.sim), dlambda = rep(dlambda, m.sim),
skew = rep(skew, m.sim), shape = rep(shape, m.sim), n = rep(n, m.sim), seed = sseed)
z = .custzdist(custom.dist, zmatrix, m.sim, n)
}
if(startMethod == "unconditional"){
z = rbind(matrix(0, nrow = m, ncol = m.sim), z)
} else{
z = rbind(matrix(tail(fit@fit$z, m), nrow = m, ncol = m.sim), z)
}
# create the presample information
if(!is.na(presigma[1])){
presigma = as.vector(presigma)
if(length(presigma)<m) stop(paste("\nugarchsim-->error: presigma must be of length ", m, sep=""))
}
if(!is.na(prereturns[1])){
prereturns = as.vector(prereturns)
if(length(prereturns)<m) stop(paste("\nugarchsim-->error: prereturns must be of length ", m, sep=""))
}
if(!is.na(preresiduals[1])){
preresiduals = as.vector(preresiduals)
if(length(preresiduals)<m) stop(paste("\nugarchsim-->error: preresiduals must be of length ", m, sep=""))
preres = matrix(preresiduals, nrow = m)
}
if(is.na(presigma[1])){
if(startMethod[1] == "unconditional"){
persist = persistence(fit)
hEst = sqrt(omega/abs(1-persist))
presigma = as.numeric(rep(hEst, times = m))}
else{
presigma = tail(sigma, m)
}
}
if(is.na(prereturns[1])){
if(startMethod[1] == "unconditional"){
prereturns = as.numeric(rep(uncmean(fit), times = m))
}
else{
prereturns = tail(data, m)
}
}
# input vectors/matrices
h = matrix(c(presigma * presigma, rep(0, n)), nrow = n + m, ncol = m.sim)
x = matrix(c(prereturns, rep(0, n)), nrow = n + m, ncol = m.sim)
constm = matrix(mu, nrow = n + m, ncol = m.sim)
# outpus matrices
sigmaSim = matrix(0, ncol = m.sim, nrow = n.sim)
seriesSim = matrix(0, ncol = m.sim, nrow = n.sim)
residSim = matrix(0, ncol = m.sim, nrow = n.sim)
if(is.na(preresiduals[1])){
if(startMethod[1] == "unconditional"){
preres = matrix( z[1:m, 1:m.sim] * presigma, nrow = m, ncol = m.sim )
} else{
preres = matrix(tail(resids, m), nrow = m, ncol = m.sim)
}
}
ngrd = which(preres<0, arr.ind = TRUE)
tmpr = matrix(0, nrow = dim(preres)[1], ncol = m.sim)
if(dim(ngrd)[1]>0) tmpr[ngrd] = preres[ngrd]*preres[ngrd]
nres = rbind(tmpr, matrix(0, nrow = n, ncol = m.sim))
res = rbind(preres, matrix(0, nrow = n, ncol = m.sim))
nindx = matrix(1, nrow = n + m, ncol = m.sim)
nindx[which(z>0, arr.ind = TRUE)] = 0.0
garchindx = as.numeric( c(garchp, garchq, m) )
# we'll do the external regressors first for speed.
if(include.vex){
vxs = sapply(vexsim, FUN = function(x) vxreg%*%t(matrix(x, ncol = vxn)))
} else{
vxs = matrix(0, nrow = m + n, ncol = m.sim)
}
e = res * res
ans = .Call("mgjrgarchsim", garchindx = garchindx, omega = as.numeric(omega), alpha = as.numeric(alpha),
beta = as.numeric(beta), gm = as.numeric(gm), h = h, z = z, res = res, e = e, nres = nres,
nindx = nindx, vxs = vxs, PACKAGE = "rgarch")
# we have to re-check the matrix argument since when n.sim = 1 (i.e. 1 row), it returns a vector!
# (THIS IS SO ANNOYING!!!!)
sigmaSim = matrix(sqrt( ans$h[(n.start + m + 1):(n+m), ] ), ncol = m.sim)
residSim = matrix(ans$res[(n.start + m + 1):(n+m), ], ncol = m.sim)
if(include.mex){
mxs = sapply(mexsim, FUN = function(x) mxreg%*%t(matrix(x, ncol = mxn)))
} else{
mxs = 0
}
if(garchInMean){
imh = inmean*(sigmaSim^im)
} else{
imh = 0
}
constm = constm + mxs + imh
armaindx = c(armap, armaq, m, n)
if(arfima){
## if(constant) constm[,i] = constm[,i]*(1-sum(ar))
for(i in 1:m.sim){
tmp = .arfimaxsim(armap, armaq, constm[1:n, i], ar, ma, darfima, x, ans$res[1:(n+armaq), i], T = n)
seriesSim[,i] = tmp$s[(n.start + armaq + 1):(n + armaq)]
}
} else{
# if(constant) constm = constm * ( 1 - sum(ar) )
tmp = .Call("marmaxsim", armaindx = armaindx, ar = ar, ma = ma, mu = constm, x = x, res = ans$res,
PACKAGE = "rgarch")
seriesSim = matrix(tmp$x[(n.start + m + 1):(n+m), ], ncol = m.sim)
}
smodel = fit@model
smodel$dates = fit@fit$dates
sim = list(sigmaSim = sigmaSim, seriesSim = seriesSim, residSim = residSim)
sim$series = data
sim$sigma = sigma
sol = new("gjrGARCHsim",
simulation = sim,
model = smodel,
seed = as.integer(sseed))
return(sol)
}
.gjrgarchpath2 = function(spec, n.sim = 1000, n.start = 0, m.sim = 1,
presigma = NA, prereturns = NA, preresiduals = NA, rseed = NA,
custom.dist = list(name = NA, distfit = NA), mexsimdata = NULL,
vexsimdata = NULL)
{
# some checks
if(is.na(rseed[1])){
sseed = as.integer(runif(1,0,as.integer(Sys.time())))
} else{
if(length(rseed) != m.sim) sseed = as.integer(rseed[1]) else sseed = rseed[1:m.sim]
}
# Enlarge Series:
n = n.sim + n.start
armaOrder = spec@mean.model$armaOrder
armap = armaOrder[1]
armaq = armaOrder[2]
include.mean=spec@mean.model$include.mean
garchOrder = spec@variance.model$garchOrder
garchp = garchOrder[1]
garchq = garchOrder[2]
include.mex = spec@mean.model$include.mex
mxn = spec@mean.model$mxn
N = 0
if(include.mex) {
mexdata = matrix(spec@mean.model$mexdata, ncol=mxn)
N = dim(mexdata)[1]
} else { mexdata = NULL }
include.vex = spec@variance.model$include.vex
vxn = spec@variance.model$vxn
if(include.vex) {
vexdata = matrix(spec@variance.model$vexdata, ncol=vxn)
N = dim(vexdata)[1]
} else { vexdata = NULL }
garchInMean = spec@mean.model$garchInMean
inMeanType = spec@mean.model$inMeanType
im = inMeanType/2
distribution = spec@distribution.model$distribution
pos = spec@optimization.model$pos.matrix
pars = unlist(spec@optimization.model$fixed.pars)
# perform checks on supplied parameters
exp.pars = rownames(spec@optimization.model$modelmatrix)
spec.pars = names(pars)
tmatch = match(sort(spec.pars), sort(exp.pars))
if(any(is.na(tmatch))){
cat("\nugarchpath-->error: supplied parameters do not match expected parameters\n")
cat(paste("epxected parameters:", sort(exp.pars), "\n", sep=""))
cat(paste("supplied parameters:", sort(spec.pars), "\n", sep=""))
stop("\n")
}
omega = 0.00001
mu = ar = ma = inmean = darfima = mxreg = alpha = beta = gm = vxreg = skew = shape = dlambda = 0
if(pos[1,3]) mu = pars[pos[1,1]:pos[1,2]]
if(pos[2,3]) ar = pars[pos[2,1]:pos[2,2]]
if(pos[3,3]) ma = pars[pos[3,1]:pos[3,2]]
if(pos[4,3]) inmean = pars[pos[4,1]:pos[4,2]]
if(pos[5,3]) darfima = pars[pos[5,1]:pos[5,2]]
if(pos[6,3]) mxreg = pars[pos[6,1]:pos[6,2]]
if(pos[7,3]) omega = pars[pos[7,1]:pos[7,2]]
if(pos[8,3]) alpha = pars[pos[8,1]:pos[8,2]]
if(pos[9,3]) gm = pars[pos[9,1]:pos[9,2]]
if(pos[10,3]) beta = pars[pos[10,1]:pos[10,2]]
if(pos[11,3]) vxreg = pars[pos[11,1]:pos[11,2]]
if(pos[12,3]) dlambda = pars[pos[12,1]:pos[12,2]]
if(pos[13,3]) skew = pars[pos[13,1]:pos[13,2]]
if(pos[14,3]) shape = pars[pos[14,1]:pos[14,2]]
arfima = spec@mean.model$arfima
m = spec@optimization.model$maxOrder2
modelnames = spec@optimization.model$modelnames
include.skew = spec@distribution.model$include.skew
include.shape = spec@distribution.model$include.shape
# check if necessary the external regressor forecasts provided first
xreg = .simregressors(include.mex, mxn, mexsimdata, mexdata, include.vex,
vxn, vexsimdata, vexdata, N, pars, n, m.sim, m)
if(mxn>0) mxreg = matrix(xreg$mxreg, ncol = mxn)
mexsim = xreg$mexsimlist
vxreg = xreg$vxreg
vexsim = xreg$vexsimlist
persist = persistence(pars = pars, distribution = distribution, model = "gjrGARCH")
if(persist >= 1)
cat(paste("\nugarchpath->warning: persitence :", round(persist,5), sep=""))
# Random Samples from the Distribution
#dst = match(distribution, c("norm","snorm","std","sstd","ged","sged","nig","ghyp","jsu"))
if(length(sseed) == 1){
zmatrix = data.frame(dist = distribution, dlambda = dlambda,
skew = skew, shape = shape, n = n * m.sim, seed = sseed[1])
z = .custzdist(custom.dist, zmatrix, m.sim, n)
} else{
zmatrix = data.frame(dist = rep(distribution, m.sim), dlambda = rep(dlambda, m.sim),
skew = rep(skew, m.sim), shape = rep(shape, m.sim), n = rep(n, m.sim), seed = sseed)
z = .custzdist(custom.dist, zmatrix, m.sim, n)
}
z = rbind(matrix(0, nrow = m, ncol = m.sim), z)
# create the presample information
if(!is.na(presigma[1])){
presigma = as.vector(presigma)
if(length(presigma)<m) stop(paste("\nugarchsim-->error: presigma must be of length ", m, sep=""))
}
if(!is.na(prereturns[1])){
prereturns = as.vector(prereturns)
if(length(prereturns)<m) stop(paste("\nugarchsim-->error: prereturns must be of length ", m, sep=""))
}
if(!is.na(preresiduals[1])){
preresiduals = as.vector(preresiduals)
if(length(preresiduals)<m) stop(paste("\nugarchsim-->error: preresiduals must be of length ", m, sep=""))
preres = matrix(preresiduals, nrow = m)
}
if(is.na(presigma[1])){
hEst = (omega/abs(1-persist))^(1/2)
presigma = as.numeric(rep(hEst, times = m))
}
if(is.na(prereturns[1])){
prereturns = as.numeric(rep(uncmean(spec), times = m))
}
# input vectors/matrices
h = matrix(c(presigma * presigma, rep(0, n)), nrow = n + m, ncol = m.sim)
x = matrix(c(prereturns, rep(0, n)), nrow = n + m, ncol = m.sim)
constm = matrix(mu, nrow = n + m, ncol = m.sim)
# outpus matrices
sigmaSim = matrix(0, ncol = m.sim, nrow = n.sim)
seriesSim = matrix(0, ncol = m.sim, nrow = n.sim)
residSim = matrix(0, ncol = m.sim, nrow = n.sim)
if(is.na(preresiduals[1])){
preres = matrix( z[1:m, 1:m.sim] * presigma, nrow = m, ncol = m.sim )
} else{
preres = matrix( preresiduals, nrow = m, ncol = m.sim )
}
ngrd = which(preres<0, arr.ind = TRUE)
tmpr = matrix(0, nrow = dim(preres)[1], ncol = m.sim)
if(dim(ngrd)[1]>0) tmpr[ngrd] = preres[ngrd]*preres[ngrd]
nres = rbind(tmpr, matrix(0, nrow = n, ncol = m.sim))
res = rbind(preres, matrix(0, nrow = n, ncol = m.sim))
nindx = matrix(1, nrow = n + m, ncol = m.sim)
nindx[which(z>0, arr.ind = TRUE)] = 0.0
garchindx = as.numeric( c(garchp, garchq, m, n) )
# we'll do the external regressors first for speed.
if(include.vex){
vxs = sapply(vexsim, FUN = function(x) vxreg%*%t(matrix(x, ncol = vxn)))
} else{
vxs = matrix(0, nrow = m + n, ncol = m.sim)
}
e = res * res
ans = .Call("mgjrgarchsim", garchindx = garchindx, omega = as.numeric(omega), alpha = as.numeric(alpha),
beta = as.numeric(beta), gm = as.numeric(gm), h = h, z = z, res = res, e = e, nres = nres,
nindx = nindx, vxs = vxs, PACKAGE = "rgarch", DUP = FALSE)
sigmaSim = matrix(sqrt( ans$h[(n.start + m + 1):(n+m), ] ), ncol = m.sim)
residSim = matrix(ans$res[(n.start + m + 1):(n+m), ], ncol = m.sim)
if(include.mex){
mxs = sapply(mexsim, FUN = function(x) mxreg%*%t(matrix(x, ncol = mxn)))
} else{
mxs = 0
}
if(garchInMean){
imh = inmean*(sigmaSim^im)
} else{
imh = 0
}
constm = constm + mxs + imh
armaindx = c(armap, armaq, m, n)
if(arfima){
## if(constant) constm[,i] = constm[,i]*(1-sum(ar))
for(i in 1:m.sim){
tmp = .arfimaxsim(armap, armaq, constm[1:n, i], ar, ma, darfima, x, ans$res[1:(n+armaq), i], T = n)
seriesSim[,i] = tmp$s[(n.start + armaq + 1):(n + armaq)]
}
} else{
# if(constant) constm = constm * ( 1 - sum(ar) )
tmp = .Call("marmaxsim", armaindx = armaindx, ar = ar, ma = ma, mu = constm, x = x, res = ans$res)
seriesSim = matrix(tmp$x[(n.start + m + 1):(n+m), ], ncol = m.sim)
}
path = list(sigmaSim = sigmaSim, seriesSim = seriesSim, residSim = residSim)
sol = new("gjrGARCHpath",
path = path,
seed = as.integer(sseed))
return(sol)
}
Try the rgarch package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.