R/rgarch-sgarch.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 sGARCH spec
#---------------------------------------------------------------------------------
.sgarchspec = 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
	im = inmeanType/2
	# inmeanType(1) means sd
	# inmeanType(2) means var
	#if(garchInMean && (inmeanType!=1 || inmeanType!=2)) stop("valid inmeanType values are 1 and 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(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 = 13)
	colnames(pos.matrix) = c("start", "stop", "include")
	rownames(pos.matrix) = c("mu", "ar", "ma", "inmean", "darfima", "mxreg", "omega", 
			"alpha", "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(garchq > 0){pos.matrix[9,1:3] = c(pos, pos+garchq-1, 1)
		pos = max(pos.matrix[1:9,2])+1}
	if(include.vex){pos.matrix[10,1:3] = c(pos, pos+vxn-1, 1)
		pos = max(pos.matrix[1:10,2])+1}
	if(include.dlambda){pos.matrix[11,1:3] = c(pos, pos, 1)
		pos = max(pos.matrix[1:11,2])+1}
	if(include.skew){pos.matrix[12,1:3] = c(pos, pos, 1)
		pos = max(pos.matrix[1:12,2])+1}
	if(include.shape){pos.matrix[13,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 - modelmatrix[,2]
	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("sGARCHspec",
			mean.model = mmodel,
			variance.model = vmodel,
			distribution.model = dmodel,
			optimization.model = omodel)
	return(ans)
}




#---------------------------------------------------------------------------------
# SECTION sGARCH fit
#---------------------------------------------------------------------------------
.sgarchfit = 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
	# create a temporary environment to store values (deleted at end of function)
	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("trace", trace, envir = tempenvir)
	
	assign("fit.control", fit.control, 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
	zdata = data/dscale
	assign("dscale", dscale, envir = tempenvir)
	# 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 = .sgarchcon
		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 = .sgarchLLH, 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 = "sGARCH", f = .sgarchLLH, 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 sGARCHfit
	ans = new("sGARCHfit",
			fit = fit,
			model = model)
	rm(tempenvir)
	return(ans)
}

#---------------------------------------------------------------------------------
# SECTION sGARCH LLH
#---------------------------------------------------------------------------------
.sgarchLLH = 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
	omega = mu = ar = ma = inmean = darfima = mxreg = alpha = beta = vxreg = skew = shape = dlambda = 0
	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]]
	if(pos[3,3] == 1) ma = pars[pos[3,1]:pos[3,2]]
	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[6,3] == 1) mxreg = pars[pos[6,1]:pos[6,2]]
	if(pos[8,3] == 1) alpha = pars[pos[8,1]:pos[8,2]]
	if(pos[9,3] == 1) beta = pars[pos[9,1]:pos[9,2]]
	if(pos[10,3] == 1) vxreg = pars[pos[10,1]:pos[10,2]]
	if(pos[11,3] == 1) dlambda = pars[pos[11,1]:pos[11,2]]
	if(pos[12,3] == 1) skew = pars[pos[12,1]:pos[12,2]]
	if(pos[13,3] == 1) shape = pars[pos[13,1]:pos[13,2]]
	
	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
	
	# sgarch persistence value
	kappa = 1
	persist = (sum(alpha) + sum(beta))
	
	# unconditional sigma value
	mvar = mean(res*res)
	if(pos[10,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 + omega
	} else{
		if(omega.calculate){
			omega = mvar * (1 - sum(alpha) - sum(beta)) - mv
			hEst = mvar
			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. Stationarity conditions not valid for
	# solnp sover which implements them internally as constraints.
	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 [sgarch.c] via wrapper

	ans = .sgarchfitC(armaOrder, garchOrder, mu, ar, ma, inmean, mxreg, mexdata,
		omega, hEst, alpha, beta, 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[(m):T]
	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 sGARCH filter
#---------------------------------------------------------------------------------
.sgarchfilter = function(spec, data, out.sample = 0, n.old = NULL)
{
	# n.old is optional and indicates the length of the original dataseries (in
	# cases when this represents a dataseries augmented by newer data). The reason
	# for using this is so that the old and new datasets agree since the original
	# recursion uses the sum of the residuals to start the recursion and therefore
	# is influenced by new data. For a small augmentation the values converge after
	# x periods, but it is sometimes preferable to have this option so that there is
	# no forward looking information contaminating the study.
	.garchenv = environment()
	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]
	if(!is.null(n.old)) Nx = n.old else Nx = length(origdata)
	
	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("\nugarchfilter-->error: parameters do not match expected specification length")
	if(is.na(all(match(modelnames, parnames), 1:length(modelnames)))) {
		stop("\nugarchfilter-->error: 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 = vxreg = skew = shape = dlambda = 0
	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]]
	if(pos[3,3] == 1) ma = pars[pos[3,1]:pos[3,2]]
	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[6,3] == 1) mxreg = pars[pos[6,1]:pos[6,2]]
	if(pos[8,3] == 1) alpha = pars[pos[8,1]:pos[8,2]]
	if(pos[9,3] == 1) beta = pars[pos[9,1]:pos[9,2]]
	if(pos[10,3] == 1) vxreg = pars[pos[10,1]:pos[10,2]]
	if(pos[11,3] == 1) dlambda = pars[pos[11,1]:pos[11,2]]
	if(pos[12,3] == 1) skew = pars[pos[12,1]:pos[12,2]]
	if(pos[13,3] == 1) shape = pars[pos[13,1]:pos[13,2]]
	
	# sgarch persistence value
	kappa = 1
	persist = (sum(alpha) + sum(beta))
	#if(persist>1) stop("\nPersistence of model is > 1...try again!\n")
	
	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)
	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)
		res2 = rx2$res
		# unconditional sigma value
		mvar = mean(res2*res2)
	} else{
		mvar = mean(res*res)
	}

	#hEst = omega + persistence*mvar
	if(pos[10,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 + omega
	} else{
		omega = mvar * (1 - sum(alpha) - sum(beta)) - mv
		hEst = mvar
	}
	
	ans = .sgarchfitC(armaOrder, garchOrder, mu, ar, ma, inmean, mxreg, mexdata,
			omega, hEst, alpha, beta, vxreg, vexdata, dlambda, skew, shape, dist, data, zrf, 
			N, res, garchenv = .garchenv)
	
	filter = list()
	filter$z = ans$z
	filter$sigma = sqrt(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)
	rm(.garchenv)
	sol = new("sGARCHfilter",
			filter = filter,
			model = model)
	return(sol)
}

