R/rdcc-dcc.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.
##
#################################################################################


.dccspec = function(uspec, VAR = FALSE, VAR.opt = list(robust = FALSE, lag = 1, lag.max = NULL, 
				lag.criterion = c("AIC", "HQ", "SC", "FPE"), external.regressors = NULL, 
				robust.control = list("gamma" = 0.25, "delta" = 0.01, "nc" = 10, "ns" = 500)), 
		dccOrder = c(1,1), distribution = c("mvnorm", "mvt", "mvlaplace"),
		start.pars = list(), fixed.pars = list())
{
	if(is.null(distribution)) distribution = "mvnorm"
	distribution = distribution[1]
	valid.distributions = c("mvnorm", "mvt", "mvlaplace")
	if(!any(distribution == valid.distributions)) stop("\nInvalid Distribution Choice\n", call. = FALSE)
	
	
	if(is.null(dccOrder)) dccOrder =c(1,1)
	if(length(dccOrder)!=2) stop("\nInvalid dccOrder\n", call. = FALSE)
	dccOrder = as.integer(dccOrder)
	
	maxdccOrder = max(dccOrder)
	maxgarchOrder =  max( sapply(uspec@spec, FUN = function(x) max(max(x@mean.model$armaOrder), max(x@variance.model$garchOrder) ) ) )
	maxOrder = max( maxdccOrder, maxgarchOrder)
	
	shape.distributions = c("mvt")
	skew.distributions  = NULL
	include.skew = include.shape = FALSE
	
	# we will have to change this eventually to accomodate vector
	# skewness and shape parameters for more interesting distributions
	if(any(distribution == shape.distributions)) include.shape = TRUE
	if(any(distribution == skew.distributions)) include.skew = TRUE
	
	mspec = list()
	mspec$dccOrder = dccOrder
	mspec$start.pars = start.pars
	mspec$fixed.pars = fixed.pars
	mspec$distribution = distribution
	mspec$include.shape = include.shape
	mspec$include.skew = include.skew
	
	modelnames = c(
			if(dccOrder[1] > 0) paste("dcca",  1:dccOrder[1], sep = ""),
			if(dccOrder[2] > 0) paste("dccb",  1:dccOrder[2], sep = ""),
			if(include.skew)  paste("skew"),
			if(include.shape) paste("shape"))
	
	
	pos = 1
	pos.matrix = matrix(0, ncol = 4, nrow = 4)
	colnames(pos.matrix) = c("start", "stop", "include", "npars")
	rownames(pos.matrix) = c("dcca","dccb","skew","shape")
	
	if(dccOrder[1] > 0){
		pos.matrix[1,1:3] = c(pos, pos+dccOrder[1]-1, 1)
		pos = max(pos.matrix[1:1, 2]) + 1
	}
	if(dccOrder[2] > 0){
		pos.matrix[2,1:3] = c(pos, pos+dccOrder[2]-1, 1)
		pos = max(pos.matrix[1:2, 2]) + 1
	}
	if(include.skew)
	{
		pos.matrix[3,1:3] = c(pos, pos, 1)
		pos = max(pos.matrix[1:3, 2]) + 1
	}
	if(include.shape)
	{
		pos.matrix[4,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
	}
	
	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,
			maxgarchOrder = maxgarchOrder, maxdccOrder = maxdccOrder, 
			pos.matrix = pos.matrix)
	
	mspec$optimization.model = omodel
	
	if( !is(uspec, "uGARCHmultispec") ) stop("\ndccspec-->error: uspec must be a uGARCHmultispec object")
	
	mspec$varmodel = list()
	if( VAR ){
		
		mspec$varmodel$VAR = TRUE
		mspec$varmodel$robust = VAR.opt$robust	
		mspec$varmodel$lag = VAR.opt$lag
		mspec$varmodel$lag.max = VAR.opt$lag.max
		mspec$varmodel$lag.criterion = VAR.opt$lag.criterion
		mspec$varmodel$external.regressors = VAR.opt$external.regressors
		mspec$varmodel$robust.control = VAR.opt$robust.control		
		if( !is.null(mspec$varmodel$external.regressors) ){
			mspec$varmodel$mxn = dim( as.matrix(VAR.opt$external.regressors) )[2]
		} else{
			mspec$varmodel$mxn = 0
		}
		# loop and reconstruct the spec to exclude a mean specification
		nl = length(uspec@spec)
		tmpspec = vector(mode ="list", length = nl)
		for(i in 1:nl){
			tmpspec[[i]] = ugarchspec(mean.model = list(include.mean = FALSE, armaOrder = c(0,0),
							external.regressors = NULL, arfima = FALSE,  garchInMean = FALSE),
					variance.model = list(model = uspec@spec[[i]]@variance.model$model,
							submodel = uspec@spec[[i]]@variance.model$submodel,
							garchOrder = uspec@spec[[i]]@variance.model$garchOrder,
							external.regressors = uspec@spec[[i]]@variance.model$vexdata,
							variance.targeting = !uspec@spec[[i]]@variance.model$include.omega),
					distribution.model = uspec@spec[[i]]@distribution.model$distribution,
					start.pars = uspec@spec[[i]]@optimization.model$start.pars,
					fixed.pars = uspec@spec[[i]]@optimization.model$fixed.pars)
		}
		uspec = multispec( tmpspec )
	} else{
		mspec$varmodel$VAR = FALSE
		mspec$varmodel$lag = 1
		mspec$varmodel$lag.max = 1
		mspec$varmodel$lag.criterion = "HQ"
		mspec$varmodel$external.regressors = NULL
		mspec$varmodel$mxn = 0
	}
	
	
	ans = new("DCCspec",
			mspec = mspec,
			uspec = uspec)
	return(ans)
}

