R/gogarch-fit.R

In rgarch: Flexible GARCH modelling in R

#################################################################################
##
##   R package rgarch by Alexios Ghalanos Copyright (C) 2008, 2009, 2010, 2011
##   This file is part of the R package rgarch.
##
##   The R package rgarch is free software: you can redistribute it and/or modify
##   it under the terms of the GNU General Public License as published by
##   the Free Software Foundation, either version 3 of the License, or
##   (at your option) any later version.
##
##   The R package rgarch is distributed in the hope that it will be useful,
##   but WITHOUT ANY WARRANTY; without even the implied warranty of
##   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
##   GNU General Public License for more details.
##
#################################################################################


#-----------------------------------------------------------------------------------------------
# distribution based methods

#-----------------------------------------------------------------------------------------------
# manig
.gogarchfit.nig = function(spec, data, out.sample = 0, solver = "solnp",
		fit.control = list(stationarity = 1), solver.control = list(), parallel = FALSE, 
		parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2), VAR.fit = NULL,...)
{
	tic = Sys.time()
	
	if( is.null( colnames(data) ) ) cnames = paste("Asset_", 1:k, sep = "") else cnames = colnames(data)
	xdata = .extractmdata(data)
	if(!is.numeric(out.sample)) stop("\ngogarchfit-->error: out.sample must be numeric\n")
	if(as.numeric(out.sample) < 0) stop("\ngogarchfit-->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("\ngogarchfit-->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
	
	armodel = vrmodel = NULL
	if( spec@mspec$mmodel$demean == "VAR"){
		if( spec@mspec$mmodel$mxn>0 ){
			ex = spec@mspec$mmodel$external.regressors
			if( dim(ex)[1] < dim(data)[1]  ) stop("\ngogarchfit-->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(data)[1]
			K = dim(data)[2]
			mu = VAR.fit$xfitted
			vrmodel = VAR.fit
			vrmodel$mxn = mxn
		} else {
			cat("\nestimating VAR model...")
			ic = spec@mspec$mmodel$lag.criterion[1]
			if( !is.null(spec@mspec$mmodel$lag.max) ){
				mp = .varxselect(y = data, lag.max = spec@mspec$mmodel$lag.max, exogen = ex)$selection
				mp = as.numeric( mp[which(substr(names(mp), 1, 2) == substr(ic, 1, 2))] )
			} else {
				mp = spec@mspec$mmodel$lag
			}
			p = mp
			vrmodel = varxfilter(X = data, p = mp, exogen = ex, Bcoef = NULL, robust = spec@mspec$mmodel$robust, 
					gamma = spec@mspec$mmodel$robust.control$gamma, delta = spec@mspec$mmodel$robust.control$delta, 
					nc = spec@mspec$mmodel$robust.control$nc, ns = spec@mspec$mmodel$robust.control$ns)	
			cat("done!\n")
			zdata = vrmodel$xresiduals
			# should be N - p
			N = dim(data)[1]
			K = dim(data)[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{
		# mu-ARX or just muX
		if( spec@mspec$mmodel$mxn>0 ){
			ex = spec@mspec$mmodel$external.regressors
			if( dim(ex)[1] < dim(data)[1]  ) stop("\ngogarchfit-->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
		}
		p = mp = spec@mspec$mmodel$lag
		arspec = arfimaspec(mean.model = list(armaOrder = c(p, 0), include.mean = TRUE, 
						arfima = FALSE, external.regressors = ex), distribution.model = "norm")
		cat("\nestimating univariate mean model...")
		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" ){
				if(!exists("mclapply")){
					require('multicore')
				}
				arfit = mclapply(1:k, FUN = function(i) arfimafit(spec = arspec, 
									data = data[,i], out.sample = 0, solver = solver, 
									solver.control = solver.control, fit.control = fit.control), 
						mc.cores = parallel.control$cores)
			} else{
				if(!exists("sfLapply")){
					require('snowfall')
				}
				sfInit(parallel=TRUE, cpus = parallel.control$cores)
				sfExport("arspec", "data", "solver", "solver.control", "fit.control",local = TRUE)
				arfit = sfLapply(1:k, fun = function(i) rgarch::arfimafit(spec = arspec, 
									data = data[,i], out.sample = 0, solver = solver, 
									solver.control = solver.control, fit.control = fit.control))
				sfStop()
			}
		} else{
			arfit = lapply(1:k, FUN = function(i) arfimafit(spec = arspec, 
								data = data[,i], out.sample = 0, solver = solver, 
								solver.control = solver.control, fit.control = fit.control))
		}
		cat("done!\n")
		zdata = sapply(arfit, FUN = function(x) residuals(x))
		if( p > 0 ) zdata = zdata[-c(1:p), , drop = FALSE]
		# should be N - p
		N = dim(data)[1]
		K = dim(data)[2]
		mu = sapply(arfit, FUN = function(x) fitted(x))
		if( p > 0 ) mu = mu[-c(1:p), , drop = FALSE]
		armodel$lag = mp
		armodel$mxn = mxn
		armodel$arCoef = arCoef = sapply(arfit, FUN = function(x) coef(x))
		# last coefficient is the sigma which we ignore

		if(out.sample>0){
			if( parallel ){
				if( parallel.control$pkg == "multicore" ){
					if(!exists("mclapply")){
						require('multicore')
					}
					arfit2 = mclapply(1:k, FUN = function(i){
								specx = arfimaspec(mean.model = list(armaOrder = c(p, 0), include.mean = TRUE, 
												arfima = FALSE, external.regressors = orex), distribution.model = "norm",
										fixed.pars = as.list(armodel$arCoef[,i]))
								ans = arfimafilter(spec = specx, data = origdata[,i], out.sample = 0)
								return(ans)}, 
							mc.cores = parallel.control$cores)
				} else{
					if(!exists("sfLapply")){
						require('snowfall')
					}
					sfInit(parallel=TRUE, cpus = parallel.control$cores)
					sfExport("origdata", "p", "orex", "arCoef", local = TRUE)
					arfit2 = sfLapply(1:k, fun = function(i) {
								specx = rgarch:::arfimaspec(mean.model = list(armaOrder = c(p, 0), include.mean = TRUE, 
												arfima = FALSE, external.regressors = orex), distribution.model = "norm",
										fixed.pars = as.list(arCoef[,i]))
								ans = rgarch::arfimafilter(spec = specx, data = origdata[,i], out.sample = 0)
								return(ans)})
					sfStop()
				}
			} else{
				arfit2 = lapply(1:k, FUN = function(i){
							specx = arfimaspec(mean.model = list(armaOrder = c(p, 0), include.mean = TRUE, 
											arfima = FALSE, external.regressors = orex), distribution.model = "norm",
									fixed.pars = as.list(armodel$arCoef[,i]))
							ans = arfimafilter(spec = specx, data = origdata[,i], out.sample = 0)
							return(ans)})
			}
			zorigdata = sapply(arfit2, FUN = function(x) residuals(x))
			if( p > 0 ) zorigdata = zorigdata[-c(1:p), , drop = FALSE]
		} else{
			zorigdata = zdata
		}
	}
	
	iidret = .makeiid(zdata, method = spec@mspec$ica, ica.fix = spec@mspec$ica.fix, ...)
	A = iidret$A
	Y = iidret$Y
	W = iidret$W

	speclist = multispec( replicate(n = K, spec@uspec) )
	
	Y = zorigdata %*% W
	
	fitlist = multifit(multispec = speclist, data = Y, out.sample = out.sample, solver = solver, 
			solver.control = solver.control, fit.control = fit.control, parallel = parallel, parallel.control = parallel.control)
	if(out.sample>0){
		filterlist = multifilter(multifitORspec = fitlist, data = Y, out.sample = 0, parallel = parallel, parallel.control = parallel.control)
	} else{
		filterlist = list()
	}

	# check all converged:
	# we return the fitted list to the global environment in case there is some non convergence
	# for debugging...
	#.fitlist <<- fitlist
	converge = sapply(fitlist@fit, FUN = function(x) x@fit$convergence)
	if(any(converge == 1)){
		# it's useless returning any more info since we have an ICA factor matrix for which the
		# user can do very little other than probably run it again hoping that the new ICA matrix
		# will return a 'nicer' set (perhaps return the random number seed used to initialize ICA).
		stop("\ngogarchdfit-->warning: convergence problem in univariate fit...", call. = FALSE)
	}
	tmp = garchica1.nig(fitlist, A, K, N - p, mu)
	
	tmp$out.sample = out.sample
	tmp$origdata = origdata
	tmp$external.regressors = orex
	tmp$mxn = spec@mspec$mmodel$mxn
	tmp$mdist = "manig"
	tmp$icamethod = spec@mspec$ica
	tmp$meanmodel = spec@mspec$mmodel$model
	tmp$Y = Y
	tmp$W = W
	tmp$A = A
	tmp$K = iidret$K
	tmp$VARmodel = vrmodel
	tmp$ARmodel = armodel
	tmp$residuals = zorigdata
	tmp$filterlist = filterlist
	tmp$timer = Sys.time() - tic
	
	ans = new("goGARCHfit",
			mfit = tmp,
			ufit = fitlist)
	return(ans)
}

.gogarchfilter.nig = function(fit, data, out.sample, n.old, parallel = FALSE, 
		parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2), ...)
{
	
	if( is.null( colnames(data) ) ) cnames = paste("Asset_", 1:k, sep = "") else cnames = colnames(data)
	xdata = .extractmdata(data)
	if(!is.numeric(out.sample)) stop("\ngogarchfilter-->error: out.sample must be numeric\n")
	if(as.numeric(out.sample) < 0) stop("\ngogarchfilter-->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("\ngogarchfilter-->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
	
	
	armodel = vrmodel = NULL
	if( fit@mfit$meanmodel == "VAR" ){
		if( fit@mfit$mxn>0 ){
			ex = fit@mfit$external.regressors
			if( dim(ex)[1] < dim(data)[1]  ) stop("\ngogarchfilter-->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
		}
		cat("\nfiltering VAR model...")
		
		vrmodel = varxfilter(X = data, p = fit@mfit$VARmodel$lag, exogen = ex, Bcoef = fit@mfit$VARmodel$Bcoef)
		cat("done!\n")
		zdata = vrmodel$xresiduals
		# should be N - p
		N = dim(data)[1]
		K = dim(data)[2]
		mu = vrmodel$xfitted
		p = fit@mfit$VARmodel$lag
		# if we have out of sample data we need to filter
		if(out.sample > 0){
			Bcoef = vrmodel$Bcoef
			vrmodel2 = varxfilter(X = origdata, p = fit@mfit$VARmodel$lag, exogen = orex, Bcoef = Bcoef)
			zorigdata = vrmodel2$xresiduals			
		} else{
			zorigdata = zdata
		}
	} else{
		arCoef = fit@mfit$ARmodel$arCoef
		if( fit@mfit$mxn>0 ){
			ex = fit@mfit$external.regressors
			if( dim(ex)[1] < dim(data)[1]  ) stop("\ngogarchfilter-->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
		}
		p = fit@mfit$VARmodel$lag
		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" ){
				if(!exists("mclapply")){
					require('multicore')
				}
				armodel = mclapply(1:k, FUN = function(i){
							specx = arfimaspec(mean.model = list(armaOrder = c(p, 0), include.mean = TRUE, 
											arfima = FALSE, external.regressors = ex), distribution.model = "norm",
									fixed.pars = as.list(arCoef[,i]))
							ans = arfimafilter(spec = specx, data = data[,i], out.sample = 0)
							return(ans)}, 
						mc.cores = parallel.control$cores)
			} else{
				if(!exists("sfLapply")){
					require('snowfall')
				}
				sfInit(parallel=TRUE, cpus = parallel.control$cores)
				sfExport("data", "p", "ex", "arCoef", local = TRUE)
				armodel = sfLapply(1:k, fun = function(i) {
							specx = rgarch:::arfimaspec(mean.model = list(armaOrder = c(p, 0), include.mean = TRUE, 
											arfima = FALSE, external.regressors = ex), distribution.model = "norm",
									fixed.pars = as.list(arCoef[,i]))
							ans = rgarch::arfimafilter(spec = specx, data = data[,i], out.sample = 0)
							return(ans)})
				sfStop()
			}
		} else{
			armodel = lapply(1:k, FUN = function(i){
						specx = arfimaspec(mean.model = list(armaOrder = c(p, 0), include.mean = TRUE, 
										arfima = FALSE, external.regressors = ex), distribution.model = "norm",
								fixed.pars = as.list(arCoef[,i]))
						ans = arfimafit(spec = specx, data = data[,i], out.sample = 0)
						return(ans)})
		}
		cat("done!\n")
		zdata = sapply(armodel, FUN = function(x) residuals(x))
		if( p > 0 ) zdata = zdata[-c(1:p), , drop = FALSE]
		# should be N - p
		N = dim(data)[1]
		K = dim(data)[2]
		mu = sapply(armodel, FUN = function(x) fitted(x))
		if( p > 0 ) mu = mu[-c(1:p), , drop = FALSE]		
		# if we have out of sample data we need to filter
		if(out.sample > 0){
			if( parallel ){
				if( parallel.control$pkg == "multicore" ){
					if(!exists("mclapply")){
						require('multicore')
					}
					armodel2 = mclapply(1:k, FUN = function(i){
								specx = arfimaspec(mean.model = list(armaOrder = c(p, 0), include.mean = TRUE, 
												arfima = FALSE, external.regressors = orex), distribution.model = "norm",
										fixed.pars = as.list(arCoef[,i]))
								ans = arfimafilter(spec = specx, data = origdata[,i], out.sample = 0)
								return(ans)}, 
							mc.cores = parallel.control$cores)
				} else{
					if(!exists("sfLapply")){
						require('snowfall')
					}
					sfInit(parallel=TRUE, cpus = parallel.control$cores)
					sfExport("origdata", "p", "orex", "arCoef", local = TRUE)
					armodel2 = sfLapply(1:k, fun = function(i) {
								specx = rgarch:::arfimaspec(mean.model = list(armaOrder = c(p, 0), include.mean = TRUE, 
												arfima = FALSE, external.regressors = orex), distribution.model = "norm",
										fixed.pars = as.list(arCoef[,i]))
								ans = rgarch::arfimafilter(spec = specx, data = origdata[,i], out.sample = 0)
								return(ans)})
					sfStop()
				}
			} else{
				armodel2 = lapply(1:k, FUN = function(i){
							specx = arfimaspec(mean.model = list(armaOrder = c(p, 0), include.mean = TRUE, 
											arfima = FALSE, external.regressors = orex), distribution.model = "norm",
									fixed.pars = as.list(arCoef[,i]))
							ans = arfimafit(spec = specx, data = origdata[,i], out.sample = 0)
							return(ans)})
			}			
			zorigdata = sapply(armodel2, FUN = function(x) residuals(x))
			if( p > 0 ) zorigdata = zorigdata[-c(1:p), , drop = FALSE]
		} else{
			zorigdata = zdata
		}
	}
	A = fit@mfit$A
	W = fit@mfit$W
	#invU = solve(K)%*%W
	#U = solve(invU)
	# should we update the covariance matrix with new information?
	Y = zorigdata %*% W
	
	filterlist = multifilter(multifitORspec = fit@ufit, data = Y, out.sample = out.sample, n.old = NULL, 
			parallel = parallel, parallel.control = parallel.control)
	tmp = garchica2.nig(filterlist, A, K, N - p, mu)
	
	tmp$out.sample = out.sample
	tmp$origdata = origdata
	tmp$external.regressors = orex
	tmp$mxn = fit@mfit$mxn
	tmp$mdist = "manig"
	tmp$icamethod = fit@mfit$icamethod
	tmp$meanmodel = fit@mfit$meanmodel
	tmp$Y = Y
	tmp$W = W
	tmp$A = A
	tmp$K = fit@mfit$K
	tmp$VARmodel = vrmodel
	tmp$ARmodel = fit@mfit$ARmodel
	tmp$residuals = zorigdata
	ans = new("goGARCHfilter",
			mfilter = tmp,
			ufilter = filterlist)
	return(ans)
}

.gogarchforecast.nig = function(fit, n.ahead, n.roll, external.forecasts,  parallel = FALSE, 
		parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2), ...)
{
	ns = fit@mfit$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")
	
	m = dim(fit@mfit$origdata)[2]
	
	# 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
	}
	vrmodel = fit@mfit$VARmodel
	if( !is.null(vrmodel) ){
		mxn = vrmodel$mxn
		if( mxn > 0 ){
			if( is.null( 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("\ngogarchforecast-->error: wrong number of external regressors!...", call. = FALSE)
				if( dim(mregfor)[1] < (n.roll + n.ahead) ) stop("\ngogarchforecast-->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("\ngogarchforecast-->error: n.ahead must be equal to 1 when using n.roll\n")
		
		if( n.ahead == 1 && (n.roll > ns) ) stop("\ngogarchforecast-->error: n.roll greater than out.sample!", call. = FALSE)
		#VARf = .varpredict(fit, n.ahead, n.roll, mregfor)$Mu
		Mu = varxforecast(X = fit@mfit$origdata, Bcoef = vrmodel$Bcoef, p = vrmodel$lag, out.sample = ns, 
				n.ahead, n.roll, mregfor)
	} else{
		# arfima model forecast
			arCoef = fit@mfit$ARmodel$arCoef
			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" ){
				if(!exists("mclapply")){
					require('multicore')
				}
				arfx = mclapply(1:m, FUN = function(i){
							specx = arfimaspec(mean.model = list(armaOrder = c(fit@mfit$ARmodel$lag, 0), include.mean = TRUE, 
											arfima = FALSE, external.regressors = fit@mfit$external.regressors), distribution.model = "norm",
									fixed.pars = as.list(arCoef[,i]))
							ans = arfimaforecast(fitORspec = specx, data = fit@mfit$origdata[,i], out.sample = ns,
									n.ahead = n.ahead, n.roll = n.roll, external.forecasts = list(mregfor = mregfor))
							return(ans)}, 
						mc.cores = parallel.control$cores)
			} else{
				if(!exists("sfLapply")){
					require('snowfall')
				}
				orgd  = fit@mfit$origdata
				p = fit@mfit$ARmodel$lag
				orex = fit@mfit$external.regressors
				sfInit(parallel=TRUE, cpus = parallel.control$cores)
				sfExport("orgd", "p", "orex", "arCoef", "mregfor", "ns", "n.ahead", "n.roll", local = TRUE)
				arfx = sfLapply(1:m, fun = function(i){
							specx = rgarch::arfimaspec(mean.model = list(armaOrder = c(p, 0), include.mean = TRUE, 
											arfima = FALSE, external.regressors = orex), distribution.model = "norm",
									fixed.pars = as.list(arCoef[,i]))
							ans = rgarch::arfimaforecast(fitORspec = specx, data = orgd[,i], out.sample = ns,
									n.ahead = n.ahead, n.roll = n.roll, external.forecasts = list(mregfor = mregfor))
							return(ans)})
				sfStop()
			}
		} else{
			arfx = lapply(1:m, FUN = function(i){
						specx = arfimaspec(mean.model = list(armaOrder = c(fit@mfit$ARmodel$lag, 0), include.mean = TRUE, 
										arfima = FALSE, external.regressors = fit@mfit$external.regressors), distribution.model = "norm",
								fixed.pars = as.list(arCoef[,i]))
						ans = arfimaforecast(fitORspec = specx, data = fit@mfit$origdata[,i], out.sample = ns,
								n.ahead = n.ahead, n.roll = n.roll, external.forecasts = list(mregfor = mregfor))
						return(ans)})
		}
		Mu = matrix(NA, ncol = m, nrow = tf)
		for(i in 1:m) Mu[,i] = sapply(arfx[[i]]@forecast$forecasts, FUN = function(x) x[,1])
	}
	
	fitlist = fit@ufit
	# we augmented n.roll before, now subtract
	forclist = multiforecast(multifitORspec = fitlist, n.ahead = n.ahead, n.roll = n.roll, parallel = parallel, 
			parallel.control = parallel.control, ...)
	A = fit@mfit$A
	tmp = garchica3.nig(forclist, A, m, n.ahead, n.roll+1, Mu)
	tmp$mdist = "manig"
	tmp$A = A
	tmp$K = fit@mfit$K
	tmp$W = fit@mfit$W
	ans = new("goGARCHforecast",
			mforecast = tmp,
			uforecast = forclist)
	return(ans)
}

.gogarchsim.nig = function(fit, n.sim, n.start, m.sim, startMethod, prereturns, mexsimdata, rseed, 
		parallel = FALSE, parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2), ...)
{
	Data = head(fit@mfit$origdata, dim(fit@mfit$origdata)[1] - fit@mfit$out.sample)
	K = dim(Data)[2]
	np = dim(Data)[1]
	A = fit@mfit$A
	simlist = vector(mode = "list", length = K)
	m = fit@ufit@fit[[1]]@model$maxOrder
	if(is.null(startMethod[1])) startMethod = "unconditional" else startMethod = startMethod[1]
	if( !is.null(rseed) ){
		if( !is.matrix(rseed) ){
			xseed = matrix(NA, ncol = K, nrow = m.sim)
			if( length(rseed) != m.sim ) stop("\ngogarchsim--error: rseed must have m.sim rows!\n")
			xseed[,1] = rseed
			for(i in 2:K) xseed[, i] = as.integer(runif(1*m.sim,0,65000))
		} else{
			if(dim(rseed)[2] != K) stop("\ngogarchsim--error: rseed must have N (= n.assets) columns!\n")
			if(dim(rseed)[1] != m.sim) stop("\ngogarchsim--error: rseed must have m.sim rows!\n")
		}
	} else{
		xseed = matrix( as.integer( runif( 1*K, 0, as.integer(Sys.time()) ) ), ncol = K )
	}
	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" ){
			if(!exists("mclapply")){
				require('multicore')
			}
			simlist = mclapply(1:K, FUN = function(i) ugarchsim(fit@ufit@fit[[i]], n.sim = n.sim + n.start, 
								n.start = 0, m.sim = m.sim, startMethod = startMethod, rseed = xseed[,i]), 
					mc.cores = parallel.control$cores)
		} else{
			if(!exists("sfLapply")){
				require('snowfall')
			}
			sfInit(parallel=TRUE, cpus = parallel.control$cores)
			sfExport("fit", "n.sim", "n.start", "m.sim", "startMethod", "xseed",local = TRUE)
			simlist = sfLapply(as.list(1:K), fun = function(i) rgarch::ugarchsim(fit@ufit@fit[[i]], 
								n.sim = n.sim + n.start, n.start = 0, m.sim = m.sim, startMethod = startMethod, 
								rseed = xseed[,i]))
			if(fit@mfit$meanmodel != "VAR") sfStop()
		}
	} else{
		for(i in 1:K){
			# increment by 1 the rseed for each individual asset so that they do not get the same
			# random numbers!!!
			simlist[[i]] = ugarchsim(fit@ufit@fit[[i]], n.sim = n.sim + n.start, n.start = 0, 
					m.sim = m.sim, startMethod = startMethod, rseed = xseed[,i])
		}
	}
	
	sigmalist = vector(mode = "list", length = m.sim)
	retlist = vector(mode = "list", length = m.sim)
	residlistx = residlist = vector(mode = "list", length = m.sim)
	zlist = vector(mode = "list", length = m.sim)
	xlist = vector(mode = "list", length = m.sim)
	# we create the dlist object for use later if convolution is called.
	dlist = vector(mode = "list", length = m.sim)
	rho =  matrix( sapply(fit@ufit@fit, FUN = function(x) rep(x@fit$coef["skew"],  n.sim) ), nrow = n.sim )
	zeta = matrix( sapply(fit@ufit@fit, FUN = function(x) rep(x@fit$coef["shape"], n.sim) ), nrow = n.sim )
	
	znigparams = vector(mode = "list", length = K)
	ynigparams = vector(mode = "list", length = K)
	# move from the location-scale parametrization of the z innovations ~(0,1,rho,chi) to the
	# alpha-beta parametrization for Y (the factors = sigma*z) ~ (0, sigma, skew, shape)==(mu, delta, alpha, beta)
	
	for(i in 1:m.sim){
		sigmalist[[i]] = matrix(sapply(simlist, FUN = function(x) tail(x@simulation$sigmaSim[,i], n.sim + n.start)), ncol = K)
		residlist[[i]] = matrix(sapply(simlist, FUN = function(x) tail(x@simulation$residSim[,i], n.sim + n.start)), ncol = K)
		zlist[[i]] = matrix(residlist[[i]]/sigmalist[[i]], ncol = K)
		residlistx[[i]] = matrix(t(apply(residlist[[i]], 1, FUN = function(x) t(A)%*%x)), ncol = K)
		residlist[[i]] = tail(residlistx[[i]], n.sim)
		sigmalist[[i]] = tail(sigmalist[[i]], n.sim)
		zlist[[i]] = tail(zlist[[i]], n.sim)
	}
	

	znigparams = lapply(1:K, FUN = function(i) .nigtransform(zeta = zeta[,i], rho = rho[,i]))	
	yniglist = vector(mode = "list", length = m.sim)
	z.D = array(data = NA, dim = c(n.sim, K, 4))
	
	for(i in 1:K){
		z.D[1:(n.sim), i, 1] = znigparams[[i]][, 1]
		z.D[1:(n.sim), i, 2] = znigparams[[i]][, 2]
		z.D[1:(n.sim), i, 3] = znigparams[[i]][, 3]
		z.D[1:(n.sim), i, 4] = znigparams[[i]][, 4]
	}
	for(j in 1:m.sim){		
		for(i in 1:K){
			ynigparams[[i]] = .nigscale2(mu = 0, sigma = as.numeric(sigmalist[[j]][,i]), 
					nigalpha = znigparams[[i]][,1], nigbeta = znigparams[[i]][,2], 
					nigdelta = znigparams[[i]][,3], nigmu = znigparams[[i]][,4])
		}
		yniglist[[j]] = ynigparams
	}
	
	
	if( fit@mfit$meanmodel == "VAR" ){
		# mexsimdata check
		mxn = fit@mfit$mxn
		if(mxn > 0 ){
			if( !is.null(mexsimdata) ){
				if( !is.list(mexsimdata) | length(mexsimdata) != m.sim ) stop("\ngogarchsim-->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("\ngogarchsim-->error: mexsimdata list submatrix must be a matrix...", call. = FALSE)
					if( dim(mexsimdata[[j]])[2] != mxn ) stop("\ngogarchsim-->error: wrong number of external regressors in list submatrix!...", call. = FALSE)
					if( dim(mexsimdata[[j]])[1] < (n.sim + n.start) ) stop("\ngogarchsim-->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[[j]] = matrix(0, ncol = mxn, nrow = (n.sim + n.start))
			}
		} else{
			mexsimdata = NULL
		}
		
		vrmodel = fit@mfit$VARmodel
		p = as.numeric(vrmodel$lag)
		if(!is.na(prereturns)[1]){
			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] != K) stop("\nwrong number of cols for prereturns matrix")
		} else{
			prereturns = as.matrix( tail(Data, p) )
		}
		if( n.sim == 1 && n.start == 0 ){
			# SPECIAL ROUTINE TO SPEED UP 1-AHEAD ESTIMATION
			tmp = varxsimXX(X = Data, vrmodel$Bcoef, p = vrmodel$lag, m.sim = m.sim, prereturns, resids = t(sapply(residlistx, 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]
		} else{
			if( parallel ){
				if( parallel.control$pkg == "multicore" ){
					xlist = mclapply(as.list(1:m.sim), FUN = function(i) varxsim(X = Data, vrmodel$Bcoef, p = vrmodel$lag, 
										n.sim = n.sim, n.start = n.start, prereturns = prereturns, resids = residlistx[[i]],
										mexsimdata = mexsimdata[[i]]), mc.cores = parallel.control$cores)
				} else{
					# when the GARCH model is without mean equation, simX == simulated residuals
					sfExport("Data", "vrmodel", "n.sim", "n.start", "prereturns", "residlistx", "mexsimdata", local = TRUE)
					xlist = sfLapply(as.list(1:m.sim), fun = function(i) rgarch:::varxsim(X = Data, vrmodel$Bcoef, p = vrmodel$lag, 
										n.sim = n.sim, n.start = n.start, prereturns = prereturns, resids = residlistx[[i]],
										mexsimdata = mexsimdata[[i]]))
					sfStop()
				}
			} else{
				xlist = lapply(as.list(1:m.sim), FUN = function(i) varxsim(X = Data, vrmodel$Bcoef, p = vrmodel$lag, 
									n.sim = n.sim, n.start = n.start, prereturns = prereturns, resids = residlistx[[i]], 
									mexsimdata = mexsimdata[[i]]))
			}
		}
	} else{
		N = dim(fit@mfit$origdata)[1]
		out.sample = fit@mfit$out.sample
		mxn = fit@mfit$mxn
		if(mxn > 0 ){
			if( !is.null(mexsimdata) ){
				if( !is.list(mexsimdata) | length(mexsimdata) != m.sim ) stop("\ngogarchsim-->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("\ngogarchsim-->error: mexsimdata list submatrix must be a matrix...", call. = FALSE)
					if( dim(mexsimdata[[j]])[2] != mxn ) stop("\ngogarchsim-->error: wrong number of external regressors in list submatrix!...", call. = FALSE)
					if( dim(mexsimdata[[j]])[1] < (n.sim + n.start) ) stop("\ngogarchsim-->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[[j]] = matrix(0, ncol = mxn, nrow = (n.sim + n.start))
			}
		} else{
			mexsimdata = NULL
		}
		
		p = fit@mfit$ARmodel$lag
		if(!is.na(prereturns)[1]){
			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] != K) stop("\nwrong number of cols for prereturns matrix")
		} else{
			prereturns = as.matrix( tail(Data, p) )
		}
		
		arCoef = fit@mfit$ARmodel$arCoef
		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" ){
				if(!exists("mclapply")){
					require('multicore')
				}
				arxsim = mclapply(1:K, FUN = function(i){
							specx = arfimaspec(mean.model = list(armaOrder = c(p, 0), include.mean = TRUE, 
											arfima = FALSE, external.regressors = fit@mfit$external.regressors[1:(N - out.sample), ,drop = FALSE]), 
									distribution.model = "norm",fixed.pars = as.list(arCoef[,i]))
							ans = arfimapath(spec = specx, n.sim = n.sim, n.start = n.start, m.sim = m.sim, prereturns = prereturns[,i],
									preresiduals = tail(fit@mfit$residuals[1:(np-p),i], p), custom.dist = list(name = "sample", 
											distfit = matrix(sapply(residlistx, FUN = function(x) x[,i]), ncol = m.sim, nrow = n.sim,
													byrow = TRUE), type = "res"), mexsimdata = mexsimdata)
							return(ans)}, mc.cores = parallel.control$cores)
			} else{
				if(!exists("sfLapply")){
					require('snowfall')
				}
				if(!is.null(fit@mfit$external.regressors)){
					ex = fit@mfit$external.regressors[1:(N - out.sample), , drop = FALSE]
				} else{
					ex = fit@mfit$external.regressors
				}
				sfInit(parallel=TRUE, cpus = parallel.control$cores)
				tres = tail(fit@mfit$residuals[1:(np-p), ], p)
				sfExport("p", "ex", "arCoef", "mexsimdata", "residlistx", "n.start", "tres", 
						"n.sim", "m.sim", "prereturns", local = TRUE)
				arxsim = sfLapply(1:K, fun = function(i){
							specx = rgarch::arfimaspec(mean.model = list(armaOrder = c(p, 0), include.mean = TRUE, 
											arfima = FALSE, external.regressors = ex), distribution.model = "norm",
									fixed.pars = as.list(arCoef[,i]))
							ans = rgarch::arfimapath(spec = specx, n.sim = n.sim, n.start = n.start, m.sim = m.sim, prereturns = prereturns[,i],
									preresiduals = tres[,i], custom.dist = list(name = "sample", 
											distfit = matrix(sapply(residlistx, FUN = function(x) x[,i]), ncol = m.sim, nrow = n.sim,
													byrow = TRUE), type = "res"), mexsimdata = mexsimdata)
							return(ans)})
				sfStop()
			}
		} else{
			arxsim = lapply(1:K, FUN = function(i){
						specx = arfimaspec(mean.model = list(armaOrder = c(p, 0), include.mean = TRUE, 
										arfima = FALSE, external.regressors = fit@mfit$external.regressors[1:(N - out.sample), ,drop = FALSE]), 
								distribution.model = "norm",fixed.pars = as.list(arCoef[,i]))
						ans = arfimapath(spec = specx, n.sim = n.sim, n.start = n.start, m.sim = m.sim, prereturns = prereturns[,i],
								preresiduals = tail(fit@mfit$residuals[1:(np-p),i], p), custom.dist = list(name = "sample", 
										distfit = matrix(sapply(residlistx, FUN = function(x) x[,i]), ncol = m.sim, nrow = n.sim,
												byrow = TRUE), type = "res"), mexsimdata = mexsimdata)
						return(ans)})
		}
		for(i in 1:m.sim) xlist[[i]] = matrix(sapply(arxsim, FUN = function(x) x@path$seriesSim[,i]), ncol = K, nrow = n.sim, byrow = TRUE)
	}
	
	tmp = list()
	tmp$distribution = fit@mfit$distr$model
	
	tmp$model = "nig"
	tmp$mdist = "manig"
	tmp$z.alphabeta = znigparams
	tmp$y.alphabeta = yniglist
	tmp$z.D = z.D
	tmp$A = A
	tmp$K = fit@mfit$K
	tmp$W = fit@mfit$W
	tmp$sigmaSim = sigmalist
	tmp$residSim = residlist
	tmp$seriesSim = xlist
	tmp$zSim = zlist
	tmp$rseed = xseed
	ans = new("goGARCHsim",
			msim = tmp,
			usim = simlist)
	return(ans)

}


#-----------------------------------------------------------------------------------------------
# mvnorm
.gogarchfit.norm = function(spec, data, out.sample = 0, solver = "solnp",
		fit.control = list(stationarity = 1), solver.control = list(), parallel = FALSE, 
		parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2), VAR.fit = NULL,...)
{
	tic = Sys.time()
	
	if( is.null( colnames(data) ) ) cnames = paste("Asset_", 1:k, sep = "") else cnames = colnames(data)
	xdata = .extractmdata(data)
	if(!is.numeric(out.sample)) stop("\ngogarchfit-->error: out.sample must be numeric\n")
	if(as.numeric(out.sample) < 0) stop("\ngogarchfit-->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("\ngogarchfit-->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
	
	armodel = vrmodel = NULL
	if( spec@mspec$mmodel$demean == "VAR" ){
		if( spec@mspec$mmodel$mxn>0 ){
			ex = spec@mspec$mmodel$external.regressors
			if( dim(ex)[1] < dim(data)[1]  ) stop("\ngogarchfit-->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(data)[1]
			K = dim(data)[2]
			mu = VAR.fit$xfitted
			vrmodel = VAR.fit
			vrmodel$mxn = mxn
		} else {
			cat("\nestimating VAR model...")
			ic = spec@mspec$mmodel$lag.criterion[1]
			if( !is.null(spec@mspec$mmodel$lag.max) ){
				mp = .varxselect(y = data, lag.max = spec@mspec$mmodel$lag.max, exogen = ex)$selection
				mp = as.numeric( mp[which(substr(names(mp), 1, 2) == substr(ic, 1, 2))] )
			} else {
				mp = spec@mspec$mmodel$lag
			}
			p = mp
			vrmodel = varxfilter(X = data, p = mp, exogen = ex, Bcoef = NULL, robust = spec@mspec$mmodel$robust, 
					gamma = spec@mspec$mmodel$robust.control$gamma, delta = spec@mspec$mmodel$robust.control$delta, 
					nc = spec@mspec$mmodel$robust.control$nc, ns = spec@mspec$mmodel$robust.control$ns)	
			cat("done!\n")
			zdata = vrmodel$xresiduals
			# should be N - p
			N = dim(data)[1]
			K = dim(data)[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{
# mu-ARX or just muX
		if( spec@mspec$mmodel$mxn>0 ){
			ex = spec@mspec$mmodel$external.regressors
			if( dim(ex)[1] < dim(data)[1]  ) stop("\ngogarchfit-->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
		}
		p = mp = spec@mspec$mmodel$lag
		arspec = arfimaspec(mean.model = list(armaOrder = c(p, 0), include.mean = TRUE, 
						arfima = FALSE, external.regressors = ex), distribution.model = "norm")
		cat("\nestimating univariate mean model...")
		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" ){
				if(!exists("mclapply")){
					require('multicore')
				}
				arfit = mclapply(1:k, FUN = function(i) arfimafit(spec = arspec, 
									data = data[,i], out.sample = 0, solver = solver, 
									solver.control = solver.control, fit.control = fit.control), 
						mc.cores = parallel.control$cores)
			} else{
				if(!exists("sfLapply")){
					require('snowfall')
				}
				sfInit(parallel=TRUE, cpus = parallel.control$cores)
				sfExport("arspec", "data", "solver", "solver.control", "fit.control",local = TRUE)
				arfit = sfLapply(1:k, fun = function(i) rgarch::arfimafit(spec = arspec, 
									data = data[,i], out.sample = 0, solver = solver, 
									solver.control = solver.control, fit.control = fit.control))
				sfStop()
			}
		} else{
			arfit = lapply(1:k, FUN = function(i) arfimafit(spec = arspec, 
								data = data[,i], out.sample = 0, solver = solver, 
								solver.control = solver.control, fit.control = fit.control))
		}
		cat("done!\n")
		zdata = sapply(arfit, FUN = function(x) residuals(x))
		if( p > 0 ) zdata = zdata[-c(1:p), , drop = FALSE]
		# should be N - p
		N = dim(data)[1]
		K = dim(data)[2]
		mu = sapply(arfit, FUN = function(x) fitted(x))
		if( p > 0 ) mu = mu[-c(1:p), , drop = FALSE]
		armodel$lag = mp
		armodel$mxn = mxn
		armodel$arCoef = arCoef = sapply(arfit, FUN = function(x) coef(x))
		# last coefficient is the sigma which we ignore
		
		if(out.sample>0){
			if( parallel ){
				if( parallel.control$pkg == "multicore" ){
					if(!exists("mclapply")){
						require('multicore')
					}
					arfit2 = mclapply(1:k, FUN = function(i){
								specx = arfimaspec(mean.model = list(armaOrder = c(p, 0), include.mean = TRUE, 
												arfima = FALSE, external.regressors = orex), distribution.model = "norm",
										fixed.pars = as.list(armodel$arCoef[,i]))
								ans = arfimafilter(spec = specx, data = origdata[,i], out.sample = 0)
								return(ans)}, 
							mc.cores = parallel.control$cores)
				} else{
					if(!exists("sfLapply")){
						require('snowfall')
					}
					sfInit(parallel=TRUE, cpus = parallel.control$cores)
					sfExport("origdata", "p", "orex", "arCoef", local = TRUE)
					arfit2 = sfLapply(1:k, fun = function(i) {
								specx = rgarch:::arfimaspec(mean.model = list(armaOrder = c(p, 0), include.mean = TRUE, 
												arfima = FALSE, external.regressors = orex), distribution.model = "norm",
										fixed.pars = as.list(arCoef[,i]))
								ans = rgarch::arfimafilter(spec = specx, data = origdata[,i], out.sample = 0)
								return(ans)})
					sfStop()
				}
			} else{
				arfit2 = lapply(1:k, FUN = function(i){
							specx = arfimaspec(mean.model = list(armaOrder = c(p, 0), include.mean = TRUE, 
											arfima = FALSE, external.regressors = orex), distribution.model = "norm",
									fixed.pars = as.list(armodel$arCoef[,i]))
							ans = arfimafilter(spec = specx, data = origdata[,i], out.sample = 0)
							return(ans)})
			}
			zorigdata = sapply(arfit2, FUN = function(x) residuals(x))
			if( p > 0 ) zorigdata = zorigdata[-c(1:p), , drop = FALSE]
		} else{
			zorigdata = zdata
		}
	}
	# iidret = .makeiid(zdata, method = spec@mspec$ica, spec@mspec$ica.fix, gfun  ="tanh")
	iidret = .makeiid(zdata, method = spec@mspec$ica, ica.fix = spec@mspec$ica.fix, ...)
	A = iidret$A
	Y = iidret$Y
	W = iidret$W
	
	speclist = multispec( replicate(n = K, spec@uspec) )
	
	Y = zorigdata %*% W
	
	fitlist = multifit(multispec = speclist, data = Y, out.sample = out.sample, solver = solver, 
			solver.control = solver.control, fit.control = fit.control, parallel = parallel, parallel.control = parallel.control)
	if(out.sample>0){
		filterlist = multifilter(multifitORspec = fitlist, data = Y, out.sample = 0, parallel = parallel, parallel.control = parallel.control)
	} else{
		filterlist = list()
	}
	
	# check all converged:
	# we return the fitted list to the global environment in case there is some non convergence
	# for debugging...
	#.fitlist <<- fitlist
	converge = sapply(fitlist@fit, FUN = function(x) x@fit$convergence)
	if(any(converge == 1)){
		# it's useless returning any more info since we have an ICA factor matrix for which the
		# user can do very little other than probably run it again hoping that the new ICA matrix
		# will return a 'nicer' set (perhaps return the random number seed used to initialize ICA).
		stop("\ngogarchdfit-->warning: convergence problem in univariate fit...", call. = FALSE)
	}
	tmp = garchica1.norm(fitlist, A, K, N - p, mu)
	
	tmp$out.sample = out.sample
	tmp$origdata = origdata
	tmp$external.regressors = orex
	tmp$mxn = spec@mspec$mmodel$mxn
	tmp$mdist = "mvnorm"
	tmp$icamethod = spec@mspec$ica
	tmp$meanmodel = spec@mspec$mmodel$model
	tmp$Y = Y
	tmp$W = iidret$W
	tmp$A = A
	tmp$K = iidret$K
	tmp$VARmodel = vrmodel
	tmp$ARmodel = armodel
	tmp$residuals = zorigdata
	tmp$filterlist = filterlist
	tmp$timer = Sys.time() - tic
	
	ans = new("goGARCHfit",

			mfit = tmp,
			ufit = fitlist)
	return(ans)
}

.gogarchfilter.norm = function(fit, data, out.sample, n.old, parallel = FALSE, parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2),...)
{
	if( is.null( colnames(data) ) ) cnames = paste("Asset_", 1:k, sep = "") else cnames = colnames(data)
	xdata = .extractmdata(data)
	if(!is.numeric(out.sample)) stop("\ngogarchfilter-->error: out.sample must be numeric\n")
	if(as.numeric(out.sample) < 0) stop("\ngogarchfilter-->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("\ngogarchfilter-->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
	
	
	armodel = vrmodel = NULL
	if( fit@mfit$meanmodel == "VAR" ){
		if( fit@mfit$mxn>0 ){
			ex = fit@mfit$external.regressors
			if( dim(ex)[1] < dim(data)[1]  ) stop("\ngogarchfilter-->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
		}
		cat("\nfiltering VAR model...")
		vrmodel = varxfilter(X = data, p = fit@mfit$VARmodel$lag, exogen = ex, Bcoef = fit@mfit$VARmodel$Bcoef)
		cat("done!\n")
		zdata = vrmodel$xresiduals
		# should be N - p
		N = dim(data)[1]
		K = dim(data)[2]
		mu = vrmodel$xfitted
		p = fit@mfit$VARmodel$lag
		# if we have out of sample data we need to filter
		if(out.sample > 0){
			Bcoef = vrmodel$Bcoef
			vrmodel2 = varxfilter(X = origdata, p = fit@mfit$VARmodel$lag, exogen = orex, Bcoef = Bcoef)
			zorigdata = vrmodel2$xresiduals			
		} else{
			zorigdata = zdata
		}
	} else{
		arCoef = fit@mfit$ARmodel$arCoef
		if( fit@mfit$mxn>0 ){
			ex = fit@mfit$external.regressors
			if( dim(ex)[1] < dim(data)[1]  ) stop("\ngogarchfilter-->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
		}
		p = fit@mfit$VARmodel$lag
		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" ){
				if(!exists("mclapply")){
					require('multicore')
				}
				armodel = mclapply(1:k, FUN = function(i){
							specx = arfimaspec(mean.model = list(armaOrder = c(p, 0), include.mean = TRUE, 
											arfima = FALSE, external.regressors = ex), distribution.model = "norm",
									fixed.pars = as.list(arCoef[,i]))
							ans = arfimafilter(spec = specx, data = data[,i], out.sample = 0)
							return(ans)}, 
						mc.cores = parallel.control$cores)
			} else{
				if(!exists("sfLapply")){
					require('snowfall')
				}
				sfInit(parallel=TRUE, cpus = parallel.control$cores)
				sfExport("data", "p", "ex", "arCoef", local = TRUE)
				armodel = sfLapply(1:k, fun = function(i) {
							specx = rgarch:::arfimaspec(mean.model = list(armaOrder = c(p, 0), include.mean = TRUE, 
											arfima = FALSE, external.regressors = ex), distribution.model = "norm",
									fixed.pars = as.list(arCoef[,i]))
							ans = rgarch::arfimafilter(spec = specx, data = data[,i], out.sample = 0)
							return(ans)})
				sfStop()
			}
		} else{
			armodel = lapply(1:k, FUN = function(i){
						specx = arfimaspec(mean.model = list(armaOrder = c(p, 0), include.mean = TRUE, 
										arfima = FALSE, external.regressors = ex), distribution.model = "norm",
								fixed.pars = as.list(arCoef[,i]))
						ans = arfimafit(spec = specx, data = data[,i], out.sample = 0)
						return(ans)})
		}
		cat("done!\n")
		zdata = sapply(armodel, FUN = function(x) residuals(x))
		if( p > 0 ) zdata = zdata[-c(1:p), , drop = FALSE]
		# should be N - p
		N = dim(data)[1]
		K = dim(data)[2]
		mu = sapply(armodel, FUN = function(x) fitted(x))
		if( p > 0 ) mu = mu[-c(1:p), , drop = FALSE]		
		# if we have out of sample data we need to filter
		if(out.sample > 0){
			if( parallel ){
				if( parallel.control$pkg == "multicore" ){
					if(!exists("mclapply")){
						require('multicore')
					}
					armodel2 = mclapply(1:k, FUN = function(i){
								specx = arfimaspec(mean.model = list(armaOrder = c(p, 0), include.mean = TRUE, 
												arfima = FALSE, external.regressors = orex), distribution.model = "norm",
										fixed.pars = as.list(arCoef[,i]))
								ans = arfimafilter(spec = specx, data = origdata[,i], out.sample = 0)
								return(ans)}, 
							mc.cores = parallel.control$cores)
				} else{
					if(!exists("sfLapply")){
						require('snowfall')
					}
					sfInit(parallel=TRUE, cpus = parallel.control$cores)
					sfExport("origdata", "p", "orex", "arCoef", local = TRUE)
					armodel2 = sfLapply(1:k, fun = function(i) {
								specx = rgarch:::arfimaspec(mean.model = list(armaOrder = c(p, 0), include.mean = TRUE, 
												arfima = FALSE, external.regressors = orex), distribution.model = "norm",
										fixed.pars = as.list(arCoef[,i]))
								ans = rgarch::arfimafilter(spec = specx, data = origdata[,i], out.sample = 0)
								return(ans)})
					sfStop()
				}
			} else{
				armodel2 = lapply(1:k, FUN = function(i){
							specx = arfimaspec(mean.model = list(armaOrder = c(p, 0), include.mean = TRUE, 
											arfima = FALSE, external.regressors = orex), distribution.model = "norm",
									fixed.pars = as.list(arCoef[,i]))
							ans = arfimafit(spec = specx, data = origdata[,i], out.sample = 0)
							return(ans)})
			}			
			zorigdata = sapply(armodel2, FUN = function(x) residuals(x))
			if( p > 0 ) zorigdata = zorigdata[-c(1:p), , drop = FALSE]
		} else{
			zorigdata = zdata
		}
	}
	A = fit@mfit$A
	W = fit@mfit$W
	#invU = solve(K)%*%W
	#U = solve(invU)
	# should we update the covariance matrix with new information?
	Y = zorigdata %*% W
	
	filterlist = multifilter(multifitORspec = fit@ufit, data = Y, out.sample = out.sample, n.old = NULL, 
			parallel = parallel, parallel.control = parallel.control)
	tmp = garchica2.norm(filterlist, A, K, N - p, mu)
	
	tmp$out.sample = out.sample
	tmp$origdata = origdata
	tmp$external.regressors = orex
	tmp$mxn = fit@mfit$mxn
	tmp$mdist = "mvnorm"
	tmp$icamethod = fit@mfit$icamethod
	tmp$meanmodel = fit@mfit$meanmodel
	tmp$Y = Y
	tmp$W = W
	tmp$A = A
	tmp$K = fit@mfit$K
	tmp$VARmodel = vrmodel
	tmp$ARmodel = fit@mfit$ARmodel
	tmp$residuals = zorigdata
	
	ans = new("goGARCHfilter",
			mfilter = tmp,
			ufilter = filterlist)
	return(ans)
}


.gogarchforecast.norm = function(fit, n.ahead, n.roll, external.forecasts, parallel = FALSE, 
		parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2), ...)
{
	ns = fit@mfit$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")
	
	m = dim(fit@mfit$origdata)[2]
	
	# 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
	}
	vrmodel = fit@mfit$VARmodel
	if( !is.null(vrmodel) ){
		mxn = vrmodel$mxn
		if( mxn > 0 ){
			if( is.null( 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("\ngogarchforecast-->error: wrong number of external regressors!...", call. = FALSE)
				if( dim(mregfor)[1] < (n.roll + n.ahead) ) stop("\ngogarchforecast-->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("\ngogarchforecast-->error: n.ahead must be equal to 1 when using n.roll\n")
		
		if( n.ahead == 1 && (n.roll > ns) ) stop("\ngogarchforecast-->error: n.roll greater than out.sample!", call. = FALSE)
		#VARf = .varpredict(fit, n.ahead, n.roll, mregfor)$Mu
		Mu = varxforecast(X = fit@mfit$origdata, Bcoef = vrmodel$Bcoef, p = vrmodel$lag, out.sample = ns, 
				n.ahead, n.roll, mregfor)
	} else{
		# arfima model forecast
		arCoef = fit@mfit$ARmodel$arCoef
		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" ){
				if(!exists("mclapply")){
					require('multicore')
				}
				arfx = mclapply(1:m, FUN = function(i){
							specx = arfimaspec(mean.model = list(armaOrder = c(fit@mfit$ARmodel$lag, 0), include.mean = TRUE, 
											arfima = FALSE, external.regressors = fit@mfit$external.regressors), distribution.model = "norm",
									fixed.pars = as.list(arCoef[,i]))
							ans = arfimaforecast(fitORspec = specx, data = fit@mfit$origdata[,i], out.sample = ns,
									n.ahead = n.ahead, n.roll = n.roll, external.forecasts = list(mregfor = mregfor))
							return(ans)}, 
						mc.cores = parallel.control$cores)
			} else{
				if(!exists("sfLapply")){
					require('snowfall')
				}
				orgd  = fit@mfit$origdata
				p = fit@mfit$ARmodel$lag
				orex = fit@mfit$external.regressors
				sfInit(parallel=TRUE, cpus = parallel.control$cores)
				sfExport("orgd", "p", "orex", "arCoef", "mregfor", "ns", "n.ahead", "n.roll", local = TRUE)
				arfx = sfLapply(1:m, fun = function(i){
							specx = rgarch::arfimaspec(mean.model = list(armaOrder = c(p, 0), include.mean = TRUE, 
											arfima = FALSE, external.regressors = orex), distribution.model = "norm",
									fixed.pars = as.list(arCoef[,i]))
							ans = rgarch::arfimaforecast(fitORspec = specx, data = orgd[,i], out.sample = ns,
									n.ahead = n.ahead, n.roll = n.roll, external.forecasts = list(mregfor = mregfor))
							return(ans)})
				sfStop()
			}
		} else{
			arfx = lapply(1:m, FUN = function(i){
						specx = arfimaspec(mean.model = list(armaOrder = c(fit@mfit$ARmodel$lag, 0), include.mean = TRUE, 
										arfima = FALSE, external.regressors = fit@mfit$external.regressors), distribution.model = "norm",
								fixed.pars = as.list(arCoef[,i]))
						ans = arfimaforecast(fitORspec = specx, data = fit@mfit$origdata[,i], out.sample = ns,
								n.ahead = n.ahead, n.roll = n.roll, external.forecasts = list(mregfor = mregfor))
						return(ans)})
		}
		Mu = matrix(NA, ncol = m, nrow = tf)
		for(i in 1:m) Mu[,i] = sapply(arfx[[i]]@forecast$forecasts, FUN = function(x) x[,1])
	}
	
	fitlist = fit@ufit
	# we augmented n.roll before, now subtract
	forclist = multiforecast(multifitORspec = fitlist, n.ahead = n.ahead, n.roll = n.roll, parallel = parallel, 
			parallel.control = parallel.control, ...)
	A = fit@mfit$A
	tmp = garchica3.norm(forclist, A, m, n.ahead, n.roll+1, Mu)
	tmp$mdist = "mvnorm"
	tmp$A = A
	tmp$K = fit@mfit$K
	tmp$W = fit@mfit$W
	ans = new("goGARCHforecast",
			mforecast = tmp,
			uforecast = forclist)
	return(ans)
}

.gogarchsim.norm = function(fit, n.sim, n.start, m.sim, startMethod, prereturns, mexsimdata, rseed, 
		parallel, parallel.control, ...)
{
	Data = head(fit@mfit$origdata, dim(fit@mfit$origdata)[1] - fit@mfit$out.sample)
	K = dim(Data)[2]
	np = dim(Data)[1]
	A = fit@mfit$A
	simlist = vector(mode = "list", length = K)
	m = fit@ufit@fit[[1]]@model$maxOrder
	if(is.null(startMethod[1])) startMethod = "unconditional" else startMethod = startMethod[1]
	
	if( !is.null(rseed) ){
		if( !is.matrix(rseed) ){
			xseed = matrix(NA, ncol = K, nrow = m.sim)
			if( length(rseed) != m.sim ) stop("\ngogarchsim--error: rseed must have m.sim rows!\n")
			xseed[,1] = rseed
			for(i in 2:K) xseed[, i] = as.integer(runif(1*m.sim,0,65000))
		} else{
			if(dim(rseed)[2] != K) stop("\ngogarchsim--error: rseed must have N (= n.assets) columns!\n")
			if(dim(rseed)[1] != m.sim) stop("\ngogarchsim--error: rseed must have m.sim rows!\n")
		}
	} else{
		xseed = matrix( as.integer( runif( 1*K, 0, as.integer(Sys.time()) ) ), ncol = K )
	}
	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" ){
			if(!exists("mclapply")){
				require('multicore')
			}
			simlist = mclapply(1:K, FUN = function(i) ugarchsim(fit@ufit@fit[[i]], n.sim = n.sim + n.start, 
								n.start = 0, m.sim = m.sim, startMethod = startMethod, rseed = xseed[,i]), 
					mc.cores = parallel.control$cores)
		} else{
			if(!exists("sfLapply")){
				require('snowfall')
			}
			sfInit(parallel=TRUE, cpus = parallel.control$cores)
			sfExport("fit", "n.sim", "n.start", "m.sim", "startMethod", "xseed",local = TRUE)
			simlist = sfLapply(as.list(1:K), fun = function(i) rgarch::ugarchsim(fit@ufit@fit[[i]], 
								n.sim = n.sim + n.start, n.start = 0, m.sim = m.sim, startMethod = startMethod, 
								rseed = xseed[,i]))
			if(fit@mfit$meanmodel != "VAR") sfStop()
		}
	} else{
		for(i in 1:K){
			# increment by 1 the rseed for each individual asset so that they do not get the same
			# random numbers!!!
			simlist[[i]] = ugarchsim(fit@ufit@fit[[i]], n.sim = n.sim + n.start, n.start = 0, 
					m.sim = m.sim, startMethod = startMethod, rseed = xseed[,i])
		}
	}
	
	sigmalist = vector(mode = "list", length = m.sim)
	retlist = vector(mode = "list", length = m.sim)
	residlistx = residlist = vector(mode = "list", length = m.sim)
	zlist = vector(mode = "list", length = m.sim)
	xlist = vector(mode = "list", length = m.sim)
	# we create the dlist object for use later if convolution is called.
	dlist = vector(mode = "list", length = m.sim)
	#rho =  sapply(fit@ufit@fit, FUN = function(x) rep(x@fit$coef["skew"],  n.sim) )
	#zeta = sapply(fit@ufit@fit, FUN = function(x) rep(x@fit$coef["shape"], n.sim) )
	
	znigparams = vector(mode = "list", length = K)
	ynigparams = vector(mode = "list", length = K)
	# move from the location-scale parametrization of the z innovations ~(0,1,rho,chi) to the
	# alpha-beta parametrization for Y (the factors = sigma*z) ~ (0, sigma, skew, shape)==(mu, delta, alpha, beta)
	
	for(i in 1:m.sim){
		sigmalist[[i]] = matrix(sapply(simlist, FUN = function(x) tail(x@simulation$sigmaSim[,i], n.sim + n.start)), ncol = K)
		residlist[[i]] = matrix(sapply(simlist, FUN = function(x) tail(x@simulation$residSim[,i], n.sim + n.start)), ncol = K)
		zlist[[i]] = matrix(residlist[[i]]/sigmalist[[i]], ncol = K)
		residlistx[[i]] = matrix(t(apply(residlist[[i]], 1, FUN = function(x) t(A)%*%x)), ncol = K)
		residlist[[i]] = tail(residlistx[[i]], n.sim)
		sigmalist[[i]] = tail(sigmalist[[i]], n.sim)
		zlist[[i]] = tail(zlist[[i]], n.sim)
	}
	
	znigparams = ynigparams = vector(mode = "list", length = K)
	
	for(i in 1:K){
		ynigparams[[i]] = znigparams[[i]] = lapply(sigmalist, FUN = function(x) as.matrix(data.frame(sigma = x[,i], mu = rep(0, n.sim))))
	}
	
	#x.R[i,]   = mu[i,] + t(A)%*%diag(sigmas[i,])%*%zres[i,]
	if( fit@mfit$meanmodel == "VAR" ){
		# mexsimdata check
		mxn = fit@mfit$mxn
		if(mxn > 0 ){
			if( !is.null(mexsimdata) ){
				if( !is.list(mexsimdata) | length(mexsimdata) != m.sim ) stop("\ngogarchsim-->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("\ngogarchsim-->error: mexsimdata list submatrix must be a matrix...", call. = FALSE)
					if( dim(mexsimdata[[j]])[2] != mxn ) stop("\ngogarchsim-->error: wrong number of external regressors in list submatrix!...", call. = FALSE)
					if( dim(mexsimdata[[j]])[1] < (n.sim + n.start) ) stop("\ngogarchsim-->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[[j]] = matrix(0, ncol = mxn, nrow = (n.sim + n.start))
			}
		} else{
			mexsimdata = NULL
		}
		
		vrmodel = fit@mfit$VARmodel
		p = as.numeric(vrmodel$lag)
		if(!is.na(prereturns)[1]){
			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] != K) stop("\nwrong number of cols for prereturns matrix")
		} else{
			prereturns = as.matrix( tail(Data, p) )
		}
		
		if( n.sim == 1 && n.start == 0 ){
			# SPECIAL ROUTINE TO SPEED UP 1-AHEAD ESTIMATION
			tmp = varxsimXX(X = Data, vrmodel$Bcoef, p = vrmodel$lag, m.sim = m.sim, prereturns, resids = t(sapply(residlistx, 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]
		} else{
			if( parallel ){
				if( parallel.control$pkg == "multicore" ){
					xlist = mclapply(as.list(1:m.sim), FUN = function(i) varxsim(X = Data, vrmodel$Bcoef, p = vrmodel$lag, 
										n.sim = n.sim, n.start = n.start, prereturns = prereturns, resids = residlistx[[i]],
										mexsimdata = mexsimdata[[i]]), mc.cores = parallel.control$cores)
				} else{
					# when the GARCH model is without mean equation, simX == simulated residuals
					sfExport("Data", "vrmodel", "n.sim", "n.start", "prereturns", "residlistx", "mexsimdata", local = TRUE)
					xlist = sfLapply(as.list(1:m.sim), fun = function(i) rgarch:::varxsim(X = Data, vrmodel$Bcoef, p = vrmodel$lag, 
										n.sim = n.sim, n.start = n.start, prereturns = prereturns, resids = residlistx[[i]],
										mexsimdata = mexsimdata[[i]]))
					sfStop()
				}
			} else{
				xlist = lapply(as.list(1:m.sim), FUN = function(i) varxsim(X = Data, vrmodel$Bcoef, p = vrmodel$lag, 
									n.sim = n.sim, n.start = n.start, prereturns = prereturns, resids = residlistx[[i]], 
									mexsimdata = mexsimdata[[i]]))
			}
		}
	} else{
		N = dim(fit@mfit$origdata)[1]
		out.sample = fit@mfit$out.sample
		mxn = fit@mfit$mxn
		if(mxn > 0 ){
			if( !is.null(mexsimdata) ){
				if( !is.list(mexsimdata) | length(mexsimdata) != m.sim ) stop("\ngogarchsim-->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("\ngogarchsim-->error: mexsimdata list submatrix must be a matrix...", call. = FALSE)
					if( dim(mexsimdata[[j]])[2] != mxn ) stop("\ngogarchsim-->error: wrong number of external regressors in list submatrix!...", call. = FALSE)
					if( dim(mexsimdata[[j]])[1] < (n.sim + n.start) ) stop("\ngogarchsim-->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[[j]] = matrix(0, ncol = mxn, nrow = (n.sim + n.start))
			}
		} else{
			mexsimdata = NULL
		}
		
		p = fit@mfit$ARmodel$lag
		if(!is.na(prereturns)[1]){
			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] != K) stop("\nwrong number of cols for prereturns matrix")
		} else{
			prereturns = as.matrix( tail(Data, p) )
		}
		
		arCoef = fit@mfit$ARmodel$arCoef
		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" ){
				if(!exists("mclapply")){
					require('multicore')
				}
				arxsim = mclapply(1:K, FUN = function(i){
							specx = arfimaspec(mean.model = list(armaOrder = c(p, 0), include.mean = TRUE, 
											arfima = FALSE, external.regressors = fit@mfit$external.regressors[1:(N - out.sample), ,drop = FALSE]), 
									distribution.model = "norm",fixed.pars = as.list(arCoef[,i]))
							ans = arfimapath(spec = specx, n.sim = n.sim, n.start = n.start, m.sim = m.sim, prereturns = prereturns[,i],
									preresiduals = tail(fit@mfit$residuals[1:(np-p),i], p), custom.dist = list(name = "sample", 
											distfit = matrix(sapply(residlistx, FUN = function(x) x[,i]), ncol = m.sim, nrow = n.sim,
													byrow = TRUE), type = "res"), mexsimdata = mexsimdata)
							return(ans)}, mc.cores = parallel.control$cores)
			} else{
				if(!exists("sfLapply")){
					require('snowfall')
				}
				if(!is.null(fit@mfit$external.regressors)){
					ex = fit@mfit$external.regressors[1:(N - out.sample), , drop = FALSE]
				} else{
					ex = fit@mfit$external.regressors
				}
				sfInit(parallel=TRUE, cpus = parallel.control$cores)
				tres = tail(fit@mfit$residuals[1:(np-p),], p)
				sfExport("p", "ex", "arCoef", "mexsimdata", "residlistx", "n.start", "tres", 
						"n.sim", "m.sim", "prereturns", local = TRUE)
				arxsim = sfLapply(1:K, fun = function(i){
							specx = rgarch::arfimaspec(mean.model = list(armaOrder = c(p, 0), include.mean = TRUE, 
											arfima = FALSE, external.regressors = ex), distribution.model = "norm",
									fixed.pars = as.list(arCoef[,i]))
							ans = rgarch::arfimapath(spec = specx, n.sim = n.sim, n.start = n.start, m.sim = m.sim, prereturns = prereturns[,i],
									preresiduals = tres[,i], custom.dist = list(name = "sample", 
											distfit = matrix(sapply(residlistx, FUN = function(x) x[,i]), ncol = m.sim, nrow = n.sim,
													byrow = TRUE), type = "res"), mexsimdata = mexsimdata)
							return(ans)})
				sfStop()
			}
		} else{
			arxsim = lapply(1:K, FUN = function(i){
						specx = arfimaspec(mean.model = list(armaOrder = c(p, 0), include.mean = TRUE, 
										arfima = FALSE, external.regressors = fit@mfit$external.regressors[1:(N - out.sample), ,drop = FALSE]), 
								distribution.model = "norm",fixed.pars = as.list(arCoef[,i]))
						ans = arfimapath(spec = specx, n.sim = n.sim, n.start = n.start, m.sim = m.sim, prereturns = prereturns[,i],
								preresiduals = tail(fit@mfit$residuals[1:(np-p),i], p), custom.dist = list(name = "sample", 
										distfit = matrix(sapply(residlistx, FUN = function(x) x[,i]), ncol = m.sim, nrow = n.sim,
												byrow = TRUE), type = "res"), mexsimdata = mexsimdata)
						return(ans)})
		}
		for(i in 1:m.sim) xlist[[i]] = matrix(sapply(arxsim, FUN = function(x) x@path$seriesSim[,i]), ncol = K, nrow = n.sim, byrow = TRUE)
	}
	
	tmp = list()
	tmp$distribution = fit@mfit$distr$model
	
	tmp$model = "norm"
	tmp$mdist = "mvnorm"
	tmp$z.alphabeta = znigparams
	tmp$y.alphabeta = ynigparams
	tmp$z.D = NULL
	tmp$A = A
	tmp$K = fit@mfit$K
	tmp$W = fit@mfit$W
	tmp$sigmaSim = sigmalist
	tmp$residSim = residlist
	tmp$seriesSim = xlist
	tmp$zSim = zlist
	tmp$rseed = xseed
	
	ans = new("goGARCHsim",
			msim = tmp,
			usim = simlist)
	return(ans)
}

# magh
.gogarchfit.ghyp = function(spec, data, out.sample = 0, solver = "solnp",
		fit.control = list(stationarity = 1), solver.control = list(), parallel = FALSE, 
		parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2), VAR.fit = NULL, ...)
{
	tic = Sys.time()
	
	if( is.null( colnames(data) ) ) cnames = paste("Asset_", 1:k, sep = "") else cnames = colnames(data)
	xdata = .extractmdata(data)
	if(!is.numeric(out.sample)) stop("\ngogarchfit-->error: out.sample must be numeric\n")
	if(as.numeric(out.sample) < 0) stop("\ngogarchfit-->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("\ngogarchfit-->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
	
	armodel = vrmodel = NULL
	if( spec@mspec$mmodel$demean == "VAR" ){
		if( spec@mspec$mmodel$mxn>0 ){
			ex = spec@mspec$mmodel$external.regressors
			if( dim(ex)[1] < dim(data)[1]  ) stop("\ngogarchfit-->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(data)[1]
			K = dim(data)[2]
			mu = VAR.fit$xfitted
			vrmodel = VAR.fit
			vrmodel$mxn = mxn
		} else {
			cat("\nestimating VAR model...")
			ic = spec@mspec$mmodel$lag.criterion[1]
			if( !is.null(spec@mspec$mmodel$lag.max) ){
				mp = .varxselect(y = data, lag.max = spec@mspec$mmodel$lag.max, exogen = ex)$selection
				mp = as.numeric( mp[which(substr(names(mp), 1, 2) == substr(ic, 1, 2))] )
			} else {
				mp = spec@mspec$mmodel$lag
			}
			p = mp
			vrmodel = varxfilter(X = data, p = mp, exogen = ex, Bcoef = NULL, robust = spec@mspec$mmodel$robust, 
					gamma = spec@mspec$mmodel$robust.control$gamma, delta = spec@mspec$mmodel$robust.control$delta, 
					nc = spec@mspec$mmodel$robust.control$nc, ns = spec@mspec$mmodel$robust.control$ns)	
			cat("done!\n")
			zdata = vrmodel$xresiduals
			# should be N - p
			N = dim(data)[1]
			K = dim(data)[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{
# mu-ARX or just muX
		if( spec@mspec$mmodel$mxn>0 ){
			ex = spec@mspec$mmodel$external.regressors
			if( dim(ex)[1] < dim(data)[1]  ) stop("\ngogarchfit-->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
		}
		p = mp = spec@mspec$mmodel$lag
		arspec = arfimaspec(mean.model = list(armaOrder = c(p, 0), include.mean = TRUE, 
						arfima = FALSE, external.regressors = ex), distribution.model = "norm")
		cat("\nestimating univariate mean model...")
		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" ){
				if(!exists("mclapply")){
					require('multicore')
				}
				arfit = mclapply(1:k, FUN = function(i) arfimafit(spec = arspec, 
									data = data[,i], out.sample = 0, solver = solver, 
									solver.control = solver.control, fit.control = fit.control), 
						mc.cores = parallel.control$cores)
			} else{
				if(!exists("sfLapply")){
					require('snowfall')
				}
				sfInit(parallel=TRUE, cpus = parallel.control$cores)
				sfExport("arspec", "data", "solver", "solver.control", "fit.control",local = TRUE)
				arfit = sfLapply(1:k, fun = function(i) rgarch::arfimafit(spec = arspec, 
									data = data[,i], out.sample = 0, solver = solver, 
									solver.control = solver.control, fit.control = fit.control))
				sfStop()
			}
		} else{
			arfit = lapply(1:k, FUN = function(i) arfimafit(spec = arspec, 
								data = data[,i], out.sample = 0, solver = solver, 
								solver.control = solver.control, fit.control = fit.control))
		}
		cat("done!\n")
		zdata = sapply(arfit, FUN = function(x) residuals(x))
		if( p > 0 ) zdata = zdata[-c(1:p), , drop = FALSE]
		# should be N - p
		N = dim(data)[1]
		K = dim(data)[2]
		mu = sapply(arfit, FUN = function(x) fitted(x))
		if( p > 0 ) mu = mu[-c(1:p), , drop = FALSE]
		armodel$lag = mp
		armodel$mxn = mxn
		armodel$arCoef = arCoef = sapply(arfit, FUN = function(x) coef(x))
		# last coefficient is the sigma which we ignore
		
		if(out.sample>0){
			if( parallel ){
				if( parallel.control$pkg == "multicore" ){
					if(!exists("mclapply")){
						require('multicore')
					}
					arfit2 = mclapply(1:k, FUN = function(i){
								specx = arfimaspec(mean.model = list(armaOrder = c(p, 0), include.mean = TRUE, 
												arfima = FALSE, external.regressors = orex), distribution.model = "norm",
										fixed.pars = as.list(armodel$arCoef[,i]))
								ans = arfimafilter(spec = specx, data = origdata[,i], out.sample = 0)
								return(ans)}, 
							mc.cores = parallel.control$cores)
				} else{
					if(!exists("sfLapply")){
						require('snowfall')
					}
					sfInit(parallel=TRUE, cpus = parallel.control$cores)
					sfExport("origdata", "p", "orex", "arCoef", local = TRUE)
					arfit2 = sfLapply(1:k, fun = function(i) {
								specx = rgarch:::arfimaspec(mean.model = list(armaOrder = c(p, 0), include.mean = TRUE, 
												arfima = FALSE, external.regressors = orex), distribution.model = "norm",
										fixed.pars = as.list(arCoef[,i]))
								ans = rgarch::arfimafilter(spec = specx, data = origdata[,i], out.sample = 0)
								return(ans)})
					sfStop()
				}
			} else{
				arfit2 = lapply(1:k, FUN = function(i){
							specx = arfimaspec(mean.model = list(armaOrder = c(p, 0), include.mean = TRUE, 
											arfima = FALSE, external.regressors = orex), distribution.model = "norm",
									fixed.pars = as.list(armodel$arCoef[,i]))
							ans = arfimafilter(spec = specx, data = origdata[,i], out.sample = 0)
							return(ans)})
			}
			zorigdata = sapply(arfit2, FUN = function(x) residuals(x))
			if( p > 0 ) zorigdata = zorigdata[-c(1:p), , drop = FALSE]
		} else{
			zorigdata = zdata
		}
	}
	
	iidret = .makeiid(zdata, method = spec@mspec$ica, ica.fix = spec@mspec$ica.fix, ...)
	A = iidret$A
	Y = iidret$Y
	W = iidret$W
	
	speclist = multispec( replicate(n = K, spec@uspec) )
	
	Y = zorigdata %*% W
	
	fitlist = multifit(multispec = speclist, data = Y, out.sample = out.sample, solver = solver, 
			solver.control = solver.control, fit.control = fit.control, parallel = parallel, parallel.control = parallel.control)
	if(out.sample>0){
		filterlist = multifilter(multifitORspec = fitlist, data = Y, out.sample = 0, parallel = parallel, parallel.control = parallel.control)
	} else{
		filterlist = list()
	}
	
	# check all converged:
	# we return the fitted list to the global environment in case there is some non convergence
	# for debugging...
	#.fitlist <<- fitlist
	converge = sapply(fitlist@fit, FUN = function(x) x@fit$convergence)
	if(any(converge == 1)){
		# it's useless returning any more info since we have an ICA factor matrix for which the
		# user can do very little other than probably run it again hoping that the new ICA matrix
		# will return a 'nicer' set (perhaps return the random number seed used to initialize ICA).
		stop("\ngogarchdfit-->warning: convergence problem in univariate fit...", call. = FALSE)
	}
	tmp = garchica1.ghyp(fitlist, A, K, N - p, mu)
	
	tmp$out.sample = out.sample
	tmp$origdata = origdata
	tmp$external.regressors = orex
	tmp$mxn = spec@mspec$mmodel$mxn
	tmp$mdist = "magh"
	tmp$icamethod = spec@mspec$ica
	tmp$meanmodel = spec@mspec$mmodel$model	
	tmp$Y = Y
	tmp$W = iidret$W
	tmp$A = A
	tmp$K = iidret$K
	tmp$VARmodel = vrmodel
	tmp$ARmodel = armodel
	tmp$residuals = zorigdata
	tmp$filterlist = filterlist
	tmp$timer = Sys.time() - tic
	
	ans = new("goGARCHfit",
			mfit = tmp,
			ufit = fitlist)
	return(ans)
}
# need to re-write it like the DCC method
.gogarchfilter.ghyp = function(fit, data, out.sample, n.old, parallel = FALSE, 
		parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2),...)
{
	if( is.null( colnames(data) ) ) cnames = paste("Asset_", 1:k, sep = "") else cnames = colnames(data)
	xdata = .extractmdata(data)
	if(!is.numeric(out.sample)) stop("\ngogarchfilter-->error: out.sample must be numeric\n")
	if(as.numeric(out.sample) < 0) stop("\ngogarchfilter-->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("\ngogarchfilter-->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
	
	
	armodel = vrmodel = NULL
	if( fit@mfit$meanmodel == "VAR" ){
		if( fit@mfit$mxn>0 ){
			ex = fit@mfit$external.regressors
			if( dim(ex)[1] < dim(data)[1]  ) stop("\ngogarchfilter-->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
		}
		cat("\nfiltering VAR model...")
		vrmodel = varxfilter(X = data, p = fit@mfit$VARmodel$lag, exogen = ex, Bcoef = fit@mfit$VARmodel$Bcoef)
		cat("done!\n")
		zdata = vrmodel$xresiduals
		# should be N - p
		N = dim(data)[1]
		K = dim(data)[2]
		mu = vrmodel$xfitted
		p = fit@mfit$VARmodel$lag
		# if we have out of sample data we need to filter
		if(out.sample > 0){
			Bcoef = vrmodel$Bcoef
			vrmodel2 = varxfilter(X = origdata, p = fit@mfit$VARmodel$lag, exogen = orex, Bcoef = Bcoef)
			zorigdata = vrmodel2$xresiduals			
		} else{
			zorigdata = zdata
		}
	} else{
		arCoef = fit@mfit$ARmodel$arCoef
		if( fit@mfit$mxn>0 ){
			ex = fit@mfit$external.regressors
			if( dim(ex)[1] < dim(data)[1]  ) stop("\ngogarchfilter-->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
		}
		p = fit@mfit$VARmodel$lag
		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" ){
				if(!exists("mclapply")){
					require('multicore')
				}
				armodel = mclapply(1:k, FUN = function(i){
							specx = arfimaspec(mean.model = list(armaOrder = c(p, 0), include.mean = TRUE, 
											arfima = FALSE, external.regressors = ex), distribution.model = "norm",
									fixed.pars = as.list(arCoef[,i]))
							ans = arfimafilter(spec = specx, data = data[,i], out.sample = 0)
							return(ans)}, 
						mc.cores = parallel.control$cores)
			} else{
				if(!exists("sfLapply")){
					require('snowfall')
				}
				sfInit(parallel=TRUE, cpus = parallel.control$cores)
				sfExport("data", "p", "ex", "arCoef", local = TRUE)
				armodel = sfLapply(1:k, fun = function(i) {
							specx = rgarch:::arfimaspec(mean.model = list(armaOrder = c(p, 0), include.mean = TRUE, 
											arfima = FALSE, external.regressors = ex), distribution.model = "norm",
									fixed.pars = as.list(arCoef[,i]))
							ans = rgarch::arfimafilter(spec = specx, data = data[,i], out.sample = 0)
							return(ans)})
				sfStop()
			}
		} else{
			armodel = lapply(1:k, FUN = function(i){
						specx = arfimaspec(mean.model = list(armaOrder = c(p, 0), include.mean = TRUE, 
										arfima = FALSE, external.regressors = ex), distribution.model = "norm",
								fixed.pars = as.list(arCoef[,i]))
						ans = arfimafit(spec = specx, data = data[,i], out.sample = 0)
						return(ans)})
		}
		cat("done!\n")
		zdata = sapply(armodel, FUN = function(x) residuals(x))
		if( p > 0 ) zdata = zdata[-c(1:p), , drop = FALSE]
		# should be N - p
		N = dim(data)[1]
		K = dim(data)[2]
		mu = sapply(armodel, FUN = function(x) fitted(x))
		if( p > 0 ) mu = mu[-c(1:p), , drop = FALSE]		
		# if we have out of sample data we need to filter
		if(out.sample > 0){
			if( parallel ){
				if( parallel.control$pkg == "multicore" ){
					if(!exists("mclapply")){
						require('multicore')
					}
					armodel2 = mclapply(1:k, FUN = function(i){
								specx = arfimaspec(mean.model = list(armaOrder = c(p, 0), include.mean = TRUE, 
												arfima = FALSE, external.regressors = orex), distribution.model = "norm",
										fixed.pars = as.list(arCoef[,i]))
								ans = arfimafilter(spec = specx, data = origdata[,i], out.sample = 0)
								return(ans)}, 
							mc.cores = parallel.control$cores)
				} else{
					if(!exists("sfLapply")){
						require('snowfall')
					}
					sfInit(parallel=TRUE, cpus = parallel.control$cores)
					sfExport("origdata", "p", "orex", "arCoef", local = TRUE)
					armodel2 = sfLapply(1:k, fun = function(i) {
								specx = rgarch:::arfimaspec(mean.model = list(armaOrder = c(p, 0), include.mean = TRUE, 
												arfima = FALSE, external.regressors = orex), distribution.model = "norm",
										fixed.pars = as.list(arCoef[,i]))
								ans = rgarch::arfimafilter(spec = specx, data = origdata[,i], out.sample = 0)
								return(ans)})
					sfStop()
				}
			} else{
				armodel2 = lapply(1:k, FUN = function(i){
							specx = arfimaspec(mean.model = list(armaOrder = c(p, 0), include.mean = TRUE, 
											arfima = FALSE, external.regressors = orex), distribution.model = "norm",
									fixed.pars = as.list(arCoef[,i]))
							ans = arfimafit(spec = specx, data = origdata[,i], out.sample = 0)
							return(ans)})
			}			
			zorigdata = sapply(armodel2, FUN = function(x) residuals(x))
			if( p > 0 ) zorigdata = zorigdata[-c(1:p), , drop = FALSE]
		} else{
			zorigdata = zdata
		}
	}
	A = fit@mfit$A
	W = fit@mfit$W
	#invU = solve(K)%*%W
	#U = solve(invU)
	# should we update the covariance matrix with new information?
	Y = zorigdata %*% W
	
	filterlist = multifilter(multifitORspec = fit@ufit, data = Y, out.sample = out.sample, n.old = NULL, 
			parallel = parallel, parallel.control = parallel.control)
	tmp = garchica2.ghyp(filterlist, A, K, N - p, mu)
	
	tmp$out.sample = out.sample
	tmp$origdata = origdata
	tmp$external.regressors = orex
	tmp$mxn = fit@mfit$mxn
	tmp$mdist = "magh"
	tmp$icamethod = fit@mfit$icamethod
	tmp$meanmodel = fit@mfit$meanmodel
	tmp$Y = Y
	tmp$W = W
	tmp$A = A
	tmp$K = fit@mfit$K
	tmp$VARmodel = vrmodel
	tmp$ARmodel = fit@mfit$ARmodel
	tmp$residuals = zorigdata
	
	ans = new("goGARCHfilter",
			mfilter = tmp,
			ufilter = filterlist)
	return(ans)
}

.gogarchforecast.ghyp = function(fit, n.ahead, n.roll, external.forecasts, parallel = FALSE, 
		parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2), ...)
{
	ns = fit@mfit$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")
	
	m = dim(fit@mfit$origdata)[2]
	
	# 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
	}
	vrmodel = fit@mfit$VARmodel
	if( !is.null(vrmodel) ){
		mxn = vrmodel$mxn
		if( mxn > 0 ){
			if( is.null( 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("\ngogarchforecast-->error: wrong number of external regressors!...", call. = FALSE)
				if( dim(mregfor)[1] < (n.roll + n.ahead) ) stop("\ngogarchforecast-->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("\ngogarchforecast-->error: n.ahead must be equal to 1 when using n.roll\n")
		
		if( n.ahead == 1 && (n.roll > ns) ) stop("\ngogarchforecast-->error: n.roll greater than out.sample!", call. = FALSE)
		#VARf = .varpredict(fit, n.ahead, n.roll, mregfor)$Mu
		Mu = varxforecast(X = fit@mfit$origdata, Bcoef = vrmodel$Bcoef, p = vrmodel$lag, out.sample = ns, 
				n.ahead, n.roll, mregfor)
	} else{
		# arfima model forecast
		arCoef = fit@mfit$ARmodel$arCoef
		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" ){
				if(!exists("mclapply")){
					require('multicore')
				}
				arfx = mclapply(1:m, FUN = function(i){
							specx = arfimaspec(mean.model = list(armaOrder = c(fit@mfit$ARmodel$lag, 0), include.mean = TRUE, 
											arfima = FALSE, external.regressors = fit@mfit$external.regressors), distribution.model = "norm",
									fixed.pars = as.list(arCoef[,i]))
							ans = arfimaforecast(fitORspec = specx, data = fit@mfit$origdata[,i], out.sample = ns,
									n.ahead = n.ahead, n.roll = n.roll, external.forecasts = list(mregfor = mregfor))
							return(ans)}, 
						mc.cores = parallel.control$cores)
			} else{
				if(!exists("sfLapply")){
					require('snowfall')
				}
				orgd  = fit@mfit$origdata
				p = fit@mfit$ARmodel$lag
				orex = fit@mfit$external.regressors
				sfInit(parallel=TRUE, cpus = parallel.control$cores)
				sfExport("orgd", "p", "orex", "arCoef", "mregfor", "ns", "n.ahead", "n.roll", local = TRUE)
				arfx = sfLapply(1:m, fun = function(i){
							specx = rgarch::arfimaspec(mean.model = list(armaOrder = c(p, 0), include.mean = TRUE, 
											arfima = FALSE, external.regressors = orex), distribution.model = "norm",
									fixed.pars = as.list(arCoef[,i]))
							ans = rgarch::arfimaforecast(fitORspec = specx, data = orgd[,i], out.sample = ns,
									n.ahead = n.ahead, n.roll = n.roll, external.forecasts = list(mregfor = mregfor))
							return(ans)})
				sfStop()
			}
		} else{
			arfx = lapply(1:m, FUN = function(i){
						specx = arfimaspec(mean.model = list(armaOrder = c(fit@mfit$ARmodel$lag, 0), include.mean = TRUE, 
										arfima = FALSE, external.regressors = fit@mfit$external.regressors), distribution.model = "norm",
								fixed.pars = as.list(arCoef[,i]))
						ans = arfimaforecast(fitORspec = specx, data = fit@mfit$origdata[,i], out.sample = ns,
								n.ahead = n.ahead, n.roll = n.roll, external.forecasts = list(mregfor = mregfor))
						return(ans)})
		}
		Mu = matrix(NA, ncol = m, nrow = tf)
		for(i in 1:m) Mu[,i] = sapply(arfx[[i]]@forecast$forecasts, FUN = function(x) x[,1])
	}
	
	fitlist = fit@ufit
	# we augmented n.roll before, now subtract
	forclist = multiforecast(multifitORspec = fitlist, n.ahead = n.ahead, n.roll = n.roll, parallel = parallel, 
			parallel.control = parallel.control, ...)
	A = fit@mfit$A
	tmp = garchica3.ghyp(forclist, A, m, n.ahead, n.roll+1, Mu)
	tmp$mdist = "magh"
	tmp$A = A
	tmp$K = fit@mfit$K
	tmp$W = fit@mfit$W
	ans = new("goGARCHforecast",
			mforecast = tmp,
			uforecast = forclist)
	return(ans)
}

.gogarchsim.ghyp = function(fit, n.sim, n.start, m.sim, startMethod, prereturns, mexsimdata, rseed,
		parallel = FALSE, parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2), ...)
{
	Data = head(fit@mfit$origdata, dim(fit@mfit$origdata)[1] - fit@mfit$out.sample)
	K = dim(Data)[2]
	np = dim(Data)[1]
	A = fit@mfit$A
	simlist = vector(mode = "list", length = K)
	m = fit@ufit@fit[[1]]@model$maxOrder
	if(is.null(startMethod[1])) startMethod = "unconditional" else startMethod = startMethod[1]
	
	if( !is.null(rseed) ){
		if( !is.matrix(rseed) ){
			xseed = matrix(NA, ncol = K, nrow = m.sim)
			if( length(rseed) != m.sim ) stop("\ngogarchsim--error: rseed must have m.sim rows!\n")
			xseed[,1] = rseed
			for(i in 2:K) xseed[, i] = as.integer(runif(1*m.sim,0,65000))
		} else{
			if(dim(rseed)[2] != K) stop("\ngogarchsim--error: rseed must have N (= n.assets) columns!\n")
			if(dim(rseed)[1] != m.sim) stop("\ngogarchsim--error: rseed must have m.sim rows!\n")
		}
	} else{
		xseed = matrix( as.integer( runif( 1*K, 0, as.integer(Sys.time()) ) ), ncol = K )
	}
	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" ){
			if(!exists("mclapply")){
				require('multicore')
			}
			simlist = mclapply(1:K, FUN = function(i) ugarchsim(fit@ufit@fit[[i]], n.sim = n.sim + n.start, 
								n.start = 0, m.sim = m.sim, startMethod = startMethod, rseed = xseed[,i]), 
					mc.cores = parallel.control$cores)
		} else{
			if(!exists("sfLapply")){
				require('snowfall')
			}
			sfInit(parallel=TRUE, cpus = parallel.control$cores)
			sfExport("fit", "n.sim", "n.start", "m.sim", "startMethod", "xseed",local = TRUE)
			simlist = sfLapply(as.list(1:K), fun = function(i) rgarch::ugarchsim(fit@ufit@fit[[i]], 
								n.sim = n.sim + n.start, n.start = 0, m.sim = m.sim, startMethod = startMethod, 
								rseed = xseed[,i]))
			if(fit@mfit$meanmodel != "VAR") sfStop()
		}
	} else{
		for(i in 1:K){
			# increment by 1 the rseed for each individual asset so that they do not get the same
			# random numbers!!!
			simlist[[i]] = ugarchsim(fit@ufit@fit[[i]], n.sim = n.sim + n.start, n.start = 0, 
					m.sim = m.sim, startMethod = startMethod, rseed = xseed[,i])
		}
	}
	
	sigmalist = vector(mode = "list", length = m.sim)
	retlist = vector(mode = "list", length = m.sim)
	residlistx = residlist = vector(mode = "list", length = m.sim)
	zlist = vector(mode = "list", length = m.sim)
	xlist = vector(mode = "list", length = m.sim)
	# we create the dlist object for use later if convolution is called.
	dlist = vector(mode = "list", length = m.sim)
	rho =  matrix( sapply(fit@ufit@fit, FUN = function(x) rep(x@fit$coef["skew"],  n.sim) ), nrow = n.sim )
	zeta = matrix( sapply(fit@ufit@fit, FUN = function(x) rep(x@fit$coef["shape"], n.sim) ), nrow = n.sim )
	dlambda = matrix( sapply( fit@ufit@fit, FUN = function(x) rep(x@fit$coef["dlambda"], n.sim) ), nrow = n.sim )
	
	if(!is.matrix(rho)) rho = matrix(rho, ncol = K)
	if(!is.matrix(zeta)) zeta = matrix(zeta, ncol = K)
	if(!is.matrix(dlambda)) dlambda = matrix(dlambda, ncol = K)
	
	zghypparams = vector(mode = "list", length = K)
	yghypparams = vector(mode = "list", length = K)
	# move from the location-scale parametrization of the z innovations ~(0,1,rho,chi) to the
	# alpha-beta parametrization for Y (the factors = sigma*z) ~ (0, sigma, skew, shape)==(mu, delta, alpha, beta)
	
	for(i in 1:m.sim){
		sigmalist[[i]] = matrix(sapply(simlist, FUN = function(x) tail(x@simulation$sigmaSim[,i], n.sim + n.start)), ncol = K)
		residlist[[i]] = matrix(sapply(simlist, FUN = function(x) tail(x@simulation$residSim[,i], n.sim + n.start)), ncol = K)
		zlist[[i]] = matrix(residlist[[i]]/sigmalist[[i]], ncol = K)
		residlistx[[i]] = matrix(t(apply(residlist[[i]], 1, FUN = function(x) t(A)%*%x)), ncol = K)
		residlist[[i]] = tail(residlistx[[i]], n.sim)
		sigmalist[[i]] = tail(sigmalist[[i]], n.sim)
		zlist[[i]] = tail(zlist[[i]], n.sim)
	}
		
	
	zghypparams = lapply(1:K, FUN = function(i) .ghyptransform(zeta = zeta[,i], rho = rho[,i], lambda = dlambda[,i]))	
	yghyplist = vector(mode = "list", length = m.sim)
	z.D = array(data = NA, dim = c(n.sim, K, 5))
	for(i in 1:K){
		z.D[1:n.sim, i, 1] = zghypparams[[i]][, 1]
		z.D[1:n.sim, i, 2] = zghypparams[[i]][, 2]
		z.D[1:n.sim, i, 3] = zghypparams[[i]][, 3]
		z.D[1:n.sim, i, 4] = zghypparams[[i]][, 4]
		z.D[1:n.sim, i, 5] = zghypparams[[i]][, 5]
	}
	
	for(j in 1:m.sim){		
		for(i in 1:K){
			yghypparams[[i]] = .ghypscale2(mu = 0, sigma = as.numeric(sigmalist[[j]][,i]), 
					ghalpha = zghypparams[[i]][,1], ghbeta = zghypparams[[i]][,2], 
					ghdelta = zghypparams[[i]][,3], ghmu = zghypparams[[i]][,4])
			yghypparams[[i]] = cbind(yghypparams[[i]], rep(dlambda[1,i], n.sim))
		}
		yghyplist[[j]] = yghypparams
	}
	
	
	if( fit@mfit$meanmodel == "VAR" ){
		# mexsimdata check
		mxn = fit@mfit$mxn
		if(mxn > 0 ){
			if( !is.null(mexsimdata) ){
				if( !is.list(mexsimdata) | length(mexsimdata) != m.sim ) stop("\ngogarchsim-->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("\ngogarchsim-->error: mexsimdata list submatrix must be a matrix...", call. = FALSE)
					if( dim(mexsimdata[[j]])[2] != mxn ) stop("\ngogarchsim-->error: wrong number of external regressors in list submatrix!...", call. = FALSE)
					if( dim(mexsimdata[[j]])[1] < (n.sim + n.start) ) stop("\ngogarchsim-->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[[j]] = matrix(0, ncol = mxn, nrow = (n.sim + n.start))
			}
		} else{
			mexsimdata = NULL
		}
		
		vrmodel = fit@mfit$VARmodel
		p = as.numeric(vrmodel$lag)
		if(!is.na(prereturns)[1]){
			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] != K) stop("\nwrong number of cols for prereturns matrix")
		} else{
			prereturns = as.matrix( tail(Data, p) )
		}
		if( n.sim == 1 && n.start == 0 ){
			# SPECIAL ROUTINE TO SPEED UP 1-AHEAD ESTIMATION
			tmp = varxsimXX(X = Data, vrmodel$Bcoef, p = vrmodel$lag, m.sim = m.sim, prereturns, resids = t(sapply(residlistx, 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]
		} else{
			if( parallel ){
				if( parallel.control$pkg == "multicore" ){
					xlist = mclapply(as.list(1:m.sim), FUN = function(i) varxsim(X = Data, vrmodel$Bcoef, p = vrmodel$lag, 
										n.sim = n.sim, n.start = n.start, prereturns = prereturns, resids = residlistx[[i]],
										mexsimdata = mexsimdata[[i]]), mc.cores = parallel.control$cores)
				} else{
					# when the GARCH model is without mean equation, simX == simulated residuals
					sfExport("Data", "vrmodel", "n.sim", "n.start", "prereturns", "residlistx", "mexsimdata", local = TRUE)
					xlist = sfLapply(as.list(1:m.sim), fun = function(i) rgarch:::varxsim(X = Data, vrmodel$Bcoef, p = vrmodel$lag, 
										n.sim = n.sim, n.start = n.start, prereturns = prereturns, resids = residlistx[[i]],
										mexsimdata = mexsimdata[[i]]))
					sfStop()
				}
			} else{
				xlist = lapply(as.list(1:m.sim), FUN = function(i) varxsim(X = Data, vrmodel$Bcoef, p = vrmodel$lag, 
									n.sim = n.sim, n.start = n.start, prereturns = prereturns, resids = residlistx[[i]], 
									mexsimdata = mexsimdata[[i]]))
			}
		}
	} else{
		N = dim(fit@mfit$origdata)[1]
		out.sample = fit@mfit$out.sample
		mxn = fit@mfit$mxn
		if(mxn > 0 ){
			if( !is.null(mexsimdata) ){
				if( !is.list(mexsimdata) | length(mexsimdata) != m.sim ) stop("\ngogarchsim-->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("\ngogarchsim-->error: mexsimdata list submatrix must be a matrix...", call. = FALSE)
					if( dim(mexsimdata[[j]])[2] != mxn ) stop("\ngogarchsim-->error: wrong number of external regressors in list submatrix!...", call. = FALSE)
					if( dim(mexsimdata[[j]])[1] < (n.sim + n.start) ) stop("\ngogarchsim-->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[[j]] = matrix(0, ncol = mxn, nrow = (n.sim + n.start))
			}
		} else{
			mexsimdata = NULL
		}
		
		p = fit@mfit$ARmodel$lag
		if(!is.na(prereturns)[1]){
			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] != K) stop("\nwrong number of cols for prereturns matrix")
		} else{
			prereturns = as.matrix( tail(Data, p) )
		}
		
		arCoef = fit@mfit$ARmodel$arCoef
		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" ){
				if(!exists("mclapply")){
					require('multicore')
				}
				arxsim = mclapply(1:K, FUN = function(i){
							specx = arfimaspec(mean.model = list(armaOrder = c(p, 0), include.mean = TRUE, 
											arfima = FALSE, external.regressors = fit@mfit$external.regressors[1:(N - out.sample), ,drop = FALSE]), 
									distribution.model = "norm",fixed.pars = as.list(arCoef[,i]))
							ans = arfimapath(spec = specx, n.sim = n.sim, n.start = n.start, m.sim = m.sim, prereturns = prereturns[,i],
									preresiduals = tail(fit@mfit$residuals[1:(np-p),i], p), custom.dist = list(name = "sample", 
											distfit = matrix(sapply(residlistx, FUN = function(x) x[,i]), ncol = m.sim, nrow = n.sim,
													byrow = TRUE), type = "res"), mexsimdata = mexsimdata)
							return(ans)}, mc.cores = parallel.control$cores)
			} else{
				if(!exists("sfLapply")){
					require('snowfall')
				}
				if(!is.null(fit@mfit$external.regressors)){
					ex = fit@mfit$external.regressors[1:(N - out.sample), , drop = FALSE]
				} else{
					ex = fit@mfit$external.regressors
				}
				sfInit(parallel=TRUE, cpus = parallel.control$cores)
				tres = tail(fit@mfit$residuals[1:(np-p), ], p)
				sfExport("p", "ex", "arCoef", "mexsimdata", "residlistx", "n.start", "tres", 
						"n.sim", "m.sim", "prereturns", local = TRUE)
				arxsim = sfLapply(1:K, fun = function(i){
							specx = rgarch::arfimaspec(mean.model = list(armaOrder = c(p, 0), include.mean = TRUE, 
											arfima = FALSE, external.regressors = ex), distribution.model = "norm",
									fixed.pars = as.list(arCoef[,i]))
							ans = rgarch::arfimapath(spec = specx, n.sim = n.sim, n.start = n.start, m.sim = m.sim, prereturns = prereturns[,i],
									preresiduals = tres[,i], custom.dist = list(name = "sample", 
											distfit = matrix(sapply(residlistx, FUN = function(x) x[,i]), ncol = m.sim, nrow = n.sim,
													byrow = TRUE), type = "res"), mexsimdata = mexsimdata)
							return(ans)})
				sfStop()
			}
		} else{
			arxsim = lapply(1:K, FUN = function(i){
						specx = arfimaspec(mean.model = list(armaOrder = c(p, 0), include.mean = TRUE, 
										arfima = FALSE, external.regressors = fit@mfit$external.regressors[1:(N - out.sample), ,drop = FALSE]), 
								distribution.model = "norm",fixed.pars = as.list(arCoef[,i]))
						ans = arfimapath(spec = specx, n.sim = n.sim, n.start = n.start, m.sim = m.sim, prereturns = prereturns[,i],
								preresiduals = tail(fit@mfit$residuals[1:(np-p),i], p), custom.dist = list(name = "sample", 
										distfit = matrix(sapply(residlistx, FUN = function(x) x[,i]), ncol = m.sim, nrow = n.sim,
												byrow = TRUE), type = "res"), mexsimdata = mexsimdata)
						return(ans)})
		}
		for(i in 1:m.sim) xlist[[i]] = matrix(sapply(arxsim, FUN = function(x) x@path$seriesSim[,i]), ncol = K, nrow = n.sim, byrow = TRUE)
	}
	
	
	tmp = list()
	tmp$distribution = fit@mfit$distr$model
	
	tmp$model = "ghyp"
	tmp$mdist = "magh"
	tmp$z.alphabeta = zghypparams
	tmp$y.alphabeta = yghyplist
	tmp$z.D = z.D
	tmp$A = A
	tmp$K = fit@mfit$K
	tmp$W = fit@mfit$W
	tmp$sigmaSim = sigmalist
	tmp$residSim = residlist
	tmp$seriesSim = xlist
	tmp$zSim = zlist
	tmp$rseed = xseed
	

	ans = new("goGARCHsim",
			msim = tmp,
			usim = simlist)
	return(ans)
	
}
#-----------------------------------------------------------------------------------------------
# secondary functions

#-----------------------------------------------------------------------------------------------
# manig
garchica1.nig = function(fit, A, m, n, mu)
{
	sigmas 	= matrix(NA, ncol = m, nrow = n)
	rho		= matrix(NA, ncol = m, nrow = n)
	zeta	= matrix(NA, ncol = m, nrow = n)
	resids 	= matrix(NA, ncol = m, nrow = n)
	zres 	= matrix(NA, ncol = m, nrow = n)
	
	sigmas	= sapply(fit@fit, FUN = function(x) x@fit$sigma, simplify = TRUE)	
	resids	= sapply(fit@fit, FUN = function(x) x@fit$residuals, simplify = TRUE)
	zeta	= sapply(fit@fit, FUN = function(x) as.matrix(rep(coef(x)["shape"], n)), simplify = TRUE)
	rho		= sapply(fit@fit, FUN = function(x) as.matrix(rep(coef(x)["skew"],  n)), simplify = TRUE)
	zres	= sapply(fit@fit, FUN = function(x) x@fit$z, simplify = TRUE)
	
	znigparams = vector(mode = "list", length = m)
	ynigparams = vector(mode = "list", length = m)
	
	
	# move from the location-scale parametrization of the z innovations ~(0,1,rho,chi) to the
	# alpha-beta parametrization for Y (the factors = sigma*z) ~ (0, sigma, skew, shape)==(mu, delta, alpha, beta)
	#
	# we should really return a scaled and non scaled version of nigparams since we will scale
	# later...
	
	for(i in 1:m){
		tmp = .nigtransform(zeta = zeta[,i], rho = rho[,i])
		# znig = the standardized innovations (0,1)
		# ynig = the rescaled residual distribution (0, sigma)
		znigparams[[i]] = tmp
		ynigparams[[i]] = .nigscale2(mu = 0, sigma = as.numeric(sigmas[,i]), nigalpha = tmp[,1], nigbeta = tmp[,2],
				nigdelta = tmp[,3], nigmu = tmp[,4])
		colnames(znigparams[[i]]) = colnames(ynigparams[[i]]) = c("alpha", "beta", "delta", "mu")
	}
	
	distr = list()
	distr$model = "nig"
	distr$z.alphabeta  	= znigparams
	distr$y.alphabeta 	= ynigparams
	
	z.D = array(data = NA, dim = c(n, m, 4))
	# z.D...-> page1: alpha, page2: beta, page3: delta, page4: mu
	for(i in 1:m)
	{
		z.D[1:n, i, 1] = znigparams[[i]][, 1]
		z.D[1:n, i, 2] = znigparams[[i]][, 2]
		z.D[1:n, i, 3] = znigparams[[i]][, 3]
		z.D[1:n, i, 4] = znigparams[[i]][, 4]
	}
	
	distr$z.D = z.D
	
	#y.expskew = vector(mode = "list", length = m)
	#y.expskew = sapply(ynigparams, FUN = function(y) apply(y, 1, FUN = function(x) .nigskew(x[1], x[2], x[3], x[4])))
	
	#y.expexkurt = vector(mode = "list", length = m)
	#y.expexkurt = sapply(ynigparams, FUN = function(y) apply(y, 1, FUN = function(x) .nigexkurt(x[1], x[2], x[3], x[4])))
	
	x.M = mu
	tmp = .gogarchllh.independent(t(A), fit)
	ans = list(sigmas = sigmas, resids = resids, x.M = x.M, llh = tmp$llh, likelihoods = tmp$likelihoods, distr = distr)
	return(ans)
}

garchica2.nig = function(filter, A, m, n, mu)
{
	sigmas 	= matrix(NA, ncol = m, nrow = n)
	rho		= matrix(NA, ncol = m, nrow = n)
	zeta	= matrix(NA, ncol = m, nrow = n)
	resids 	= matrix(NA, ncol = m, nrow = n)
	zres 	= matrix(NA, ncol = m, nrow = n)
	
	sigmas	= sapply(filter@filter, FUN = function(x) x@filter$sigma, simplify = TRUE)
	resids	= sapply(filter@filter, FUN = function(x) x@filter$residuals, simplify = TRUE)
	zeta	= sapply(filter@filter, FUN = function(x) as.matrix(rep(coef(x)["shape"], n)), simplify = TRUE)
	rho		= sapply(filter@filter, FUN = function(x) as.matrix(rep(coef(x)["skew"],  n)), simplify = TRUE)
	zres	= sapply(filter@filter, FUN = function(x) x@filter$z, simplify = TRUE)
	
	if(!is.matrix(sigmas)) sigmas = matrix(sigmas, ncol = m)
	if(!is.matrix(rho)) rho = matrix(rho, ncol = m)
	if(!is.matrix(zeta)) zeta = matrix(zeta, ncol = m)
	if(!is.matrix(resids)) resids = matrix(resids, ncol = m)
	if(!is.matrix(zres)) zres = matrix(zres, ncol = m)
	
	znigparams = vector(mode = "list", length = m)
	ynigparams = vector(mode = "list", length = m)
	# move from the location-scale parametrization of the z innovations ~(0,1,rho,chi) to the
	# alpha-beta parametrization for Y (the factors = sigma*z) ~ (0, sigma, skew, shape)==(mu, delta, alpha, beta)
	
	for(i in 1:m){
		# transform from (\zeta, \rho) to (\alpha, \beta, \delta, \mu)
		tmp = .nigtransform(zeta = zeta[,i], rho = rho[,i])
		# znig = the standardized innovations distribution 	: (0,1)
		# ynig = the rescaled residuals distribution 		: (0, sigma)
		znigparams[[i]] = tmp
		ynigparams[[i]] = .nigscale2(mu = 0, sigma = as.numeric(sigmas[,i]), nigalpha = tmp[,1], nigbeta = tmp[,2],
				nigdelta = tmp[,3], nigmu = tmp[,4])
		colnames(znigparams[[i]]) = colnames(ynigparams[[i]]) = c("alpha", "beta", "delta", "mu")
	}
	
	distr = list()
	distr$model = "nig"
	distr$z.alphabeta  	= znigparams
	distr$y.alphabeta 	= ynigparams
	
	z.D = array(data = NA, dim = c(n, m, 4))
	# z.D...-> page1: alpha, page2: beta, page3: delta, page4: mu
	for(i in 1:m)
	{
		z.D[1:n, i, 1] = znigparams[[i]][, 1]
		z.D[1:n, i, 2] = znigparams[[i]][, 2]
		z.D[1:n, i, 3] = znigparams[[i]][, 3]
		z.D[1:n, i, 4] = znigparams[[i]][, 4]
	}
	
	distr$z.D = z.D
	
	x.M = mu

	ans = list(sigmas = sigmas, resids = resids, x.M = x.M, x.R = NA, llh = NA, likelihoods = NA, distr = distr)
	return(ans)
}

# for simplicity and uncertainty of the forward-n density we only allow n.roll with n.ahead = 1
garchica3.nig = function(forclist, A, m, n.ahead = 1, n.roll, mu)
{
	nr = n.ahead * n.roll
	sigmas 	= matrix(NA, nrow = n.ahead*n.roll, ncol = m)
	rho		= matrix(NA, nrow = n.ahead*n.roll, ncol = m)
	zeta	= matrix(NA, nrow = n.ahead*n.roll, ncol = m)
	
	# continue here
	sigmas = sapply(forclist@forecast, FUN = function(x) sapply(x@forecast$forecasts, FUN = function(y) y[,1], simplify = TRUE))
	rho =  sapply(forclist@forecast, FUN = function(x) rep(x@filter$pars["skew"], nr) )
	zeta = sapply(forclist@forecast, FUN = function(x) rep(x@filter$pars["shape"],nr) )
	if(!is.matrix(sigmas)) sigmas = matrix(sigmas, ncol = m)
	if(!is.matrix(rho)) rho = matrix(rho, ncol = m)
	if(!is.matrix(zeta)) zeta = matrix(zeta, ncol = m)
	
	znigparams = vector(mode = "list", length = m)
	ynigparams = vector(mode = "list", length = m)
	# move from the location-scale parametrization of the z innovations ~(0,1,rho,chi) to the
	# alpha-beta parametrization for Y (the factors = sigma*z) ~ (0, sigma, skew, shape)==(mu, delta, alpha, beta)
	
	for(i in 1:m){
		znigparams[[i]] = .nigtransform(zeta = zeta[,i], rho = rho[,i])
		ynigparams[[i]] = .nigscale2(mu = 0, sigma = as.numeric(sigmas[,i]), nigalpha = znigparams[[i]][,1], 
				nigbeta = znigparams[[i]][,2], nigdelta = znigparams[[i]][,3], nigmu = znigparams[[i]][,4])
		# rownames (forecast values)
	}
	
	distr = list()
	distr$model = "nig"
	distr$z.alphabeta  	= znigparams
	distr$y.alphabeta 	= ynigparams
	
	z.D = array(data = NA, dim = c(nr, m, 4))
	# z.D...-> page1: alpha, page2: beta, page3: delta, page4: mu
	for(i in 1:m)
	{
		z.D[1:nr, i, 1] = znigparams[[i]][, 1]
		z.D[1:nr, i, 2] = znigparams[[i]][, 2]
		z.D[1:nr, i, 3] = znigparams[[i]][, 3]
		z.D[1:nr, i, 4] = znigparams[[i]][, 4]
	}
	
	distr$z.D = z.D
	x.M = mu
	#y.expskew = vector(mode = "list", length = m)
	#y.expskew = sapply(ynigparams, FUN = function(y) apply(y, 1, FUN = function(x) .nigskew(x[1], x[2], x[3], x[4])))
	
	#y.expexkurt = vector(mode = "list", length = m)
	#y.expexkurt = sapply(ynigparams, FUN = function(y) apply(y, 1, FUN = function(x) .nigexkurt(x[1], x[2], x[3], x[4])))
	
	#tmp = .goacdllh.independent(t(A), fit)
	ans = list(sigmas = sigmas, resids = NA, x.M = x.M, x.R = NA, llh = NA, likelihoods = NA,  distr = distr)
	return(ans)
}


# magh
garchica1.ghyp = function(fit, A, m, n, mu)
{
	sigmas 	= matrix(NA, ncol = m, nrow = n)
	rho		= matrix(NA, ncol = m, nrow = n)
	dlambda	= matrix(NA, ncol = m, nrow = n)
	zeta	= matrix(NA, ncol = m, nrow = n)
	resids 	= matrix(NA, ncol = m, nrow = n)
	zres 	= matrix(NA, ncol = m, nrow = n)
	
	sigmas	= sapply(fit@fit, FUN = function(x) x@fit$sigma, simplify = TRUE)	
	resids	= sapply(fit@fit, FUN = function(x) x@fit$residuals, simplify = TRUE)
	zeta	= sapply(fit@fit, FUN = function(x) as.matrix(rep(coef(x)["shape"], n)), simplify = TRUE)
	rho		= sapply(fit@fit, FUN = function(x) as.matrix(rep(coef(x)["skew"],  n)), simplify = TRUE)
	dlambda	= sapply(fit@fit, FUN = function(x) as.matrix(rep(coef(x)["dlambda"],  n)), simplify = TRUE)
	zres	= sapply(fit@fit, FUN = function(x) x@fit$z, simplify = TRUE)
	
	zghypparams = vector(mode = "list", length = m)
	yghypparams = vector(mode = "list", length = m)
	
	
	# move from the location-scale parametrization of the z innovations ~(0,1,rho,chi) to the
	# alpha-beta parametrization for Y (the factors = sigma*z) ~ (0, sigma, skew, shape)==(mu, delta, alpha, beta)
	#
	# we should really return a scaled and non scaled version of nigparams since we will scale
	# later...
	
	for(i in 1:m){
		tmp = .ghyptransform(zeta = zeta[,i], rho = rho[,i], lambda = dlambda[1,i])
		# znig = the standardized innovations (0,1)
		# ynig = the rescaled residual distribution (0, sigma)
		zghypparams[[i]] = tmp
		yghypparams[[i]] = .ghypscale2(mu = 0, sigma = as.numeric(sigmas[,i]), ghalpha = tmp[,1], 
				ghbeta = tmp[,2], ghdelta = tmp[,3], ghmu = tmp[,4])
		yghypparams[[i]] = cbind(yghypparams[[i]], rep(dlambda[1,i], n))
		colnames(zghypparams[[i]]) = c("alpha", "beta", "delta", "mu", "lambda")
		colnames(yghypparams[[i]]) = c("alpha", "beta", "delta", "mu", "lambda")
	}
	
	distr = list()
	distr$model = "ghyp"
	distr$z.alphabeta  	= zghypparams
	distr$y.alphabeta 	= yghypparams
	
	z.D = array(data = NA, dim = c(n, m, 5))
	# z.D...-> page1: alpha, page2: beta, page3: delta, page4: mu
	for(i in 1:m)
	{
		z.D[1:n, i, 1] = zghypparams[[i]][, 1]
		z.D[1:n, i, 2] = zghypparams[[i]][, 2]
		z.D[1:n, i, 3] = zghypparams[[i]][, 3]
		z.D[1:n, i, 4] = zghypparams[[i]][, 4]
		z.D[1:n, i, 5] = zghypparams[[i]][, 5]
	}
	
	distr$z.D = z.D
	
	#y.expskew = vector(mode = "list", length = m)
	#y.expskew = sapply(ynigparams, FUN = function(y) apply(y, 1, FUN = function(x) .nigskew(x[1], x[2], x[3], x[4])))
	
	#y.expexkurt = vector(mode = "list", length = m)
	#y.expexkurt = sapply(ynigparams, FUN = function(y) apply(y, 1, FUN = function(x) .nigexkurt(x[1], x[2], x[3], x[4])))
	
	x.M = mu
	tmp = .gogarchllh.independent(t(A), fit)
	ans = list(sigmas = sigmas, resids = resids, x.M = x.M, llh = tmp$llh, likelihoods = tmp$likelihoods, distr = distr)
	return(ans)
}

garchica2.ghyp = function(filter, A, m, n, mu)
{
	sigmas 	= matrix(NA, ncol = m, nrow = n)
	rho		= matrix(NA, ncol = m, nrow = n)
	dlambda	= matrix(NA, ncol = m, nrow = n)
	zeta	= matrix(NA, ncol = m, nrow = n)
	resids 	= matrix(NA, ncol = m, nrow = n)
	zres 	= matrix(NA, ncol = m, nrow = n)
	
	sigmas	= sapply(filter@filter, FUN = function(x) x@filter$sigma, simplify = TRUE)
	resids	= sapply(filter@filter, FUN = function(x) x@filter$residuals, simplify = TRUE)
	zeta	= sapply(filter@filter, FUN = function(x) as.matrix(rep(coef(x)["shape"], n)), simplify = TRUE)
	rho		= sapply(filter@filter, FUN = function(x) as.matrix(rep(coef(x)["skew"],  n)), simplify = TRUE)
	dlambda	= sapply(filter@filter, FUN = function(x) as.matrix(rep(coef(x)["dlambda"],  n)), simplify = TRUE)
	zres	= sapply(filter@filter, FUN = function(x) x@filter$z, simplify = TRUE)
	
	if(!is.matrix(sigmas)) sigmas = matrix(sigmas, ncol = m)
	if(!is.matrix(rho)) rho = matrix(rho, ncol = m)
	if(!is.matrix(zeta)) zeta = matrix(zeta, ncol = m)
	if(!is.matrix(resids)) resids = matrix(resids, ncol = m)
	if(!is.matrix(zres)) zres = matrix(zres, ncol = m)
	if(!is.matrix(dlambda)) dlambda = matrix(dlambda, ncol = m)
	
	zghypparams = vector(mode = "list", length = m)
	yghypparams = vector(mode = "list", length = m)
	# move from the location-scale parametrization of the z innovations ~(0,1,rho,chi) to the
	# alpha-beta parametrization for Y (the factors = sigma*z) ~ (0, sigma, skew, shape)==(mu, delta, alpha, beta)
	
	for(i in 1:m){
		# transform from (\zeta, \rho) to (\alpha, \beta, \delta, \mu)
		tmp = .ghyptransform(zeta = zeta[,i], rho = rho[,i], lambda = dlambda[,i])
		# znig = the standardized innovations distribution 	: (0,1)
		# ynig = the rescaled residuals distribution 		: (0, sigma)
		zghypparams[[i]] = tmp
		yghypparams[[i]] = .ghypscale2(mu = 0, sigma = as.numeric(sigmas[,i]), ghalpha = tmp[,1], 
				ghbeta = tmp[,2], ghdelta = tmp[,3], ghmu = tmp[,4])
		yghypparams[[i]] = cbind(yghypparams[[i]], rep(dlambda[1,i], n))
		colnames(zghypparams[[i]]) = c("alpha", "beta", "delta", "mu", "lambda")
		colnames(yghypparams[[i]]) = c("alpha", "beta", "delta", "mu", "lambda")
	}
	
	distr = list()
	distr$model = "nig"
	distr$z.alphabeta  	= zghypparams
	distr$y.alphabeta 	= yghypparams
	
	z.D = array(data = NA, dim = c(n, m, 5))
	# z.D...-> page1: alpha, page2: beta, page3: delta, page4: mu
	for(i in 1:m)
	{
		z.D[1:n, i, 1] = zghypparams[[i]][, 1]
		z.D[1:n, i, 2] = zghypparams[[i]][, 2]
		z.D[1:n, i, 3] = zghypparams[[i]][, 3]
		z.D[1:n, i, 4] = zghypparams[[i]][, 4]
		z.D[1:n, i, 5] = zghypparams[[i]][, 5]
	}
	
	distr$z.D = z.D
	
	x.M = mu
	
	ans = list(sigmas = sigmas, resids = resids, x.M = x.M, x.R = NA, llh = NA, likelihoods = NA, distr = distr)
	return(ans)
}

# for simplicity and uncertainty of the forward-n density we only allow n.roll with n.ahead = 1
garchica3.ghyp = function(forclist, A, m, n.ahead = 1, n.roll, mu)
{
	nr = n.ahead * n.roll
	sigmas 	= matrix(NA, nrow = n.ahead*n.roll, ncol = m)
	rho	= matrix(NA, nrow = n.ahead*n.roll, ncol = m)
	zeta	= matrix(NA, nrow = n.ahead*n.roll, ncol = m)
	dlambda	= matrix(NA, nrow = n.ahead*n.roll, ncol = m)
	
	# continue here
	sigmas = sapply(forclist@forecast, FUN = function(x) sapply(x@forecast$forecasts, FUN = function(y) y[,1], simplify = TRUE))
	rho =  sapply(forclist@forecast, FUN = function(x) rep(x@filter$pars["skew"], nr) )
	zeta = sapply(forclist@forecast, FUN = function(x) rep(x@filter$pars["shape"],nr) )
	dlambda = sapply(forclist@forecast, FUN = function(x) rep(x@filter$pars["dlambda"],nr) )
	
	if(!is.matrix(sigmas)) sigmas = matrix(sigmas, ncol = m)
	if(!is.matrix(rho)) rho = matrix(rho, ncol = m)
	if(!is.matrix(zeta)) zeta = matrix(zeta, ncol = m)
	if(!is.matrix(dlambda)) dlambda = matrix(dlambda, ncol = m)
	
	zghypparams = vector(mode = "list", length = m)
	yghypparams = vector(mode = "list", length = m)
	# move from the location-scale parametrization of the z innovations ~(0,1,rho,chi) to the
	# alpha-beta parametrization for Y (the factors = sigma*z) ~ (0, sigma, skew, shape)==(mu, delta, alpha, beta)
	
	for(i in 1:m){
		zghypparams[[i]] = .ghyptransform(zeta = zeta[,i], rho = rho[,i], lambda = dlambda[1,i])
		yghypparams[[i]] = .ghypscale2(mu = 0, sigma = as.numeric(sigmas[,i]), ghalpha = zghypparams[[i]][,1], 
				ghbeta = zghypparams[[i]][,2], ghdelta = zghypparams[[i]][,3], ghmu = zghypparams[[i]][,4])
		yghypparams[[i]] = cbind(yghypparams[[i]], rep(dlambda[1,i], nr))
	}
	
	distr = list()
	distr$model = "nig"
	distr$z.alphabeta  	= zghypparams
	distr$y.alphabeta 	= yghypparams
	
	z.D = array(data = NA, dim = c(nr, m, 5))
	# z.D...-> page1: alpha, page2: beta, page3: delta, page4: mu
	for(i in 1:m)
	{
		z.D[1:nr, i, 1] = zghypparams[[i]][, 1]
		z.D[1:nr, i, 2] = zghypparams[[i]][, 2]
		z.D[1:nr, i, 3] = zghypparams[[i]][, 3]
		z.D[1:nr, i, 4] = zghypparams[[i]][, 4]
		z.D[1:nr, i, 5] = zghypparams[[i]][, 5]
	}
	
	distr$z.D = z.D
	x.M = mu
	#y.expskew = vector(mode = "list", length = m)
	#y.expskew = sapply(ynigparams, FUN = function(y) apply(y, 1, FUN = function(x) .nigskew(x[1], x[2], x[3], x[4])))
	
	#y.expexkurt = vector(mode = "list", length = m)
	#y.expexkurt = sapply(ynigparams, FUN = function(y) apply(y, 1, FUN = function(x) .nigexkurt(x[1], x[2], x[3], x[4])))
	
	#tmp = .goacdllh.independent(t(A), fit)
	ans = list(sigmas = sigmas, resids = NA, x.M = x.M, x.R = NA, llh = NA, likelihoods = NA,  distr = distr)
	return(ans)
}

#-----------------------------------------------------------------------------------------------
# mvnorm
garchica1.norm = function(fit, A, m, n, mu)
{
	sigmas 	= matrix(NA, ncol = m, nrow = n)
	resids 	= matrix(NA, ncol = m, nrow = n)
	zres 	= matrix(NA, ncol = m, nrow = n)
	
	sigmas	= sapply(fit@fit, FUN = function(x) x@fit$sigma, simplify = TRUE)	
	resids	= sapply(fit@fit, FUN = function(x) x@fit$residuals, simplify = TRUE)
	zres	= sapply(fit@fit, FUN = function(x) x@fit$z, simplify = TRUE)
	
	if(!is.matrix(sigmas)) sigmas = matrix(sigmas, ncol = m)
	if(!is.matrix(resids)) resids = matrix(resids, ncol = m)
	if(!is.matrix(zres)) zres = matrix(zres, ncol = m)
	
	znigparams = vector(mode = "list", length = m)
	ynigparams = vector(mode = "list", length = m)
	
	for(i in 1:m){
		znigparams[[i]] = sigmas[,i]
		ynigparams[[i]] = znigparams[[i]]
	}
	
	distr = list()
	distr$model = "norm"
	distr$z.alphabeta  	= znigparams
	distr$y.alphabeta 	= ynigparams
	
	z.D = array(data = NA, dim = c(n, m, 1))
	# z.D...-> page1: sigma, page2: mu
	for(i in 1:m)
	{
		z.D[1:n, i, 1] = znigparams[[i]]
	}
	
	distr$z.D = z.D
	
	x.M = mu
	tmp = .gogarchllh.independent(t(A), fit)
	ans = list(sigmas = sigmas, resids = resids, x.M = x.M, llh = tmp$llh, likelihoods = tmp$likelihoods, distr = distr)
	return(ans)
}

garchica2.norm = function(filter, A, m, n, mu)
{
	sigmas 	= matrix(NA, ncol = m, nrow = n)
	resids 	= matrix(NA, ncol = m, nrow = n)
	zres 	= matrix(NA, ncol = m, nrow = n)
	
	sigmas	= sapply(filter@filter, FUN = function(x) x@filter$sigma, simplify = TRUE)
	resids	= sapply(filter@filter, FUN = function(x) x@filter$residuals, simplify = TRUE)
	zres	= sapply(filter@filter, FUN = function(x) x@filter$z, simplify = TRUE)
	
	if(!is.matrix(sigmas)) sigmas = matrix(sigmas, ncol = m)
	if(!is.matrix(resids)) resids = matrix(resids, ncol = m)
	if(!is.matrix(zres)) zres = matrix(zres, ncol = m)
	
	znigparams = vector(mode = "list", length = m)
	ynigparams = vector(mode = "list", length = m)
	
	for(i in 1:m){
		znigparams[[i]] = cbind(sigmas[,i], rep(0, n))
		ynigparams[[i]] = znigparams[[i]]
		colnames(znigparams[[i]]) = colnames(ynigparams[[i]]) = c("sigma", "mu")
	}
	
	distr = list()
	distr$model = "nig"
	distr$z.alphabeta  	= znigparams
	distr$y.alphabeta 	= ynigparams
	
	z.D = array(data = NA, dim = c(n, m, 2))
	# z.D...-> page1: sigma, page2: mu
	for(i in 1:m)
	{
		z.D[1:n, i, 1] = znigparams[[i]][, 1]
		z.D[1:n, i, 2] = znigparams[[i]][, 2]
	}
	
	distr$z.D = z.D
	
	x.M = mu
	
	ans = list(sigmas = sigmas, resids = resids, x.M = x.M, x.R = NA, llh = NA, likelihoods = NA, distr = distr)
	return(ans)
}

# for simplicity and uncertainty of the forward-n density we only allow n.roll with n.ahead = 1
garchica3.norm = function(forclist, A, m, n.ahead = 1, n.roll, mu)
{
	nr = n.ahead * n.roll
	sigmas 	= matrix(NA, nrow = n.ahead*n.roll, ncol = m)
	
	# continue here
	sigmas = sapply(forclist@forecast, FUN = function(x) sapply(x@forecast$forecasts, FUN = function(y) y[,1], simplify = TRUE))
	if(!is.matrix(sigmas)) sigmas = matrix(sigmas, ncol = m)
	
	znigparams = ynigparams = vector(mode = "list", length = m)
	
	for(i in 1:m){
		znigparams[[i]] = cbind(sigmas[,i], rep(0, nr))
		ynigparams[[i]] = znigparams[[i]]
		colnames(znigparams[[i]]) = colnames(ynigparams[[i]]) = c("sigma", "mu")
	}
	
	distr = list()
	distr$model = "norm"
	distr$z.alphabeta  	= znigparams
	distr$y.alphabeta 	= ynigparams
	
	z.D = array(data = NA, dim = c(nr, m, 4))
	# z.D...-> page1: alpha, page2: beta, page3: delta, page4: mu
	for(i in 1:m)
	{
		z.D[1:nr, i, 1] = znigparams[[i]][, 1]
		z.D[1:nr, i, 2] = znigparams[[i]][, 2]
	}
	
	distr$z.D = z.D
	x.M = mu
	ans = list(sigmas = sigmas, resids = NA, x.M = x.M, x.R = NA, llh = NA, likelihoods = NA,  distr = distr)
	return(ans)
}


# The likelihood of the affine linear model with independent margins is simply the
# sum of the univariate/independent log-likelihoods plus a term for the mixing matrix A
.gogarchllh.independent = function(Z, fit)
{
	# from the density transformation theorem (see for example Schmidt Technical Report P.14)
	lik = sapply(fit@fit, FUN = function(x) -x@fit$log.likelihoods)
	n = dim(lik)[1]
	#U = Z %*% fit@mfit$K
	dt = log(abs(det(solve(Z))))
	likelihoods = apply(lik, 1, FUN = function(x) sum(x))
	llh = n*dt + sum(likelihoods)
	return(list(llh = llh, likelihoods = likelihoods))
}



# make changes to subdivisions and rel.tol
.safeintegrate = function(fun, lower, upper, ...)
{
	ans = try(integrate(fun, lower, upper, ...))
	if(inherits(ans, "try-error")){
		ans = integrate(fun, lower, upper, ...)
	}
	return(ans)
}

# quadrature based numerical moments from the convoluted density
.dcintmoments = function(z, y, sm1 = NULL, fix.sm1 = TRUE, subdivisions = 200, 
		rel.tol = .Machine$double.eps^0.5, n.approx = 100, error.reltol = 0.05, ...)
{
	# convoluted density arising from numerical inversion of characteristic function of multivariate affine distribution
	dmad = .makeDNew(x = z, dx = y, h = diff(z[1:10])[1], standM = "sum")
	if(!is.null(sm1) && fix.sm1){
		m1 = sm1
	} else{
		fm1 = function(x) x * dmad(x)
			m1 = integrate(fm1, -Inf, Inf, subdivisions = 150, rel.tol=.Machine$double.eps^0.25)$value
			if(!is.null(sm1)){
				if( ((m1-sm1)/m1)>error.reltol || ((m1-sm1)/m1)<(-error.reltol) ) {
					m1 = try(integrate(fm1, -Inf, Inf, subdivisions = 200, rel.tol=.Machine$double.eps^0.5)$value, silent = TRUE)
					if(inherits(m1, "try-error")) {
						m1 = sm1
					} else{
						if( ((m1-sm1)/m1)>error.reltol || ((m1-sm1)/m1)<(-error.reltol) ) m1 = sm1
					}
				}
			}
	}
	fm2 = function(x) ((x-m1)^2) * dmad(x)
	m2 = sqrt(integrate(fm2, -Inf, Inf, subdivisions = subdivisions, rel.tol = rel.tol, stop.on.error = FALSE)$value)
	fm3 =  function(x) ((x-m1)^3) * dmad(x)
	m3 = integrate(fm3, -Inf, Inf, subdivisions = subdivisions, rel.tol = rel.tol, stop.on.error = FALSE)$value
	fm4 =  function(x) ((x-m1)^4) * dmad(x)
	m4 = integrate(fm4, -Inf, Inf, subdivisions = subdivisions, rel.tol = rel.tol, stop.on.error = FALSE)$value
	return(c(mean = m1, sd = m2, m3 = m3, m4 = m4))
}

.dfun = function(z, y, fn, lower = -Inf, upper = Inf, subdivisions = 200, rel.tol = .Machine$double.eps^0.5)
{
	# convoluted density arising from numerical inversion of characteristic function of multivariate affine distribution
	dmad = .makeDNew(x = z, dx = y, h = diff(z[1:10])[1], standM = "sum")
	ff = eval(parse(text = paste("function(x) ", fn, "*dmad(x)", sep = "")))
	ans = integrate(ff, lower, upper, subdivisions = subdivisions, rel.tol = rel.tol)$value
	return(ans)
}


.simulate.constant = function(resid, Mu, n, k)
{
	ans = matrix(NA, ncol = k, nrow = n)
	for(i in 1:n){
		ans[i,] = Mu + resid[i,]
	}
	return(ans)
}