#---------------------------------------------------------------------------------
# SECTION sGARCH forecast
#---------------------------------------------------------------------------------
.sgarchforecast = function(fitORspec, data = NULL, n.ahead = 10, n.roll = 0, out.sample = 0, 
		external.forecasts = list(mregfor = NULL, vregfor = NULL), ...)
{
	# Recipe for a GARCH forecast:
	# 1. split forecast into arma and garch forecasts.
	# 2. take into account garch-in-mean, arfima, exogenous regressors
	# 3. check kappa and persistence
	# 4. allow for rolling forecast based on filtering out.sample data from fit stage
	
	# n.roll signifies whether to filter/roll the sigma else will use unconditional
	# i.e. n.ahead=10 with n.roll=1 means that the first forecast is based on the
	# previous value whereas all subsequent forecasts are based on the unconditional
	# expectation formula
	# with n.roll = n means that while n.ahead < n.roll we filter the data using the coef
	# of the garch model to obtain filtered sigma estimates which are used to roll
	# the forecast. This requires that during the fit process the out.sample option
	# was used and that (n.roll) < out.sample (otherwise we revert to the unconditional
	# expectation formula for the long run sigma.
	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("\nugarchforecast-->error: 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 = 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){
		beta = pars[pos[9,1]:pos[9,2]]
		garchq = garchOrder[2]
	}
	if(pos[11,3] == 1) dlambda = pars[pos[11,1]:pos[11,2]]
	if(pos[12,3] == 1) skew = pars[pos[12,1]:pos[12,2]]
	if(pos[13,3] == 1) shape = pars[pos[13,1]:pos[13,2]]
	
	
	# filter data (check external regressor data - must equal length of origData)
	fcreq = ifelse(ns >= (n.ahead+n.roll), n.ahead+n.roll, ns)
	fspec = ugarchspec(variance.model = list(model = "sGARCH", 
					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 = .sgarchfilter(data = tmp[,1,drop=FALSE], spec = fspec, 
			n.old = length(fit@fit$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 = .nsgarchforecast(N = np, n.ahead, include.vex, vxfi, vxreg, vxn, h, 
				omega, alpha, beta, garchp, garchq, epsx, z, data = x, armaOrder, 
				mu, ar, ma, arfima, inmean, im, include.mex, mxreg, mxfi, 
				darfima = darfima)
		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("sGARCHforecast",
			forecast = fcst,
			model = model,
			flt)
	return(ans)
}

.nsgarchforecast = function(N, n.ahead, include.vex, vxf, vxreg, vxn, h, omega,
		alpha, beta, garchp, garchq, epsx, z, data, armaOrder, mu, ar, ma, 
		arfima, inmean, im, include.mex, mxreg, mxf, darfima)
{
	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){
				s = h[N + i - j]^2
			} else{ 
				s = epsx[N + i - j]^2
			}
			h[N+i] = h[N+i] + alpha[j] * s
		}
		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
.sgarchforecast2 = function(fitORspec, data = NULL, n.ahead = 10, n.roll = 0, out.sample = 0, 
		external.forecasts = list(mregfor = NULL, vregfor = NULL), ...)
{
	# first we filter the data to get the results:
	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 = 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){
		beta = pars[pos[9,1]:pos[9,2]]
		garchq=garchOrder[2]
	}
	if(pos[11,3] == 1) dlambda = pars[pos[11,1]:pos[11,2]]
	if(pos[12,3] == 1) skew = pars[pos[12,1]:pos[12,2]]
	if(pos[13,3] == 1) shape = pars[pos[13,1]:pos[13,2]]
	
	
	# filter data (check external regressor data - must equal length of origData)
	fcreq = ifelse(ns >= (n.ahead+n.roll), n.ahead+n.roll, ns)
	fspec = ugarchspec(variance.model = list(model = "sGARCH", 
					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 = .sgarchfilter(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 = .nsgarchforecast(N = np, n.ahead, include.vex, vxfi, vxreg, vxn, h, 
				omega, alpha, beta, garchp, garchq, epsx, z, data = x, armaOrder, 
				mu, ar, ma, arfima, inmean, im, include.mex, mxreg, mxfi, 
				darfima = darfima)	
		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("sGARCHforecast",
			forecast = fcst,
			model = model,
			flt)
	return(ans)
	
}
#---------------------------------------------------------------------------------
# SECTION sGARCH simulate
#---------------------------------------------------------------------------------
.sgarchsim = 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 = .sgarchsim2(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 = .sgarchsim1(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 )
}
.sgarchsim1 = 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:
	# need to allow for arfima case:
	arfima = fit@model$arfima
	armap = fit@model$armaOrder[1]
	armaq = fit@model$armaOrder[2]
	if(arfima) n.start = max(armap+armaq, n.start)
	
	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
	#if(im==1) im = 0.5 else im = 1
	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 = 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){
		beta = pars[pos[9,1]:pos[9,2]]
		garchq=garchOrder[2]
	}
	if(pos[11,3] == 1) dlambda = pars[pos[11,1]:pos[11,2]]
	if(pos[12,3] == 1) skew = pars[pos[12,1]:pos[12,2]]
	if(pos[13,3] == 1) shape = pars[pos[13,1]:pos[13,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
	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"){
			kappa = 1
			persist = (sum(alpha)*kappa + sum(beta))
			hEst = (omega/abs(1-persist))^(1/2)
			presigma = as.numeric(rep(hEst, m))}
		else{
			presigma  = tail(sigma, m)
		}
	}
	if(is.na(prereturns[1])){
		if(startMethod[1] == "unconditional"){
			prereturns = as.numeric(rep(uncmean(fit), m))
		}
		else{
			prereturns = tail(data, m)
		}
	}
	
	# input vectors/matrices
	h = c(presigma^2, 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)
			}
		}
		res = c(preres, rep(0, n))
		
		ans1 = .sgarchsimC(garchp, garchq, incvex = vxn, omega, alpha, beta, vxreg, 
				h, as.numeric(z[,i]), res, 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)]
		}
	}
	smodel = fit@model
	smodel$dates = fit@fit$dates
	sim = list(sigmaSim = sigmaSim, seriesSim = seriesSim, residSim = residSim)
	sim$series = data
	sim$sigma = sigma
	sol = new("sGARCHsim",
			simulation = sim,
			model = smodel,
			seed = as.integer(sseed))
	return(sol)
}
.sgarchsim2 = 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:
	# need to allow for arfima case:
	arfima = fit@model$arfima
	armap = fit@model$armaOrder[1]
	armaq = fit@model$armaOrder[2]
	if(arfima) n.start = max(armap+armaq, n.start)
	
	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
	#if(im==1) im = 0.5 else im = 1
	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 = 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){
		beta = pars[pos[9,1]:pos[9,2]]
		garchq=garchOrder[2]
	}
	if(pos[11,3] == 1) dlambda = pars[pos[11,1]:pos[11,2]]
	if(pos[12,3] == 1) skew = pars[pos[12,1]:pos[12,2]]
	if(pos[13,3] == 1) shape = pars[pos[13,1]:pos[13,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
	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"){
			kappa = 1
			persist = (sum(alpha)*kappa + sum(beta))
			hEst = (omega/abs(1-persist))^(1/2)
			presigma = as.numeric(rep(hEst, m))}
		else{
			presigma  = tail(sigma, m)
		}
	}
	if(is.na(prereturns[1])){
		if(startMethod[1] == "unconditional"){
			prereturns = as.numeric(rep(uncmean(fit), 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)
		}
	}
	res =  rbind(preres, matrix(0, nrow = n, ncol = m.sim))
	garchindx = as.numeric( c(garchp, garchq, m) )
	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("msgarchsim", garchindx = garchindx, omega = as.numeric(omega), alpha = as.numeric(alpha), 
			beta = as.numeric(beta), h = h, z = z, res = res, e = e, vxs = vxs, PACKAGE = "rgarch")
	
	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("sGARCHsim",
			simulation = sim,
			model = smodel,
			seed = as.integer(sseed))
	return(sol)
}