dccfit.norm = function(spec, data, out.sample = 0, solver = "solnp", solver.control = list(), 
		fit.control = list(eval.se = TRUE, stationarity = TRUE, scale = TRUE), 
		parallel = FALSE, parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2), 
		fit = NULL, VAR.fit = NULL, ...)
{
	# we only do 2-stage estimation
	# first get the univariate fits
	tic = Sys.time()
	ufit.control = list()
	if(is.null(fit.control$stationarity)){
		ufit.control$stationarity = TRUE 
	} else {
		ufit.control$stationarity = fit.control$stationarity
		fit.control$stationarity = NULL
	}
	if(is.null(fit.control$scale)){
		ufit.control$scale = TRUE 
	} else{
		ufit.control$scale = fit.control$scale
		fit.control$scale = NULL
	}
		
	if(is.null(fit.control$eval.se)) fit.control$eval.se = TRUE
	##################
	# Parallel Execution Prelim Checks (we need them because of the hessian estimation)
	if( parallel ){
		os = .Platform$OS.type
		if(is.null(parallel.control$pkg)){
			if( os == "windows" ) parallel.control$pkg = "snowfall" else parallel.control$pkg = "multicore"
			if( is.null(parallel.control$cores) ) parallel.control$cores = 2
		} else{
			mtype = match(tolower(parallel.control$pkg[1]), c("multicore", "snowfall"))
			if(is.na(mtype)) stop("\nParallel Package type not recognized in parallel.control\n")
			parallel.control$pkg = tolower(parallel.control$pkg[1])
			if( os == "windows" && parallel.control$pkg == "multicore" ) stop("\nmulticore not supported on windows O/S\n")
			if( is.null(parallel.control$cores) ) parallel.control$cores = 2 else parallel.control$cores = as.integer(parallel.control$cores[1])
		}
	}
	
	if( is.null( colnames(data) ) ) cnames = paste("Asset_", 1:k, sep = "") else cnames = colnames(data)
	xdata = .extractmdata(data)
	if(!is.numeric(out.sample)) stop("\ndccfit-->error: out.sample must be numeric\n")
	if(as.numeric(out.sample) < 0) stop("\ndccfit-->error: out.sample must be positive\n")
	n.start = round(out.sample, 0)
	
	n = dim(xdata$data)[1]
	k = dim(xdata$data)[2]
		
	if( (n-n.start) < 100) stop("\ndccfit-->error: function requires at least 100 data\n points to run\n")
	data 		= xdata$data[1:(n-n.start), , drop = FALSE]
	dates 		= xdata$pos[1:(n-n.start), drop = FALSE]
	origdata  	= xdata$data
	origdates 	= xdata$pos
	dformat 	= xdata$dformat
	
	
	vrmodel = NULL
	if( spec@mspec$varmodel$VAR ){
		if(spec@mspec$varmodel$mxn>0){
			ex = spec@mspec$varmodel$external.regressors
			if( dim(ex)[1] < dim(data)[1]  ) stop("\ndccfit-->error: external regressor matrix has less points than data!...", call. = FALSE)
			orex = ex
			ex = ex[1:(n - n.start), , drop = FALSE]
			mxn = dim(orex)[2]
		} else{
			orex = NULL
			ex = NULL
			mxn = 0
		}
		if( !is.null(VAR.fit) ){
			zdata = VAR.fit$xresiduals
			# should be N - p
			p = mp = VAR.fit$lag
			N = dim(zdata)[1]
			K = dim(zdata)[2]
			mu = VAR.fit$xfitted
			vrmodel = VAR.fit
			#spec@mspec$optimization.model$maxOrder = max( mp, spec@mspec$optimization.model$maxdccOrder )
			#spec@mspec$optimization.model$maxgarchOrder  = mp
		} else{

			cat("\nestimating VAR model...")

			ic = spec@mspec$varmodel$lag.criterion[1]
			if( !is.null(spec@mspec$varmodel$lag.max) ){
				mp = .varxselect(y = data, lag.max = spec@mspec$varmodel$lag.max, exogen = ex)$selection
				mp = as.numeric( mp[which(substr(names(mp), 1, 2) == substr(ic, 1, 2))] )
			} else {
				mp = spec@mspec$varmodel$lag
			}
			vrmodel = varxfilter(X = data, p = mp, exogen = ex, Bcoef = NULL, robust = spec@mspec$varmodel$robust, 
					gamma = spec@mspec$varmodel$robust.control$gamma, delta = spec@mspec$varmodel$robust.control$delta, 
					nc = spec@mspec$varmodel$robust.control$nc, ns = spec@mspec$varmodel$robust.control$ns)
			cat("done!\n")
			zdata = vrmodel$xresiduals
			# should be N - p
			N = dim(zdata)[1]
			K = dim(zdata)[2]
			mu = vrmodel$xfitted
			vrmodel$lag = mp
			vrmodel$mxn = mxn
			#spec@mspec$optimization.model$maxOrder = max( mp, spec@mspec$optimization.model$maxdccOrder )
			#spec@mspec$optimization.model$maxgarchOrder  = mp
		}
		# if we have out of sample data we need to filter
		if(out.sample>0){
			Bcoef = vrmodel$Bcoef
			vrmodel2 = varxfilter(X = origdata, p = mp, exogen = orex, Bcoef = Bcoef)
			zorigdata = vrmodel2$xresiduals			
		} else{
			zorigdata = zdata
		}
	} else{
		zdata = data
		zorigdata = origdata
		orex = NULL
		ex = NULL
	}
	
	if( is(spec@uspec, "uGARCHmultispec") ){
		speclist = spec@uspec
	} else{
		stop("\ndccfit-->error: DCCspec must contain multispec object with fixed parameters")
	}
	maxgarchOrder = spec@mspec$optimization.model$maxgarchOrder
	
	tmpidx  = .dccparidx( speclist )
	paridx  = tmpidx$idx
	paridxi = tmpidx$idxi
	paridxa = tmpidx$idxa
	if( !is.null(fit) && is(fit, "uGARCHmultifit") ){
		fitlist = fit
	} else{
		fitlist = multifit(multispec = speclist, data = zorigdata, out.sample = n.start, solver = solver, 
				solver.control = solver.control, fit.control = ufit.control, parallel = parallel, parallel.control = parallel.control)
		assign(".fitlist", fitlist, envir = .GlobalEnv)
	}
	garchpars = as.numeric( unlist(coef(fitlist) ) )
	
	converge = sapply(fitlist@fit, FUN = function(x) x@fit$convergence)
	if( any( converge == 1 ) ){
		pr = which(converge != 1)
		cat("\nNon-Converged:\n")
		print(pr)
		cat("\ndccfit-->error: convergence problem in univariate fit...")
		cat("\n...returning uGARCHmultifit object instead...check and resubmit...")
		return( fitlist )
	}
	
	# create a temporary environment to store values (deleted at end of function)
	tempenvir = new.env(hash = TRUE, parent = .GlobalEnv)
	
	assign("parallel", parallel, envir = tempenvir)
	assign("parallel.control", parallel.control, envir = tempenvir)
	assign("eval.se", fit.control$eval.se, envir = tempenvir)
	
	assign("paridx", paridx, envir = tempenvir)
	assign("paridxi", paridxi, envir = tempenvir)
	assign("paridxa", paridxa, envir = tempenvir)
	
	assign("speclist", speclist, envir = tempenvir)
	assign("fitlist", fitlist, envir = tempenvir)
	
	omodel = spec@mspec$optimization.model
	
	sig = sigma(fitlist)
	res = residuals(fitlist)
	stdresid = res/sig
	if( maxgarchOrder > 0 )  stdresid = stdresid[-(1:maxgarchOrder), , drop = FALSE]
	
	Qbar = cov(stdresid)
	H = sig^2
	assign("omodel", omodel, envir = tempenvir)
	assign("mspec", spec@mspec, envir = tempenvir)
	assign("stdresid", stdresid,  envir = tempenvir)
	assign("Qbar", Qbar,  envir = tempenvir)
	assign("solver", solver,  envir = tempenvir)
	assign("fit.control", fit.control, envir = tempenvir)
	assign("H", H, envir = tempenvir)
	
	Ifn = .dcccon
	ILB = 0.0001
	IUB = 0.9998
	if( solver == "solnp" ) fit.control$stationarity = FALSE else fit.control$stationarity = TRUE
	assign("fit.control", fit.control, envir = tempenvir)
	
	# DCC parameters
	ipars = .dccstart(data = zdata, garchenv = tempenvir)
	pars  = ipars$pars
	LB 	= ipars$modelLB
	UB 	= ipars$modelUB
	assign("npar", length(pars), envir = tempenvir)
	
	if( any( LB == UB ) )
	{
		fixid = which(LB == UB)
		if( length(fixid) == length(pars) )
		{
			if( !fit.control$eval.se )
			{
				# if all parameters are fixed an no standard erros are to
				# be calculated then we return a ugarchfilter object
				cat("\nrdccfit-->warning: all parameters fixed...returning dccfilter object instead\n")
				return(dccfilter(spec = spec, data = zdata, out.sample = out.sample, 
								parallel = parallel, parallel.control = parallel.control, VAR.fit = VAR.fit))
			} else{
				# if all parameters are fixed but we require standard errors, we
				# skip the solver
				use.solver 	= FALSE
				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 	= TRUE
			fixed 		= pars[fixid]
			pars		= pars[-fixid]
			LB 			= LB[-fixid]
			UB 			= UB[-fixid]
		} 
	} else{
		# no fixed parameters, all normal
		fixed = NULL
		fixid = NULL
		use.solver = TRUE
	}
	assign("fixed", fixed, envir = tempenvir)
	assign("fixid", fixid, envir = tempenvir)	
	assign(".llh", 1, envir = tempenvir)
	
	# get
	if( use.solver )
	{
		solution = .dccsolver(solver, pars, fun = normal.dccLLH1, Ifn, ILB, 
		IUB, gr = NULL, hessian = NULL, control = solver.control, 
		LB = LB, UB = UB, data = zdata, returnType = "llh", garchenv = tempenvir)
		sol 	= solution$sol
		hess 	= solution$hess
		opt.pars = sol$par
		names(opt.pars) = names(pars)
		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()
	model$vrmodel = vrmodel
	model$asset.names = cnames
	model$distribution = spec@mspec$distribution
	model$include.skew = spec@mspec$include.skew
	model$include.shape = spec@mspec$include.shape
	model$dccOrder = spec@mspec$dccOrder
	model$out.sample = out.sample
	model$pos.matrix = spec@mspec$optimization.model$pos.matrix
	model$model.matrix = spec@mspec$optimization.model$modelmatrix
	fit$spec = spec
	# check convergence else write message/return
	if( convergence == 0 ){
		if( !is.null(fixed) ){
			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 = .dccmakefitmodel(garchmodel = "dccnorm", f = normal.dccLLH2, 
				pars = opt.pars, garchpars, zdata, garchenv = tempenvir,
				timer = 0, message = sol$message, fname = "normal.dccLLH2")
		fit$dates = dates
		fit$date.format = dformat
		fit$origdata = origdata
		fit$origdates = origdates
		fit$timer	  = Sys.time() - tic
	} else{
		fit$message = sol$message
		fit$convergence = 1
	}
	fit$model = model
	# make model list to return some usefule information which
	
	ans = new("DCCfit",
			mfit = fit,
			ufit = fitlist)
	rm(tempenvir)
	return(ans)
}

dccfit.student = function(spec, data, out.sample = 0, solver = "solnp", solver.control = list(), 
		fit.control = list(eval.se = TRUE, stationarity = TRUE, scale = FALSE), parallel = FALSE, 
		parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2), fit = NULL, VAR.fit = NULL, ...)
{
	# we only do 2-stage estimation
	# first get the univariate fits
	tic = Sys.time()
	ufit.control = list()
	if(is.null(fit.control$stationarity)){
		ufit.control$stationarity = TRUE 
	} else {
		ufit.control$stationarity = fit.control$stationarity
		fit.control$stationarity = NULL
	}
	if(is.null(fit.control$scale)){
		ufit.control$scale = TRUE 
	} else{
		ufit.control$scale = fit.control$scale
		fit.control$scale = NULL
	}
	
	if(is.null(fit.control$eval.se)) 		fit.control$eval.se = TRUE
	##################
	# Parallel Execution Prelim Checks (we need them because of the hessian estimation)
	if( parallel ){
		os = .Platform$OS.type
		if(is.null(parallel.control$pkg)){
			if( os == "windows" ) parallel.control$pkg = "snowfall" else parallel.control$pkg = "multicore"
			if( is.null(parallel.control$cores) ) parallel.control$cores = 2
		} else{
			mtype = match(tolower(parallel.control$pkg[1]), c("multicore", "snowfall"))
			if(is.na(mtype)) stop("\nParallel Package type not recognized in parallel.control\n")
			parallel.control$pkg = tolower(parallel.control$pkg[1])
			if( os == "windows" && parallel.control$pkg == "multicore" ) stop("\nmulticore not supported on windows O/S\n")
			if( is.null(parallel.control$cores) ) parallel.control$cores = 2 else parallel.control$cores = as.integer(parallel.control$cores[1])
		}
	}
	
	if( is.null( colnames(data) ) ) cnames = paste("Asset_", 1:k, sep = "") else cnames = colnames(data)
	xdata = .extractmdata(data)
	if(!is.numeric(out.sample)) stop("\ndccfit-->error: out.sample must be numeric\n")
	if(as.numeric(out.sample) < 0) stop("\ndccfit-->error: out.sample must be positive\n")
	n.start = round(out.sample, 0)
	
	n = dim(xdata$data)[1]
	k = dim(xdata$data)[2]
	
	if( (n-n.start) < 100) stop("\nrdcc-->error: function requires at least 100 data\n points to run\n")
	data 		= xdata$data[1:(n-n.start), , drop = FALSE]
	dates 		= xdata$pos[1:(n-n.start), drop = FALSE]
	origdata  	= xdata$data
	origdates 	= xdata$pos
	dformat 	= xdata$dformat
	
	vrmodel = NULL
	if( spec@mspec$varmodel$VAR ){
		if(spec@mspec$varmodel$mxn>0){
			ex = spec@mspec$varmodel$external.regressors
			if( dim(ex)[1] < dim(data)[1]  ) stop("\ndccfit-->error: external regressor matrix has less points than data!...", call. = FALSE)
			orex = ex
			ex = ex[1:(n - n.start), , drop = FALSE]
			mxn = dim(orex)[2]
		} else{
			orex = NULL
			ex = NULL
			mxn = 0
		}
		if( !is.null(VAR.fit) ){
			zdata = VAR.fit$xresiduals
			# should be N - p
			p = mp = VAR.fit$lag
			N = dim(zdata)[1]
			K = dim(zdata)[2]
			mu = VAR.fit$xfitted
			vrmodel = VAR.fit
			#spec@mspec$optimization.model$maxOrder = max( mp, spec@mspec$optimization.model$maxdccOrder )
			#spec@mspec$optimization.model$maxgarchOrder  = mp
		} else{
			
			cat("\nestimating VAR model...")
			
			ic = spec@mspec$varmodel$lag.criterion[1]
			if( !is.null(spec@mspec$varmodel$lag.max) ){
				mp = .varxselect(y = data, lag.max = spec@mspec$varmodel$lag.max, exogen = ex)$selection
				mp = as.numeric( mp[which(substr(names(mp), 1, 2) == substr(ic, 1, 2))] )
			} else {
				mp = spec@mspec$varmodel$lag
			}
			vrmodel = varxfilter(X = data, p = mp, exogen = ex, Bcoef = NULL, robust = spec@mspec$varmodel$robust, 
					gamma = spec@mspec$varmodel$robust.control$gamma, delta = spec@mspec$varmodel$robust.control$delta, 
					nc = spec@mspec$varmodel$robust.control$nc, ns = spec@mspec$varmodel$robust.control$ns)
			cat("done!\n")
			zdata = vrmodel$xresiduals
			# should be N - p
			N = dim(zdata)[1]
			K = dim(zdata)[2]
			mu = vrmodel$xfitted
			vrmodel$lag = mp
			vrmodel$mxn = mxn
			#spec@mspec$optimization.model$maxOrder = max( mp, spec@mspec$optimization.model$maxdccOrder )
			#spec@mspec$optimization.model$maxgarchOrder  = mp
		}
		# if we have out of sample data we need to filter
		if(out.sample>0){
			Bcoef = vrmodel$Bcoef
			vrmodel2 = varxfilter(X = origdata, p = mp, exogen = orex, Bcoef = Bcoef)
			zorigdata = vrmodel2$xresiduals			
		} else{
			zorigdata = zdata
		}
	} else{
		zdata = data
		zorigdata = origdata
		orex = NULL
		ex = NULL
	}
	
	if( is(spec@uspec, "uGARCHmultispec") ){
		speclist = spec@uspec
	} else{
		stop("\ndccfit-->error: DCCspec must contain multispec object with fixed parameters")
	}	
	maxgarchOrder = spec@mspec$optimization.model$maxgarchOrder

	
	if( !is.null(fit) && is(fit, "uGARCHmultifit") ){
		fitlist = fit
	} else{
		fitlist = multifit(multispec = speclist, data = zorigdata, out.sample = n.start, solver = solver, 
				solver.control = solver.control, fit.control = ufit.control, parallel = parallel, parallel.control = parallel.control)
		assign(".fitlist", fitlist, envir = .GlobalEnv)
	}
	garchpars = as.numeric( unlist(coef(fitlist) ) )
	
	
	tmpidx  = .dccparidx( speclist )
	paridx  = tmpidx$idx
	paridxi = tmpidx$idxi
	paridxa = tmpidx$idxa
	converge = sapply(fitlist@fit, FUN = function(x) x@fit$convergence)
	if( any( converge == 1 ) ){
		pr = which(converge != 1)
		cat("\nNon-Converged:\n")
		print(pr)
		cat("\ndccfit-->error: convergence problem in univariate fit...")
		cat("\n...returning uGARCHmultifit object instead...check and resubmit...")
		return( fitlist )
	}
	
	# create a temporary environment to store values (deleted at end of function)
	tempenvir = new.env(hash = TRUE, parent = .GlobalEnv)
	
	assign("parallel", parallel, envir = tempenvir)
	assign("parallel.control", parallel.control, envir = tempenvir)
	assign("eval.se", fit.control$eval.se, envir = tempenvir)
	
	assign("paridx", paridx, envir = tempenvir)
	assign("paridxi", paridxi, envir = tempenvir)
	assign("paridxa", paridxa, envir = tempenvir)
	
	assign("speclist", speclist, envir = tempenvir)
	assign("fitlist", fitlist, envir = tempenvir)
	
	omodel = spec@mspec$optimization.model
	
	sig = sigma(fitlist)
	res = residuals(fitlist)
	stdresid = res/sig
	if(maxgarchOrder > 0 )  stdresid = stdresid[-(1:maxgarchOrder), ]
	
	Qbar = cov(stdresid)
	H = sig^2
	assign("omodel", omodel, envir = tempenvir)
	assign("mspec", spec@mspec, envir = tempenvir)
	assign("stdresid", stdresid,  envir = tempenvir)
	assign("Qbar", Qbar,  envir = tempenvir)
	assign("solver", solver,  envir = tempenvir)
	assign("fit.control", fit.control, envir = tempenvir)
	assign("H", H, envir = tempenvir)
	
	Ifn = .dcccon
	ILB = 0.0001
	IUB = 0.9998
	if(solver == "solnp") fit.control$stationarity = FALSE else fit.control$stationarity = TRUE
	assign("fit.control", fit.control, envir = tempenvir)
	
	# DCC parameters
	ipars = .dccstart(data = zdata, garchenv = tempenvir)
	pars  = ipars$pars
	LB 	= ipars$modelLB
	UB 	= ipars$modelUB
	assign("npar", length(pars), envir = tempenvir)
	
	if( any( LB == UB ) )
	{
		fixid = which(LB == UB)
		if( length(fixid) == length(pars) )
		{
			if( !fit.control$eval.se )
			{
				# if all parameters are fixed an no standard erros are to
				# be calculated then we return a ugarchfilter object
				cat("\nrdccfit-->warning: all parameters fixed...returning dccfilter object instead\n")
				return(dccfilter(spec = spec, data = zdata, out.sample = out.sample, 
								parallel = parallel, parallel.control = parallel.control, VAR.fit = VAR.fit))
			} else{
				# if all parameters are fixed but we require standard errors, we
				# skip the solver
				use.solver 	= FALSE
				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 	= TRUE
			fixed 		= pars[fixid]
			pars		= pars[-fixid]
			LB 			= LB[-fixid]
			UB 			= UB[-fixid]
		} 
	} else{
		# no fixed parameters, all normal
		fixed = NULL
		fixid = NULL
		use.solver = TRUE
	}
	assign("fixed", fixed, envir = tempenvir)
	assign("fixid", fixid, envir = tempenvir)	
	assign(".llh", 1, envir = tempenvir)
	
	# get
	if( use.solver )
	{
		solution = .dccsolver(solver, pars, fun = student.dccLLH1, Ifn, ILB, 
				IUB, gr = NULL, hessian = NULL, control = solver.control, 
				LB = LB, UB = UB, data = zdata, returnType = "llh", garchenv = tempenvir)
		sol 	= solution$sol
		hess 	= solution$hess
		opt.pars = sol$par
		names(opt.pars) = names(pars)
		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()
	model$vrmodel = vrmodel
	model$asset.names = cnames
	model$distribution = spec@mspec$distribution
	model$include.skew = spec@mspec$include.skew
	model$include.shape = spec@mspec$include.shape
	model$dccOrder = spec@mspec$dccOrder
	model$out.sample = out.sample
	model$pos.matrix = spec@mspec$optimization.model$pos.matrix
	model$model.matrix = spec@mspec$optimization.model$modelmatrix
	fit$spec = spec
	# check convergence else write message/return
	if( convergence == 0 ){
		if( !is.null(fixed) ){
			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 = .dccmakefitmodel(garchmodel = "dccstudent", f = student.dccLLH2, 
				pars = opt.pars, garchpars, zdata, garchenv = tempenvir,
				timer = 0, message = sol$message, fname = "student.dccLLH2")
		fit$dates = dates
		fit$date.format = dformat
		fit$origdata = origdata
		fit$origdates = origdates
		fit$timer	  = Sys.time() - tic
	} else{
		fit$message = sol$message
		fit$convergence = 1
	}
	fit$model = model
	# make model list to return some usefule information which
	
	ans = new("DCCfit",
			mfit = fit,
			ufit = fitlist)
	rm(tempenvir)
	return(ans)
}

dccfit.laplace = function(spec, data, out.sample = 0, solver = "solnp", solver.control = list(), 
		fit.control = list(eval.se = TRUE, stationarity = TRUE, scale = FALSE), 
		parallel = FALSE, parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2), fit = NULL, VAR.fit = NULL, ...)
{
	# we only do 2-stage estimation
	# first get the univariate fits
	tic = Sys.time()
	ufit.control = list()
	if(is.null(fit.control$stationarity)){
		ufit.control$stationarity = TRUE 
	} else {
		ufit.control$stationarity = fit.control$stationarity
		fit.control$stationarity = NULL
	}
	if(is.null(fit.control$scale)){
		ufit.control$scale = TRUE 
	} else{
		ufit.control$scale = fit.control$scale
		fit.control$scale = NULL
	}
	
	if(is.null(fit.control$eval.se)) 		fit.control$eval.se = TRUE
	##################
	# Parallel Execution Prelim Checks (we need them because of the hessian estimation)
	if( parallel ){
		os = .Platform$OS.type
		if(is.null(parallel.control$pkg)){
			if( os == "windows" ) parallel.control$pkg = "snowfall" else parallel.control$pkg = "multicore"
			if( is.null(parallel.control$cores) ) parallel.control$cores = 2
		} else{
			mtype = match(tolower(parallel.control$pkg[1]), c("multicore", "snowfall"))
			if(is.na(mtype)) stop("\nParallel Package type not recognized in parallel.control\n")
			parallel.control$pkg = tolower(parallel.control$pkg[1])
			if( os == "windows" && parallel.control$pkg == "multicore" ) stop("\nmulticore not supported on windows O/S\n")
			if( is.null(parallel.control$cores) ) parallel.control$cores = 2 else parallel.control$cores = as.integer(parallel.control$cores[1])
		}
	}
	
	if( is.null( colnames(data) ) ) cnames = paste("Asset_", 1:k, sep = "") else cnames = colnames(data)
	xdata = .extractmdata(data)
	if(!is.numeric(out.sample)) stop("\ndccfit-->error: out.sample must be numeric\n")
	if(as.numeric(out.sample) < 0) stop("\ndccfit-->error: out.sample must be positive\n")
	n.start = round(out.sample, 0)
	
	n = dim(xdata$data)[1]
	k = dim(xdata$data)[2]
	
	if( (n-n.start) < 100) stop("\nrdcc-->error: function requires at least 100 data\n points to run\n")
	data 		= xdata$data[1:(n-n.start), , drop = FALSE]
	dates 		= xdata$pos[1:(n-n.start), drop = FALSE]
	origdata  	= xdata$data
	origdates 	= xdata$pos
	dformat 	= xdata$dformat
	
	vrmodel = NULL
	if( spec@mspec$varmodel$VAR ){
		if(spec@mspec$varmodel$mxn>0){
			ex = spec@mspec$varmodel$external.regressors
			if( dim(ex)[1] < dim(data)[1]  ) stop("\ndccfit-->error: external regressor matrix has less points than data!...", call. = FALSE)
			orex = ex
			ex = ex[1:(n - n.start), , drop = FALSE]
			mxn = dim(orex)[2]
		} else{
			orex = NULL
			ex = NULL
			mxn = 0
		}
		if( !is.null(VAR.fit) ){
			zdata = VAR.fit$xresiduals
			# should be N - p
			p = mp = VAR.fit$lag
			N = dim(zdata)[1]
			K = dim(zdata)[2]
			mu = VAR.fit$xfitted
			vrmodel = VAR.fit
			
			#spec@mspec$optimization.model$maxOrder = max( mp, spec@mspec$optimization.model$maxdccOrder )
			#spec@mspec$optimization.model$maxgarchOrder  = mp
		} else{
			
			cat("\nestimating VAR model...")
			
			ic = spec@mspec$varmodel$lag.criterion[1]
			if( !is.null(spec@mspec$varmodel$lag.max) ){
				mp = .varxselect(y = data, lag.max = spec@mspec$varmodel$lag.max, exogen = ex)$selection
				mp = as.numeric( mp[which(substr(names(mp), 1, 2) == substr(ic, 1, 2))] )
			} else {
				mp = spec@mspec$varmodel$lag
			}
			vrmodel = varxfilter(X = data, p = mp, exogen = ex, Bcoef = NULL, robust = spec@mspec$varmodel$robust, 
					gamma = spec@mspec$varmodel$robust.control$gamma, delta = spec@mspec$varmodel$robust.control$delta, 
					nc = spec@mspec$varmodel$robust.control$nc, ns = spec@mspec$varmodel$robust.control$ns)
			cat("done!\n")
			zdata = vrmodel$xresiduals
			# should be N - p
			N = dim(zdata)[1]
			K = dim(zdata)[2]
			mu = vrmodel$xfitted
			vrmodel$lag = mp
			vrmodel$mxn = mxn
			#spec@mspec$optimization.model$maxOrder = max( mp, spec@mspec$optimization.model$maxdccOrder )
			#spec@mspec$optimization.model$maxgarchOrder  = mp
		}
		# if we have out of sample data we need to filter
		if(out.sample>0){
			Bcoef = vrmodel$Bcoef
			vrmodel2 = varxfilter(X = origdata, p = mp, exogen = orex, Bcoef = Bcoef)
			zorigdata = vrmodel2$xresiduals			
		} else{
			zorigdata = zdata
		}
	} else{
		zdata = data
		zorigdata = origdata
		orex = NULL
		ex = NULL
	}
	
	if( is(spec@uspec, "uGARCHmultispec") ){
		speclist = spec@uspec
	} else{
		stop("\ndccfit-->error: DCCspec must contain multispec object with fixed parameters")
	}
	maxgarchOrder = spec@mspec$optimization.model$maxgarchOrder
	
	tmpidx  = .dccparidx( speclist )
	paridx  = tmpidx$idx
	paridxi = tmpidx$idxi
	paridxa = tmpidx$idxa
	
	if( !is.null(fit) && is(fit, "uGARCHmultifit") ){
		fitlist = fit
	} else{
		fitlist = multifit(multispec = speclist, data = zorigdata, out.sample = n.start, solver = solver, 
				solver.control = solver.control, fit.control = ufit.control, parallel = parallel, parallel.control = parallel.control)
		assign(".fitlist", fitlist, envir = .GlobalEnv)
	}
	garchpars = as.numeric( unlist(coef(fitlist) ) )
	
	converge = sapply(fitlist@fit, FUN = function(x) x@fit$convergence)
	if( any( converge == 1 ) ){
		pr = which(converge != 1)
		cat("\nNon-Converged:\n")
		print(pr)
		cat("\ndccfit-->error: convergence problem in univariate fit...")
		cat("\n...returning uGARCHmultifit object instead...check and resubmit...")
		return( fitlist )
	}
	
	# create a temporary environment to store values (deleted at end of function)
	tempenvir = new.env(hash = TRUE, parent = .GlobalEnv)
	
	assign("parallel", parallel, envir = tempenvir)
	assign("parallel.control", parallel.control, envir = tempenvir)
	assign("eval.se", fit.control$eval.se, envir = tempenvir)
	
	assign("paridx", paridx, envir = tempenvir)
	assign("paridxi", paridxi, envir = tempenvir)
	assign("paridxa", paridxa, envir = tempenvir)
	
	assign("speclist", speclist, envir = tempenvir)
	assign("fitlist", fitlist, envir = tempenvir)
	
	omodel = spec@mspec$optimization.model
	
	sig = sigma(fitlist)
	res = residuals(fitlist)
	stdresid = res/sig
	if( maxgarchOrder > 0 )  stdresid = stdresid[-(1:maxgarchOrder), , drop = FALSE]
	
	Qbar = cov(stdresid)
	H = sig^2
	assign("omodel", omodel, envir = tempenvir)
	assign("mspec", spec@mspec, envir = tempenvir)
	assign("stdresid", stdresid,  envir = tempenvir)
	assign("Qbar", Qbar,  envir = tempenvir)
	assign("solver", solver,  envir = tempenvir)
	assign("fit.control", fit.control, envir = tempenvir)
	assign("H", H, envir = tempenvir)
	
	Ifn = .dcccon
	ILB = 0.0001
	IUB = 0.9998
	if(solver == "solnp") fit.control$stationarity = FALSE else fit.control$stationarity = TRUE
	assign("fit.control", fit.control, envir = tempenvir)
	
	# DCC parameters
	ipars = .dccstart(data = zdata, garchenv = tempenvir)
	pars  = ipars$pars
	LB 	= ipars$modelLB
	UB 	= ipars$modelUB
	assign("npar", length(pars), envir = tempenvir)
	
	if( any( LB == UB ) )
	{
		fixid = which(LB == UB)
		if( length(fixid) == length(pars) )
		{
			if( !fit.control$eval.se )
			{
				# if all parameters are fixed an no standard erros are to
				# be calculated then we return a ugarchfilter object
				cat("\nrdccfit-->warning: all parameters fixed...returning dccfilter object instead\n")
				return(dccfilter(spec = spec, data = zdata, out.sample = out.sample, 
								parallel = parallel, parallel.control = parallel.control, VAR.fit = VAR.fit))
			} else{
				# if all parameters are fixed but we require standard errors, we
				# skip the solver
				use.solver 	= FALSE
				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 	= TRUE
			fixed 		= pars[fixid]
			pars		= pars[-fixid]
			LB 			= LB[-fixid]
			UB 			= UB[-fixid]
		} 
	} else{
		# no fixed parameters, all normal
		fixed = NULL
		fixid = NULL
		use.solver = TRUE
	}
	assign("fixed", fixed, envir = tempenvir)
	assign("fixid", fixid, envir = tempenvir)	
	assign(".llh", 1, envir = tempenvir)
	
	# get
	if( use.solver )
	{
		solution = .dccsolver(solver, pars, fun = laplace.dccLLH1, Ifn, ILB, 
				IUB, gr = NULL, hessian = NULL, control = solver.control, 
				LB = LB, UB = UB, data = zdata, returnType = "llh", garchenv = tempenvir)
		sol 	= solution$sol
		hess 	= solution$hess
		opt.pars = sol$par
		names(opt.pars) = names(pars)
		convergence = sol$convergence
		# normal.dccLLH2(opt.pars, garchpars, data, returnType = "all", garchenv = tempenvir)$llh
	} else{
		hess  = NULL
		timer = Sys.time()-tic
		opt.pars = NULL
		convergence = 0
		sol = list()
		sol$message = "all parameters fixed"
	}
	
	fit = list()
	model = list()
	model$vrmodel = vrmodel
	model$asset.names = cnames
	model$distribution = spec@mspec$distribution
	model$include.skew = spec@mspec$include.skew
	model$include.shape = spec@mspec$include.shape
	model$dccOrder = spec@mspec$dccOrder
	model$out.sample = out.sample
	model$pos.matrix = spec@mspec$optimization.model$pos.matrix
	model$model.matrix = spec@mspec$optimization.model$modelmatrix
	fit$spec = spec
	# check convergence else write message/return
	if( convergence == 0 ){
		if( !is.null(fixed) ){
			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 = .dccmakefitmodel(garchmodel = "dcclaplace", f = laplace.dccLLH2, 
				pars = opt.pars, garchpars, zdata, garchenv = tempenvir,
				timer = 0, message = sol$message, fname = "laplace.dccLLH2")
		fit$dates = dates
		fit$date.format = dformat
		fit$origdata = origdata
		fit$origdates = origdates
		fit$timer	  = Sys.time() - tic
	} else{
		fit$message = sol$message
		fit$convergence = 1
	}
	fit$model = model
	# make model list to return some usefule information which
	
	ans = new("DCCfit",
			mfit = fit,
			ufit = fitlist)
	rm(tempenvir)
	return(ans)
}

dccfilter.norm = function(spec, data, out.sample = 0, filter.control = list(n.old = NULL), 
		parallel = FALSE, parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2), VAR.fit = NULL, ...)
{
	tic = Sys.time()
	
	if( parallel ){
		os = .Platform$OS.type
		if(is.null(parallel.control$pkg)){
			if( os == "windows" ) parallel.control$pkg = "snowfall" else parallel.control$pkg = "multicore"
			if( is.null(parallel.control$cores) ) parallel.control$cores = 2
		} else{
			mtype = match(tolower(parallel.control$pkg[1]), c("multicore", "snowfall"))
			if(is.na(mtype)) stop("\nParallel Package type not recognized in parallel.control\n")
			parallel.control$pkg = tolower(parallel.control$pkg[1])
			if( os == "windows" && parallel.control$pkg == "multicore" ) stop("\nmulticore not supported on windows O/S\n")
			if( is.null(parallel.control$cores) ) parallel.control$cores = 2 else parallel.control$cores = as.integer(parallel.control$cores[1])
		}
	}
	n.old = filter.control$n.old
	
	if( is.null( colnames(data) ) ) cnames = paste("Asset_", 1:k, sep = "") else cnames = colnames(data)
	xdata = .extractmdata(data)
	if(!is.numeric(out.sample)) stop("\ndccfilter-->error: out.sample must be numeric\n")
	if(as.numeric(out.sample) < 0) stop("\ndccfilter-->error: out.sample must be positive\n")
	n.start = round(out.sample, 0)
	
	n = dim(xdata$data)[1]
	k = dim(xdata$data)[2]
	
	if( (n - n.start) < 100) stop("\ndccfilter-->error: function requires at least 100 data\n points to run\n")
	data 		= xdata$data[1:(n - n.start), , drop = FALSE]
	dates 		= xdata$pos[1:(n - n.start)]
	origdata  	= xdata$data
	origdates 	= xdata$pos
	dformat 	= xdata$dformat
	
	vrmodel = NULL
	if( spec@mspec$varmodel$VAR ){
		if(spec@mspec$varmodel$mxn>0){
			ex = spec@mspec$varmodel$external.regressors
			if( dim(ex)[1] < dim(data)[1]  ) stop("\ndccfilter-->error: external regressor matrix has less points than data!...", call. = FALSE)
			orex = ex
			ex = ex[1:(n - n.start), , drop = FALSE]
			mxn = dim(orex)[2]
		} else{
			orex = NULL
			ex = NULL
			mxn = 0
		}
		if( !is.null(VAR.fit) ){
			zdata = VAR.fit$xresiduals
			# should be N - p
			p = mp = VAR.fit$lag
			N = dim(zdata)[1]
			K = dim(zdata)[2]
			mu = VAR.fit$xfitted
			vrmodel = VAR.fit
			
			
			if(out.sample>0){
				Bcoef = VAR.fit$Bcoef
				vrmodel2 = varxfilter(X = origdata, p = mp, exogen = orex, Bcoef = Bcoef)
				zorigdata = vrmodel2$xresiduals			
			} else{
				zorigdata = zdata
			}
			
		} else{
			
			cat("\nestimating VAR model...")
			
			ic = spec@mspec$varmodel$lag.criterion[1]
			if( !is.null(spec@mspec$varmodel$lag.max) ){
				mp = .varxselect(y = data, lag.max = spec@mspec$varmodel$lag.max, exogen = ex)$selection
				mp = as.numeric( mp[which(substr(names(mp), 1, 2) == substr(ic, 1, 2))] )
			} else {
				mp = spec@mspec$varmodel$lag
			}
			vrmodel = varxfilter(X = data, p = mp, exogen = ex, Bcoef = NULL, robust = spec@mspec$varmodel$robust, 
					gamma = spec@mspec$varmodel$robust.control$gamma, delta = spec@mspec$varmodel$robust.control$delta, 
					nc = spec@mspec$varmodel$robust.control$nc, ns = spec@mspec$varmodel$robust.control$ns)
			cat("done!\n")
			zdata = vrmodel$xresiduals
			# should be N - p
			N = dim(zdata)[1]
			K = dim(zdata)[2]
			mu = vrmodel$xfitted
			vrmodel$lag = mp
			vrmodel$mxn = mxn
			# if we have out of sample data we need to filter
			if(out.sample>0){
				Bcoef = vrmodel$Bcoef
				vrmodel2 = varxfilter(X = origdata, p = mp, exogen = orex, Bcoef = Bcoef)
				zorigdata = vrmodel2$xresiduals			
			} else{
				zorigdata = zdata
			}
		}
	} else{
		zdata = data
		zorigdata = origdata
		orex = NULL
		ex = NULL
	}
	
	if( is(spec@uspec, "uGARCHmultispec") ){
		speclist = spec@uspec
	} else{
		stop("\ndccfilter-->error: DCCspec must contain multispec object with fixed parameters")
	}
	maxgarchOrder = spec@mspec$optimization.model$maxgarchOrder
	
	filterlist = multifilter(multifitORspec = spec@uspec, data = as.data.frame(zorigdata), out.sample = out.sample, 
			parallel = parallel, parallel.control = parallel.control, n.old = n.old)
	
	
	tmpidx  = .dccparidx( spec@uspec )
	paridx  = tmpidx$idx
	paridxi = tmpidx$idxi
	paridxa = tmpidx$idxa
	
	garchpars = as.numeric( unlist(coef(filterlist) ) )
	
	tempenvir = new.env(hash = TRUE, parent = .GlobalEnv)
	fit.control = list(stationarity = TRUE, eval.se = FALSE)
	
	assign("fixed", NULL, envir = tempenvir)
	assign("fixid", NULL, envir = tempenvir)	
	assign("parallel", parallel, envir = tempenvir)
	assign("parallel.control", parallel.control, envir = tempenvir)
	assign("eval.se", FALSE, envir = tempenvir)
	assign("mspec", spec@mspec, envir = tempenvir)
	
	assign("paridx", paridx, envir = tempenvir)
	assign("paridxi", paridxi, envir = tempenvir)
	assign("paridxa", paridxa, envir = tempenvir)
	
	assign("speclist", speclist, envir = tempenvir)
	assign("filterlist", filterlist, envir = tempenvir)
	assign("fit.control", fit.control, envir = tempenvir)
	
	omodel = spec@mspec$optimization.model
	
	sig = sigma(filterlist)
	res = residuals(filterlist)
	stdresid = res/sig
	if(maxgarchOrder > 0 )  stdresid = stdresid[-(1:maxgarchOrder), ]
	
	Qbar = cov(stdresid)
	H = sig^2
	assign("omodel", omodel, envir = tempenvir)
	assign("mspec", spec@mspec, envir = tempenvir)
	assign("stdresid", stdresid,  envir = tempenvir)
	assign("Qbar", Qbar,  envir = tempenvir)
	assign("H", H, envir = tempenvir)
	dccpars = as.numeric( unlist(spec@mspec$optimization.model$fixed.pars) )
	
	assign("npar", length(dccpars), envir = tempenvir)
	
	filt = .dccmakefiltermodel(garchmodel = "dccnorm", f = normal.dccLLH2, 
			pars = dccpars, garchpars, zdata, garchenv = tempenvir,
			timer = 0, message = 0, fname = "normal.dccLLH2")
	
	model = list()
	model$vrmodel = vrmodel
	model$asset.names = cnames
	model$distribution = spec@mspec$distribution
	model$include.skew = spec@mspec$include.skew
	model$include.shape = spec@mspec$include.shape
	model$dccOrder = spec@mspec$dccOrder
	model$out.sample = out.sample
	
	filt$dates = dates
	filt$date.format = dformat
	filt$origdata = origdata
	filt$origdates = origdates
	filt$model = model
	filt$timer = Sys.time() - tic
	
	ans = new("DCCfilter",
			mfilter = filt,
			ufilter = filterlist)
	rm(tempenvir)
	return(ans)
	
}

dccfilter.student = function(spec, data, out.sample = 0, filter.control = list(n.old = NULL), 
		parallel = FALSE, parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2), VAR.fit = NULL, ...)
{
	tic = Sys.time()
	
	if( parallel ){
		os = .Platform$OS.type
		if(is.null(parallel.control$pkg)){
			if( os == "windows" ) parallel.control$pkg = "snowfall" else parallel.control$pkg = "multicore"
			if( is.null(parallel.control$cores) ) parallel.control$cores = 2
		} else{
			mtype = match(tolower(parallel.control$pkg[1]), c("multicore", "snowfall"))
			if(is.na(mtype)) stop("\nParallel Package type not recognized in parallel.control\n")
			parallel.control$pkg = tolower(parallel.control$pkg[1])
			if( os == "windows" && parallel.control$pkg == "multicore" ) stop("\nmulticore not supported on windows O/S\n")
			if( is.null(parallel.control$cores) ) parallel.control$cores = 2 else parallel.control$cores = as.integer(parallel.control$cores[1])
		}
	}
	n.old = filter.control$n.old
	
	if( is.null( colnames(data) ) ) cnames = paste("Asset_", 1:k, sep = "") else cnames = colnames(data)
	xdata = .extractmdata(data)
	if(!is.numeric(out.sample)) stop("\ndccfilter-->error: out.sample must be numeric\n")
	if(as.numeric(out.sample) < 0) stop("\ndccfilter-->error: out.sample must be positive\n")
	n.start = round(out.sample, 0)
	
	n = dim(xdata$data)[1]
	k = dim(xdata$data)[2]
	
	if( (n - n.start) < 100) stop("\ndccfilter-->error: function requires at least 100 data\n points to run\n")
	data 		= xdata$data[1:(n - n.start), , drop = FALSE]
	dates 		= xdata$pos[1:(n - n.start)]
	origdata  	= xdata$data
	origdates 	= xdata$pos
	dformat 	= xdata$dformat
	
	vrmodel = NULL
	if( spec@mspec$varmodel$VAR ){
		if(spec@mspec$varmodel$mxn>0){
			ex = spec@mspec$varmodel$external.regressors
			if( dim(ex)[1] < dim(data)[1]  ) stop("\ndccfilter-->error: external regressor matrix has less points than data!...", call. = FALSE)
			orex = ex
			ex = ex[1:(n - n.start), , drop = FALSE]
			mxn = dim(orex)[2]
		} else{
			orex = NULL
			ex = NULL
			mxn = 0
		}
		if( !is.null(VAR.fit) ){
			zdata = VAR.fit$xresiduals
			# should be N - p
			p = mp = VAR.fit$lag
			N = dim(zdata)[1]
			K = dim(zdata)[2]
			mu = VAR.fit$xfitted
			vrmodel = VAR.fit
			
			
			if(out.sample>0){
				Bcoef = VAR.fit$Bcoef
				vrmodel2 = varxfilter(X = origdata, p = mp, exogen = orex, Bcoef = Bcoef)
				zorigdata = vrmodel2$xresiduals			
			} else{
				zorigdata = zdata
			}
			
		} else{
			
			cat("\nestimating VAR model...")
			
			ic = spec@mspec$varmodel$lag.criterion[1]
			if( !is.null(spec@mspec$varmodel$lag.max) ){
				mp = .varxselect(y = data, lag.max = spec@mspec$varmodel$lag.max, exogen = ex)$selection
				mp = as.numeric( mp[which(substr(names(mp), 1, 2) == substr(ic, 1, 2))] )
			} else {
				mp = spec@mspec$varmodel$lag
			}
			vrmodel = varxfilter(X = data, p = mp, exogen = ex, Bcoef = NULL, robust = spec@mspec$varmodel$robust, 
					gamma = spec@mspec$varmodel$robust.control$gamma, delta = spec@mspec$varmodel$robust.control$delta, 
					nc = spec@mspec$varmodel$robust.control$nc, ns = spec@mspec$varmodel$robust.control$ns)
			cat("done!\n")
			zdata = vrmodel$xresiduals
			# should be N - p
			N = dim(zdata)[1]
			K = dim(zdata)[2]
			mu = vrmodel$xfitted
			vrmodel$lag = mp
			vrmodel$mxn = mxn
			# if we have out of sample data we need to filter
			if(out.sample>0){
				Bcoef = vrmodel$Bcoef
				vrmodel2 = varxfilter(X = origdata, p = mp, exogen = orex, Bcoef = Bcoef)
				zorigdata = vrmodel2$xresiduals			
			} else{
				zorigdata = zdata
			}
		}
	} else{
		zdata = data
		zorigdata = origdata
		orex = NULL
		ex = NULL
	}
	
	if( is(spec@uspec, "uGARCHmultispec") ){
		speclist = spec@uspec
	} else{
		stop("\ndccfilter-->error: DCCspec must contain multispec object with fixed parameters")
	}
	maxgarchOrder = spec@mspec$optimization.model$maxgarchOrder
	
	filterlist = multifilter(multifitORspec = spec@uspec, data = as.data.frame(zorigdata), out.sample = out.sample, 
			parallel = parallel, parallel.control = parallel.control, n.old = n.old)
	
	
	tmpidx  = .dccparidx( spec@uspec )
	paridx  = tmpidx$idx
	paridxi = tmpidx$idxi
	paridxa = tmpidx$idxa
	
	garchpars = as.numeric( unlist(coef(filterlist) ) )
	
	tempenvir = new.env(hash = TRUE, parent = .GlobalEnv)
	fit.control = list(stationarity = TRUE, eval.se = FALSE)
	
	assign("fixed", NULL, envir = tempenvir)
	assign("fixid", NULL, envir = tempenvir)	
	assign("parallel", parallel, envir = tempenvir)
	assign("parallel.control", parallel.control, envir = tempenvir)
	assign("eval.se", FALSE, envir = tempenvir)
	assign("mspec", spec@mspec, envir = tempenvir)
	
	assign("paridx", paridx, envir = tempenvir)
	assign("paridxi", paridxi, envir = tempenvir)
	assign("paridxa", paridxa, envir = tempenvir)
	
	assign("speclist", speclist, envir = tempenvir)
	assign("filterlist", filterlist, envir = tempenvir)
	assign("fit.control", fit.control, envir = tempenvir)
	
	omodel = spec@mspec$optimization.model
	
	sig = sigma(filterlist)
	res = residuals(filterlist)
	stdresid = res/sig
	if(maxgarchOrder > 0 )  stdresid = stdresid[-(1:maxgarchOrder), ]
	
	Qbar = cov(stdresid)
	H = sig^2
	assign("omodel", omodel, envir = tempenvir)
	assign("mspec", spec@mspec, envir = tempenvir)
	assign("stdresid", stdresid,  envir = tempenvir)
	assign("Qbar", Qbar,  envir = tempenvir)
	assign("H", H, envir = tempenvir)
	dccpars = as.numeric( unlist(spec@mspec$optimization.model$fixed.pars) )
	
	assign("npar", length(dccpars), envir = tempenvir)
	
	filt = .dccmakefiltermodel(garchmodel = "dccstudent", f = student.dccLLH2, 
			pars = dccpars, garchpars, zdata, garchenv = tempenvir,
			timer = 0, message = 0, fname = "student.dccLLH2")
	
	model = list()
	model$vrmodel = vrmodel
	model$asset.names = cnames
	model$distribution = spec@mspec$distribution
	model$include.skew = spec@mspec$include.skew
	model$include.shape = spec@mspec$include.shape
	model$dccOrder = spec@mspec$dccOrder
	model$out.sample = out.sample
	
	filt$dates = dates
	filt$date.format = dformat
	filt$origdata = origdata
	filt$origdates = origdates
	filt$model = model
	filt$timer	  = Sys.time() - tic
	
	ans = new("DCCfilter",
			mfilter = filt,
			ufilter = filterlist)
	rm(tempenvir)
	return(ans)
	
}

dccfilter.laplace = function(spec, data, out.sample = 0, filter.control = list(n.old = NULL), 
		parallel = FALSE, parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2), VAR.fit = NULL, ...)
{
	tic = Sys.time()
	
	if( parallel ){
		os = .Platform$OS.type
		if(is.null(parallel.control$pkg)){
			if( os == "windows" ) parallel.control$pkg = "snowfall" else parallel.control$pkg = "multicore"
			if( is.null(parallel.control$cores) ) parallel.control$cores = 2
		} else{
			mtype = match(tolower(parallel.control$pkg[1]), c("multicore", "snowfall"))
			if(is.na(mtype)) stop("\nParallel Package type not recognized in parallel.control\n")
			parallel.control$pkg = tolower(parallel.control$pkg[1])
			if( os == "windows" && parallel.control$pkg == "multicore" ) stop("\nmulticore not supported on windows O/S\n")
			if( is.null(parallel.control$cores) ) parallel.control$cores = 2 else parallel.control$cores = as.integer(parallel.control$cores[1])
		}
	}
	n.old = filter.control$n.old
	
	if( is.null( colnames(data) ) ) cnames = paste("Asset_", 1:k, sep = "") else cnames = colnames(data)
	xdata = .extractmdata(data)
	if(!is.numeric(out.sample)) stop("\ndccfilter-->error: out.sample must be numeric\n")
	if(as.numeric(out.sample) < 0) stop("\ndccfilter-->error: out.sample must be positive\n")
	n.start = round(out.sample, 0)
	
	n = dim(xdata$data)[1]
	k = dim(xdata$data)[2]
	
	if( (n - n.start) < 100) stop("\ndccfilter-->error: function requires at least 100 data\n points to run\n")
	data 		= xdata$data[1:(n - n.start), , drop = FALSE]
	dates 		= xdata$pos[1:(n - n.start)]
	origdata  	= xdata$data
	origdates 	= xdata$pos
	dformat 	= xdata$dformat
	
	vrmodel = NULL
	if( spec@mspec$varmodel$VAR ){
		if(spec@mspec$varmodel$mxn>0){
			ex = spec@mspec$varmodel$external.regressors
			if( dim(ex)[1] < dim(data)[1]  ) stop("\ndccfilter-->error: external regressor matrix has less points than data!...", call. = FALSE)
			orex = ex
			ex = ex[1:(n - n.start), , drop = FALSE]
			mxn = dim(orex)[2]
		} else{
			orex = NULL
			ex = NULL
			mxn = 0
		}
		if( !is.null(VAR.fit) ){
			zdata = VAR.fit$xresiduals
			# should be N - p
			p = mp = VAR.fit$lag
			N = dim(zdata)[1]
			K = dim(zdata)[2]
			mu = VAR.fit$xfitted
			vrmodel = VAR.fit
			
			
			if(out.sample>0){
				Bcoef = VAR.fit$Bcoef
				vrmodel2 = varxfilter(X = origdata, p = mp, exogen = orex, Bcoef = Bcoef)
				zorigdata = vrmodel2$xresiduals			
			} else{
				zorigdata = zdata
			}
			
		} else{
			
			cat("\nestimating VAR model...")
			
			ic = spec@mspec$varmodel$lag.criterion[1]
			if( !is.null(spec@mspec$varmodel$lag.max) ){
				mp = .varxselect(y = data, lag.max = spec@mspec$varmodel$lag.max, exogen = ex)$selection
				mp = as.numeric( mp[which(substr(names(mp), 1, 2) == substr(ic, 1, 2))] )
			} else {
				mp = spec@mspec$varmodel$lag
			}
			vrmodel = varxfilter(X = data, p = mp, exogen = ex, Bcoef = NULL, robust = spec@mspec$varmodel$robust, 
					gamma = spec@mspec$varmodel$robust.control$gamma, delta = spec@mspec$varmodel$robust.control$delta, 
					nc = spec@mspec$varmodel$robust.control$nc, ns = spec@mspec$varmodel$robust.control$ns)
			cat("done!\n")
			zdata = vrmodel$xresiduals
			# should be N - p
			N = dim(zdata)[1]
			K = dim(zdata)[2]
			mu = vrmodel$xfitted
			vrmodel$lag = mp
			vrmodel$mxn = mxn
			# if we have out of sample data we need to filter
			if(out.sample>0){
				Bcoef = vrmodel$Bcoef
				vrmodel2 = varxfilter(X = origdata, p = mp, exogen = orex, Bcoef = Bcoef)
				zorigdata = vrmodel2$xresiduals			
			} else{
				zorigdata = zdata
			}
		}
	} else{
		zdata = data
		zorigdata = origdata
		orex = NULL
		ex = NULL
	}
	
	if( is(spec@uspec, "uGARCHmultispec") ){
		speclist = spec@uspec
	} else{
		stop("\ndccfilter-->error: DCCspec must contain multispec object with fixed parameters")
	}
	maxgarchOrder = spec@mspec$optimization.model$maxgarchOrder
	
	filterlist = multifilter(multifitORspec = spec@uspec, data = as.data.frame(zorigdata), out.sample = out.sample, 
			parallel = parallel, parallel.control = parallel.control, n.old = n.old)
	
	
	tmpidx  = .dccparidx( spec@uspec )
	paridx  = tmpidx$idx
	paridxi = tmpidx$idxi
	paridxa = tmpidx$idxa
	
	garchpars = as.numeric( unlist(coef(filterlist) ) )
	
	tempenvir = new.env(hash = TRUE, parent = .GlobalEnv)
	fit.control = list(stationarity = TRUE, eval.se = FALSE)
	
	assign("fixed", NULL, envir = tempenvir)
	assign("fixid", NULL, envir = tempenvir)	
	assign("parallel", parallel, envir = tempenvir)
	assign("parallel.control", parallel.control, envir = tempenvir)
	assign("eval.se", FALSE, envir = tempenvir)
	assign("mspec", spec@mspec, envir = tempenvir)
	
	assign("paridx", paridx, envir = tempenvir)
	assign("paridxi", paridxi, envir = tempenvir)
	assign("paridxa", paridxa, envir = tempenvir)
	
	assign("speclist", speclist, envir = tempenvir)
	assign("filterlist", filterlist, envir = tempenvir)
	assign("fit.control", fit.control, envir = tempenvir)
	
	omodel = spec@mspec$optimization.model
	
	sig = sigma(filterlist)
	res = residuals(filterlist)
	stdresid = res/sig
	if(maxgarchOrder > 0 )  stdresid = stdresid[-(1:maxgarchOrder), ]
	
	Qbar = cov(stdresid)
	H = sig^2
	assign("omodel", omodel, envir = tempenvir)
	assign("mspec", spec@mspec, envir = tempenvir)
	assign("stdresid", stdresid,  envir = tempenvir)
	assign("Qbar", Qbar,  envir = tempenvir)
	assign("H", H, envir = tempenvir)
	dccpars = as.numeric( unlist(spec@mspec$optimization.model$fixed.pars) )
	
	assign("npar", length(dccpars), envir = tempenvir)
	
	filt = .dccmakefiltermodel(garchmodel = "dcclaplace", f = laplace.dccLLH2, 
			pars = dccpars, garchpars, zdata, garchenv = tempenvir,
			timer = 0, message = 0, fname = "laplace.dccLLH2")
	
	model = list()
	model$vrmodel = vrmodel
	model$asset.names = cnames
	model$distribution = spec@mspec$distribution
	model$include.skew = spec@mspec$include.skew
	model$include.shape = spec@mspec$include.shape
	model$dccOrder = spec@mspec$dccOrder
	model$out.sample = out.sample
	
	filt$dates = dates
	filt$date.format = dformat
	filt$origdata  = origdata
	filt$origdates = origdates
	filt$model = model
	filt$timer = Sys.time() - tic
	
	ans = new("DCCfilter",
			mfilter = filt,
			ufilter = filterlist)
	rm(tempenvir)
	return(ans)
	
}

.dccforecast = function(fit, n.ahead = 1, n.roll = 0, external.forecasts = list(mregfor = NULL, vregfor = NULL), 
		parallel = FALSE, parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2), ...)
{
	# checks for out.sample wrt n.roll
	
	ns = fit@mfit$model$out.sample
	if( n.roll > ns ) stop("n.roll must not be greater than out.sample!")
	if(n.roll>1 && n.ahead>1) stop("\ngogarchforecast-->error: n.ahead must be equal to 1 when using n.roll\n")
	
	# checks for external forecasts
	tf = n.ahead + n.roll
	if( !is.null( external.forecasts$mregfor ) ){
		mregfor = external.forecasts$mregfor
		if( !is.matrix(mregfor) ) stop("\nmregfor must be a matrix.")
		if( dim(mregfor)[1] < tf ) stop("\nmregfor must have at least n.ahead + n.roll observations to be used")
		mregfor = mregfor[1:tf, , drop = FALSE]
	} else{
		mregfor = NULL
	}
	if( !is.null( external.forecasts$vregfor ) ){
		if( !is.matrix(vregfor) ) stop("\nvregfor must be a matrix.")
		if( dim(vregfor)[1] < tf ) stop("\nvregfor must have at least n.ahead + n.roll observations to be used")
		vregfor = vregfor[1:tf, , drop = FALSE]
	}
	vrmodel = fit@mfit$model$vrmodel
	if( !is.null(vrmodel) ){
		mxn = vrmodel$mxn
		if( mxn > 0 ){
			if( is.null(external.forecasts$mregfor ) ){
				warning("\nExternal Regressor Forecasts Matrix NULL...setting to zero...\n")
				mregfor = matrix(0, ncol = mxn, nrow = (n.roll + n.ahead) )
			} else{
				mxn = vrmodel$mxn
				if( dim(mregfor)[2] != mxn ) stop("\ndccforecast-->error: wrong number of external regressors!...", call. = FALSE)
				if( dim(mregfor)[1] < (n.roll + n.ahead) ) stop("\ndccforecast-->error: external regressor matrix has less points than requested forecast length (1+n.roll) x n.ahead!...", call. = FALSE)
			}
		} else{
			mregfor = NULL
		}
		if(n.roll>1 && n.ahead>1) stop("\ndccforecast-->error: n.ahead must be equal to 1 when using n.roll\n")
		
		if( n.ahead == 1 && (n.roll > ns) ) stop("\ndccforecast-->error: n.roll greater than out.sample!", call. = FALSE)
		#VARf = .varpredict(fit, n.ahead, n.roll, mregfor)$Mu
		VARf = varxforecast(X = fit@mfit$origdata, Bcoef = vrmodel$Bcoef, p = vrmodel$lag, out.sample = ns, 
				n.ahead, n.roll, mregfor)
	} else{
		VARf = NULL
	}
	exf = external.forecasts
	if(  !is.null(fit@mfit$model$vrmodel) ){
		exf$mregfor = NULL
	}
	ans = .dccforecastm(fit, n.ahead = n.ahead, n.roll = n.roll, external.forecasts = exf, 
			parallel = parallel, parallel.control = parallel.control, ...)
	
	model = list()
	model$n.roll = n.roll
	model$n.ahead = n.ahead
	model$asset.names = fit@mfit$model$asset.names
	model$distribution = fit@mfit$model$distribution
	model$include.skew = fit@mfit$model$include.skew
	model$include.shape = fit@mfit$model$include.shape
	model$dccOrder = fit@mfit$model$dccOrder
	model$out.sample = ns
	model$dccpars = coef(fit, type = "dcc")
	model$garchpars = coef(fit, type = "garch")
	model$origdata = fit@mfit$origdata
	model$origdates = fit@mfit$origdates
	
	# keep last 50 points of R and H for plotting
	model$R = last( rcor(fit), min(length(fit@mfit$R), 50) )
	model$H = last( rcov(fit), min(length(fit@mfit$R), 50) )
	model$fitted = tail(fitted(fit), min(length(fit@mfit$R), 50) )
	model$sigma = tail(sigma(fit), min(length(fit@mfit$R), 50) )
	mforecast = list( H = ans$H, R = ans$R, Q = ans$Q, Rbar = ans$Rbar, VARf = VARf )
	mforecast$model = model
	uforecast = ans$ufor
	mforecast$shape = mforecast$skew = NULL
	if( mforecast$model$include.shape ) mforecast$shape = rshape(fit)
	if( mforecast$model$include.skew ) mforecast$skew = rskew(fit)
	
	ans = new("DCCforecast",
			mforecast = mforecast,
			uforecast = uforecast)
	
	return( ans )
}


dccsim1.norm = function(fitORspec, n.sim = 1000, n.start = 0, m.sim = 1, startMethod = c("unconditional", "sample"), 
		presigma = NULL, preresiduals = NULL, prereturns = NULL, preR = NULL, preH = NULL, rseed = NULL, mexsimdata = NULL, 
		vexsimdata = NULL, parallel = FALSE, parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2), VAR.fit = NULL, ...)
{
	fit = fitORspec
	Data = head(fit@mfit$origdata, dim(fit@mfit$origdata)[1] - fit@mfit$model$out.sample)
	m = dim(Data)[2]
	n = dim(Data)[1]
	
	dccOrder = fit@mfit$model$dccOrder
	mo = max( dccOrder )
	
	if( is.null(rseed) ){
		rseed = as.integer(runif(1, 1, Sys.time())) 
	} else {
		if(length(rseed) == 1) rseed = as.integer(rseed[1]) else rseed = as.integer( rseed[1:m.sim] )
	}
	
	# we use the GED distribution  with nu = 1 which corresponds to the Laplace case.
	if(length(rseed) == 1){
		tmp = rnorm(m * (n.sim + n.start + mo) * m.sim, 0, 1)
		z = array(tmp,  dim = c(n.sim + n.start + mo, m, m.sim))
	} else{
		z = array(NA, dim = c(n.sim + n.start + mo, m, m.sim))
		for(i in 1:m.sim){
			set.seed( rseed[i] )
			z[,,i] = matrix(rnorm(m * (n.sim + n.start + mo), 0, 1), nrow = n.sim + n.start + mo, ncol = m)
		}
	}
	# ok now to expand rseed
	if(length(rseed) == 1){
		rseed = c(rseed, as.integer(runif(m.sim, 1, Sys.time()))) 
	}
	if( is.null(preR) ){
		if( startMethod[1] == "sample" ){
			preR = last(rcor(fit))[,,1]
			diag(preR) = 1
		} else{
			preR = cor(as.matrix(Data))
		}
	} else{
		chk = dcc.symcheck(preR, m, d = rep(1, m))
	}
	
	preQ = preR
		
	uncv = sapply(fit@ufit@fit, FUN = function(x) uncvariance(x))
	
	
	if( is.null(preH) ){
		preH = diag(sqrt(uncv)) %*% preR %*% diag(sqrt(uncv))
	} else{
		chk = dcc.symcheck(preH, m, d = NULL)
	}
	if( !is.null(presigma) ){
		if( !is.matrix(presigma) ) stop("\ndccsim-->error: presigma must be a matrix.")
		if( dim(presigma)[2] != m ) stop("\ndccsim-->error: wrong column dimension for presigma.")
	}
	
	if( !is.null(preresiduals) ){
		if( !is.matrix(preresiduals) ) stop("\ndccsim-->error: preresiduals must be a matrix.")
		if( dim(preresiduals)[2] != m ) stop("\ndccsim-->error: wrong column dimension for preresiduals.")
	}
	
	if( !is.null(prereturns) ){
		if( !is.matrix(prereturns) ) stop("\ndccsim-->error: prereturns must be a matrix.")
		if( dim(prereturns)[2] != m ) stop("\ndccsim-->error: wrong column dimension for prereturns.")
	}
	
	garchpars	= coef(fit, type = "garch")
	dccpars		= coef(fit, type = "dcc")
	dcca = dccpars[1:dccOrder[1]]
	dccb = dccpars[(dccOrder[1]+1):sum(dccOrder[1:2])]
	
	simRes = simX = simR = simQ = simH = simSeries = vector(mode = "list", length = m.sim)
	
	if( parallel ){
		os = .Platform$OS.type
		if(is.null(parallel.control$pkg)){
			if( os == "windows" ) parallel.control$pkg = "snowfall" else parallel.control$pkg = "multicore"
			if( is.null(parallel.control$cores) ) parallel.control$cores = 2
		} else{
			mtype = match(tolower(parallel.control$pkg[1]), c("multicore", "snowfall"))
			if(is.na(mtype)) stop("\nParallel Package type not recognized in parallel.control\n")
			parallel.control$pkg = tolower(parallel.control$pkg[1])
			if( os == "windows" && parallel.control$pkg == "multicore" ) stop("\nmulticore not supported on windows O/S\n")
			if( is.null(parallel.control$cores) ) parallel.control$cores = 2 else parallel.control$cores = as.integer(parallel.control$cores[1])
		}
		if( parallel.control$pkg == "multicore" ){
			simH = vector(mode = "list", length = m.sim)
			simX = vector(mode = "list", length = m.sim)
			mtmp = mclapply(as.list(1:m.sim), FUN = function(j)
						.dccsimf(Z = z[,,j], Qbar = preQ, Rbar = preR, dcca, dccb, mo, n.sim, n.start, m, rseed[j]), 
					mc.cores = parallel.control$cores)
			# need to pass m.sim and matrix of simulated z's to ugarchsim (speedup)
			simR = lapply(mtmp, FUN = function(x) if(is.matrix(x$R)) array(x$R, dim = c(m, m, n.sim)) else last(x$R, n.sim))
			simQ = lapply(mtmp, FUN = function(x) if(is.matrix(x$Q)) array(x$Q, dim = c(m, m, n.sim)) else last(x$Q, n.sim))
			stdresid = vector(mode = "list", length = m)
			for(i in 1:m) stdresid[[i]] = sapply(mtmp, FUN = function(x) x$stdresid[,i])
			xtmp = mclapply(as.list(1:m), FUN = function(j){
						maxx = fit@ufit@fit[[j]]@model$maxOrder2;
						htmp = ugarchsim(fit@ufit@fit[[j]], n.sim = n.sim + n.start, n.start = 0, m.sim = m.sim,
								startMethod = startMethod[1], custom.dist = list(name = "sample", 
										distfit = matrix(stdresid[[j]][-(1:mo), ], ncol = m.sim)),
								presigma = if( is.null(presigma) ) NA else tail(presigma[,j], maxx), 
								preresiduals = if( is.null(preresiduals) ) NA else tail(preresiduals[,j], maxx), 
								prereturns = if( is.null(prereturns) | !is.null(fit@mfit$model$vrmodel) ) NA else tail(prereturns[,j], maxx),
								mexsimdata = if( is.null(fit@mfit$model$vrmodel) ) mexsimdata[[j]] else NULL, vexsimdata = vexsimdata[[j]] )
						h = matrix(tail(htmp@simulation$sigmaSim^2, n.sim), nrow = n.sim)
						x = matrix(htmp@simulation$seriesSim,  nrow = n.sim + n.start)
						return(list(h = h, x = x))}, mc.cores = parallel.control$cores)
			H = array(NA, dim = c(n.sim, m, m.sim))
			tmpH = array(NA, dim = c(m, m, n.sim))
			for(i in 1:n.sim) H[i,,] = t(sapply(xtmp, FUN = function(x) as.numeric(x$h[i,])))
			
			for(i in 1:m.sim){
				for(j in 1:n.sim){
					tmpH[ , , j] = diag(sqrt( H[j, , i]) ) %*% simR[[i]][ , , j] %*% diag(sqrt( H[j, , i] ) )
				}
				simH[[i]] = tmpH
			}
			simxX = array(NA, dim = c(n.sim, m, m.sim))
			for(i in 1:n.sim) simxX[i,,] = t(sapply(xtmp, FUN = function(x) as.numeric(x$x[i,])))
			simX = vector(mode = "list", length = m.sim)
			for(i in 1:m.sim) simX[[i]] = matrix(simxX[,,i], nrow = n.sim)
		} else{
			simH = vector(mode = "list", length = m.sim)
			simX = vector(mode = "list", length = m.sim)
			sfInit(parallel = TRUE, cpus = parallel.control$cores)
			sfExport("z", "preQ", "preR", "dcca", "dccb", "mo", "n.sim", "n.start", "m", "rseed", local = TRUE)
			mtmp = sfLapply(as.list(1:m.sim), fun = function(j)
						rgarch:::.dccsimf(Z = z[,,j], Qbar = preQ, Rbar = preR, dcca, dccb, mo, n.sim, n.start, m, rseed[j]))
			
			# need to pass m.sim and matrix of simulated z's to ugarchsim (speedup)
			simR = lapply(mtmp, FUN = function(x) if(is.matrix(x$R)) array(x$R, dim = c(m, m, n.sim)) else last(x$R, n.sim))
			simQ = lapply(mtmp, FUN = function(x) if(is.matrix(x$Q)) array(x$Q, dim = c(m, m, n.sim)) else last(x$Q, n.sim))
			stdresid = vector(mode = "list", length = m)
			for(i in 1:m) stdresid[[i]] = sapply(mtmp, FUN = function(x) x$stdresid[,i])			
			sfExport("fit", "n.sim", "n.start", "m.sim", "startMethod", "stdresid", "presigma", 
					"preresiduals", "prereturns", "mexsimdata", "vexsimdata", local = TRUE)
			xtmp = sfLapply(as.list(1:m), fun = function(j){
						maxx = fit@ufit@fit[[j]]@model$maxOrder2;
						htmp = rgarch::ugarchsim(fit@ufit@fit[[j]], n.sim = n.sim + n.start, n.start = 0, m.sim = m.sim,
								startMethod = startMethod[1], custom.dist = list(name = "sample", 
										distfit = matrix(stdresid[[j]][-(1:mo), ], ncol = m.sim)),
								presigma = if( is.null(presigma) ) NA else tail(presigma[,j], maxx), 
								preresiduals = if( is.null(preresiduals) ) NA else tail(preresiduals[,j], maxx), 
								prereturns = if( is.null(prereturns) | !is.null(fit@mfit$model$vrmodel) ) NA else tail(prereturns[,j], maxx),
								mexsimdata = if( is.null(fit@mfit$model$vrmodel) ) mexsimdata[[j]] else NULL, vexsimdata = vexsimdata[[j]] )
						h = matrix(tail(htmp@simulation$sigmaSim^2, n.sim), nrow = n.sim)
						x = matrix(htmp@simulation$seriesSim,  nrow = n.sim + n.start)
						return(list(h = h, x = x))})
			if( is.null(fit@mfit$model$vrmodel) ) sfStop()
			H = array(NA, dim = c(n.sim, m, m.sim))
			tmpH = array(NA, dim = c(m, m, n.sim))
			for(i in 1:n.sim) H[i,,] = t(sapply(xtmp, FUN = function(x) as.numeric(x$h[i,])))
			
			for(i in 1:m.sim){
				for(j in 1:n.sim){
					tmpH[ , , j] = diag(sqrt( H[j, , i]) ) %*% simR[[i]][ , , j] %*% diag(sqrt( H[j, , i] ) )
				}
				simH[[i]] = tmpH
			}
			simxX = array(NA, dim = c(n.sim, m, m.sim))
			for(i in 1:n.sim) simxX[i,,] = t(sapply(xtmp, FUN = function(x) as.numeric(x$x[i,])))
			simX = vector(mode = "list", length = m.sim)
			for(i in 1:m.sim) simX[[i]] = matrix(simxX[,,i], nrow = n.sim)
		}
	} else{
		simH = vector(mode = "list", length = m.sim)
		simX = vector(mode = "list", length = m.sim)
		mtmp = lapply(as.list(1:m.sim), FUN = function(j) .dccsimf(Z = z[,,j], Qbar = preQ, Rbar = preR, dcca, dccb, mo, n.sim, n.start, m, rseed[j]))
		# need to pass m.sim and matrix of simulated z's to ugarchsim (speedup)
		simR = lapply(mtmp, FUN = function(x) if(is.matrix(x$R)) array(x$R, dim = c(m, m, n.sim)) else last(x$R, n.sim))
		simQ = lapply(mtmp, FUN = function(x) if(is.matrix(x$Q)) array(x$Q, dim = c(m, m, n.sim)) else last(x$Q, n.sim))
		stdresid = vector(mode = "list", length = m)
		for(i in 1:m) stdresid[[i]] = sapply(mtmp, FUN = function(x) x$stdresid[,i])			
		xtmp = lapply(as.list(1:m), FUN = function(j){
					maxx = fit@ufit@fit[[j]]@model$maxOrder2;
					htmp = ugarchsim(fit@ufit@fit[[j]], n.sim = n.sim + n.start, n.start = 0, m.sim = m.sim,
							startMethod = startMethod[1], custom.dist = list(name = "sample", 
									distfit = matrix(stdresid[[j]][-(1:mo), ], ncol = m.sim)),
							presigma = if( is.null(presigma) ) NA else tail(presigma[,j], maxx), 
							preresiduals = if( is.null(preresiduals) ) NA else tail(preresiduals[,j], maxx), 
							prereturns = if( is.null(prereturns) | !is.null(fit@mfit$model$vrmodel) ) NA else tail(prereturns[,j], maxx),
							mexsimdata = if( is.null(fit@mfit$model$vrmodel) ) mexsimdata[[j]] else NULL, vexsimdata = vexsimdata[[j]] )
					h = matrix(tail(htmp@simulation$sigmaSim^2, n.sim), nrow = n.sim)
					x = matrix(htmp@simulation$seriesSim,  nrow = n.sim + n.start)
					return(list(h = h, x = x))})
		H = array(NA, dim = c(n.sim, m, m.sim))
		tmpH = array(NA, dim = c(m, m, n.sim))
		for(i in 1:n.sim) H[i,,] = t(sapply(xtmp, FUN = function(x) as.numeric(x$h[i,])))	
		for(i in 1:m.sim){
			for(j in 1:n.sim){
				tmpH[ , , j] = diag(sqrt( H[j, , i]) ) %*% simR[[i]][ , , j] %*% diag(sqrt( H[j, , i] ) )
			}
			simH[[i]] = tmpH
		}
		simxX = array(NA, dim = c(n.sim, m, m.sim))
		for(i in 1:n.sim) simxX[i,,] = t(sapply(xtmp, FUN = function(x) as.numeric(x$x[i,])))
		simX = vector(mode = "list", length = m.sim)
		for(i in 1:m.sim) simX[[i]] = matrix(simxX[,,i], nrow = n.sim)
	}
	
	if( !is.null(fit@mfit$model$vrmodel) ){
		# mexsimdata check
		mxn = fit@mfit$model$vrmodel$mxn
		if(mxn > 0 ){
			if( !is.null(mexsimdata) ){
				if( !is.list(mexsimdata) | length(mexsimdata) != m.sim ) stop("\ndccsim-->error: mexsimdata must be a list of length equal to m.sim...", call. = FALSE)
				for(j in 1:m.sim){
					if( !is.matrix(mexsimdata[[j]]) ) stop("\ndccsim-->error: mexsimdata list submatrix must be a matrix...", call. = FALSE)
					if( dim(mexsimdata[[j]])[2] != mxn ) stop("\ndccsim-->error: wrong number of external regressors in list submatrix!...", call. = FALSE)
					if( dim(mexsimdata[[j]])[1] < (n.sim + n.start) ) stop("\ndccsim-->error: external regressor list submatrix has less points than requested simulation length (n.sim + n.start)...", call. = FALSE)
				}
			} else{
				mexsimdata = vector(mode = "list", length = m.sim)
				for(j in 1:m.sim) mexsimdata[[i]] = matrix(0, ncol = mxn, nrow = (n.sim + n.start))
			}
		} else{
			mexsimdata = NULL
		}
		
		vrmodel = fit@mfit$model$vrmodel
		p = as.numeric(vrmodel$lag)
		
		if(startMethod == "sample"){
			if( !is.null(prereturns) ){
				if(!is.matrix(prereturns)) stop("\nprereturns must be a matrix\n")
				if(dim(prereturns)[1] != p) stop("\nwrong number of rows for prereturns matrix")
				if(dim(prereturns)[2] != m) stop("\nwrong number of cols for prereturns matrix")
				prereturns = matrix(prereturns, nrow = p, ncol = m)
			} else{
				prereturns = as.matrix( tail(Data, p) )
			}
		} else{
			prereturns = matrix(0, nrow = p, ncol = m, byrow = TRUE)
		}
		xlist = vector(mode = "list", length = m.sim)
		if( n.sim == 1 && n.start == 0 ){
			# SPECIAL ROUTINE TO SPEED UP 1-AHEAD ESTIMATION
			tmp = varxsimXX(X = Data, vrmodel$Bcoef, p, m.sim = m.sim, prereturns, resids = t(sapply(simX, FUN = function(x) x)), if(!is.null(mexsimdata)) t(sapply(mexsimdata, FUN = function(x) x)) else NULL)
			for(i in 1:m.sim) xlist[[i]] = tmp[i, , drop = FALSE]
			simRes = simX
			simX = xlist
		} else{
			if( parallel ){
				if( parallel.control$pkg == "multicore" ){
					xlist = mclapply(as.list(1:m.sim), FUN = function(i) varxsim(X = Data, vrmodel$Bcoef, p, n.sim, 
										n.start, prereturns, resids = simX[[i]], mexsimdata[[i]]), 
							mc.cores = parallel.control$cores)
				} else{
					# when the GARCH model is without mean equation, simX == simulated residuals
					sfExport("Data", "vrmodel", "p", "n.sim", "n.start", "prereturns", "simX", "mexsimdata", local = TRUE)
					xlist = sfLapply(as.list(1:m.sim), fun = function(i) rgarch:::varxsim(X = Data, vrmodel$Bcoef, p, 
										n.sim, n.start, prereturns, resids = simX[[i]], mexsimdata[[i]]))
					sfStop()
				}
			} else{
				xlist = lapply(as.list(1:m.sim), FUN = function(i) varxsim(X = Data, vrmodel$Bcoef, p, 
									n.sim, n.start, prereturns, resids = simX[[i]], mexsimdata[[i]]))
			}
			simRes = simX
			simX = xlist
		}
	} else{
		# reshape
		for(j in 1:m.sim) simX[[j]] = tail(simX[[j]], n.sim)
	}
	
	msim = list()
	msim$varmodel = fit@mfit$vrmodel
	msim$simH = simH
	msim$simR = simR
	msim$simQ = simQ
	msim$simX = simX
	msim$simRes = simRes
	msim$Z = z
	msim$rseed = rseed
	msim$model$Data = Data
	msim$model$garchpars = garchpars
	msim$model$dccpars = dccpars
	msim$model$n.sim = n.sim
	msim$model$m.sim = m.sim
	msim$model$n.start = n.start
	msim$model$startMethod = startMethod[1]
	msim$model$dccOrder = dccOrder
	msim$model$distribution = "mvnorm"
	ans = new("DCCsim",
			msim = msim)
	
	return( ans )
	
}

dccsim2.norm = function(fitORspec, n.sim = 1000, n.start = 0, m.sim = 1, startMethod = c("unconditional", 
				"sample"), presigma = NULL, preresiduals = NULL, 
		prereturns = NULL, preR = NULL, preH = NULL, rseed = NULL, mexsimdata = NULL, vexsimdata = NULL, 
		parallel = FALSE, parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2), VAR.fit = NULL, ...)
{
	spec = fitORspec
	
	# check that VAR.fit is included if used
	if( spec@mspec$varmodel$VAR ){
		if( is.null(VAR.fit) ) stop("\ndccsim-->error: VAR.fit must not be NULL for VAR method when calling dccsim using spec!", call. = FALSE)
	}
	
	dccOrder = spec@mspec$dccOrder
	mo = max( dccOrder )
	m = length(spec@uspec@spec)
	if( is.null(rseed) ){
		rseed = as.integer(runif(1, 1, Sys.time())) 
	} else {
		if(length(rseed) == 1) rseed = as.integer(rseed[1]) else rseed = as.integer( rseed[1:m.sim] )
	}
	
	# we use the GED distribution  with nu = 1 which corresponds to the Laplace case.
	if(length(rseed) == 1){
		tmp = rnorm(m * (n.sim + n.start + mo) * m.sim, 0, 1)
		z = array(tmp,  dim = c(n.sim + n.start + mo, m, m.sim))
	} else{
		z = array(NA, dim = c(n.sim + n.start + mo, m, m.sim))
		for(i in 1:m.sim){
			set.seed( rseed[i] )
			z[,,i] = matrix(rnorm(m * (n.sim + n.start + mo), 0, 1), nrow = n.sim + n.start + mo, ncol = m)
		}
	}
	if(length(rseed) == 1){
		rseed = c(rseed, as.integer(runif(m.sim, 1, Sys.time()))) 
	}
	if( is.null(preR) ){
		stop("\ndccsim-->error: preR cannot be NULL when method uses only DCC specification!")
	} else{
		chk = dcc.symcheck(preR, m, d = rep(1, m))
	}
	
	preQ = preR
	
	# uncvariance has uGARCHfit AND uGARCHspec dispatch methods
	uncv = sapply(spec@uspec@spec, FUN = function(x) uncvariance(x))
	
	
	if( is.null(preH) ){
		preH = diag(sqrt(uncv)) %*% preR %*% diag(sqrt(uncv))
	} else{
		chk = dcc.symcheck(preH, m, d = NULL)
	}
	if( !is.null(presigma) ){
		if( !is.matrix(presigma) ) stop("\ndccsim-->error: presigma must be a matrix.")
		if( dim(presigma)[2] != m ) stop("\ndccsim-->error: wrong column dimension for presigma.")
	}
	
	if( !is.null(preresiduals) ){
		if( !is.matrix(preresiduals) ) stop("\ndccsim-->error: preresiduals must be a matrix.")
		if( dim(preresiduals)[2] != m ) stop("\ndccsim-->error: wrong column dimension for preresiduals.")
	}
	
	if( !is.null(prereturns) ){
		if( !is.matrix(prereturns) ) stop("\ndccsim-->error: prereturns must be a matrix.")
		if( dim(prereturns)[2] != m ) stop("\ndccsim-->error: wrong column dimension for prereturns.")
	}
	
	garchpars	= sapply(spec@uspec@spec, FUN = function(x) unlist(x@optimization.model$fixed.pars))
	dccpars		= unlist(spec@mspec$optimization.model$fixed.pars)
	dcca = dccpars[1:dccOrder[1]]
	dccb = dccpars[(dccOrder[1]+1):sum(dccOrder[1:2])]
	
	simRes = simX = simR = simQ = simH = simSeries = vector(mode = "list", length = m.sim)
	
	if( parallel ){
		os = .Platform$OS.type
		if(is.null(parallel.control$pkg)){
			if( os == "windows" ) parallel.control$pkg = "snowfall" else parallel.control$pkg = "multicore"
			if( is.null(parallel.control$cores) ) parallel.control$cores = 2
		} else{
			mtype = match(tolower(parallel.control$pkg[1]), c("multicore", "snowfall"))
			if(is.na(mtype)) stop("\nParallel Package type not recognized in parallel.control\n")
			parallel.control$pkg = tolower(parallel.control$pkg[1])
			if( os == "windows" && parallel.control$pkg == "multicore" ) stop("\nmulticore not supported on windows O/S\n")
			if( is.null(parallel.control$cores) ) parallel.control$cores = 2 else parallel.control$cores = as.integer(parallel.control$cores[1])
		}
		if( parallel.control$pkg == "multicore" ){
			simH = vector(mode = "list", length = m.sim)
			simX = vector(mode = "list", length = m.sim)
			mtmp = mclapply(as.list(1:m.sim), FUN = function(j)
						.dccsimf(Z = z[,,j], Qbar = preQ, Rbar = preR, dcca, dccb, mo, n.sim, n.start, m, rseed[j]), 
					mc.cores = parallel.control$cores)
			simR = lapply(mtmp, FUN = function(x) if(is.matrix(x$R)) array(x$R, dim = c(m, m, n.sim)) else last(x$R, n.sim))
			simQ = lapply(mtmp, FUN = function(x) if(is.matrix(x$Q)) array(x$Q, dim = c(m, m, n.sim)) else last(x$Q, n.sim))
			stdresid = vector(mode = "list", length = m)
			for(i in 1:m) stdresid[[i]] = sapply(mtmp, FUN = function(x) x$stdresid[,i])			
			xtmp = mclapply(as.list(1:m), FUN = function(j){
						maxx = spec@uspec@spec[[j]]@optimization.model$maxOrder2;
						htmp = ugarchpath(spec@uspec@spec[[j]], n.sim = n.sim + n.start, n.start = 0, m.sim = m.sim,
								custom.dist = list(name = "sample", distfit = matrix(stdresid[[j]][-(1:mo), ], ncol = m.sim)),
								presigma = if( is.null(presigma) ) NA else tail(presigma[,j], maxx), 
								preresiduals = if( is.null(preresiduals) ) NA else tail(preresiduals[,j], maxx), 
								prereturns = if( is.null(prereturns) | spec@mspec$varmodel$VAR ) NA else tail(prereturns[,j], maxx),
								mexsimdata = if( !spec@mspec$varmodel$VAR ) mexsimdata[[j]] else NULL, vexsimdata = vexsimdata[[j]] )
						h = matrix(tail(htmp@path$sigmaSim^2, n.sim), nrow = n.sim)
						x = matrix(htmp@path$seriesSim,  nrow = n.sim + n.start)
						return(list(h = h, x = x))}, mc.cores = parallel.control$cores)
			H = array(NA, dim = c(n.sim, m, m.sim))
			tmpH = array(NA, dim = c(m, m, n.sim))
			for(i in 1:n.sim) H[i,,] = t(sapply(xtmp, FUN = function(x) as.numeric(x$h[i,])))
			
			for(i in 1:m.sim){
				for(j in 1:n.sim){
					tmpH[ , , j] = diag(sqrt( H[j, , i]) ) %*% simR[[i]][ , , j] %*% diag(sqrt( H[j, , i] ) )
				}
				simH[[i]] = tmpH
			}
			simxX = array(NA, dim = c(n.sim, m, m.sim))
			for(i in 1:n.sim) simxX[i,,] = t(sapply(xtmp, FUN = function(x) as.numeric(x$x[i,])))
			simX = vector(mode = "list", length = m.sim)
			for(i in 1:m.sim) simX[[i]] = matrix(simxX[,,i], nrow = n.sim)
		} else{
			sfInit(parallel = TRUE, cpus = parallel.control$cores)
			sfExport("spec", "z", "preQ", "preR", "dcca", "dccb", "mo", "n.sim", "n.start", 
					"m", "prereturns", "preresiduals", "presigma", "mexsimdata", "vexsimdata", 
					"rseed", local = TRUE)
			simH = vector(mode = "list", length = m.sim)
			simX = vector(mode = "list", length = m.sim)
			mtmp = sfLapply(as.list(1:m.sim), fun = function(j) .dccsimf(Z = z[,,j], Qbar = preQ, Rbar = preR, dcca, dccb, mo, n.sim, n.start, m, rseed[j]))
			simR = lapply(mtmp, FUN = function(x) if(is.matrix(x$R)) array(x$R, dim = c(m, m, n.sim)) else last(x$R, n.sim))
			simQ = lapply(mtmp, FUN = function(x) if(is.matrix(x$Q)) array(x$Q, dim = c(m, m, n.sim)) else last(x$Q, n.sim))
			stdresid = vector(mode = "list", length = m)
			for(i in 1:m) stdresid[[i]] = sapply(mtmp, FUN = function(x) x$stdresid[,i])			
			xtmp = sfLapply(as.list(1:m), fun = function(j){
						maxx = spec@uspec@spec[[j]]@optimization.model$maxOrder2;
						htmp = rgarch::ugarchpath(spec@uspec@spec[[j]], n.sim = n.sim + n.start, n.start = 0, m.sim = m.sim,
								custom.dist = list(name = "sample", distfit = matrix(stdresid[[j]][-(1:mo), ], ncol = m.sim)),
								presigma = if( is.null(presigma) ) NA else tail(presigma[,j], maxx), 
								preresiduals = if( is.null(preresiduals) ) NA else tail(preresiduals[,j], maxx), 
								prereturns = if( is.null(prereturns) | spec@mspec$varmodel$VAR ) NA else tail(prereturns[,j], maxx),
								mexsimdata = if( !spec@mspec$varmodel$VAR ) mexsimdata[[j]] else NULL, vexsimdata = vexsimdata[[j]] )
						h = matrix(tail(htmp@path$sigmaSim^2, n.sim), nrow = n.sim)
						x = matrix(htmp@path$seriesSim,  nrow = n.sim + n.start)
						return(list(h = h, x = x))})
			if( !spec@mspec$varmodel$VAR ) sfStop()
			H = array(NA, dim = c(n.sim, m, m.sim))
			tmpH = array(NA, dim = c(m, m, n.sim))
			for(i in 1:n.sim) H[i,,] = t(sapply(xtmp, FUN = function(x) as.numeric(x$h[i,])))
			
			for(i in 1:m.sim){
				for(j in 1:n.sim){
					tmpH[ , , j] = diag(sqrt( H[j, , i]) ) %*% simR[[i]][ , , j] %*% diag(sqrt( H[j, , i] ) )
				}
				simH[[i]] = tmpH
			}
			simxX = array(NA, dim = c(n.sim, m, m.sim))
			for(i in 1:n.sim) simxX[i,,] = t(sapply(xtmp, FUN = function(x) as.numeric(x$x[i,])))
			simX = vector(mode = "list", length = m.sim)
			for(i in 1:m.sim) simX[[i]] = matrix(simxX[,,i], nrow = n.sim)
		}
	} else{
		simH = vector(mode = "list", length = m.sim)
		simX = vector(mode = "list", length = m.sim)
		mtmp = lapply(as.list(1:m.sim), FUN = function(j) .dccsimf(Z = z[,,j], Qbar = preQ, Rbar = preR, dcca, dccb, mo, n.sim, n.start, m, rseed[j]))
		simR = lapply(mtmp, FUN = function(x) if(is.matrix(x$R)) array(x$R, dim = c(m, m, n.sim)) else last(x$R, n.sim))
		simQ = lapply(mtmp, FUN = function(x) if(is.matrix(x$Q)) array(x$Q, dim = c(m, m, n.sim)) else last(x$Q, n.sim))
		stdresid = vector(mode = "list", length = m)
		for(i in 1:m) stdresid[[i]] = sapply(mtmp, FUN = function(x) x$stdresid[,i])			
		xtmp = lapply(as.list(1:m), FUN = function(j){
					maxx = spec@uspec@spec[[j]]@optimization.model$maxOrder2;
					htmp = ugarchpath(spec@uspec@spec[[j]], n.sim = n.sim + n.start, n.start = 0, m.sim = m.sim,
							custom.dist = list(name = "sample", distfit = matrix(stdresid[[j]][-(1:mo), ], ncol = m.sim)),
							presigma = if( is.null(presigma) ) NA else tail(presigma[,j], maxx), 
							preresiduals = if( is.null(preresiduals) ) NA else tail(preresiduals[,j], maxx), 
							prereturns = if( is.null(prereturns) | spec@mspec$varmodel$VAR ) NA else tail(prereturns[,j], maxx),
							mexsimdata = if( !spec@mspec$varmodel$VAR ) mexsimdata[[j]] else NULL, vexsimdata = vexsimdata[[j]] )
					h = matrix(tail(htmp@path$sigmaSim^2, n.sim), nrow = n.sim)
					x = matrix(htmp@path$seriesSim,  nrow = n.sim + n.start)
					return(list(h = h, x = x))})
		H = array(NA, dim = c(n.sim, m, m.sim))
		tmpH = array(NA, dim = c(m, m, n.sim))
		for(i in 1:n.sim) H[i,,] = t(sapply(xtmp, FUN = function(x) as.numeric(x$h[i,])))
		
		for(i in 1:m.sim){
			for(j in 1:n.sim){
				tmpH[ , , j] = diag(sqrt( H[j, , i]) ) %*% simR[[i]][ , , j] %*% diag(sqrt( H[j, , i] ) )
			}
			simH[[i]] = tmpH
		}
		simxX = array(NA, dim = c(n.sim, m, m.sim))
		for(i in 1:n.sim) simxX[i,,] = t(sapply(xtmp, FUN = function(x) as.numeric(x$x[i,])))
		simX = vector(mode = "list", length = m.sim)
		for(i in 1:m.sim) simX[[i]] = matrix(simxX[,,i], nrow = n.sim)
	}
	if( spec@mspec$varmodel$VAR ){
		# mexsimdata check
		mxn = VAR.fit$mxn
		Data = VAR.fit$xfitted
		if(mxn > 0 ){
			if( !is.null(mexsimdata) ){
				if( !is.list(mexsimdata) | length(mexsimdata) != m.sim ) stop("\ndccsim-->error: mexsimdata must be a list of length equal to m.sim...", call. = FALSE)
				for(j in 1:m.sim){
					if( !is.matrix(mexsimdata[[j]]) ) stop("\ndccsim-->error: mexsimdata list submatrix must be a matrix...", call. = FALSE)
					if( dim(mexsimdata[[j]])[2] != mxn ) stop("\ndccsim-->error: wrong number of external regressors in list submatrix!...", call. = FALSE)
					if( dim(mexsimdata[[j]])[1] < (n.sim + n.start) ) stop("\ndccsim-->error: external regressor list submatrix has less points than requested simulation length (n.sim + n.start)...", call. = FALSE)
				}
			} else{
				mexsimdata = vector(mode = "list", length = m.sim)
				for(j in 1:m.sim) mexsimdata[[i]] = matrix(0, ncol = mxn, nrow = (n.sim + n.start))
			}
		} else{
			mexsimdata = NULL
		}
		
		p = VAR.fit$lag
		
		if( !is.null(prereturns) ){
			if(!is.matrix(prereturns)) stop("\nprereturns must be a matrix\n")
			if(dim(prereturns)[1] != p) stop("\nwrong number of rows for prereturns matrix")
			if(dim(prereturns)[2] != m) stop("\nwrong number of cols for prereturns matrix")
			prereturns = matrix(prereturns, nrow = p, ncol = m)
		} else{
			prereturns = matrix(0, nrow = p, ncol = m, byrow = TRUE)
		}
		xlist = vector(mode = "list", length = m.sim)
		if( n.sim == 1 && n.start == 0 ){
			# SPECIAL ROUTINE TO SPEED UP 1-AHEAD ESTIMATION
			tmp = varxsimXX(X = Data, VAR.fit$Bcoef, p, m.sim = m.sim, prereturns, resids = t(sapply(simX, FUN = function(x) x)), if(!is.null(mexsimdata)) t(sapply(mexsimdata, FUN = function(x) x)) else NULL)
			for(i in 1:m.sim) xlist[[i]] = tmp[i, , drop = FALSE]
			simRes = simX
			simX = xlist
		} else{
			if( parallel ){
				if( parallel.control$pkg == "multicore" ){
					xlist = mclapply(as.list(1:m.sim), FUN = function(i) varxsim(X = Data, VAR.fit$Bcoef, p, n.sim, 
										n.start, prereturns, resids = simX[[i]], mexsimdata[[i]]), 
							mc.cores = parallel.control$cores)
				} else{
					# when the GARCH model is without mean equation, simX == simulated residuals
					sfExport("Data", "VAR.fit", "p", "n.sim", "n.start", "prereturns", "simX", "mexsimdata", local = TRUE)
					xlist = sfLapply(as.list(1:m.sim), fun = function(i) rgarch:::varxsim(X = Data, VAR.fit$Bcoef, p, 
										n.sim, n.start, prereturns, resids = simX[[i]], mexsimdata[[i]]))
					sfStop()
				}
			} else{
				xlist = lapply(as.list(1:m.sim), FUN = function(i) varxsim(X = Data, VAR.fit$Bcoef, p, 
									n.sim, n.start, prereturns, resids = simX[[i]], mexsimdata[[i]]))
			}
			simRes = simX
			simX = xlist
		}
	} else{
		# reshape
		for(j in 1:m.sim) simX[[j]] = tail(simX[[j]], n.sim)
	}
	
	msim = list()
	msim$simH = simH
	msim$simR = simR
	msim$simQ = simQ
	msim$simX = simX
	msim$simRes = simRes
	msim$Z = z
	msim$rseed = rseed
	msim$model$Data = NULL
	msim$model$garchpars = garchpars
	msim$model$dccpars = dccpars
	msim$model$n.sim = n.sim
	msim$model$m.sim = m.sim
	msim$model$n.start = n.start
	msim$model$startMethod = "unconditional"
	msim$model$dccOrder = dccOrder
	msim$model$distribution = "mvnorm"
	ans = new("DCCsim",
			msim = msim)
	
	return( ans )
	
}

dccsim1.student = function(fitORspec, n.sim = 1000, n.start = 0, m.sim = 1, startMethod = c("unconditional", "sample"), 
		presigma = NULL, preresiduals = NULL, prereturns = NULL, preR = NULL, preH = NULL, rseed = NULL, mexsimdata = NULL, 
		vexsimdata = NULL, parallel = FALSE, parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2), VAR.fit = NULL, ...)
{
	fit = fitORspec
	Data = head(fit@mfit$origdata, dim(fit@mfit$origdata)[1] - fit@mfit$model$out.sample)
	m = dim(Data)[2]
	n = dim(Data)[1]
	
	dccOrder = fit@mfit$model$dccOrder
	shape = rshape(fit)
	mo = max( dccOrder )
	
	if( is.null(rseed) ){
		rseed = as.integer(runif(1, 1, Sys.time())) 
	} else {
		if(length(rseed) == 1) rseed = as.integer(rseed[1]) else rseed = as.integer( rseed[1:m.sim] )
	}
	
	# we use the GED distribution  with nu = 1 which corresponds to the Laplace case.
	if(length(rseed) == 1){
		tmp = rstd(m * (n.sim + n.start + mo) * m.sim, 0, 1, nu = shape)
		z = array(tmp,  dim = c(n.sim + n.start + mo, m, m.sim))
	} else{
		z = array(NA, dim = c(n.sim + n.start + mo, m, m.sim))
		for(i in 1:m.sim){
			set.seed( rseed[i] )
			z[,,i] = matrix(rstd(m * (n.sim + n.start + mo), 0, 1, nu = shape), nrow = n.sim + n.start + mo, ncol = m)
		}
	}
	if(length(rseed) == 1){
		rseed = c(rseed, as.integer(runif(m.sim, 1, Sys.time()))) 
	}
	if( is.null(preR) ){
		if( startMethod[1] == "sample" ){
			preR = last(rcor(fit))[,,1]
			diag(preR) = 1
		} else{
			preR = cor(as.matrix(Data))
		}
	} else{
		chk = dcc.symcheck(preR, m, d = rep(1, m))
	}
	
	preQ = preR
	
	uncv = sapply(fit@ufit@fit, FUN = function(x) uncvariance(x))
	
	
	if( is.null(preH) ){
		preH = diag(sqrt(uncv)) %*% preR %*% diag(sqrt(uncv))
	} else{
		chk = dcc.symcheck(preH, m, d = NULL)
	}
	
	if( !is.null(presigma) ){
		if( !is.matrix(presigma) ) stop("\ndccsim-->error: presigma must be a matrix.")
		if( dim(presigma)[2] != m ) stop("\ndccsim-->error: wrong column dimension for presigma.")
	}
	
	if( !is.null(preresiduals) ){
		if( !is.matrix(preresiduals) ) stop("\ndccsim-->error: preresiduals must be a matrix.")
		if( dim(preresiduals)[2] != m ) stop("\ndccsim-->error: wrong column dimension for preresiduals.")
	}
	
	if( !is.null(prereturns) ){
		if( !is.matrix(prereturns) ) stop("\ndccsim-->error: prereturns must be a matrix.")
		if( dim(prereturns)[2] != m ) stop("\ndccsim-->error: wrong column dimension for prereturns.")
	}
	
	garchpars	= coef(fit, type = "garch")
	dccpars		= coef(fit, type = "dcc")
	dcca = dccpars[1:dccOrder[1]]
	dccb = dccpars[(dccOrder[1]+1):sum(dccOrder[1:2])]
	simRes = simX = simR = simQ = simH = simSeries = vector(mode = "list", length = m.sim)
	
	if( parallel ){
		os = .Platform$OS.type
		if(is.null(parallel.control$pkg)){
			if( os == "windows" ) parallel.control$pkg = "snowfall" else parallel.control$pkg = "multicore"
			if( is.null(parallel.control$cores) ) parallel.control$cores = 2
		} else{
			mtype = match(tolower(parallel.control$pkg[1]), c("multicore", "snowfall"))
			if(is.na(mtype)) stop("\nParallel Package type not recognized in parallel.control\n")
			parallel.control$pkg = tolower(parallel.control$pkg[1])
			if( os == "windows" && parallel.control$pkg == "multicore" ) stop("\nmulticore not supported on windows O/S\n")
			if( is.null(parallel.control$cores) ) parallel.control$cores = 2 else parallel.control$cores = as.integer(parallel.control$cores[1])
		}
		if( parallel.control$pkg == "multicore" ){
			simH = vector(mode = "list", length = m.sim)
			simX = vector(mode = "list", length = m.sim)
			mtmp = mclapply(as.list(1:m.sim), FUN = function(j)
						.dccsimf(Z = z[,,j], Qbar = preQ, Rbar = preR, dcca, dccb, mo, n.sim, n.start, m, rseed[j]), 
					mc.cores = parallel.control$cores)
			# need to pass m.sim and matrix of simulated z's to ugarchsim (speedup)
			simR = lapply(mtmp, FUN = function(x) if(is.matrix(x$R)) array(x$R, dim = c(m, m, n.sim)) else last(x$R, n.sim))
			simQ = lapply(mtmp, FUN = function(x) if(is.matrix(x$Q)) array(x$Q, dim = c(m, m, n.sim)) else last(x$Q, n.sim))
			stdresid = vector(mode = "list", length = m)
			for(i in 1:m) stdresid[[i]] = sapply(mtmp, FUN = function(x) x$stdresid[,i])			
			xtmp = mclapply(as.list(1:m), FUN = function(j){
						maxx = fit@ufit@fit[[j]]@model$maxOrder2;
						htmp = ugarchsim(fit@ufit@fit[[j]], n.sim = n.sim + n.start, n.start = 0, m.sim = m.sim,
								startMethod = startMethod[1], custom.dist = list(name = "sample", 
										distfit = matrix(stdresid[[j]][-(1:mo), ], ncol = m.sim)),
								presigma = if( is.null(presigma) ) NA else tail(presigma[,j], maxx), 
								preresiduals = if( is.null(preresiduals) ) NA else tail(preresiduals[,j], maxx), 
								prereturns = if( is.null(prereturns) | !is.null(fit@mfit$model$vrmodel) ) NA else tail(prereturns[,j], maxx),
								mexsimdata = if( is.null(fit@mfit$model$vrmodel) ) mexsimdata[[j]] else NULL, vexsimdata = vexsimdata[[j]] )
						h = matrix(tail(htmp@simulation$sigmaSim^2, n.sim), nrow = n.sim)
						x = matrix(htmp@simulation$seriesSim,  nrow = n.sim + n.start)
						return(list(h = h, x = x))}, mc.cores = parallel.control$cores)
			H = array(NA, dim = c(n.sim, m, m.sim))
			tmpH = array(NA, dim = c(m, m, n.sim))
			for(i in 1:n.sim) H[i,,] = t(sapply(xtmp, FUN = function(x) as.numeric(x$h[i,])))
			
			for(i in 1:m.sim){
				for(j in 1:n.sim){
					tmpH[ , , j] = diag(sqrt( H[j, , i]) ) %*% simR[[i]][ , , j] %*% diag(sqrt( H[j, , i] ) )
				}
				simH[[i]] = tmpH
			}
			simxX = array(NA, dim = c(n.sim, m, m.sim))
			for(i in 1:n.sim) simxX[i,,] = t(sapply(xtmp, FUN = function(x) as.numeric(x$x[i,])))
			simX = vector(mode = "list", length = m.sim)
			for(i in 1:m.sim) simX[[i]] = matrix(simxX[,,i], nrow = n.sim)
		} else{
			simH = vector(mode = "list", length = m.sim)
			simX = vector(mode = "list", length = m.sim)
			sfInit(parallel = TRUE, cpus = parallel.control$cores)
			sfExport("z", "preQ", "preR", "dcca", "dccb", "mo", "n.sim", "n.start", "m", "rseed", local = TRUE)
			mtmp = sfLapply(as.list(1:m.sim), fun = function(j)
						rgarch:::.dccsimf(Z = z[,,j], Qbar = preQ, Rbar = preR, dcca, dccb, mo, n.sim, n.start, m, rseed[j]))
			
			# need to pass m.sim and matrix of simulated z's to ugarchsim (speedup)
			simR = lapply(mtmp, FUN = function(x) if(is.matrix(x$R)) array(x$R, dim = c(m, m, n.sim)) else last(x$R, n.sim))
			simQ = lapply(mtmp, FUN = function(x) if(is.matrix(x$Q)) array(x$Q, dim = c(m, m, n.sim)) else last(x$Q, n.sim))
			stdresid = vector(mode = "list", length = m)
			for(i in 1:m) stdresid[[i]] = sapply(mtmp, FUN = function(x) x$stdresid[,i])			
			sfExport("fit", "n.sim", "n.start", "m.sim", "startMethod", "stdresid", "presigma", 
					"preresiduals", "prereturns", "mexsimdata", "vexsimdata", local = TRUE)
			xtmp = sfLapply(as.list(1:m), fun = function(j){
						maxx = fit@ufit@fit[[j]]@model$maxOrder2;
						htmp = rgarch::ugarchsim(fit@ufit@fit[[j]], n.sim = n.sim + n.start, n.start = 0, m.sim = m.sim,
								startMethod = startMethod[1], custom.dist = list(name = "sample", 
										distfit = matrix(stdresid[[j]][-(1:mo), ], ncol = m.sim)),
								presigma = if( is.null(presigma) ) NA else tail(presigma[,j], maxx), 
								preresiduals = if( is.null(preresiduals) ) NA else tail(preresiduals[,j], maxx), 
								prereturns = if( is.null(prereturns) | !is.null(fit@mfit$model$vrmodel) ) NA else tail(prereturns[,j], maxx),
								mexsimdata = if( is.null(fit@mfit$model$vrmodel) ) mexsimdata[[j]] else NULL, vexsimdata = vexsimdata[[j]] )
						h = matrix(tail(htmp@simulation$sigmaSim^2, n.sim), nrow = n.sim)
						x = matrix(htmp@simulation$seriesSim,  nrow = n.sim + n.start)
						return(list(h = h, x = x))})
			if( is.null(fit@mfit$model$vrmodel) ) sfStop()
			H = array(NA, dim = c(n.sim, m, m.sim))
			tmpH = array(NA, dim = c(m, m, n.sim))
			for(i in 1:n.sim) H[i,,] = t(sapply(xtmp, FUN = function(x) as.numeric(x$h[i,])))
			
			for(i in 1:m.sim){
				for(j in 1:n.sim){
					tmpH[ , , j] = diag(sqrt( H[j, , i]) ) %*% simR[[i]][ , , j] %*% diag(sqrt( H[j, , i] ) )
				}
				simH[[i]] = tmpH
			}
			simxX = array(NA, dim = c(n.sim, m, m.sim))
			for(i in 1:n.sim) simxX[i,,] = t(sapply(xtmp, FUN = function(x) as.numeric(x$x[i,])))
			simX = vector(mode = "list", length = m.sim)
			for(i in 1:m.sim) simX[[i]] = matrix(simxX[,,i], nrow = n.sim)
		}
	} else{
		simH = vector(mode = "list", length = m.sim)
		simX = vector(mode = "list", length = m.sim)
		mtmp = lapply(as.list(1:m.sim), FUN = function(j) .dccsimf(Z = z[,,j], Qbar = preQ, Rbar = preR, dcca, dccb, mo, n.sim, n.start, m, rseed[j]))
		# need to pass m.sim and matrix of simulated z's to ugarchsim (speedup)
		simR = lapply(mtmp, FUN = function(x) if(is.matrix(x$R)) array(x$R, dim = c(m, m, n.sim)) else last(x$R, n.sim))
		simQ = lapply(mtmp, FUN = function(x) if(is.matrix(x$Q)) array(x$Q, dim = c(m, m, n.sim)) else last(x$Q, n.sim))
		stdresid = vector(mode = "list", length = m)
		for(i in 1:m) stdresid[[i]] = sapply(mtmp, FUN = function(x) x$stdresid[,i])			
		xtmp = lapply(as.list(1:m), FUN = function(j){
					maxx = fit@ufit@fit[[j]]@model$maxOrder2;
					htmp = ugarchsim(fit@ufit@fit[[j]], n.sim = n.sim + n.start, n.start = 0, m.sim = m.sim,
							startMethod = startMethod[1], custom.dist = list(name = "sample", 
									distfit = matrix(stdresid[[j]][-(1:mo), ], ncol = m.sim)),
							presigma = if( is.null(presigma) ) NA else tail(presigma[,j], maxx), 
							preresiduals = if( is.null(preresiduals) ) NA else tail(preresiduals[,j], maxx), 
							prereturns = if( is.null(prereturns) | !is.null(fit@mfit$model$vrmodel) ) NA else tail(prereturns[,j], maxx),
							mexsimdata = if( is.null(fit@mfit$model$vrmodel) ) mexsimdata[[j]] else NULL, vexsimdata = vexsimdata[[j]] )
					h = matrix(tail(htmp@simulation$sigmaSim^2, n.sim), nrow = n.sim)
					x = matrix(htmp@simulation$seriesSim,  nrow = n.sim + n.start)
					return(list(h = h, x = x))})
		H = array(NA, dim = c(n.sim, m, m.sim))
		tmpH = array(NA, dim = c(m, m, n.sim))
		for(i in 1:n.sim) H[i,,] = t(sapply(xtmp, FUN = function(x) as.numeric(x$h[i,])))	
		for(i in 1:m.sim){
			for(j in 1:n.sim){
				tmpH[ , , j] = diag(sqrt( H[j, , i]) ) %*% simR[[i]][ , , j] %*% diag(sqrt( H[j, , i] ) )
			}
			simH[[i]] = tmpH
		}
		simxX = array(NA, dim = c(n.sim, m, m.sim))
		for(i in 1:n.sim) simxX[i,,] = t(sapply(xtmp, FUN = function(x) as.numeric(x$x[i,])))
		simX = vector(mode = "list", length = m.sim)
		for(i in 1:m.sim) simX[[i]] = matrix(simxX[,,i], nrow = n.sim)
	}
	
	if( !is.null(fit@mfit$model$vrmodel) ){
		# mexsimdata check
		mxn = fit@mfit$model$vrmodel$mxn
		if(mxn > 0 ){
			if( !is.null(mexsimdata) ){
				if( !is.list(mexsimdata) | length(mexsimdata) != m.sim ) stop("\ndccsim-->error: mexsimdata must be a list of length equal to m.sim...", call. = FALSE)
				for(j in 1:m.sim){
					if( !is.matrix(mexsimdata[[j]]) ) stop("\ndccsim-->error: mexsimdata list submatrix must be a matrix...", call. = FALSE)
					if( dim(mexsimdata[[j]])[2] != mxn ) stop("\ndccsim-->error: wrong number of external regressors in list submatrix!...", call. = FALSE)
					if( dim(mexsimdata[[j]])[1] < (n.sim + n.start) ) stop("\ndccsim-->error: external regressor list submatrix has less points than requested simulation length (n.sim + n.start)...", call. = FALSE)
				}
			} else{
				mexsimdata = vector(mode = "list", length = m.sim)
				for(j in 1:m.sim) mexsimdata[[i]] = matrix(0, ncol = mxn, nrow = (n.sim + n.start))
			}
		} else{
			mexsimdata = NULL
		}
		
		vrmodel = fit@mfit$model$vrmodel
		p = as.numeric(vrmodel$lag)
		
		if(startMethod == "sample"){
			if( !is.null(prereturns) ){
				if(!is.matrix(prereturns)) stop("\nprereturns must be a matrix\n")
				if(dim(prereturns)[1] != p) stop("\nwrong number of rows for prereturns matrix")
				if(dim(prereturns)[2] != m) stop("\nwrong number of cols for prereturns matrix")
				prereturns = matrix(prereturns, nrow = p, ncol = m)
			} else{
				prereturns = as.matrix( tail(Data, p) )
			}
		} else{
			prereturns = matrix(0, nrow = p, ncol = m, byrow = TRUE)
		}
		xlist = vector(mode = "list", length = m.sim)
		if( n.sim == 1 && n.start == 0 ){
			# SPECIAL ROUTINE TO SPEED UP 1-AHEAD ESTIMATION
			tmp = varxsimXX(X = Data, vrmodel$Bcoef, p, m.sim = m.sim, prereturns, resids = t(sapply(simX, FUN = function(x) x)), if(!is.null(mexsimdata)) t(sapply(mexsimdata, FUN = function(x) x)) else NULL)
			for(i in 1:m.sim) xlist[[i]] = tmp[i, , drop = FALSE]
			simRes = simX
			simX = xlist
		} else{
			if( parallel ){
				if( parallel.control$pkg == "multicore" ){
					xlist = mclapply(as.list(1:m.sim), FUN = function(i) varxsim(X = Data, vrmodel$Bcoef, p, n.sim, 
										n.start, prereturns, resids = simX[[i]], mexsimdata[[i]]), 
							mc.cores = parallel.control$cores)
				} else{
					# when the GARCH model is without mean equation, simX == simulated residuals
					sfExport("Data", "vrmodel", "p", "n.sim", "n.start", "prereturns", "simX", "mexsimdata", local = TRUE)
					xlist = sfLapply(as.list(1:m.sim), fun = function(i) rgarch:::varxsim(X = Data, vrmodel$Bcoef, p, 
										n.sim, n.start, prereturns, resids = simX[[i]], mexsimdata[[i]]))
					sfStop()
				}
			} else{
				xlist = lapply(as.list(1:m.sim), FUN = function(i) varxsim(X = Data, vrmodel$Bcoef, p, 
									n.sim, n.start, prereturns, resids = simX[[i]], mexsimdata[[i]]))
			}
			simRes = simX
			simX = xlist
		}
	} else{
		# reshape
		for(j in 1:m.sim) simX[[j]] = tail(simX[[j]], n.sim)
	}
	
	msim = list()
	msim$simH = simH
	msim$simR = simR
	msim$simQ = simQ
	msim$simX = simX
	msim$simRes = simRes
	msim$Z = z
	msim$rseed = rseed
	msim$model$Data = Data
	msim$model$garchpars = garchpars
	msim$model$dccpars = dccpars
	msim$model$n.sim = n.sim
	msim$model$m.sim = m.sim
	msim$model$n.start = n.start
	msim$model$startMethod = startMethod[1]
	msim$model$dccOrder = dccOrder
	msim$model$distribution = "mvt"
	ans = new("DCCsim",
			msim = msim)
	
	return( ans )
	
}