#---------------------------------------------------------------------------------
# SECTION sGARCH path
#---------------------------------------------------------------------------------
.sgarchpath = 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 = .sgarchpath2(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 = .sgarchpath1(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 )
}

.sgarchpath1 = 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("expected parameters: \n")
		cat(paste(sort(exp.pars), " ", sep="")) 
		cat("\nsupplied parameters: \n")
		cat(paste(sort(spec.pars), " ", sep="")) 
		stop("\n", call. = FALSE)
	}
	
	omega = 0.00001
	mu = ar = ma = inmean = darfima = mxreg = alpha = beta = vxreg = skew = shape = dlambda = 0
	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]]
	if(pos[3,3] == 1) ma = pars[pos[3,1]:pos[3,2]]
	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[6,3] == 1) mxreg = pars[pos[6,1]:pos[6,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]]
	if(pos[9,3] == 1) beta = pars[pos[9,1]:pos[9,2]]
	if(pos[10,3] == 1) vxreg = pars[pos[10,1]:pos[10,2]]
	if(pos[11,3] == 1) dlambda = pars[pos[11,1]:pos[11,2]]
	if(pos[12,3] == 1) skew = pars[pos[12,1]:pos[12,2]]
	if(pos[13,3] == 1) shape = pars[pos[13,1]:pos[13,2]]
	
	arfima = spec@mean.model$arfima
	m = spec@optimization.model$maxOrder2
	modelnames = spec@optimization.model$modelnames
	
	include.dlambda = spec@distribution.model$include.dlambda
	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
	
	kappa = 1
	persist = (sum(alpha)*kappa + sum(beta))
	if(persist >= 1) 
		cat(paste("\nugarchpath->warning: persitence :", round(persist,5), sep=""))
	
	# Random Samples from the Distribution
	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^2, 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))
		
		ans1 = .sgarchsimC(garchp, garchq, incvex = vxn, omega, alpha, beta, vxreg, 
				h, as.numeric(z[,i]), res, 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("sGARCHpath",
			path = path,
			seed = as.integer(sseed))
	return(sol)
}