dccsim2.student = function(fitORspec, n.sim = 1000, n.start = 0, m.sim = 1, startMethod = c("unconditional", 
				"sample"), presigma = NULL, preresiduals = NULL, 
		prereturns = NULL, preR = NULL, preH = NULL, rseed = NULL, mexsimdata = NULL, 
		vexsimdata = NULL, parallel = FALSE, parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2), VAR.fit = NULL, ...)
{
	spec = fitORspec
	
	# check that VAR.fit is included if used
	if( spec@mspec$varmodel$VAR ){
		if( is.null(VAR.fit) ) stop("\ndccsim-->error: VAR.fit must not be NULL for VAR method when calling dccsim using spec!", call. = FALSE)
	}
	dccOrder = spec@mspec$dccOrder
	mo = max( dccOrder )
	m = length(spec@uspec@spec)
	shape = spec@mspec$optimization.model$fixed.pars["shape"]
	if( is.null(rseed) ){
		rseed = as.integer(runif(1, 1, Sys.time())) 
	} else {
		if(length(rseed) == 1) rseed = as.integer(rseed[1]) else rseed = as.integer( rseed[1:m.sim] )
	}
	
	# we use the GED distribution  with nu = 1 which corresponds to the Laplace case.
	if(length(rseed) == 1){
		tmp = rstd(m * (n.sim + n.start + mo) * m.sim, 0, 1, nu = shape)
		z = array(tmp,  dim = c(n.sim + n.start + mo, m, m.sim))
	} else{
		z = array(NA, dim = c(n.sim + n.start + mo, m, m.sim))
		for(i in 1:m.sim){
			set.seed( rseed[i] )
			z[,,i] = matrix(rstd(m * (n.sim + n.start + mo), 0, 1, nu = shape), nrow = n.sim + n.start + mo, ncol = m)
		}
	}
	if(length(rseed) == 1){
		rseed = c(rseed, as.integer(runif(m.sim, 1, Sys.time()))) 
	}
	if( is.null(preR) ){
		stop("\ndccsim-->error: preR cannot be NULL when method uses only DCC specification!\n")
	} else{
		chk = dcc.symcheck(preR, m, d = rep(1, m))
	}
	preQ = preR
	
	# uncvariance has uGARCHfit AND uGARCHspec dispatch methods
	uncv = sapply(spec@uspec@spec, FUN = function(x) uncvariance(x))
	
	
	if( is.null(preH) ){
		preH = diag(sqrt(uncv)) %*% preR %*% diag(sqrt(uncv))
	} else{
		chk = dcc.symcheck(preH, m, d = NULL)
	}
	if( !is.null(presigma) ){
		if( !is.matrix(presigma) ) stop("\ndccsim-->error: presigma must be a matrix.")
		if( dim(presigma)[2] != m ) stop("\ndccsim-->error: wrong column dimension for presigma.")
	}
	
	if( !is.null(preresiduals) ){
		if( !is.matrix(preresiduals) ) stop("\ndccsim-->error: preresiduals must be a matrix.")
		if( dim(preresiduals)[2] != m ) stop("\ndccsim-->error: wrong column dimension for preresiduals.")
	}
	
	if( !is.null(prereturns) ){
		if( !is.matrix(prereturns) ) stop("\ndccsim-->error: prereturns must be a matrix.")
		if( dim(prereturns)[2] != m ) stop("\ndccsim-->error: wrong column dimension for prereturns.")
	}
	
	garchpars	= sapply(spec@uspec@spec, FUN = function(x) unlist(x@optimization.model$fixed.pars))
	dccpars		= unlist(spec@mspec$optimization.model$fixed.pars)
	dcca = dccpars[1:dccOrder[1]]
	dccb = dccpars[(dccOrder[1]+1):sum(dccOrder[1:2])]
	simRes = simX = simR = simQ = simH = simSeries = vector(mode = "list", length = m.sim)
	
	if( parallel ){
		os = .Platform$OS.type
		if(is.null(parallel.control$pkg)){
			if( os == "windows" ) parallel.control$pkg = "snowfall" else parallel.control$pkg = "multicore"
			if( is.null(parallel.control$cores) ) parallel.control$cores = 2
		} else{
			mtype = match(tolower(parallel.control$pkg[1]), c("multicore", "snowfall"))
			if(is.na(mtype)) stop("\nParallel Package type not recognized in parallel.control\n")
			parallel.control$pkg = tolower(parallel.control$pkg[1])
			if( os == "windows" && parallel.control$pkg == "multicore" ) stop("\nmulticore not supported on windows O/S\n")
			if( is.null(parallel.control$cores) ) parallel.control$cores = 2 else parallel.control$cores = as.integer(parallel.control$cores[1])
		}
		if( parallel.control$pkg == "multicore" ){
			simH = vector(mode = "list", length = m.sim)
			simX = vector(mode = "list", length = m.sim)
			mtmp = mclapply(as.list(1:m.sim), FUN = function(j)
						.dccsimf(Z = z[,,j], Qbar = preQ, Rbar = preR, dcca, dccb, mo, n.sim, n.start, m, rseed[j]), 
					mc.cores = parallel.control$cores)
			simR = lapply(mtmp, FUN = function(x) if(is.matrix(x$R)) array(x$R, dim = c(m, m, n.sim)) else last(x$R, n.sim))
			simQ = lapply(mtmp, FUN = function(x) if(is.matrix(x$Q)) array(x$Q, dim = c(m, m, n.sim)) else last(x$Q, n.sim))
			stdresid = vector(mode = "list", length = m)
			for(i in 1:m) stdresid[[i]] = sapply(mtmp, FUN = function(x) x$stdresid[,i])			
			xtmp = mclapply(as.list(1:m), FUN = function(j){
						maxx = spec@uspec@spec[[j]]@optimization.model$maxOrder2;
						htmp = ugarchpath(spec@uspec@spec[[j]], n.sim = n.sim + n.start, n.start = 0, m.sim = m.sim,
								custom.dist = list(name = "sample", distfit = matrix(stdresid[[j]][-(1:mo), ], ncol = m.sim)),
								presigma = if( is.null(presigma) ) NA else tail(presigma[,j], maxx), 
								preresiduals = if( is.null(preresiduals) ) NA else tail(preresiduals[,j], maxx), 
								prereturns = if( is.null(prereturns) | spec@mspec$varmodel$VAR ) NA else tail(prereturns[,j], maxx),
								mexsimdata = if( !spec@mspec$varmodel$VAR ) mexsimdata[[j]] else NULL, vexsimdata = vexsimdata[[j]] )
						h = matrix(tail(htmp@path$sigmaSim^2, n.sim), nrow = n.sim)
						x = matrix(htmp@path$seriesSim,  nrow = n.sim + n.start)
						return(list(h = h, x = x))}, mc.cores = parallel.control$cores)
			H = array(NA, dim = c(n.sim, m, m.sim))
			tmpH = array(NA, dim = c(m, m, n.sim))
			for(i in 1:n.sim) H[i,,] = t(sapply(xtmp, FUN = function(x) as.numeric(x$h[i,])))
			
			for(i in 1:m.sim){
				for(j in 1:n.sim){
					tmpH[ , , j] = diag(sqrt( H[j, , i]) ) %*% simR[[i]][ , , j] %*% diag(sqrt( H[j, , i] ) )
				}
				simH[[i]] = tmpH
			}
			simxX = array(NA, dim = c(n.sim, m, m.sim))
			for(i in 1:n.sim) simxX[i,,] = t(sapply(xtmp, FUN = function(x) as.numeric(x$x[i,])))
			simX = vector(mode = "list", length = m.sim)
			for(i in 1:m.sim) simX[[i]] = matrix(simxX[,,i], nrow = n.sim)
		} else{
			sfInit(parallel = TRUE, cpus = parallel.control$cores)
			sfExport("spec", "z", "preQ", "preR", "dcca", "dccb", "mo", "n.sim", "n.start", 
					"m", "prereturns", "preresiduals", "presigma", "mexsimdata", "vexsimdata", 
					"rseed", local = TRUE)
			simH = vector(mode = "list", length = m.sim)
			simX = vector(mode = "list", length = m.sim)
			mtmp = sfLapply(as.list(1:m.sim), fun = function(j) .dccsimf(Z = z[,,j], Qbar = preQ, Rbar = preR, dcca, dccb, mo, n.sim, n.start, m, rseed[j]))
			simR = lapply(mtmp, FUN = function(x) if(is.matrix(x$R)) array(x$R, dim = c(m, m, n.sim)) else last(x$R, n.sim))
			simQ = lapply(mtmp, FUN = function(x) if(is.matrix(x$Q)) array(x$Q, dim = c(m, m, n.sim)) else last(x$Q, n.sim))
			stdresid = vector(mode = "list", length = m)
			for(i in 1:m) stdresid[[i]] = sapply(mtmp, FUN = function(x) x$stdresid[,i])			
			xtmp = sfLapply(as.list(1:m), fun = function(j){
						maxx = spec@uspec@spec[[j]]@optimization.model$maxOrder2;
						htmp = rgarch::ugarchpath(spec@uspec@spec[[j]], n.sim = n.sim + n.start, n.start = 0, m.sim = m.sim,
								custom.dist = list(name = "sample", distfit = matrix(stdresid[[j]][-(1:mo), ], ncol = m.sim)),
								presigma = if( is.null(presigma) ) NA else tail(presigma[,j], maxx), 
								preresiduals = if( is.null(preresiduals) ) NA else tail(preresiduals[,j], maxx), 
								prereturns = if( is.null(prereturns) | spec@mspec$varmodel$VAR ) NA else tail(prereturns[,j], maxx),
								mexsimdata = if( !spec@mspec$varmodel$VAR ) mexsimdata[[j]] else NULL, vexsimdata = vexsimdata[[j]] )
						h = matrix(tail(htmp@path$sigmaSim^2, n.sim), nrow = n.sim)
						x = matrix(htmp@path$seriesSim,  nrow = n.sim + n.start)
						return(list(h = h, x = x))})
			if( !spec@mspec$varmodel$VAR ) sfStop()
			H = array(NA, dim = c(n.sim, m, m.sim))
			tmpH = array(NA, dim = c(m, m, n.sim))
			for(i in 1:n.sim) H[i,,] = t(sapply(xtmp, FUN = function(x) as.numeric(x$h[i,])))
			
			for(i in 1:m.sim){
				for(j in 1:n.sim){
					tmpH[ , , j] = diag(sqrt( H[j, , i]) ) %*% simR[[i]][ , , j] %*% diag(sqrt( H[j, , i] ) )
				}
				simH[[i]] = tmpH
			}
			simxX = array(NA, dim = c(n.sim, m, m.sim))
			for(i in 1:n.sim) simxX[i,,] = t(sapply(xtmp, FUN = function(x) as.numeric(x$x[i,])))
			simX = vector(mode = "list", length = m.sim)
			for(i in 1:m.sim) simX[[i]] = matrix(simxX[,,i], nrow = n.sim)
		}
	} else{
		simH = vector(mode = "list", length = m.sim)
		simX = vector(mode = "list", length = m.sim)
		mtmp = lapply(as.list(1:m.sim), FUN = function(j) .dccsimf(Z = z[,,j], Qbar = preQ, Rbar = preR, dcca, dccb, mo, n.sim, n.start, m, rseed[j]))
		simR = lapply(mtmp, FUN = function(x) if(is.matrix(x$R)) array(x$R, dim = c(m, m, n.sim)) else last(x$R, n.sim))
		simQ = lapply(mtmp, FUN = function(x) if(is.matrix(x$Q)) array(x$Q, dim = c(m, m, n.sim)) else last(x$Q, n.sim))
		stdresid = vector(mode = "list", length = m)
		for(i in 1:m) stdresid[[i]] = sapply(mtmp, FUN = function(x) x$stdresid[,i])			
		xtmp = lapply(as.list(1:m), FUN = function(j){
					maxx = spec@uspec@spec[[j]]@optimization.model$maxOrder2;
					htmp = ugarchpath(spec@uspec@spec[[j]], n.sim = n.sim + n.start, n.start = 0, m.sim = m.sim,
							custom.dist = list(name = "sample", distfit = matrix(stdresid[[j]][-(1:mo), ], ncol = m.sim)),
							presigma = if( is.null(presigma) ) NA else tail(presigma[,j], maxx), 
							preresiduals = if( is.null(preresiduals) ) NA else tail(preresiduals[,j], maxx), 
							prereturns = if( is.null(prereturns) | spec@mspec$varmodel$VAR ) NA else tail(prereturns[,j], maxx),
							mexsimdata = if( !spec@mspec$varmodel$VAR ) mexsimdata[[j]] else NULL, vexsimdata = vexsimdata[[j]] )
					h = matrix(tail(htmp@path$sigmaSim^2, n.sim), nrow = n.sim)
					x = matrix(htmp@path$seriesSim,  nrow = n.sim + n.start)
					return(list(h = h, x = x))})
		H = array(NA, dim = c(n.sim, m, m.sim))
		tmpH = array(NA, dim = c(m, m, n.sim))
		for(i in 1:n.sim) H[i,,] = t(sapply(xtmp, FUN = function(x) as.numeric(x$h[i,])))
		
		for(i in 1:m.sim){
			for(j in 1:n.sim){
				tmpH[ , , j] = diag(sqrt( H[j, , i]) ) %*% simR[[i]][ , , j] %*% diag(sqrt( H[j, , i] ) )
			}
			simH[[i]] = tmpH
		}
		simxX = array(NA, dim = c(n.sim, m, m.sim))
		for(i in 1:n.sim) simxX[i,,] = t(sapply(xtmp, FUN = function(x) as.numeric(x$x[i,])))
		simX = vector(mode = "list", length = m.sim)
		for(i in 1:m.sim) simX[[i]] = matrix(simxX[,,i], nrow = n.sim)
	}
	if( spec@mspec$varmodel$VAR ){
		# mexsimdata check
		mxn = VAR.fit$mxn
		Data = VAR.fit$xfitted
		if(mxn > 0 ){
			if( !is.null(mexsimdata) ){
				if( !is.list(mexsimdata) | length(mexsimdata) != m.sim ) stop("\ndccsim-->error: mexsimdata must be a list of length equal to m.sim...", call. = FALSE)
				for(j in 1:m.sim){
					if( !is.matrix(mexsimdata[[j]]) ) stop("\ndccsim-->error: mexsimdata list submatrix must be a matrix...", call. = FALSE)
					if( dim(mexsimdata[[j]])[2] != mxn ) stop("\ndccsim-->error: wrong number of external regressors in list submatrix!...", call. = FALSE)
					if( dim(mexsimdata[[j]])[1] < (n.sim + n.start) ) stop("\ndccsim-->error: external regressor list submatrix has less points than requested simulation length (n.sim + n.start)...", call. = FALSE)
				}
			} else{
				mexsimdata = vector(mode = "list", length = m.sim)
				for(j in 1:m.sim) mexsimdata[[i]] = matrix(0, ncol = mxn, nrow = (n.sim + n.start))
			}
		} else{
			mexsimdata = NULL
		}
		
		p = VAR.fit$lag
		
		if( !is.null(prereturns) ){
			if(!is.matrix(prereturns)) stop("\nprereturns must be a matrix\n")
			if(dim(prereturns)[1] != p) stop("\nwrong number of rows for prereturns matrix")
			if(dim(prereturns)[2] != m) stop("\nwrong number of cols for prereturns matrix")
			prereturns = matrix(prereturns, nrow = p, ncol = m)
		} else{
			prereturns = matrix(0, nrow = p, ncol = m, byrow = TRUE)
		}
		xlist = vector(mode = "list", length = m.sim)
		if( n.sim == 1 && n.start == 0 ){
			# SPECIAL ROUTINE TO SPEED UP 1-AHEAD ESTIMATION
			tmp = varxsimXX(X = Data, VAR.fit$Bcoef, p, m.sim = m.sim, prereturns, resids = t(sapply(simX, FUN = function(x) x)), if(!is.null(mexsimdata)) t(sapply(mexsimdata, FUN = function(x) x)) else NULL)
			for(i in 1:m.sim) xlist[[i]] = tmp[i, , drop = FALSE]
			simRes = simX
			simX = xlist
		} else{
			if( parallel ){
				if( parallel.control$pkg == "multicore" ){
					xlist = mclapply(as.list(1:m.sim), FUN = function(i) varxsim(X = Data, VAR.fit$Bcoef, p, n.sim, 
										n.start, prereturns, resids = simX[[i]], mexsimdata[[i]]), 
							mc.cores = parallel.control$cores)
				} else{
					# when the GARCH model is without mean equation, simX == simulated residuals
					sfExport("Data", "VAR.fit", "p", "n.sim", "n.start", "prereturns", "simX", "mexsimdata", local = TRUE)
					xlist = sfLapply(as.list(1:m.sim), fun = function(i) rgarch:::varxsim(X = Data, VAR.fit$Bcoef, p, 
										n.sim, n.start, prereturns, resids = simX[[i]], mexsimdata[[i]]))
					sfStop()
				}
			} else{
				xlist = lapply(as.list(1:m.sim), FUN = function(i) varxsim(X = Data, VAR.fit$Bcoef, p, 
									n.sim, n.start, prereturns, resids = simX[[i]], mexsimdata[[i]]))
			}
			simRes = simX
			simX = xlist
		}
	} else{
		# reshape
		for(j in 1:m.sim) simX[[j]] = tail(simX[[j]], n.sim)
	}
	
	msim = list()
	msim$simH = simH
	msim$simR = simR
	msim$simQ = simQ
	msim$simX = simX
	msim$simRes = simRes
	msim$Z = z
	msim$rseed = rseed
	msim$model$Data = NULL
	msim$model$garchpars = garchpars
	msim$model$dccpars = dccpars
	msim$model$n.sim = n.sim
	msim$model$m.sim = m.sim
	msim$model$n.start = n.start
	msim$model$startMethod = "unconditional"
	msim$model$dccOrder = dccOrder
	msim$model$distribution = "mvt"
	ans = new("DCCsim",
			msim = msim)
	
	return( ans )
	
}

dccsim1.laplace = function(fitORspec, n.sim = 1000, n.start = 0, m.sim = 1, startMethod = c("unconditional", "sample"), 
		presigma = NULL, preresiduals = NULL, prereturns = NULL, preR = NULL, preH = NULL, rseed = NULL, mexsimdata = NULL, 
		vexsimdata = NULL, parallel = FALSE, parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2), VAR.fit = NULL, ...)
{
	fit = fitORspec
	Data = head(fit@mfit$origdata, dim(fit@mfit$origdata)[1] - fit@mfit$model$out.sample)
	m = dim(Data)[2]
	n = dim(Data)[1]
	
	dccOrder = fit@mfit$model$dccOrder
	mo = max( dccOrder )
	if( is.null(rseed) ){
		rseed = as.integer(runif(1, 1, Sys.time())) 
	} else {
		if(length(rseed) == 1) rseed = as.integer(rseed[1]) else rseed = as.integer( rseed[1:m.sim] )
	}
	
	# we use the GED distribution  with nu = 1 which corresponds to the Laplace case.
	if(length(rseed) == 1){
		tmp = rged(m * (n.sim + n.start + mo) * m.sim, 0, 1, nu = 1)
		z = array(tmp,  dim = c(n.sim + n.start + mo, m, m.sim))
	} else{
		z = array(NA, dim = c(n.sim + n.start + mo, m, m.sim))
		for(i in 1:m.sim){
			set.seed( rseed[i] )
			z[,,i] = matrix(rged(m * (n.sim + n.start + mo), 0, 1, nu = 1), nrow = n.sim + n.start + mo, ncol = m)
		}
	}
	if(length(rseed) == 1){
		rseed = c(rseed, as.integer(runif(m.sim, 1, Sys.time()))) 
	}
	if( is.null(preR) ){
		if( startMethod[1] == "sample" ){
			preR = last(rcor(fit))[,,1]
			diag(preR) = 1
		} else{
			preR = cor(as.matrix(Data))
		}
	} else{
		chk = dcc.symcheck(preR, m, d = rep(1, m))
	}
	
	preQ = preR
	
	uncv = sapply(fit@ufit@fit, FUN = function(x) uncvariance(x))
	
	
	if( is.null(preH) ){
		preH = diag(sqrt(uncv)) %*% preR %*% diag(sqrt(uncv))
	} else{
		chk = dcc.symcheck(preH, m, d = NULL)
	}
	if( !is.null(presigma) ){
		if( !is.matrix(presigma) ) stop("\ndccsim-->error: presigma must be a matrix.")
		if( dim(presigma)[2] != m ) stop("\ndccsim-->error: wrong column dimension for presigma.")
	}
	
	if( !is.null(preresiduals) ){
		if( !is.matrix(preresiduals) ) stop("\ndccsim-->error: preresiduals must be a matrix.")
		if( dim(preresiduals)[2] != m ) stop("\ndccsim-->error: wrong column dimension for preresiduals.")
	}
	
	if( !is.null(prereturns) ){
		if( !is.matrix(prereturns) ) stop("\ndccsim-->error: prereturns must be a matrix.")
		if( dim(prereturns)[2] != m ) stop("\ndccsim-->error: wrong column dimension for prereturns.")
	}
	garchpars	= coef(fit, type = "garch")
	dccpars		= coef(fit, type = "dcc")
	dcca = dccpars[1:dccOrder[1]]
	dccb = dccpars[(dccOrder[1]+1):sum(dccOrder[1:2])]
	simRes = simX = simR = simQ = simH = simSeries = vector(mode = "list", length = m.sim)
	
	if( parallel ){
		os = .Platform$OS.type
		if(is.null(parallel.control$pkg)){
			if( os == "windows" ) parallel.control$pkg = "snowfall" else parallel.control$pkg = "multicore"
			if( is.null(parallel.control$cores) ) parallel.control$cores = 2
		} else{
			mtype = match(tolower(parallel.control$pkg[1]), c("multicore", "snowfall"))
			if(is.na(mtype)) stop("\nParallel Package type not recognized in parallel.control\n")
			parallel.control$pkg = tolower(parallel.control$pkg[1])
			if( os == "windows" && parallel.control$pkg == "multicore" ) stop("\nmulticore not supported on windows O/S\n")
			if( is.null(parallel.control$cores) ) parallel.control$cores = 2 else parallel.control$cores = as.integer(parallel.control$cores[1])
		}
		if( parallel.control$pkg == "multicore" ){
			simH = vector(mode = "list", length = m.sim)
			simX = vector(mode = "list", length = m.sim)
			mtmp = mclapply(as.list(1:m.sim), FUN = function(j)
						.dccsimf(Z = z[,,j], Qbar = preQ, Rbar = preR, dcca, dccb, mo, n.sim, n.start, m, rseed[j]), 
					mc.cores = parallel.control$cores)
			# need to pass m.sim and matrix of simulated z's to ugarchsim (speedup)
			simR = lapply(mtmp, FUN = function(x) if(is.matrix(x$R)) array(x$R, dim = c(m, m, n.sim)) else last(x$R, n.sim))
			simQ = lapply(mtmp, FUN = function(x) if(is.matrix(x$Q)) array(x$Q, dim = c(m, m, n.sim)) else last(x$Q, n.sim))
			stdresid = vector(mode = "list", length = m)
			for(i in 1:m) stdresid[[i]] = sapply(mtmp, FUN = function(x) x$stdresid[,i])			
			xtmp = mclapply(as.list(1:m), FUN = function(j){
						maxx = fit@ufit@fit[[j]]@model$maxOrder2;
						htmp = ugarchsim(fit@ufit@fit[[j]], n.sim = n.sim + n.start, n.start = 0, m.sim = m.sim,
								startMethod = startMethod[1], custom.dist = list(name = "sample", 
										distfit = matrix(stdresid[[j]][-(1:mo), ], ncol = m.sim)),
								presigma = if( is.null(presigma) ) NA else tail(presigma[,j], maxx), 
								preresiduals = if( is.null(preresiduals) ) NA else tail(preresiduals[,j], maxx), 
								prereturns = if( is.null(prereturns) | !is.null(fit@mfit$model$vrmodel) ) NA else tail(prereturns[,j], maxx),
								mexsimdata = if( is.null(fit@mfit$model$vrmodel) ) mexsimdata[[j]] else NULL, vexsimdata = vexsimdata[[j]] )
						h = matrix(tail(htmp@simulation$sigmaSim^2, n.sim), nrow = n.sim)
						x = matrix(htmp@simulation$seriesSim,  nrow = n.sim + n.start)
						return(list(h = h, x = x))}, mc.cores = parallel.control$cores)
			H = array(NA, dim = c(n.sim, m, m.sim))
			tmpH = array(NA, dim = c(m, m, n.sim))
			for(i in 1:n.sim) H[i,,] = t(sapply(xtmp, FUN = function(x) as.numeric(x$h[i,])))
			
			for(i in 1:m.sim){
				for(j in 1:n.sim){
					tmpH[ , , j] = diag(sqrt( H[j, , i]) ) %*% simR[[i]][ , , j] %*% diag(sqrt( H[j, , i] ) )
				}
				simH[[i]] = tmpH
			}
			simxX = array(NA, dim = c(n.sim, m, m.sim))
			for(i in 1:n.sim) simxX[i,,] = t(sapply(xtmp, FUN = function(x) as.numeric(x$x[i,])))
			simX = vector(mode = "list", length = m.sim)
			for(i in 1:m.sim) simX[[i]] = matrix(simxX[,,i], nrow = n.sim)
		} else{
			simH = vector(mode = "list", length = m.sim)
			simX = vector(mode = "list", length = m.sim)
			sfInit(parallel = TRUE, cpus = parallel.control$cores)
			sfExport("z", "preQ", "preR", "dcca", "dccb", "mo", "n.sim", "n.start", "m", "rseed", local = TRUE)
			mtmp = sfLapply(as.list(1:m.sim), fun = function(j)
						rgarch:::.dccsimf(Z = z[,,j], Qbar = preQ, Rbar = preR, dcca, dccb, mo, n.sim, n.start, m, rseed[j]))
			
			# need to pass m.sim and matrix of simulated z's to ugarchsim (speedup)
			simR = lapply(mtmp, FUN = function(x) if(is.matrix(x$R)) array(x$R, dim = c(m, m, n.sim)) else last(x$R, n.sim))
			simQ = lapply(mtmp, FUN = function(x) if(is.matrix(x$Q)) array(x$Q, dim = c(m, m, n.sim)) else last(x$Q, n.sim))
			stdresid = vector(mode = "list", length = m)
			for(i in 1:m) stdresid[[i]] = sapply(mtmp, FUN = function(x) x$stdresid[,i])			
			sfExport("fit", "n.sim", "n.start", "m.sim", "startMethod", "stdresid", "presigma", 
					"preresiduals", "prereturns", "mexsimdata", "vexsimdata", local = TRUE)
			xtmp = sfLapply(as.list(1:m), fun = function(j){
						maxx = fit@ufit@fit[[j]]@model$maxOrder2;
						htmp = rgarch::ugarchsim(fit@ufit@fit[[j]], n.sim = n.sim + n.start, n.start = 0, m.sim = m.sim,
								startMethod = startMethod[1], custom.dist = list(name = "sample", 
										distfit = matrix(stdresid[[j]][-(1:mo), ], ncol = m.sim)),
								presigma = if( is.null(presigma) ) NA else tail(presigma[,j], maxx), 
								preresiduals = if( is.null(preresiduals) ) NA else tail(preresiduals[,j], maxx), 
								prereturns = if( is.null(prereturns) | !is.null(fit@mfit$model$vrmodel) ) NA else tail(prereturns[,j], maxx),
								mexsimdata = if( is.null(fit@mfit$model$vrmodel) ) mexsimdata[[j]] else NULL, vexsimdata = vexsimdata[[j]] )
						h = matrix(tail(htmp@simulation$sigmaSim^2, n.sim), nrow = n.sim)
						x = matrix(htmp@simulation$seriesSim,  nrow = n.sim + n.start)
						return(list(h = h, x = x))})
			if( is.null(fit@mfit$model$vrmodel) ) sfStop()
			H = array(NA, dim = c(n.sim, m, m.sim))
			tmpH = array(NA, dim = c(m, m, n.sim))
			for(i in 1:n.sim) H[i,,] = t(sapply(xtmp, FUN = function(x) as.numeric(x$h[i,])))
			
			for(i in 1:m.sim){
				for(j in 1:n.sim){
					tmpH[ , , j] = diag(sqrt( H[j, , i]) ) %*% simR[[i]][ , , j] %*% diag(sqrt( H[j, , i] ) )
				}
				simH[[i]] = tmpH
			}
			simxX = array(NA, dim = c(n.sim, m, m.sim))
			for(i in 1:n.sim) simxX[i,,] = t(sapply(xtmp, FUN = function(x) as.numeric(x$x[i,])))
			simX = vector(mode = "list", length = m.sim)
			for(i in 1:m.sim) simX[[i]] = matrix(simxX[,,i], nrow = n.sim)
		}
	} else{
		simH = vector(mode = "list", length = m.sim)
		simX = vector(mode = "list", length = m.sim)
		mtmp = lapply(as.list(1:m.sim), FUN = function(j) .dccsimf(Z = z[,,j], Qbar = preQ, Rbar = preR, dcca, dccb, mo, n.sim, n.start, m, rseed[j]))
		# need to pass m.sim and matrix of simulated z's to ugarchsim (speedup)
		simR = lapply(mtmp, FUN = function(x) if(is.matrix(x$R)) array(x$R, dim = c(m, m, n.sim)) else last(x$R, n.sim))
		simQ = lapply(mtmp, FUN = function(x) if(is.matrix(x$Q)) array(x$Q, dim = c(m, m, n.sim)) else last(x$Q, n.sim))
		stdresid = vector(mode = "list", length = m)
		for(i in 1:m) stdresid[[i]] = sapply(mtmp, FUN = function(x) x$stdresid[,i])			
		xtmp = lapply(as.list(1:m), FUN = function(j){
					maxx = fit@ufit@fit[[j]]@model$maxOrder2;
					htmp = ugarchsim(fit@ufit@fit[[j]], n.sim = n.sim + n.start, n.start = 0, m.sim = m.sim,
							startMethod = startMethod[1], custom.dist = list(name = "sample", 
									distfit = matrix(stdresid[[j]][-(1:mo), ], ncol = m.sim)),
							presigma = if( is.null(presigma) ) NA else tail(presigma[,j], maxx), 
							preresiduals = if( is.null(preresiduals) ) NA else tail(preresiduals[,j], maxx), 
							prereturns = if( is.null(prereturns) | !is.null(fit@mfit$model$vrmodel) ) NA else tail(prereturns[,j], maxx),
							mexsimdata = if( is.null(fit@mfit$model$vrmodel) ) mexsimdata[[j]] else NULL, vexsimdata = vexsimdata[[j]] )
					h = matrix(tail(htmp@simulation$sigmaSim^2, n.sim), nrow = n.sim)
					x = matrix(htmp@simulation$seriesSim,  nrow = n.sim + n.start)
					return(list(h = h, x = x))})
		H = array(NA, dim = c(n.sim, m, m.sim))
		tmpH = array(NA, dim = c(m, m, n.sim))
		for(i in 1:n.sim) H[i,,] = t(sapply(xtmp, FUN = function(x) as.numeric(x$h[i,])))	
		for(i in 1:m.sim){
			for(j in 1:n.sim){
				tmpH[ , , j] = diag(sqrt( H[j, , i]) ) %*% simR[[i]][ , , j] %*% diag(sqrt( H[j, , i] ) )
			}
			simH[[i]] = tmpH
		}
		simxX = array(NA, dim = c(n.sim, m, m.sim))
		for(i in 1:n.sim) simxX[i,,] = t(sapply(xtmp, FUN = function(x) as.numeric(x$x[i,])))
		simX = vector(mode = "list", length = m.sim)
		for(i in 1:m.sim) simX[[i]] = matrix(simxX[,,i], nrow = n.sim)
	}
	
	if( !is.null(fit@mfit$model$vrmodel) ){
		# mexsimdata check
		mxn = fit@mfit$model$vrmodel$mxn
		if(mxn > 0 ){
			if( !is.null(mexsimdata) ){
				if( !is.list(mexsimdata) | length(mexsimdata) != m.sim ) stop("\ndccsim-->error: mexsimdata must be a list of length equal to m.sim...", call. = FALSE)
				for(j in 1:m.sim){
					if( !is.matrix(mexsimdata[[j]]) ) stop("\ndccsim-->error: mexsimdata list submatrix must be a matrix...", call. = FALSE)
					if( dim(mexsimdata[[j]])[2] != mxn ) stop("\ndccsim-->error: wrong number of external regressors in list submatrix!...", call. = FALSE)
					if( dim(mexsimdata[[j]])[1] < (n.sim + n.start) ) stop("\ndccsim-->error: external regressor list submatrix has less points than requested simulation length (n.sim + n.start)...", call. = FALSE)
				}
			} else{
				mexsimdata = vector(mode = "list", length = m.sim)
				for(j in 1:m.sim) mexsimdata[[i]] = matrix(0, ncol = mxn, nrow = (n.sim + n.start))
			}
		} else{
			mexsimdata = NULL
		}
		
		vrmodel = fit@mfit$model$vrmodel
		p = as.numeric(vrmodel$lag)
		
		if(startMethod == "sample"){
			if( !is.null(prereturns) ){
				if(!is.matrix(prereturns)) stop("\nprereturns must be a matrix\n")
				if(dim(prereturns)[1] != p) stop("\nwrong number of rows for prereturns matrix")
				if(dim(prereturns)[2] != m) stop("\nwrong number of cols for prereturns matrix")
				prereturns = matrix(prereturns, nrow = p, ncol = m)
			} else{
				prereturns = as.matrix( tail(Data, p) )
			}
		} else{
			prereturns = matrix(0, nrow = p, ncol = m, byrow = TRUE)
		}
		xlist = vector(mode = "list", length = m.sim)
		if( n.sim == 1 && n.start == 0 ){
			# SPECIAL ROUTINE TO SPEED UP 1-AHEAD ESTIMATION
			tmp = varxsimXX(X = Data, vrmodel$Bcoef, p, m.sim = m.sim, prereturns, resids = t(sapply(simX, FUN = function(x) x)), if(!is.null(mexsimdata)) t(sapply(mexsimdata, FUN = function(x) x)) else NULL)
			for(i in 1:m.sim) xlist[[i]] = tmp[i, , drop = FALSE]
			simRes = simX
			simX = xlist
		} else{
			if( parallel ){
				if( parallel.control$pkg == "multicore" ){
					xlist = mclapply(as.list(1:m.sim), FUN = function(i) varxsim(X = Data, vrmodel$Bcoef, p, n.sim, 
										n.start, prereturns, resids = simX[[i]], mexsimdata[[i]]), 
							mc.cores = parallel.control$cores)
				} else{
					# when the GARCH model is without mean equation, simX == simulated residuals
					sfExport("Data", "vrmodel", "p", "n.sim", "n.start", "prereturns", "simX", "mexsimdata", local = TRUE)
					xlist = sfLapply(as.list(1:m.sim), fun = function(i) rgarch:::varxsim(X = Data, vrmodel$Bcoef, p, 
										n.sim, n.start, prereturns, resids = simX[[i]], mexsimdata[[i]]))
					sfStop()
				}
			} else{
				xlist = lapply(as.list(1:m.sim), FUN = function(i) varxsim(X = Data, vrmodel$Bcoef, p, 
									n.sim, n.start, prereturns, resids = simX[[i]], mexsimdata[[i]]))
			}
			simRes = simX
			simX = xlist
		}
	} else{
		# reshape
		for(j in 1:m.sim) simX[[j]] = tail(simX[[j]], n.sim)
	}
	
	msim = list()
	msim$simH = simH
	msim$simR = simR
	msim$simQ = simQ
	msim$simX = simX
	msim$simRes = simRes
	msim$Z = z
	msim$rseed = rseed
	msim$model$Data = Data
	msim$model$garchpars = garchpars
	msim$model$dccpars = dccpars
	msim$model$n.sim = n.sim
	msim$model$m.sim = m.sim
	msim$model$n.start = n.start
	msim$model$startMethod = startMethod[1]
	msim$model$dccOrder = dccOrder
	msim$model$distribution = "mvlaplace"
	ans = new("DCCsim",
			msim = msim)
	
	return( ans )
	
}