.sgarchpath2 = 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("expected parameters: \n")
		cat(paste(sort(exp.pars), " ", sep="")) 
		cat("\nsupplied parameters: \n")
		cat(paste(sort(spec.pars), " ", sep="")) 
		stop("\n", call. = FALSE)
	}
	
	omega = 0.00001
	mu = ar = ma = inmean = darfima = mxreg = alpha = beta = vxreg = skew = shape = dlambda = 0
	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]]
	if(pos[3,3] == 1) ma = pars[pos[3,1]:pos[3,2]]
	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[6,3] == 1) mxreg = pars[pos[6,1]:pos[6,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]]
	if(pos[9,3] == 1) beta = pars[pos[9,1]:pos[9,2]]
	if(pos[10,3] == 1) vxreg = pars[pos[10,1]:pos[10,2]]
	if(pos[11,3] == 1) dlambda = pars[pos[11,1]:pos[11,2]]
	if(pos[12,3] == 1) skew = pars[pos[12,1]:pos[12,2]]
	if(pos[13,3] == 1) shape = pars[pos[13,1]:pos[13,2]]

	arfima = spec@mean.model$arfima
	m = spec@optimization.model$maxOrder2
	modelnames = spec@optimization.model$modelnames
	
	include.dlambda = spec@distribution.model$include.dlambda
	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

	kappa = 1
	persist = (sum(alpha)*kappa + sum(beta))
	if(persist >= 1) 
		cat(paste("\nugarchpath->warning: persitence :", round(persist,5), sep=""))
	
	# Random Samples from the Distribution
	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 )
	}	
	res =  rbind(preres, matrix(0, nrow = n, ncol = m.sim))
	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("msgarchsim", garchindx = garchindx, omega = as.numeric(omega), alpha = as.numeric(alpha), 
			beta = as.numeric(beta), h = h, z = z, res = res, e = e, vxs = vxs, PACKAGE = "rgarch")
	
	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("sGARCHpath",
			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.