dccsim2.laplace = function(fitORspec, n.sim = 1000, n.start = 0, m.sim = 1,  startMethod = c("unconditional", 
				"sample"), presigma = NULL, preresiduals = NULL, 
		prereturns = NULL, preR = NULL, preH = NULL, rseed = NULL, mexsimdata = NULL, vexsimdata = NULL, 
		parallel = FALSE, parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2), VAR.fit = NULL, ...)
{
	spec = fitORspec
	
	# check that VAR.fit is included if used
	if( spec@mspec$varmodel$VAR ){
		if( is.null(VAR.fit) ) stop("\ndccsim-->error: VAR.fit must not be NULL for VAR method when calling dccsim using spec!", call. = FALSE)
	}
	
	dccOrder = spec@mspec$dccOrder
	mo = max( dccOrder )
	m = length(spec@uspec@spec)
	
	if( is.null(rseed) ){
		rseed = as.integer(runif(1, 1, Sys.time())) 
	} else {
		if(length(rseed) == 1) rseed = as.integer(rseed[1]) else rseed = as.integer( rseed[1:m.sim] )
	}
	# we use the GED distribution  with nu = 1 which corresponds to the Laplace case.
	if(length(rseed) == 1){
		tmp = rged(m * (n.sim + n.start + mo) * m.sim, 0, 1, nu = 1)
		z = array(tmp,  dim = c(n.sim + n.start + mo, m, m.sim))
	} else{
		z = array(NA, dim = c(n.sim + n.start + mo, m, m.sim))
		for(i in 1:m.sim){
			set.seed( rseed[i] )
			z[,,i] = matrix(rged(m * (n.sim + n.start + mo), 0, 1, nu = 1), nrow = n.sim + n.start + mo, ncol = m)
		}
	}
	if(length(rseed) == 1){
		rseed = c(rseed, as.integer(runif(m.sim, 1, Sys.time()))) 
	}
	if( is.null(preR) ){
		stop("\ndccsim-->error: preR cannot be NULL when method uses only DCC specification!")
	} else{
		chk = dcc.symcheck(preR, m, d = rep(1, m))
	}
	
	preQ = preR
	
	# uncvariance has uGARCHfit AND uGARCHspec dispatch methods
	uncv = sapply(spec@uspec@spec, FUN = function(x) uncvariance(x))
	
	
	if( is.null(preH) ){
		preH = diag(sqrt(uncv)) %*% preR %*% diag(sqrt(uncv))
	} else{
		chk = dcc.symcheck(preH, m, d = NULL)
	}
	if( !is.null(presigma) ){
		if( !is.matrix(presigma) ) stop("\ndccsim-->error: presigma must be a matrix.")
		if( dim(presigma)[2] != m ) stop("\ndccsim-->error: wrong column dimension for presigma.")
	}
	
	if( !is.null(preresiduals) ){
		if( !is.matrix(preresiduals) ) stop("\ndccsim-->error: preresiduals must be a matrix.")
		if( dim(preresiduals)[2] != m ) stop("\ndccsim-->error: wrong column dimension for preresiduals.")
	}
	
	if( !is.null(prereturns) ){
		if( !is.matrix(prereturns) ) stop("\ndccsim-->error: prereturns must be a matrix.")
		if( dim(prereturns)[2] != m ) stop("\ndccsim-->error: wrong column dimension for prereturns.")
	}
	
	garchpars	= sapply(spec@uspec@spec, FUN = function(x) unlist(x@optimization.model$fixed.pars))
	dccpars		= unlist(spec@mspec$optimization.model$fixed.pars)
	dcca = dccpars[1:dccOrder[1]]
	dccb = dccpars[(dccOrder[1]+1):sum(dccOrder[1:2])]
	simRes = simX = simR = simQ = simH = simSeries = vector(mode = "list", length = m.sim)
	
	if( parallel ){
		os = .Platform$OS.type
		if(is.null(parallel.control$pkg)){
			if( os == "windows" ) parallel.control$pkg = "snowfall" else parallel.control$pkg = "multicore"
			if( is.null(parallel.control$cores) ) parallel.control$cores = 2
		} else{
			mtype = match(tolower(parallel.control$pkg[1]), c("multicore", "snowfall"))
			if(is.na(mtype)) stop("\nParallel Package type not recognized in parallel.control\n")
			parallel.control$pkg = tolower(parallel.control$pkg[1])
			if( os == "windows" && parallel.control$pkg == "multicore" ) stop("\nmulticore not supported on windows O/S\n")
			if( is.null(parallel.control$cores) ) parallel.control$cores = 2 else parallel.control$cores = as.integer(parallel.control$cores[1])
		}
		if( parallel.control$pkg == "multicore" ){
			simH = vector(mode = "list", length = m.sim)
			simX = vector(mode = "list", length = m.sim)
			mtmp = mclapply(as.list(1:m.sim), FUN = function(j)
						.dccsimf(Z = z[,,j], Qbar = preQ, Rbar = preR, dcca, dccb, mo, n.sim, n.start, m, rseed[j]), 
					mc.cores = parallel.control$cores)
			simR = lapply(mtmp, FUN = function(x) if(is.matrix(x$R)) array(x$R, dim = c(m, m, n.sim)) else last(x$R, n.sim))
			simQ = lapply(mtmp, FUN = function(x) if(is.matrix(x$Q)) array(x$Q, dim = c(m, m, n.sim)) else last(x$Q, n.sim))
			stdresid = vector(mode = "list", length = m)
			for(i in 1:m) stdresid[[i]] = sapply(mtmp, FUN = function(x) x$stdresid[,i])			
			xtmp = mclapply(as.list(1:m), FUN = function(j){
						maxx = spec@uspec@spec[[j]]@optimization.model$maxOrder2;
						htmp = ugarchpath(spec@uspec@spec[[j]], n.sim = n.sim + n.start, n.start = 0, m.sim = m.sim,
								custom.dist = list(name = "sample", distfit = matrix(stdresid[[j]][-(1:mo), ], ncol = m.sim)),
								presigma = if( is.null(presigma) ) NA else tail(presigma[,j], maxx), 
								preresiduals = if( is.null(preresiduals) ) NA else tail(preresiduals[,j], maxx), 
								prereturns = if( is.null(prereturns) | spec@mspec$varmodel$VAR ) NA else tail(prereturns[,j], maxx),
								mexsimdata = if( !spec@mspec$varmodel$VAR ) mexsimdata[[j]] else NULL, vexsimdata = vexsimdata[[j]] )
						h = matrix(tail(htmp@path$sigmaSim^2, n.sim), nrow = n.sim)
						x = matrix(htmp@path$seriesSim,  nrow = n.sim + n.start)
						return(list(h = h, x = x))}, mc.cores = parallel.control$cores)
			H = array(NA, dim = c(n.sim, m, m.sim))
			tmpH = array(NA, dim = c(m, m, n.sim))
			for(i in 1:n.sim) H[i,,] = t(sapply(xtmp, FUN = function(x) as.numeric(x$h[i,])))
			
			for(i in 1:m.sim){
				for(j in 1:n.sim){
					tmpH[ , , j] = diag(sqrt( H[j, , i]) ) %*% simR[[i]][ , , j] %*% diag(sqrt( H[j, , i] ) )
				}
				simH[[i]] = tmpH
			}
			simxX = array(NA, dim = c(n.sim, m, m.sim))
			for(i in 1:n.sim) simxX[i,,] = t(sapply(xtmp, FUN = function(x) as.numeric(x$x[i,])))
			simX = vector(mode = "list", length = m.sim)
			for(i in 1:m.sim) simX[[i]] = matrix(simxX[,,i], nrow = n.sim)
		} else{
			sfInit(parallel = TRUE, cpus = parallel.control$cores)
			sfExport("spec", "z", "preQ", "preR", "dcca", "dccb", "mo", "n.sim", "n.start", 
					"m", "prereturns", "preresiduals", "presigma", "mexsimdata", "vexsimdata", 
					"rseed", local = TRUE)
			simH = vector(mode = "list", length = m.sim)
			simX = vector(mode = "list", length = m.sim)
			mtmp = sfLapply(as.list(1:m.sim), fun = function(j) .dccsimf(Z = z[,,j], Qbar = preQ, Rbar = preR, dcca, dccb, mo, n.sim, n.start, m, rseed[j]))
			simR = lapply(mtmp, FUN = function(x) if(is.matrix(x$R)) array(x$R, dim = c(m, m, n.sim)) else last(x$R, n.sim))
			simQ = lapply(mtmp, FUN = function(x) if(is.matrix(x$Q)) array(x$Q, dim = c(m, m, n.sim)) else last(x$Q, n.sim))
			stdresid = vector(mode = "list", length = m)
			for(i in 1:m) stdresid[[i]] = sapply(mtmp, FUN = function(x) x$stdresid[,i])			
			xtmp = sfLapply(as.list(1:m), fun = function(j){
						maxx = spec@uspec@spec[[j]]@optimization.model$maxOrder2;
						htmp = rgarch::ugarchpath(spec@uspec@spec[[j]], n.sim = n.sim + n.start, n.start = 0, m.sim = m.sim,
								custom.dist = list(name = "sample", distfit = matrix(stdresid[[j]][-(1:mo), ], ncol = m.sim)),
								presigma = if( is.null(presigma) ) NA else tail(presigma[,j], maxx), 
								preresiduals = if( is.null(preresiduals) ) NA else tail(preresiduals[,j], maxx), 
								prereturns = if( is.null(prereturns) | spec@mspec$varmodel$VAR ) NA else tail(prereturns[,j], maxx),
								mexsimdata = if( !spec@mspec$varmodel$VAR ) mexsimdata[[j]] else NULL, vexsimdata = vexsimdata[[j]] )
						h = matrix(tail(htmp@path$sigmaSim^2, n.sim), nrow = n.sim)
						x = matrix(htmp@path$seriesSim,  nrow = n.sim + n.start)
						return(list(h = h, x = x))})
			if( !spec@mspec$varmodel$VAR ) sfStop()
			H = array(NA, dim = c(n.sim, m, m.sim))
			tmpH = array(NA, dim = c(m, m, n.sim))
			for(i in 1:n.sim) H[i,,] = t(sapply(xtmp, FUN = function(x) as.numeric(x$h[i,])))
			
			for(i in 1:m.sim){
				for(j in 1:n.sim){
					tmpH[ , , j] = diag(sqrt( H[j, , i]) ) %*% simR[[i]][ , , j] %*% diag(sqrt( H[j, , i] ) )
				}
				simH[[i]] = tmpH
			}
			simxX = array(NA, dim = c(n.sim, m, m.sim))
			for(i in 1:n.sim) simxX[i,,] = t(sapply(xtmp, FUN = function(x) as.numeric(x$x[i,])))
			simX = vector(mode = "list", length = m.sim)
			for(i in 1:m.sim) simX[[i]] = matrix(simxX[,,i], nrow = n.sim)
		}
	} else{
		simH = vector(mode = "list", length = m.sim)
		simX = vector(mode = "list", length = m.sim)
		mtmp = lapply(as.list(1:m.sim), FUN = function(j) .dccsimf(Z = z[,,j], Qbar = preQ, Rbar = preR, dcca, dccb, mo, n.sim, n.start, m, rseed[j]))
		simR = lapply(mtmp, FUN = function(x) if(is.matrix(x$R)) array(x$R, dim = c(m, m, n.sim)) else last(x$R, n.sim))
		simQ = lapply(mtmp, FUN = function(x) if(is.matrix(x$Q)) array(x$Q, dim = c(m, m, n.sim)) else last(x$Q, n.sim))
		stdresid = vector(mode = "list", length = m)
		for(i in 1:m) stdresid[[i]] = sapply(mtmp, FUN = function(x) x$stdresid[,i])			
		xtmp = lapply(as.list(1:m), FUN = function(j){
					maxx = spec@uspec@spec[[j]]@optimization.model$maxOrder2;
					htmp = ugarchpath(spec@uspec@spec[[j]], n.sim = n.sim + n.start, n.start = 0, m.sim = m.sim,
							custom.dist = list(name = "sample", distfit = matrix(stdresid[[j]][-(1:mo), ], ncol = m.sim)),
							presigma = if( is.null(presigma) ) NA else tail(presigma[,j], maxx), 
							preresiduals = if( is.null(preresiduals) ) NA else tail(preresiduals[,j], maxx), 
							prereturns = if( is.null(prereturns) | spec@mspec$varmodel$VAR ) NA else tail(prereturns[,j], maxx),
							mexsimdata = if( !spec@mspec$varmodel$VAR ) mexsimdata[[j]] else NULL, vexsimdata = vexsimdata[[j]] )
					h = matrix(tail(htmp@path$sigmaSim^2, n.sim), nrow = n.sim)
					x = matrix(htmp@path$seriesSim,  nrow = n.sim + n.start)
					return(list(h = h, x = x))})
		H = array(NA, dim = c(n.sim, m, m.sim))
		tmpH = array(NA, dim = c(m, m, n.sim))
		for(i in 1:n.sim) H[i,,] = t(sapply(xtmp, FUN = function(x) as.numeric(x$h[i,])))
		
		for(i in 1:m.sim){
			for(j in 1:n.sim){
				tmpH[ , , j] = diag(sqrt( H[j, , i]) ) %*% simR[[i]][ , , j] %*% diag(sqrt( H[j, , i] ) )
			}
			simH[[i]] = tmpH
		}
		simxX = array(NA, dim = c(n.sim, m, m.sim))
		for(i in 1:n.sim) simxX[i,,] = t(sapply(xtmp, FUN = function(x) as.numeric(x$x[i,])))
		simX = vector(mode = "list", length = m.sim)
		for(i in 1:m.sim) simX[[i]] = matrix(simxX[,,i], nrow = n.sim)
	}
	if( spec@mspec$varmodel$VAR ){
		# mexsimdata check
		mxn = VAR.fit$mxn
		Data = VAR.fit$xfitted
		if(mxn > 0 ){
			if( !is.null(mexsimdata) ){
				if( !is.list(mexsimdata) | length(mexsimdata) != m.sim ) stop("\ndccsim-->error: mexsimdata must be a list of length equal to m.sim...", call. = FALSE)
				for(j in 1:m.sim){
					if( !is.matrix(mexsimdata[[j]]) ) stop("\ndccsim-->error: mexsimdata list submatrix must be a matrix...", call. = FALSE)
					if( dim(mexsimdata[[j]])[2] != mxn ) stop("\ndccsim-->error: wrong number of external regressors in list submatrix!...", call. = FALSE)
					if( dim(mexsimdata[[j]])[1] < (n.sim + n.start) ) stop("\ndccsim-->error: external regressor list submatrix has less points than requested simulation length (n.sim + n.start)...", call. = FALSE)
				}
			} else{
				mexsimdata = vector(mode = "list", length = m.sim)
				for(j in 1:m.sim) mexsimdata[[i]] = matrix(0, ncol = mxn, nrow = (n.sim + n.start))
			}
		} else{
			mexsimdata = NULL
		}
		
		p = VAR.fit$lag
		
		if( !is.null(prereturns) ){
			if(!is.matrix(prereturns)) stop("\nprereturns must be a matrix\n")
			if(dim(prereturns)[1] != p) stop("\nwrong number of rows for prereturns matrix")
			if(dim(prereturns)[2] != m) stop("\nwrong number of cols for prereturns matrix")
			prereturns = matrix(prereturns, nrow = p, ncol = m)
		} else{
			prereturns = matrix(0, nrow = p, ncol = m, byrow = TRUE)
		}
		xlist = vector(mode = "list", length = m.sim)
		if( n.sim == 1 && n.start == 0 ){
			# SPECIAL ROUTINE TO SPEED UP 1-AHEAD ESTIMATION
			tmp = varxsimXX(X = Data, VAR.fit$Bcoef, p, m.sim = m.sim, prereturns, resids = t(sapply(simX, FUN = function(x) x)), if(!is.null(mexsimdata)) t(sapply(mexsimdata, FUN = function(x) x)) else NULL)
			for(i in 1:m.sim) xlist[[i]] = tmp[i, , drop = FALSE]
			simRes = simX
			simX = xlist
		} else{
			if( parallel ){
				if( parallel.control$pkg == "multicore" ){
					xlist = mclapply(as.list(1:m.sim), FUN = function(i) varxsim(X = Data, VAR.fit$Bcoef, p, n.sim, 
										n.start, prereturns, resids = simX[[i]], mexsimdata[[i]]), 
							mc.cores = parallel.control$cores)
				} else{
					# when the GARCH model is without mean equation, simX == simulated residuals
					sfExport("Data", "VAR.fit", "p", "n.sim", "n.start", "prereturns", "simX", "mexsimdata", local = TRUE)
					xlist = sfLapply(as.list(1:m.sim), fun = function(i) rgarch:::varxsim(X = Data, VAR.fit$Bcoef, p, 
										n.sim, n.start, prereturns, resids = simX[[i]], mexsimdata[[i]]))
					sfStop()
				}
			} else{
				xlist = lapply(as.list(1:m.sim), FUN = function(i) varxsim(X = Data, VAR.fit$Bcoef, p, 
									n.sim, n.start, prereturns, resids = simX[[i]], mexsimdata[[i]]))
			}
			simRes = simX
			simX = xlist
		}
	} else{
		# reshape
		for(j in 1:m.sim) simX[[j]] = tail(simX[[j]], n.sim)
	}
	
	msim = list()
	msim$simH = simH
	msim$simR = simR
	msim$simQ = simQ
	msim$simX = simX
	msim$simRes = simRes
	msim$Z = z
	msim$rseed = rseed
	msim$model$Data = NULL
	msim$model$garchpars = garchpars
	msim$model$dccpars = dccpars
	msim$model$n.sim = n.sim
	msim$model$m.sim = m.sim
	msim$model$n.start = n.start
	msim$model$startMethod = "unconditional"
	msim$model$dccOrder = dccOrder
	msim$model$distribution = "mvlaplace"
	ans = new("DCCsim",
			msim = msim)
	
	return( ans )
	
}

###################################################################
# Functions for parallel evaluation based on large m.sim
.mcsim = function(fit, z, preQ, preR, dcca, dccb, mo, n.sim, n.start, m, startMethod, prereturns, preresiduals, 
		presigma, mexsimdata, vexsimdata, rseed)
{
	tmp = .dccsimf(Z = z, Qbar = preQ, Rbar = preR, dcca, dccb, mo, n.sim, n.start, m, rseed)
	simR = if(is.matrix(tmp$R)) array(tmp$R, dim = c(m, m, n.sim)) else last(tmp$R, n.sim)
	simQ = if(is.matrix(tmp$Q)) array(tmp$Q, dim = c(m, m, n.sim)) else last(tmp$Q, n.sim)
	stdresid = tmp$stdresid
	X = H = matrix(NA, nrow = n.sim + n.start, ncol = m)
	simH = array(NA, dim = c(m, m, n.sim))
	for(i in 1:m){
		maxx = fit@ufit@fit[[i]]@model$maxOrder2
		htmp = ugarchsim(fit@ufit@fit[[i]], n.sim = n.sim + n.start, n.start = 0, m.sim = 1,
				startMethod = startMethod[1], custom.dist = list(name = "sample", 
						distfit = matrix(stdresid[-(1:mo), i], ncol = 1)),
				presigma = if( is.null(presigma) ) NA else tail(presigma[,i], maxx), 
				preresiduals = if( is.null(preresiduals) ) NA else tail(preresiduals[,i], maxx), 
				prereturns = if( is.null(prereturns) | !is.null(fit@mfit$model$vrmodel) ) NA else tail(prereturns[,i], maxx),
				mexsimdata = if( is.null(fit@mfit$model$vrmodel) ) mexsimdata[[i]] else NULL, vexsimdata = vexsimdata[[i]] )
		H[, i] = htmp@simulation$sigmaSim^2
		X[, i] = htmp@simulation$seriesSim
	}
	H = tail(H, n.sim)
	for(i in 1:n.sim){
		simH[ , , i] = diag(sqrt( H[i, ]) ) %*% simR[ , , i] %*% diag(sqrt( H[i, ]) )
	}
	simX = X
	return( list(simH = simH, simX = simX, simR = simR, simQ = simQ) )
}

.mcpath = function(spec, z, preQ, preR, dcca, dccb, mo, n.sim, n.start, m, prereturns, preresiduals, 
		presigma, mexsimdata, vexsimdata, rseed)
{
	tmp = .dccsimf(Z = z, Qbar = preQ, Rbar = preR, dcca, dccb, mo, n.sim, n.start, m, rseed)
	simR = if(is.matrix(tmp$R)) array(tmp$R, dim = c(m, m, n.sim)) else last(tmp$R, n.sim)
	simQ = if(is.matrix(tmp$Q)) array(tmp$Q, dim = c(m, m, n.sim)) else last(tmp$Q, n.sim)
	stdresid = tmp$stdresid
	X = H = matrix(NA, nrow = n.sim + n.start, ncol = m)
	simH = array(NA, dim = c(m, m, n.sim))
	for(i in 1:m){
		maxx = spec@uspec@spec[[i]]@optimization.model$maxOrder2
		htmp = ugarchpath(spec@uspec@spec[[i]], n.sim = n.sim + n.start, n.start = 0, m.sim = 1,
				custom.dist = list(name = "sample", distfit = matrix(stdresid[-(1:mo), i], ncol = 1)),
				presigma = if( is.null(presigma) ) NA else tail(presigma[,i], maxx), 
				preresiduals = if( is.null(preresiduals) ) NA else tail(preresiduals[,i], maxx), 
				prereturns = if( is.null(prereturns) | spec@mspec$varmodel$VAR ) NA else tail(prereturns[,i], maxx),
				mexsimdata = if( !spec@mspec$varmodel$VAR ) mexsimdata[[i]] else NULL, vexsimdata = vexsimdata[[i]] )
		H[, i] = htmp@path$sigmaSim^2
		X[, i] = htmp@path$seriesSim
	}
	H = tail(H, n.sim)
	for(i in 1:n.sim){
		simH[ , , i] = diag(sqrt( H[i, ]) ) %*% simR[ , , i] %*% diag(sqrt( H[i, ]) )
	}
	simX = X
	return( list(simH = simH, simX = simX, simR = simR, simQ = simQ) )
}

.rolldcc = function(spec, data, n.ahead = 1, forecast.length = 50, refit.every = 25, 
		refit.window = c("recursive", "moving"), solver = "solnp", 
		fit.control = list(eval.se = TRUE), 
		parallel = FALSE, parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2), 
		solver.control = list(), save.fit = FALSE, save.wdir = NULL, trace = FALSE, ...)
{
	forecast.length = floor(forecast.length/refit.every) * refit.every	
	WD <- getwd()
	if(is.null(save.wdir)){
		if (!is.null(WD)) setwd(WD)
	} else{
		ND = save.wdir
		if (!is.null(ND)) setwd(ND)
	}
	# Setup Indexed Data for Fitting
	cnames = colnames(data)
	if( is.null(cnames) ) cnames = paste("Asset_", 1:dim(data)[2], sep = "")
	xdata = .extractmdata(data)
	data = as.matrix(data)
	dates = xdata$pos
	if(is.null(fit.control$eval.se)) fit.control$eval.se = FALSE
	T = dim(data)[1]
	startT = T - forecast.length
	# forecast window index
	fwindex = t(.embed((T-forecast.length - n.ahead + 2):T, refit.every, by = refit.every, TRUE ) )
	fitindx.start = c(fwindex[1,] - 1)
	fwindex.end = fwindex[refit.every,]
	nf = length(fwindex.end)
	fitdata = vector(mode = "list", length = nf)
	mcfit = vector(mode = "list", length = nf)
	
	#gofit = vector(mode = "list", length = nf)
	outdata = vector(mode = "list", length = nf)
	distribution = spec@mspec$distribution
	forecastlist = cf = lik = vector(mode = "list", length = nf)
	n.data = rep(0, nf)
	mi = 0
	for(i in 1:nf){
		if(refit.window[1]=="recursive"){
			fitdata = data[1:(fitindx.start[i]+refit.every),]
		} else{
			fitdata = data[(i*refit.every):(fitindx.start[i]+refit.every),]
		}
		mcfit = dccfit(spec, data = fitdata, out.sample = refit.every, solver = solver,
				fit.control = fit.control, solver.control = solver.control, parallel = parallel, 
				parallel.control = parallel.control)
		
		cf[[i]] = coef(mcfit)
		lik[[i]] = likelihood(mcfit)
		n.data[i] = dim(fitdata)[1]
		# solver.control = list(tol=1e-6, delta=1e-5)
		# mcfit = dccfit(spec, data = fitdata, out.sample = refit.every, solver = solver,
		# fit.control = fit.control, solver.control = solver.control)	
		
		forecastlist[[i]] = dccforecast(mcfit, n.ahead = n.ahead, n.roll = refit.every - 1, 
				parallel = parallel, parallel.control = parallel.control)
		if(save.fit){
			eval(parse(text = paste("mcfitroll",i,"=mcfit",sep = "")))
			eval(parse(text = paste("save(mcfitroll",i,",file='mcfitroll",i,".rda')",sep = "")))
		}
		if(trace) print(i)
		mi = mi + 1
	}
	mspec = list()
	mspec$refit.every = refit.every
	mspec$forecast.length = forecast.length
	mspec$asset.names = cnames
	mspec$dccspec = spec
	mspec$refit.window = refit.window
	mspec$coef = cf
	mspec$likelihood = lik
	mspec$n.data = n.data
	ans = new("DCCroll",
			forecast = forecastlist,
			spec = mspec)
	# return to the current directory
	setwd(WD)
	return(ans)
}
.dccparidx = function( speclist )
{
	m = length(speclist@spec)
	type = speclist@type
	idx.s = idx.e = vector(mode = "numeric", length = m)
	idxi.s = idxi.e = vector(mode = "numeric", length = m)
	idxa.s = idxa.e = vector(mode = "numeric", length = m)
	pars = lapply( speclist@spec, FUN = function(x) x@optimization.model$modelnames )
	for(i in 1:m){
		n = length( pars[[i]] )
		nx = 0
		if( any(substr(pars[[i]], 1, 2) == "mu") ){
			nx = c( nx, which( substr(pars[[i]], 1, 2) == "mu" ) )
		}
		
		if( any(substr(pars[[i]], 1, 2) == "ar") ){
			nx = c( nx, which( substr(pars[[i]], 1, 2) == "ar" ) ) 
		}
		
		if( any(substr(pars[[i]], 1, 2) == "ma") ){
			nx = c( nx, which( substr(pars[[i]], 1, 2) == "ma" ) ) 
		}
		
		if( any(substr(pars[[i]], 1, 6) == "inmean") ){
			nx = c( nx, which( substr(pars[[i]], 1, 6) == "inmean" ) ) 
		}
		
		if( any(substr(pars[[i]], 1, 7) == "darfima") ){
			nx = c( nx, which( substr(pars[[i]], 1, 7) == "darfima" ) ) 
		}
		
		if( any(substr(pars[[i]], 1, 5) == "mxreg") ){
			nx = c( nx, which( substr(pars[[i]], 1, 5) == "mxreg" ) ) 
		}
		
		if( i == 1){
			idx.s[i] = max(nx) + 1
			idx.e[i] = length(pars[[i]])
			
			idxa.s[i] = 1
			idxa.e[i] = length(pars[[i]])
		} else{
			idx.s[i] = idx.e[i - 1] + max(nx) + 1
			idx.e[i] = idx.e[i - 1] + n
			
			idxa.s[i] = idxa.e[i - 1] + 1
			idxa.e[i] = idxa.e[i - 1] + n
		}
		idxi.s[i] = max(nx) + 1
		idxi.e[i] =  n
		
		
	}
	idx = cbind(idx.s, idx.e)
	idxi = cbind(idxi.s, idxi.e)
	idxa = cbind(idxa.s, idxa.e)
	
	return( list( idx = idx, idxi = idxi, idxa = idxa ) )
}

.dccfilterspec1 = function(spec, dccpars, fitlist)
{
	uspec = spec@uspec
	m = length(spec@uspec@spec)
	newspec = vector(mode = "list", length = m)
	for(i in 1:m){
		newspec[[i]] = uspec@spec[[i]]
		newspec[[i]]@optimization.model$fixed.pars = as.list(coef(fitlist@fit[[i]]))
	}
	uspec = multispec( newspec )
	dspec = dccspec(uspec, dccOrder = spec@mspec$dccOrder, 
			distribution = spec@mspec$distribution, fixed.pars = as.list(dccpars))
	return( dspec )
}

.dccfilterspec2 = function(spec, dccpars, garchpars)
{	
	uspec = spec@uspec
	idx = .dccparidx( uspec )$idxa
	m = length(spec)
	newspec = vector(mode = "list", length = m)
	for(i in 1:m){
		newspec[[i]] = uspec[[i]]
		gnames = uspec[[i]]@optimization.model$modelnames
		gpars = garchpars[idx[i,1]:idx[i,2]]
		names(gpars) = gnames
		newspec[[i]]@optimization.model$fixed.pars = as.list(gpars)
	}
	uspec = multispec( newspec )
	dspec = dccspec(uspec, dccOrder = spec@mspec$dccOrder, 
			distribution = spec@mspec$distribution, fixed.pars = as.list(dccpars))
	return( dspec )
}

# n-ahead forecast based on the approximation method described in:
# ENGLE, R. and SHEPPARD (2001), K. Theoretical and empirical properties of
# dynamic conditional correlation multivariate garch. NBER Working Papers, No. 8554

.dccforecastm = function(fit, n.ahead = 1, n.roll = 10, external.forecasts = list(mregfor = NULL, vregfor = NULL), 
		parallel = FALSE, parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2), ...)
{
	filterlist = multifilter(multifitORspec = fit@ufit, data = fit@mfit$origdata, out.sample = 0, n.old = NULL,
			parallel = parallel, parallel.control = parallel.control)
	n.roll = n.roll + 1
	m = length(fit@ufit@fit)
	out.sample = fit@mfit$out.sample
	pos = fit@mfit$model$pos.matrix
	# we only allow n.ahead = 1 when using n.roll (for simlicity and avoidance of arrays	
	fitlist = fit@ufit
	mo = max(fit@mfit$model$dccOrder)
	# we augmented n.roll before, now subtract
	# forclist = multiforecast(multifitORspec = fitlist, n.ahead = n.ahead, n.roll = n.roll - 1)
	forclist = multiforecast(multifitORspec = fitlist, n.ahead = n.ahead, n.roll = n.roll - 1, 
			external.forecasts = external.forecasts, parallel = parallel, 
			parallel.control = parallel.control, ...)
	sig = sigma(filterlist)
	resid = residuals(filterlist)
	stdresid = resid/sig
	
	T = dim(fit@mfit$H)[3]
	pars = coef(fit, type = "dcc")
	dcca = dccb = skew = shape = 0
	if(pos[1,3] == 1) dcca = pars[pos[1,1]:pos[1,2]]
	if(pos[2,3] == 1) dccb = pars[pos[2,1]:pos[2,2]]
	
	## call .nsdccforecast with rolling window
	Rbar = Rtfor = Htfor = Qtfor = vector(mode = "list", length = n.roll)
	Qstart = last( rcor(fit, type = "Q"), mo )
	Rstart = last( rcor(fit, type = "R"), mo )
	Hstart = last( rcov(fit), mo)
	
	for(i in 1:n.roll){
		xQbar = cov(stdresid[1:(T + i - 1), ])
		xstdresids = stdresid[(T - mo + i ):(T  +  i - 1), , drop = FALSE]
		xfsig = matrix( sapply( forclist@forecast, FUN = function(x) unlist(x@forecast$forecasts[[i]]['sigma'])  ), ncol = m)
		ans = .dccforecastn(Qstart, Rstart, Hstart, xQbar, xstdresids, xfsig, dcca, dccb, n.ahead, mo)
		Rtfor[[i]] = ans$Rtfor
		Qtfor[[i]] = ans$Qtfor
		Htfor[[i]] = ans$Htfor
		Rbar[[i]] = ans$Rbar
		Qstart = last( .abind(Qstart, ans$Qtfor[, , 1]), mo )
		Rstart = last( .abind(Rstart, ans$Rtfor[, , 1]), mo )
		Hstart = last( .abind(Hstart, ans$Htfor[, , 1]), mo )
	}
	
	forc = list( H = Htfor, R = Rtfor, Q = Qtfor, Rbar = Rbar, ufor = forclist )
	
	return(forc)
}
.dccforecastn = function(Qstart, Rstart, Hstart, Qbar, stdresids, fsig, dcca, dccb, n.ahead, mo)
{
	m = dim(Qbar)[1]
	Qtfor = Rtfor = Htfor = array(NA, dim = c(m, m, n.ahead + mo))
	Qtfor[ , , 1:mo]  =  Qstart[, , 1:mo]
	Rtfor[ , , 1:mo]  =  Rstart[, , 1:mo]
	Htfor[ , , 1:mo] = Hstart[, , 1:mo]
	dccsum = sum(dcca) + sum(dccb)
	Qt_1 = (1 - dccsum) * Qbar
	for(i in 1:n.ahead){
		Qtfor[, , mo + i] = Qt_1
		if( i == 1 ){
			for(j in 1:length(dcca)){
				Qtfor[ , , mo + 1] = Qtfor[ , , mo + 1] + dcca[j] * (stdresids[(mo + 1 - j), ] %*% t(stdresids[(mo + 1 - j), ]))
			}
			for(j in 1:length(dccb)){
				Qtfor[ , , mo + 1] = Qtfor[ , , mo + 1] + dccb[j] * Qtfor[ , , mo + 1 - j]
			}
			Qtmp = diag( 1/sqrt( diag(Qtfor[ , , mo + 1]) ) , m, m)
			Rtfor[ , , mo + 1] =  Qtmp %*% Qtfor[ , , mo + 1] %*% t(Qtmp)
			Dtmp = diag(fsig[1, ], m, m)
			Htfor[ , , mo + 1] = Dtmp %*% Rtfor[ , , mo + 1] %*% Dtmp
			# now the unconditional calculations
			Qt_1star = diag( 1/sqrt( diag(Qtfor[, , mo + 1]) ) , m, m)
			ER_1 = Qt_1star %*% Qtfor[, , mo + 1] %*% t(Qt_1star)
			Qbarstar = diag( 1/sqrt( diag(Qbar) ) , m, m)
			Rbar = Qbarstar %*% Qbar %*% t(Qbarstar)
		} else{
			Rtfor[, , mo + i] = (1 - dccsum^(i - 1) ) * Rbar + dccsum^(i - 1) * ER_1
			Dtmp = diag(fsig[i, ], m, m)
			Htfor[, , mo + i] = Dtmp %*% Rtfor[, , mo + i] %*% Dtmp
			Qtfor[, , mo + i] = Qtfor[, , mo + 1]
		}
	}
	return( list( Rtfor = Rtfor[, , -(1:mo), drop = FALSE], Htfor = Htfor[, , -(1:mo), drop = FALSE], 
					Qtfor = Qtfor[, , -(1:mo), drop = FALSE], Rbar = Rbar, ER_1 = ER_1) )
		
}

.dccsimf = function(Z, Qbar, Rbar, dcca, dccb, mo, n.sim, n.start, m, rseed)
{
	n = n.sim + n.start + mo
	
	set.seed(rseed[1]+1)
	stdresid = matrix(rnorm(m * (n.sim+n.start+mo)), nrow = n.sim + n.start + mo, ncol = m)
	sumdcca = sum(dcca)
	sumdccb = sum(dccb)
	dccsumx = sumdcca + sumdccb
	ldcca = length(dcca)
	ldccb = length(dccb)
	# ToDo: Assumes that no order is zero!
	dccorderx = c(ldcca, ldccb, mo, n)
	
	res = .Call(name = "dccsim", Qbar = as.matrix(Qbar), Rbar = as.matrix(Rbar), 
			Z = as.matrix(Z), Res = as.matrix(stdresid),  dcca = as.numeric(dcca), 
			dccb = as.numeric(dccb), dccorder = as.numeric(dccorderx), dccsum = as.numeric(dccsumx))

	Q = array(NA, dim = c(m, m, n.sim + n.start + mo))
	R = array(NA, dim = c(m, m, n.sim + n.start + mo))
	for(i in 1:(n.sim + n.start + mo)){
		R[,,i] = res[[2]][[i]]
		Q[,,i] = res[[3]][[i]]
	}
	ans = list( Q = Q, R = R, stdresid = res[[3]])
	return( ans )
}
	

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.