R/gogarch-methods.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.
##
#################################################################################


#------------------------------------------------------------------------------
# spec
#------------------------------------------------------------------------------

gogarchspec = function(mean.model = list(demean = c("constant", "AR", "VAR", "robVAR"), lag = 1, lag.max = NULL, 
				lag.criterion = c("AIC", "HQ", "SC", "FPE"), external.regressors = NULL, 
				robust.control = list("gamma" = 0.25, "delta" = 0.01, "nc" = 10, "ns" = 500)), 
		variance.model = list(model = "sGARCH", garchOrder = c(1,1), submodel = NULL, variance.targeting = FALSE), 
		distribution.model = list(distribution = c("mvnorm", "manig", "magh")),
		ica = c("fastica", "pearson", "jade", "radical"), ica.fix = list(A = NULL, K = NULL), ...)
{
	UseMethod("gogarchspec")
}

setMethod("gogarchspec",  definition = .gogarchspec)

#------------------------------------------------------------------------------
# fit
#------------------------------------------------------------------------------
gogarchfit = 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, ...)
{
	UseMethod("gogarchfit")
}

setMethod("gogarchfit",  signature(spec = "goGARCHspec"), definition = .gogarchfit)


#------------------------------------------------------------------------------
# filter
#------------------------------------------------------------------------------
gogarchfilter = function(fit, data, out.sample = 0, n.old = NULL, parallel = FALSE, 
		parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2), ...)
{
	UseMethod("gogarchfilter")
}

setMethod("gogarchfilter",  signature(fit = "goGARCHfit"), definition = .gogarchfilter)


#------------------------------------------------------------------------------
# forecast
#------------------------------------------------------------------------------
gogarchforecast = function(fit, n.ahead = 10, n.roll = 0, external.forecasts = list(mregfor = NULL), 
		parallel = FALSE, parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2), ...)
{
	UseMethod("gogarchforecast")
}

setMethod("gogarchforecast",  signature(fit = "goGARCHfit"), definition = .gogarchforecast)

#------------------------------------------------------------------------------
# simulate
#------------------------------------------------------------------------------
gogarchsim = function(fit, n.sim = 1, n.start = 0, m.sim = 1, startMethod = c("unconditional", 
				"sample"), prereturns = NA, mexsimdata = NULL, rseed = NULL, parallel = FALSE, 
		parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2), ...)
{
	UseMethod("gogarchsim")
}

setMethod("gogarchsim",  signature(fit = "goGARCHfit"), definition = .gogarchsim)

#------------------------------------------------------------------------------
gogarchroll = function(spec, data, n.ahead = 1, forecast.length = 50, 
		refit.every = 25, refit.window = c("recursive", "moving"), parallel = FALSE, 
		parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2), 
		solver = "solnp", external.forecasts = list(mregfor = NULL), fit.control = list(), 
		solver.control = list(), rseed = NULL, save.fit = FALSE, save.wdir = NULL, ...)
{
	UseMethod("gogarchroll")
}

setMethod("gogarchroll",  signature(spec = "goGARCHspec"), definition = .rollgogarch)


#------------------------------------------------------------------------------
# fit show
#------------------------------------------------------------------------------
setMethod("show",
		signature(object = "goGARCHspec"),
		function(object){
			model = strsplit(class(object), "spec")
			cat(paste("\n*------------------------------*", sep = ""))
			cat(paste("\n*       GO-GARCH Spec          *", sep = ""))
			cat(paste("\n*------------------------------*", sep = ""))
			cat(paste("\n\nMean Model\t\t: ", object@mspec$mmodel$model, sep = ""))
			cat(paste("\nGARCH Model\t\t: ", object@mspec$vmodel$model, sep = ""))
			cat(paste("\nDistribution\t: ", object@mspec$distribution, sep = ""))
			cat(paste("\nICA Method\t\t: ", object@mspec$ica, sep = ""))
			cat("\n")
			invisible(object)
		})
		
setMethod("show",
		signature(object = "goGARCHfit"),
		function(object){
			model = strsplit(class(object),"fit")
			cat(paste("\n*------------------------------*", sep = ""))
			cat(paste("\n*        GO-GARCH Fit          *", sep = ""))
			cat(paste("\n*------------------------------*", sep = ""))
			cat(paste("\n\nMean Model\t: ",object@mfit$meanmodel,sep=""))
			cat(paste("\nGARCH Model\t: ",strsplit(class(object@ufit@fit[[1]]),"fit"),sep=""))
			cat(paste("\nDistribution\t: ",object@mfit$mdist,sep=""))
			cat(paste("\nICA Method\t: ",object@mfit$icamethod,sep=""))
			cat(paste("\nNo. Factors\t: ",dim(object@mfit$Y)[2], sep=""))
			cat(paste("\nLog-Likelihood : ",round(object@mfit$llh,2),sep=""))
			cat(paste("\n------------------------------------\n",sep=""))
			cat(paste("\nU (rotation matrix) : \n\n",sep=""))
			print(as.data.frame(object, type = "U"), digits = 3)
			cat(paste("\nA (mixing matrix) : \n\n",sep=""))
			print(as.data.frame(object, type = "A"), digits = 3)
			cat("\n")
			invisible(object)
		})


setMethod("show",
		signature(object = "goGARCHfilter"),
		function(object){
			model = strsplit(class(object),"filter")
			cat(paste("\n*------------------------------*", sep = ""))
			cat(paste("\n*        GO-GARCH Filter       *", sep = ""))
			cat(paste("\n*------------------------------*", sep = ""))
			cat(paste("\n\nMean Model\t: ",object@mfilter$meanmodel,sep=""))
			cat(paste("\nGARCH Model\t: ",strsplit(class(object@ufilter@filter[[1]]),"filter"),sep=""))
			cat(paste("\nDistribution\t: ",object@mfilter$mdist,sep=""))
			cat(paste("\nICA Method\t: ",object@mfilter$icamethod,sep=""))
			cat(paste("\nNo. Factors\t: ",dim(object@mfilter$Y)[2], sep=""))
			#cat(paste("\nLog-Quasi Likelihood : ",round(object@mfit$llh,2),sep=""))
			cat(paste("\n------------------------------------\n",sep=""))
			cat(paste("\nU (rotation matrix) : \n\n",sep=""))
			print(as.data.frame(object, type = "U"), digits = 3)
			cat(paste("\nA (mixing matrix) : \n\n",sep=""))
			print(as.data.frame(object, type = "A"), digits = 3)
			cat("\n")
			invisible(object)
		})

setMethod("show",
		signature(object = "goGARCHforecast"),
		function(object){
			model = strsplit(class(object),"forecast")
			cat(paste("\n*------------------------------*", sep = ""))
			cat(paste("\n*        GO-GARCH Forecast     *", sep = ""))
			cat(paste("\n*------------------------------*", sep = ""))
			cat(paste("\n\nMean Model\t: ",object@mforecast$meanmodel,sep=""))
			cat(paste("\nGARCH Model\t: ",strsplit(class(object@uforecast@forecast[[1]]),"forecast"),sep=""))
			cat(paste("\nDistribution\t: ",object@mforecast$mdist,sep=""))
			cat(paste("\nICA Method\t: ",object@mforecast$icamethod,sep=""))
			cat(paste("\nNo. Factors\t: ",dim(object@mforecast$Y)[2], sep=""))
			#cat(paste("\nLog-Quasi Likelihood : ",round(object@mfit$llh,2),sep=""))
			cat("\n\n")
			invisible(object)
		})
#------------------------------------------------------------------------------
# likelihood
#------------------------------------------------------------------------------
.gogarchllh = function(object)
{
	return(object@mfit$llh)
}

setMethod("likelihood", signature(object = "goGARCHfit"), .gogarchllh)

.gogarchllhfilter = function(object)
{
	n = dim(object@mfilter$sigma)[1]
	lik1 = sum(sapply(object@ufilter@filter,FUN = function(x) likelihood(x)))
	sumA = n*log(abs(det(solve(object@mfilter$A))))
	lik =  sumA + lik1
	return(lik)
}

setMethod("likelihood", signature(object = "goGARCHfilter"), .gogarchllhfilter)



#------------------------------------------------------------------------------
# coef, fitted, dmu, dscale, dshape, dskew, mu, sigma, skew, kurt, 
#------------------------------------------------------------------------------
.gogarchfitcoef = function(object)
{
	return(t(sapply(object@ufit@fit, FUN = function(x) coef(x))))
}

setMethod("coef", signature(object = "goGARCHfit"), .gogarchfitcoef)


.gogarchfltcoef = function(object)
{
	return(t(sapply(object@ufilter@filter, FUN = function(x) coef(x))))
}

setMethod("coef", signature(object = "goGARCHfilter"), .gogarchfltcoef)


.gogarchfitfitted = function(object)
{
	return(object@mfit$x.M)
}

setMethod("fitted", signature(object = "goGARCHfit"), .gogarchfitfitted)

.gogarchfltfitted = function(object)
{
	return(object@mfilter$x.M)
}

setMethod("fitted", signature(object = "goGARCHfilter"), .gogarchfltfitted)


.gogarchfitresids = function(object)
{
	return(object@mfit$residuals)
}

setMethod("residuals", signature(object = "goGARCHfit"), .gogarchfitresids)

.gogarchfltresids = function(object)
{
	return(object@mfilter$residuals)
}

setMethod("residuals", signature(object = "goGARCHfilter"), .gogarchfltresids)
#------------------------------------------------------------------------------

#------------------------------------------------------------------------------
as.dfgofit = function(x, row.names = NULL, optional = FALSE, type = "A", ...)
{
	validmat = c("A", "S", "U", "K", "mean", "residuals", "sigma", "z")
	if(!any(tolower(type) == tolower(validmat))) stop("Invalid type for gogarch matrix chosen.", call. = FALSE)
	ans = switch(type,
			A = .gogarchfitA(x),
			S = .gogarchfitS(x),
			U = .gogarchfitU(x),
			K = .gogarchfitK(x))
	return( as.data.frame( ans ) )
}

setMethod("as.data.frame", signature(x = "goGARCHfit"), as.dfgofit)

as.dfgofilter = function(x, row.names = NULL, optional = FALSE, type = "A", ...)
{
	validmat = c("A", "S", "U", "K")
	if(!any(tolower(type) == tolower(validmat))) stop("Invalid type for gogarch matrix chosen. Valid choices as A, S, U and K.", call. = FALSE)
	ans = switch(type,
			A = .gogarchfilterA(x),
			S = .gogarchfilterS(x),
			U = .gogarchfilterU(x),
			K = .gogarchfilterK(x))
	return( as.data.frame( ans ) )
}

setMethod("as.data.frame", signature(x = "goGARCHfilter"), as.dfgofilter)

# mixing matrix A
.gogarchfitA = function(object)
{
	return( object@mfit$A )
}

.gogarchfilterA = function(object)
{
	return( object@mfilter$A )
}

# unconditional covariance
.gogarchfitS = function(object)
{
	return( t(object@mfit$A)%*%object@mfit$A )
}

.gogarchfilterS = function(object)
{
	return( t(object@mfilter$A)%*%object@mfilter$A )
}


# rotational matrix
.gogarchfitU = function(object)
{
	A = object@mfit$A
	K = (object@mfit$K)
	if(is.null(K)) stop("ICA algorithm does not return required information to calculate U", call. = FALSE)
	return( A%*%K )
}

.gogarchfilterU = function(object)
{
	A = object@mfilter$A
	K = (object@mfilter$K)
	if(is.null(K)) stop("ICA algorithm does not return required information to calculate U", call. = FALSE)
	return( A%*%K )
}

# whitening matrix
.gogarchfitK = function(object)
{
	K = (object@mfit$K)
	if(is.null(K)) stop("ICA algorithm does not return the whitening matrix", call. = FALSE)
	return( K )
}

.gogarchfilterK = function(object)
{
	K = (object@mfilter$K)
	if(is.null(K)) stop("ICA algorithm does not return the whitening matrix", call. = FALSE)
	return( K )
}
#------------------------------------------------------------------------------

###############################################################################

#------------------------------------------------------------------------------
# cov, cor, coskew, cokurt
#------------------------------------------------------------------------------
rcov = function(object, ...)
{
	UseMethod("rcov")
}

rcor = function(object, ...)
{
	UseMethod("rcor")
}

rcoskew = function(object, ...)
{
	UseMethod("rcoskew")
}

rcokurt = function(object, ...)
{
	UseMethod("rcokurt")
}
#------------------------------------------------------------------------------
# co-moment matrices

# Covariance
.rcovgarch = function(object, type = 1)
{
	A = object@mfit$A
	sigmas = object@mfit$sigmas
	n = dim(sigmas)[1]
	m = dim(sigmas)[2]
	y.V = array(data = NA, dim = c(m, m, n))
	x.V = array(data = NA, dim = c(m, m, n))
	#y.v[1:m,1:m,1:n] = apply(sigmas, 1, FUN = function(x) diag(x)^2)
	for(i in 1:n)
	{
		y.V[,,i] =  diag(sigmas[i,]^2)
		x.V[,,i] = t(A)%*%y.V[,,i]%*%A
	}
	
	if(type == 1) return(x.V) else return(y.V)
}

setMethod("rcov", signature(object = "goGARCHfit"), .rcovgarch)


.rcovgarchfilter = function(object, type = 1)
{
	A = object@mfilter$A
	sigmas = object@mfilter$sigmas
	n = dim(sigmas)[1]
	m = dim(sigmas)[2]
	y.V = array(data = NA, dim = c(m, m, n))
	x.V = array(data = NA, dim = c(m, m, n))
	for(i in 1:n)
	{
		y.V[,,i] =  diag(sigmas[i,]^2)
		x.V[,,i] = t(A)%*%y.V[,,i]%*%A
	}
	
	if(type == 1) return(x.V) else return(y.V)
}

setMethod("rcov", signature(object = "goGARCHfilter"), .rcovgarchfilter)

.rcovgarchforecast = function(object, type = 1)
{
	A = object@mforecast$A
	sigmas = sapply(object@uforecast@forecast, FUN = function(x) sapply(x@forecast$forecasts, FUN = function(y) y[,1], simplify = TRUE))
	if(!is.matrix(sigmas)) sigmas = matrix(sigmas, ncol = dim(A)[2])
	n = dim(sigmas)[1]
	m = dim(sigmas)[2]
	y.V = array(data = NA, dim = c(m, m, n))
	x.V = array(data = NA, dim = c(m, m, n))
	for(i in 1:n)
	{
		y.V[,,i] =  diag(sigmas[i,]^2)
		x.V[,,i] = t(A)%*%y.V[,,i]%*%A
	}
	
	if(type == 1) return(x.V) else return(y.V)
}

setMethod("rcov", signature(object = "goGARCHforecast"), .rcovgarchforecast)


.rcovgarchsim = function(object, type = 1, sim = 1)
{
	m.sim = length(object@msim$seriesSim)
	if(sim > m.sim) stop("\nsim > m.sim!", call. = FALSE)
	A = object@msim$A
	sigmas = object@msim$sigmaSim[[sim]]
	if(!is.matrix(sigmas)) sigmas = matrix(sigmas, ncol = dim(A)[2])
	n = dim(sigmas)[1]
	m = dim(sigmas)[2]
	y.V = array(data = NA, dim = c(m, m, n))
	x.V = array(data = NA, dim = c(m, m, n))
	for(i in 1:n)
	{
		y.V[,,i] =  diag(sigmas[i,]^2)
		x.V[,,i] = t(A)%*%y.V[,,i]%*%A
	}
	
	if(type == 1) return(x.V) else return(y.V)
}

setMethod("rcov", signature(object = "goGARCHsim"), .rcovgarchsim)



.rcovroll = function(object, type = 1)
{
	nf = length(object@forecast)
	np = dim(object@forecast[[1]]@mforecast$A)[2]
	refit.every = object@spec$refit.every
	forecast.length = object@spec$forecast.length
	idx = matrix(NA, ncol = 2, nrow = nf)
	idx[,1] = seq(1, forecast.length, by = refit.every)
	idx[,2] = seq(refit.every, forecast.length, by = refit.every)
	forecast.length = object@spec$forecast.length
	cv = array(NA, dim = c(np, np, forecast.length))
	for(i in 1:nf) cv[,,idx[i,1]:idx[i,2]] = rcov(object@forecast[[i]], type = type)
	return(cv)
}

setMethod("rcov", signature(object = "goGARCHroll"), .rcovroll)


#------------------------------------------------------------------------------
# Correlation
.rcorgarch = function(object)
{
	A = object@mfit$A
	sigmas = object@mfit$sigmas
	n = dim(sigmas)[1]
	m = dim(sigmas)[2]
	x.C = array(data = NA, dim = c(m, m, n))
	y.V = array(data = NA, dim = c(m, m, n))
	x.V = array(data = NA, dim = c(m, m, n))
	for(i in 1:n)
	{
		y.V[,,i] =  diag(sigmas[i,]^2)
		x.V[,,i] = t(A)%*%y.V[,,i]%*%A
		# cov->cor
		x.C[,,i] = cov2cor(x.V[,,i])
		# skew, shape and mu parameters
	}
	
	return(x.C)
}

setMethod("rcor", signature(object = "goGARCHfit"), .rcorgarch)

.rcorgarchfilter = function(object)
{
	A = object@mfilter$A
	sigmas = object@mfilter$sigmas
	n = dim(sigmas)[1]
	m = dim(sigmas)[2]
	x.C = array(data = NA, dim = c(m, m, n))
	y.V = array(data = NA, dim = c(m, m, n))
	x.V = array(data = NA, dim = c(m, m, n))
	for(i in 1:n)
	{
		y.V[,,i] =  diag(sigmas[i,]^2)
		x.V[,,i] = t(A)%*%y.V[,,i]%*%A
		# cov->cor
		x.C[,,i] = cov2cor(x.V[,,i])
		# skew, shape and mu parameters
	}
	
	return(x.C)
}

setMethod("rcor", signature(object = "goGARCHfilter"), .rcorgarchfilter)

.rcorgarchforecast = function(object)
{
	A = object@mforecast$A
	sigmas = sapply(object@uforecast@forecast, FUN = function(x) sapply(x@forecast$forecasts, FUN = function(y) y[,1], simplify = TRUE))
	if(!is.matrix(sigmas)) sigmas = matrix(sigmas, ncol = dim(A)[2])
	n = dim(sigmas)[1]
	m = dim(sigmas)[2]
	x.C = array(data = NA, dim = c(m, m, n))
	y.V = array(data = NA, dim = c(m, m, n))
	x.V = array(data = NA, dim = c(m, m, n))
	for(i in 1:n)
	{
		y.V[,,i] =  diag(sigmas[i,]^2)
		x.V[,,i] = t(A)%*%y.V[,,i]%*%A
		# cov->cor
		x.C[,,i] = cov2cor(x.V[,,i])
		# skew, shape and mu parameters
	}
	
	return(x.C)
}

setMethod("rcor", signature(object = "goGARCHforecast"), .rcorgarchforecast)

.rcorgarchsim = function(object, sim = 1)
{
	m.sim = length(object@msim$seriesSim)
	if(sim > m.sim) stop("\nsim > m.sim!", call. = FALSE)
	A = object@msim$A
	sigmas = object@msim$sigmaSim[[sim]]
	if(!is.matrix(sigmas)) sigmas = matrix(sigmas, ncol = dim(A)[2])
	n = dim(sigmas)[1]
	m = dim(sigmas)[2]
	x.C = array(data = NA, dim = c(m, m, n))
	y.V = array(data = NA, dim = c(m, m, n))
	x.V = array(data = NA, dim = c(m, m, n))
	for(i in 1:n)
	{
		y.V[,,i] =  diag(sigmas[i,]^2)
		x.V[,,i] = t(A)%*%y.V[,,i]%*%A
		x.C[,,i] = cov2cor(x.V[,,i])
	}
	
	return(x.C)
}

setMethod("rcor", signature(object = "goGARCHsim"), .rcorgarchsim)


.rcorroll = function(object)
{
	nf = length(object@forecast)
	np = dim(object@forecast[[1]]@mforecast$A)[2]
	refit.every = object@spec$refit.every
	forecast.length = object@spec$forecast.length
	idx = matrix(NA, ncol = 2, nrow = nf)
	idx[,1] = seq(1, forecast.length, by = refit.every)
	idx[,2] = seq(refit.every, forecast.length, by = refit.every)
	forecast.length = object@spec$forecast.length
	cv = array(NA, dim = c(np, np, forecast.length))
	for(i in 1:nf) cr[,,idx[i,1]:idx[i,2]] = rcor(object@forecast[[i]])
	return(cr)
}

setMethod("rcor", signature(object = "goGARCHroll"), .rcorroll)

#------------------------------------------------------------------------------
# coskewness
.rcoskew = function(object, standardize = TRUE, from = 1, to = 1)
{
	mdist = object@mfit$mdist
	if(mdist == "mvnorm") stop("\nNo co-skewness for mvnorm distribution\n", call. = FALSE)
	ans = switch(mdist,
			manig = .rcoskew.nig(object, standardize, from, to),
			magh = .rcoskew.ghyp(object, standardize, from, to))
	return(ans)
}


setMethod("rcoskew", signature(object = "goGARCHfit"), .rcoskew)


.rcoskew.nig = function(object, standardize = TRUE, from = 1, to = 1)
{
	n = length(seq(from, to, by = 1))
	if( n>100 ) stop("\nMaximum from/to is 100 points due to size restrictions\n", .call = FALSE)
	A = object@mfit$A
	m = dim(A)[2]
	xs = array(NA, dim = c(m , m^2, n))
	sigmas = object@mfit$sigmas
	y.expskew = vector(mode = "list", length = m)
	y.expskew = sapply(object@mfit$distr$y.alphabeta, FUN = function(y) apply(y[from:to,,drop = FALSE], 1, FUN = function(x) .nigskew(x[1], x[2], x[3], x[4])))
	y.expskew = matrix(y.expskew, ncol = m)
	# convert to moments as since the standardized moments do not retain their geometric properties in transformation
	sigmas = sigmas[from:to,,drop = FALSE]
	yskew = y.expskew*(sigmas^3)
	for(i in 1:n){
		S = t(A)%*%diag(sigmas[i,]^2)%*%A
		Sx = sqrt(diag(S))
		tmp = .coskew.ind(yskew[i,])
		if(standardize) xs[,,i] = t(A)%*%tmp%*%kronecker(A, A)/.coskew.sigma(Sx) else xs[,,i] = t(A)%*%tmp%*%kronecker(A, A)
	}
	return(xs)
}


.rcoskew.ghyp = function(object, standardize = TRUE, from = 1, to = 1)
{
	n = length(seq(from, to, by = 1))
	if( n>100 ) stop("\nMaximum from/to is 100 points due to size restrictions\n", .call = FALSE)
	A = object@mfit$A
	m = dim(A)[2]
	xs = array(NA, dim = c(m , m^2, n))
	sigmas = object@mfit$sigmas
	y.expskew = vector(mode = "list", length = m)
	y.expskew = sapply(object@mfit$distr$y.alphabeta, FUN = function(y) apply(y[from:to,,drop = FALSE], 1, FUN = function(x) .ghypskew(x[5], x[1], x[2], x[3], x[4])))
	y.expskew = matrix(y.expskew, ncol = m)
	# convert to moments as since the standardized moments do not retain their geometric properties in transformation
	sigmas = sigmas[from:to,,drop = FALSE]
	yskew = y.expskew*(sigmas^3)
	for(i in 1:n){
		S = t(A)%*%diag(sigmas[i,]^2)%*%A
		Sx = sqrt(diag(S))
		tmp = .coskew.ind(yskew[i,])
		if(standardize) xs[,,i] = t(A)%*%tmp%*%kronecker(A, A)/.coskew.sigma(Sx) else xs[,,i] = t(A)%*%tmp%*%kronecker(A, A)
	}
	return(xs)
}

.rcoskewfilter = function(object, standardize = TRUE, from = 1, to = 1)
{
	mdist = object@mfilter$mdist
	if(mdist == "mvnorm") stop("\nNo co-skewness for mvnorm distribution\n", call. = FALSE)
	ans = switch(mdist,
			manig = .rcoskewfilter.nig(object, standardize, from, to),
			magh = .rcoskewfilter.ghyp(object, standardize, from, to))
	return(ans)
}


setMethod("rcoskew", signature(object = "goGARCHfilter"), .rcoskewfilter)

.rcoskewfilter.nig = function(object, standardize = TRUE, from = 1, to = 1)
{
	n = length(seq(from, to, by = 1))
	if( n>100 ) stop("\nMaximum from/to is 100 points due to size restrictions\n", .call = FALSE)
	A = object@mfilter$A
	m = dim(A)[2]
	xs = array(NA, dim = c(m , m^2, n))
	sigmas = object@mfilter$sigmas
	y.expskew = vector(mode = "list", length = m)
	y.expskew = sapply(object@mfilter$distr$y.alphabeta, FUN = function(y) apply(y[from:to,,drop = FALSE], 1, FUN = function(x) .nigskew(x[1], x[2], x[3], x[4])))
	y.expskew = matrix(y.expskew, ncol = m)
	sigmas = sigmas[from:to,,drop = FALSE]
	yskew = y.expskew*(sigmas^3)
	for(i in 1:n){
		S = t(A)%*%diag(sigmas[i,]^2)%*%A
		Sx = sqrt(diag(S))
		tmp = .coskew.ind(yskew[i,])
		if(standardize) xs[,,i] = t(A)%*%tmp%*%kronecker(A, A)/.coskew.sigma(Sx) else xs[,,i] = t(A)%*%tmp%*%kronecker(A, A)
	}
	return(xs)
}

.rcoskewfilter.ghyp = function(object, standardize = TRUE, from = 1, to = 1)
{
	n = length(seq(from, to, by = 1))
	if( n>100 ) stop("\nMaximum from/to is 100 points due to size restrictions\n", .call = FALSE)
	A = object@mfilter$A
	m = dim(A)[2]
	xs = array(NA, dim = c(m , m^2, n))
	sigmas = object@mfilter$sigmas
	y.expskew = vector(mode = "list", length = m)
	y.expskew = sapply(object@mfilter$distr$y.alphabeta, FUN = function(y) apply(y[from:to,,drop = FALSE], 1, FUN = function(x) .ghypskew(x[5], x[1], x[2], x[3], x[4])))
	y.expskew = matrix(y.expskew, ncol = m)
	sigmas = sigmas[from:to,,drop = FALSE]
	yskew = y.expskew*(sigmas^3)
	for(i in 1:n){
		S = t(A)%*%diag(sigmas[i,]^2)%*%A
		Sx = sqrt(diag(S))
		tmp = .coskew.ind(yskew[i,])
		if(standardize) xs[,,i] = t(A)%*%tmp%*%kronecker(A, A)/.coskew.sigma(Sx) else xs[,,i] = t(A)%*%tmp%*%kronecker(A, A)
	}
	return(xs)
}


.rcoskewforecast = function(object, standardize = TRUE, from = 1, to = 1)
{
	mdist = object@mforecast$mdist
	if(mdist == "mvnorm") stop("\nNo co-skewness for mvnorm distribution\n", call. = FALSE)
	ans = switch(mdist,
			manig = .rcoskewforecast.nig(object, standardize, from, to),
			magh = .rcoskewforecast.ghyp(object, standardize, from, to))
	return(ans)
}


setMethod("rcoskew", signature(object = "goGARCHforecast"), .rcoskewforecast)

.rcoskewforecast.nig = function(object, standardize = TRUE, from = 1, to = 1)
{
	n = length(seq(from, to, by = 1))
	if( n>100 ) stop("\nMaximum from/to is 100 points due to size restrictions\n", .call = FALSE)
	A = object@mforecast$A
	m = dim(A)[2]
	xs = array(NA, dim = c(m , m^2, n))
	sigmas = sapply(object@uforecast@forecast, FUN = function(x) sapply(x@forecast$forecasts, FUN = function(y) y[,1], simplify = TRUE))
	if(!is.matrix(sigmas)) sigmas = matrix(sigmas, ncol = dim(A)[2])
	
	y.expskew = vector(mode = "list", length = m)
	y.expskew = sapply(object@mforecast$distr$y.alphabeta, FUN = function(y) apply(y[from:to,,drop = FALSE], 1, FUN = function(x) .nigskew(x[1], x[2], x[3], x[4])))
	y.expskew = matrix(y.expskew, ncol = m)
	sigmas = sigmas[from:to,,drop = FALSE]
	yskew = y.expskew*(sigmas^3)
	for(i in 1:n){
		S = t(A)%*%diag(sigmas[i,]^2)%*%A
		Sx = sqrt(diag(S))
		tmp = .coskew.ind(yskew[i,])
		if(standardize) xs[,,i] = t(A)%*%tmp%*%kronecker(A, A)/.coskew.sigma(Sx) else xs[,,i] = t(A)%*%tmp%*%kronecker(A, A)
	}
	return(xs)
}

.rcoskewforecast.ghyp = function(object, standardize = TRUE, from = 1, to = 1)
{
	n = length(seq(from, to, by = 1))
	if( n>100 ) stop("\nMaximum from/to is 100 points due to size restrictions\n", .call = FALSE)
	A = object@mforecast$A
	m = dim(A)[2]
	xs = array(NA, dim = c(m , m^2, n))
	sigmas = sapply(object@uforecast@forecast, FUN = function(x) sapply(x@forecast$forecasts, FUN = function(y) y[,1], simplify = TRUE))
	if(!is.matrix(sigmas)) sigmas = matrix(sigmas, ncol = dim(A)[2])
	
	y.expskew = vector(mode = "list", length = m)
	y.expskew = sapply(object@mforecast$distr$y.alphabeta, FUN = function(y) apply(y[from:to,,drop = FALSE], 1, FUN = function(x) .ghypskew(x[5], x[1], x[2], x[3], x[4])))
	y.expskew = matrix(y.expskew, ncol = m)
	sigmas = sigmas[from:to,,drop = FALSE]
	yskew = y.expskew*(sigmas^3)
	for(i in 1:n){
		S = t(A)%*%diag(sigmas[i,]^2)%*%A
		Sx = sqrt(diag(S))
		tmp = .coskew.ind(yskew[i,])
		if(standardize) xs[,,i] = t(A)%*%tmp%*%kronecker(A, A)/.coskew.sigma(Sx) else xs[,,i] = t(A)%*%tmp%*%kronecker(A, A)
	}
	return(xs)
}



.rcoskewsim = function(object, standardize = TRUE, from = 1, to = 1, sim = 1)
{
	mdist = object@msim$mdist
	if(mdist == "mvnorm") stop("\nNo co-skewness for mvnorm distribution\n", call. = FALSE)
	ans = switch(mdist,
			manig = .rcoskewsim.nig(object, standardize, from, to, sim),
			magh = .rcoskewsim.ghyp(object, standardize, from, to, sim))
	return(ans)
}


setMethod("rcoskew", signature(object = "goGARCHsim"), .rcoskewsim)


.rcoskewsim.nig = function(object,  standardize = TRUE, from = 1, to = 1,  sim = 1)
{
	n = length(seq(from, to, by = 1))
	if( n>100 ) stop("\nMaximum from/to is 100 points due to size restrictions\n", .call = FALSE)
	A = object@msim$A
	m = dim(A)[2]
	xs = array(NA, dim = c(m , m^2, n))
	sigmas = object@msim$sigmaSim[[sim]]
	if(!is.matrix(sigmas)) sigmas = matrix(sigmas, ncol = dim(A)[2])
	y.expskew = vector(mode = "list", length = m)
	y.expskew = sapply(object@msim$y.alphabeta, FUN = function(y) apply(y[[sim]][from:to,,drop = FALSE], 1, FUN = function(x) .nigskew(x[1], x[2], x[3], x[4])))
	y.expskew = matrix(y.expskew, ncol = m)
	sigmas = sigmas[from:to,,drop = FALSE]
	yskew = y.expskew*(sigmas^3)
	for(i in 1:n){
		S = t(A)%*%diag(sigmas[i,]^2)%*%A
		Sx = sqrt(diag(S))
		tmp = .coskew.ind(yskew[i,])
		if(standardize) xs[,,i] = t(A)%*%tmp%*%kronecker(A, A)/.coskew.sigma(Sx) else xs[,,i] = t(A)%*%tmp%*%kronecker(A, A)
	}
	return(xs)
}

.rcoskewsim.ghyp = function(object,  standardize = TRUE, from = 1, to = 1,  sim = 1)
{
	n = length(seq(from, to, by = 1))
	if( n>100 ) stop("\nMaximum from/to is 100 points due to size restrictions\n", .call = FALSE)
	A = object@msim$A
	m = dim(A)[2]
	xs = array(NA, dim = c(m , m^2, n))
	sigmas = object@msim$sigmaSim[[sim]]
	if(!is.matrix(sigmas)) sigmas = matrix(sigmas, ncol = dim(A)[2])
	y.expskew = vector(mode = "list", length = m)
	y.expskew = sapply(object@msim$y.alphabeta, FUN = function(y) apply(y[[sim]][from:to,,drop = FALSE], 1, FUN = function(x) .ghypskew(x[5], x[1], x[2], x[3], x[4])))
	y.expskew = matrix(y.expskew, ncol = m)
	sigmas = sigmas[from:to,,drop = FALSE]
	yskew = y.expskew*(sigmas^3)
	for(i in 1:n){
		S = t(A)%*%diag(sigmas[i,]^2)%*%A
		Sx = sqrt(diag(S))
		tmp = .coskew.ind(yskew[i,])
		if(standardize) xs[,,i] = t(A)%*%tmp%*%kronecker(A, A)/.coskew.sigma(Sx) else xs[,,i] = t(A)%*%tmp%*%kronecker(A, A)
	}
	return(xs)
}

#------------------------------------------------------------------------------
# cokurtosis
.rcokurt = function(object, standardize = TRUE, from = 1, to = 1, sim = 1)
{
	mdist = object@mfit$mdist
	if(mdist == "mvnorm"){
		stop("\nNo co-kurtosis for mvnorm distribution\n", call. = FALSE)
	} else{
		cat("\nrgarch: cokurtosis under revision")
	}
	#ans = switch(mdist,
	#		manig = .rcokurt.nig(object, standardize, from, to),
	#		magh = .rcokurt.ghyp(object, standardize, from, to))
	# return(ans)
	return( 0 )
}

setMethod("rcokurt", signature(object = "goGARCHfit"), .rcokurt)

.rcokurt.nig = function(object, standardize = TRUE, from = 1, to = 1)
{
	n = length(seq(from, to, by = 1))
	if( n>100 ) stop("\nMaximum from/to is 100 points due to size restrictions\n", .call = FALSE)
	A = object@mfit$A
	m = dim(A)[2]
	if( m>20 ) stop("\nMaximum number of assets allowed is 20 due to size restrictions\n", .call = FALSE)
	xs = array(NA, dim = c(m , m^3, n))
	sigmas = object@mfit$sigmas
	y.expexkurt = vector(mode = "list", length = m)
	y.expexkurt = sapply(object@mfit$distr$y.alphabeta, FUN = function(y) apply(y[from:to,,drop = FALSE], 1, FUN = function(x) .nigexkurt(x[1], x[2], x[3], x[4])))
	y.expexkurt = matrix(y.expexkurt, ncol = m)
	sigmas = sigmas[from:to,,drop = FALSE]
	ykurt = y.expexkurt*(sigmas^4)
	
	# The indices need to be adjusted to include {iijj}  (and permutations) 
	# not only {iiii} and {jjjj}.
	# {iijj}  = sigma_i * sigma_j (proofs in thesis)
	
	for(i in 1:n){
		S = t(A)%*%diag(sigmas[i,]^2)%*%A
		Sx = sqrt(diag(S))
		tmp = .cokurt.ind(ykurt[i,])
		if(standardize) xs[,,i] = t(A)%*%tmp%*%kronecker(A, kronecker(A, A))/.cokurt.sigma(Sx) else xs[,,i] = t(A)%*%tmp%*%kronecker(A, kronecker(A, A))
	}
	return(xs)
}

.rcokurt.ghyp = function(object, standardize = TRUE, from = 1, to = 1)
{
	n = length(seq(from, to, by = 1))
	if( n>100 ) stop("\nMaximum from/to is 100 points due to size restrictions\n", .call = FALSE)
	A = object@mfit$A
	m = dim(A)[2]
	if( m>20 ) stop("\nMaximum number of assets allowed is 20 due to size restrictions\n", .call = FALSE)
	xs = array(NA, dim = c(m , m^3, n))
	sigmas = object@mfit$sigmas
	y.expexkurt = vector(mode = "list", length = m)
	y.expexkurt = sapply(object@mfit$distr$y.alphabeta, FUN = function(y) apply(y[from:to,,drop = FALSE], 1, FUN = function(x) .ghypexkurt(x[5], x[1], x[2], x[3], x[4])))
	y.expexkurt = matrix(y.expexkurt, ncol = m)
	sigmas = sigmas[from:to,,drop = FALSE]
	ykurt = y.expexkurt*(sigmas^4)
	# The indices need to be adjusted to include {iijj}  (and permutations) 
	# not only {iiii} and {jjjj}.
	# {iijj}  = sigma_i * sigma_j (proofs in thesis)
	for(i in 1:n){
		S = t(A)%*%diag(sigmas[i,]^2)%*%A
		Sx = sqrt(diag(S))
		tmp = .cokurt.ind(ykurt[i,])
		if(standardize) xs[,,i] = t(A)%*%tmp%*%kronecker(A, kronecker(A, A))/.cokurt.sigma(Sx) else xs[,,i] = t(A)%*%tmp%*%kronecker(A, kronecker(A, A))
	}
	return(xs)
}

.rcokurtfilter = function(object, standardize = TRUE, from = 1, to = 1)
{
	mdist = object@mfilter$mdist
	if(mdist == "mvnorm"){
		stop("\nNo co-kurtosis for mvnorm distribution\n", call. = FALSE)
	} else{
		cat("\nrgarch: cokurtosis under revision")
	}
	#ans = switch(mdist,
	#		manig = .rcokurtfilter.nig(object, standardize, from, to),
	#		magh = .rcokurtfilter.ghyp(object, standardize, from, to))
	#return(ans)
	return( 0 )
}

setMethod("rcokurt", signature(object = "goGARCHfilter"), .rcokurtfilter)


.rcokurtfilter.nig = function(object, standardize = TRUE, from = 1, to = 1)
{
	n = length(seq(from, to, by = 1))
	if( n>100 ) stop("\nMaximum from/to is 100 points due to size restrictions\n", .call = FALSE)
	A = object@mfilter$A
	m = dim(A)[2]
	if( m>20 ) stop("\nMaximum number of assets allowed is 20 due to size restrictions\n", .call = FALSE)
	xs = array(NA, dim = c(m , m^3, n))
	sigmas = object@mfilter$sigmas
	y.expexkurt = vector(mode = "list", length = m)
	y.expexkurt = sapply(object@mfilter$distr$y.alphabeta, FUN = function(y) apply(y[from:to,,drop = FALSE], 1, FUN = function(x) .nigexkurt(x[1], x[2], x[3], x[4])))
	y.expexkurt = matrix(y.expexkurt, ncol = m)
	sigmas = sigmas[from:to,,drop = FALSE]
	ykurt = y.expexkurt*(sigmas^4)
	# The indices need to be adjusted to include {iijj}  (and permutations) 
	# not only {iiii} and {jjjj}.
	# {iijj}  = sigma_i * sigma_j (proofs in thesis)
	for(i in 1:n){
		S = t(A)%*%diag(sigmas[i,]^2)%*%A
		Sx = sqrt(diag(S))
		tmp = .cokurt.ind(ykurt[i,])
		if(standardize) xs[,,i] = t(A)%*%tmp%*%kronecker(A, kronecker(A, A))/.cokurt.sigma(Sx) else xs[,,i] = t(A)%*%tmp%*%kronecker(A, kronecker(A, A))
	}
	return(xs)
}

.rcokurtfilter.ghyp = function(object, standardize = TRUE, from = 1, to = 1)
{
	n = length(seq(from, to, by = 1))
	if( n>100 ) stop("\nMaximum from/to is 100 points due to size restrictions\n", .call = FALSE)
	A = object@mfilter$A
	m = dim(A)[2]
	if( m>20 ) stop("\nMaximum number of assets allowed is 20 due to size restrictions\n", .call = FALSE)
	xs = array(NA, dim = c(m , m^3, n))
	sigmas = object@mfilter$sigmas
	y.expexkurt = vector(mode = "list", length = m)
	y.expexkurt = sapply(object@mfilter$distr$y.alphabeta, FUN = function(y) apply(y[from:to,,drop = FALSE], 1, FUN = function(x) .ghypexkurt(x[5], x[1], x[2], x[3], x[4])))
	y.expexkurt = matrix(y.expexkurt, ncol = m)
	sigmas = sigmas[from:to,,drop = FALSE]
	ykurt = y.expexkurt*(sigmas^4)
	# The indices need to be adjusted to include {iijj}  (and permutations) 
	# not only {iiii} and {jjjj}.
	# {iijj}  = sigma_i * sigma_j (proofs in thesis)
	for(i in 1:n){
		S = t(A)%*%diag(sigmas[i,]^2)%*%A
		Sx = sqrt(diag(S))
		tmp = .cokurt.ind(ykurt[i,])
		if(standardize) xs[,,i] = t(A)%*%tmp%*%kronecker(A, kronecker(A, A))/.cokurt.sigma(Sx) else xs[,,i] = t(A)%*%tmp%*%kronecker(A, kronecker(A, A))
	}
	return(xs)
}

.rcokurtforecast = function(object, standardize = TRUE, from = 1, to = 1)
{
	mdist = object@mforecast$mdist
	if(mdist == "mvnorm"){
		stop("\nNo co-kurtosis for mvnorm distribution\n", call. = FALSE)
	} else{
		cat("\nrgarch: cokurtosis under revision")
	}
	#ans = switch(mdist,
	#		manig = .rcokurtforecast.nig(object, standardize, from, to),
	#		magh = .rcokurtforecast.ghyp(object, standardize, from, to))
	#return(ans)
	return(0)
}

setMethod("rcokurt", signature(object = "goGARCHforecast"), .rcokurtforecast)

.rcokurtforecast.nig = function(object, standardize = TRUE, from = 1, to = 1)
{
	n = length(seq(from, to, by = 1))
	if( n>100 ) stop("\nMaximum from/to is 100 points due to size restrictions\n", .call = FALSE)
	A = object@mforecast$A
	m = dim(A)[2]
	if( m>20 ) stop("\nMaximum number of assets allowed is 20 due to size restrictions\n", .call = FALSE)
	xs = array(NA, dim = c(m , m^3, n))
	sigmas = sapply(object@uforecast@forecast, FUN = function(x) sapply(x@forecast$forecasts, FUN = function(y) y[,1], simplify = TRUE))
	if(!is.matrix(sigmas)) sigmas = matrix(sigmas, ncol = dim(A)[2])
	y.expexkurt = vector(mode = "list", length = m)
	y.expexkurt = sapply(object@mforecast$distr$y.alphabeta, FUN = function(y) apply(y[from:to,,drop = FALSE], 1, FUN = function(x) .nigexkurt(x[1], x[2], x[3], x[4])))
	y.expexkurt = matrix(y.expexkurt, ncol = m)
	sigmas = sigmas[from:to,,drop = FALSE]
	ykurt = y.expexkurt*(sigmas^4)
	# The indices need to be adjusted to include {iijj}  (and permutations) 
	# not only {iiii} and {jjjj}.
	# {iijj}  = sigma_i * sigma_j (proofs in thesis)
	for(i in 1:n){
		S = t(A)%*%diag(sigmas[i,]^2)%*%A
		Sx = sqrt(diag(S))
		tmp = .cokurt.ind(ykurt[i,])
		if(standardize) xs[,,i] = t(A)%*%tmp%*%kronecker(A, kronecker(A, A))/.cokurt.sigma(Sx) else xs[,,i] = t(A)%*%tmp%*%kronecker(A, kronecker(A, A))
	}
	return(xs)
}

.rcokurtforecast.ghyp = function(object, standardize = TRUE, from = 1, to = 1)
{
	n = length(seq(from, to, by = 1))
	if( n>100 ) stop("\nMaximum from/to is 100 points due to size restrictions\n", .call = FALSE)
	A = object@mforecast$A
	m = dim(A)[2]
	if( m>20 ) stop("\nMaximum number of assets allowed is 20 due to size restrictions\n", .call = FALSE)
	xs = array(NA, dim = c(m , m^3, n))
	sigmas = sapply(object@uforecast@forecast, FUN = function(x) sapply(x@forecast$forecasts, FUN = function(y) y[,1], simplify = TRUE))
	if(!is.matrix(sigmas)) sigmas = matrix(sigmas, ncol = dim(A)[2])
	y.expexkurt = vector(mode = "list", length = m)
	y.expexkurt = sapply(object@mforecast$distr$y.alphabeta, FUN = function(y) apply(y[from:to,,drop = FALSE], 1, FUN = function(x) .ghypexkurt(x[5], x[1], x[2], x[3], x[4])))
	y.expexkurt = matrix(y.expexkurt, ncol = m)
	sigmas = sigmas[from:to,,drop = FALSE]
	ykurt = y.expexkurt*(sigmas^4)
	# The indices need to be adjusted to include {iijj}  (and permutations) 
	# not only {iiii} and {jjjj}.
	# {iijj}  = sigma_i * sigma_j (proofs in thesis)
	for(i in 1:n){
		S = t(A)%*%diag(sigmas[i,]^2)%*%A
		Sx = sqrt(diag(S))
		tmp = .cokurt.ind(ykurt[i,])
		if(standardize) xs[,,i] = t(A)%*%tmp%*%kronecker(A, kronecker(A, A))/.cokurt.sigma(Sx) else xs[,,i] = t(A)%*%tmp%*%kronecker(A, kronecker(A, A))
	}
	return(xs)
}

.rcokurtsim = function(object, standardize = TRUE, from = 1, to = 1, sim = 1)
{
	mdist = object@msim$mdist
	if(mdist == "mvnorm"){
		stop("\nNo co-kurtosis for mvnorm distribution\n", call. = FALSE)
	} else{
		cat("\nrgarch: cokurtosis under revision")
	}
	#ans = switch(mdist,
	#		manig = .rcokurtsim.nig(object, standardize, from, to, sim),
	#		magh = .rcokurtsim.ghyp(object, standardize, from, to, sim))
	#return(ans)
	return(0)
}

setMethod("rcokurt", signature(object = "goGARCHsim"), .rcokurtsim)

.rcokurtsim.nig = function(object, standardize = TRUE, from = 1, to = 1, sim = 1)
{
	n = length(seq(from, to, by = 1))
	if( n>100 ) stop("\nMaximum from/to is 100 points due to size restrictions\n", .call = FALSE)
	A = object@msim$A
	m = dim(A)[2]
	if( m>20 ) stop("\nMaximum number of assets allowed is 20 due to size restrictions\n", .call = FALSE)
	xs = array(NA, dim = c(m , m^3, n))
	sigmas = object@msim$sigmaSim[[sim]]
	if(!is.matrix(sigmas)) sigmas = matrix(sigmas, ncol = dim(A)[2])
	y.expexkurt = vector(mode = "list", length = m)
	y.expexkurt = sapply(object@msim$y.alphabeta, FUN = function(y) apply(y[[sim]][from:to,,drop = FALSE], 1, FUN = function(x) .nigexkurt(x[1], x[2], x[3], x[4])))
	y.expexkurt = matrix(y.expexkurt, ncol = m)
	sigmas = sigmas[from:to,,drop = FALSE]
	ykurt = y.expexkurt*(sigmas^4)
	# The indices need to be adjusted to include {iijj}  (and permutations) 
	# not only {iiii} and {jjjj}.
	# {iijj}  = sigma_i * sigma_j (proofs in thesis)
	for(i in 1:n){
		S = t(A)%*%diag(sigmas[i,]^2)%*%A
		Sx = sqrt(diag(S))
		tmp = .cokurt.ind(ykurt[i,])
		if(standardize) xs[,,i] = t(A)%*%tmp%*%kronecker(A, kronecker(A, A))/.cokurt.sigma(Sx) else xs[,,i] = t(A)%*%tmp%*%kronecker(A, kronecker(A, A))
	}
	return(xs)
}

.rcokurtsim.ghyp = function(object, standardize = TRUE, from = 1, to = 1, sim = 1)
{
	n = length(seq(from, to, by = 1))
	if( n>100 ) stop("\nMaximum from/to is 100 points due to size restrictions\n", .call = FALSE)
	A = object@msim$A
	m = dim(A)[2]
	if( m>20 ) stop("\nMaximum number of assets allowed is 20 due to size restrictions\n", .call = FALSE)
	xs = array(NA, dim = c(m , m^3, n))
	sigmas = object@msim$sigmaSim[[sim]]
	if(!is.matrix(sigmas)) sigmas = matrix(sigmas, ncol = dim(A)[2])
	y.expexkurt = vector(mode = "list", length = m)
	y.expexkurt = sapply(object@msim$y.alphabeta, FUN = function(y) apply(y[[sim]][from:to,,drop = FALSE], 1, FUN = function(x) .ghypexkurt(x[5], x[1], x[2], x[3], x[4])))
	y.expexkurt = matrix(y.expexkurt, ncol = m)
	sigmas = sigmas[from:to,,drop = FALSE]
	ykurt = y.expexkurt*(sigmas^4)
	# The indices need to be adjusted to include {iijj}  (and permutations) 
	# not only {iiii} and {jjjj}.
	# {iijj}  = sigma_i * sigma_j (proofs in thesis)
	for(i in 1:n){
		S = t(A)%*%diag(sigmas[i,]^2)%*%A
		Sx = sqrt(diag(S))
		tmp = .cokurt.ind(ykurt[i,])
		if(standardize) xs[,,i] = t(A)%*%tmp%*%kronecker(A, kronecker(A, A))/.cokurt.sigma(Sx) else xs[,,i] = t(A)%*%tmp%*%kronecker(A, kronecker(A, A))
	}
	return(xs)
}

#------------------------------------------------------------------------------
###############################################################################


###############################################################################
#------------------------------------------------------------------------------
# Geometric Portfolio Moments
# 
gportmoments = function(object, ...)
{
	UseMethod("gportmoments")
}

## filter out among the various distributions:
## TODO
.gportmoments.garchfit = function(object, weights, debug = FALSE)
{
	
	A = object@mfit$A
	n = dim(object@mfit$sigmas)[1]
	m = dim(A)[1]
	if(is.vector(weights)){
		if(length(weights)!=dim(A)[2]) stop("\nIncorrect weight vector length\n", call. = FALSE)
		weights = matrix(weights, ncol = m, nrow = n, byrow = TRUE)
	} else if(is.matrix(weights)){
		n = dim(weights)[1]
		m = dim(weights)[2]
		if(m!=dim(A)[2]) stop("\nIncorrect column dimension for weights matrix\n", call. = FALSE)
		if(n!=dim(object@mfit$sigmas)[1]) stop("\nIncorrect row dimension for weights matrix\n", call. = FALSE)
	}

	if(m<20) pmom = matrix(NA, ncol = 4, nrow = n) else pmom = matrix(NA, ncol = 3, nrow = n)
	
	x.V = rcov(object)
	for(i in 1:n){
		w = weights[i,]
		pmom[i, 1] = object@mfit$x.M[i,]%*%w
		pmom[i, 2] = sqrt(w%*%x.V[,,i]%*%w)
		pmom[i, 3] = (w%*%rcoskew(object, standardize = FALSE, from = i, to = i)[,,1]%*%kronecker(w, w))/pmom[i, 2]^3
		#if(m<20) pmom[i, 4] = w%*%rcokurt(object, standardize = FALSE, from = i, to = i)[,,1]%*%kronecker(w, kronecker(w, w))/pmom[i, 2]^4
		if(debug) print(i)
	}
	return(pmom)
}

setMethod("gportmoments", signature(object = "goGARCHfit"), .gportmoments.garchfit)


.gportmoments.garchfilter = function(object, weights, debug = FALSE)
{
	A = object@mfilter$A
	n = dim(object@mfilter$sigmas)[1]
	m = dim(A)[1]
	if(is.vector(weights)){
		if(length(weights)!=dim(A)[2]) stop("\nIncorrect weight vector length\n", call. = FALSE)
		weights = matrix(weights, ncol = m, nrow = n, byrow = TRUE)
	} else if(is.matrix(weights)){
		n = dim(weights)[1]
		m = dim(weights)[2]
		if(m!=dim(A)[2]) stop("\nIncorrect column dimension for weights matrix\n", call. = FALSE)
		if(n!=dim(object@mfilter$sigmas)[1]) stop("\nIncorrect row dimension for weights matrix\n", call. = FALSE)
	}
	
	if(m<20) pmom = matrix(NA, ncol = 4, nrow = n) else pmom = matrix(NA, ncol = 3, nrow = n)
	A = object@mfilter$A
	x.V = rcov(object)
	for(i in 1:n){
		w = weights[i,]
		pmom[i, 1] = object@mfilter$x.M[i,]%*%w
		pmom[i, 2] = sqrt(w%*%x.V[,,i]%*%w)
		pmom[i, 3] = (w%*%rcoskew(object, standardize = FALSE, from = i, to = i)[,,1]%*%kronecker(w, w))/pmom[i, 2]^3
		#if(m<20) pmom[i, 4] = w%*%rcokurt(object, standardize = FALSE, from = i, to = i)[,,1]%*%kronecker(w, kronecker(w, w))/pmom[i, 2]^4
		if(debug) print(i)
	}
	return(pmom)
}

setMethod("gportmoments", signature(object = "goGARCHfilter"), .gportmoments.garchfilter)


.gportmoments.garchforecast = function(object, weights, debug = FALSE)
{
	A = object@mforecast$A
	n = dim(object@mforecast$sigmas)[1]
	m = dim(A)[1]
	if(is.vector(weights)){
		if(length(weights)!=dim(A)[2]) stop("\nIncorrect weight vector length\n", call. = FALSE)
		weights = matrix(weights, ncol = m, nrow = n, byrow = TRUE)
	} else if(is.matrix(weights)){
		n = dim(weights)[1]
		m = dim(weights)[2]
		if(m!=dim(A)[2]) stop("\nIncorrect column dimension for weights matrix\n", call. = FALSE)
		# we won't force the forecast matrix to be equal to actual forecast length as long as it
		# is less.
		if(n>dim(object@mforecast$sigmas)[1]) stop("\nIncorrect row dimension for weights matrix\n", call. = FALSE)
	}
	
	if(m<20) pmom = matrix(NA, ncol = 4, nrow = n) else pmom = matrix(NA, ncol = 3, nrow = n)
	A = object@mforecast$A
	x.V = rcov(object)
	for(i in 1:n){
		w = weights[i,]
		pmom[i, 1] = object@mforecast$x.M[i,]%*%w
		pmom[i, 2] = sqrt(w%*%x.V[,,i]%*%w)
		pmom[i, 3] = (w%*%rcoskew(object, standardize = FALSE, from = i, to = i)[,,1]%*%kronecker(w, w))/pmom[i, 2]^3
		#if(m<20) pmom[i, 4] = (w%*%rcokurt(object, standardize = FALSE, from = i, to = i)[,,1]%*%kronecker(w, kronecker(w, w)))/pmom[i, 2]^4
		if(debug) print(i)
	}
	return(pmom)
}

setMethod("gportmoments", signature(object = "goGARCHforecast"), .gportmoments.garchforecast)


.gportmoments.garchsim = function(object, weights, debug = FALSE, sim = 1)
{
	A = object@msim$A
	n = dim(object@msim$sigmas)[1]
	m = dim(A)[1]
	if(is.vector(weights)){
		if(length(weights)!=dim(A)[2]) stop("\nIncorrect weight vector length\n", call. = FALSE)
		weights = matrix(weights, ncol = m, nrow = n, byrow = TRUE)
	} else if(is.matrix(weights)){
		n = dim(weights)[1]
		m = dim(weights)[2]
		if(m!=dim(A)[2]) stop("\nIncorrect column dimension for weights matrix\n", call. = FALSE)
		# we won't force the forecast matrix to be equal to actual forecast length as long as it
		# is less.
		if(n>dim(object@msim$sigmas[[sim]])[1]) stop("\nIncorrect row dimension for weights matrix\n", call. = FALSE)
	}
	if(m<20) pmom = matrix(NA, ncol = 4, nrow = n) else pmom = matrix(NA, ncol = 3, nrow = n)
	A = object@msim$A
	x.V = rcov(object, sim = sim)
	for(i in 1:n){
		w = weights[i,]
		pmom[i, 1] = object@msim$seriesSim[[sim]][i,]%*%w
		pmom[i, 2] = sqrt(w%*%x.V[,,i]%*%w)
		pmom[i, 3] = (w%*%rcoskew(object, standardize = FALSE, from = i, to = i, sim = sim)[,,1]%*%kronecker(w, w))/pmom[i, 2]^3
		#if(m<20) pmom[i, 4] = (w%*%rcokurt(object, standardize = FALSE, from = i, to = i, sim = sim)[,,1]%*%kronecker(w, kronecker(w, w)))/pmom[i, 2]^4
		if(debug) print(i)
	}
	return(pmom)
}

setMethod("gportmoments", signature(object = "goGARCHsim"), .gportmoments.garchsim)
#------------------------------------------------------------------------------
###############################################################################


###############################################################################
#------------------------------------------------------------------------------
# Numerical Portfolio Moments (FFT based density integration)
nportmoments = function(object, ...)
{
	UseMethod("nportmoments")
}


.nportmoments = function(object, subdivisions = 200, rel.tol = .Machine$double.eps^0.5, trace = 0, ...)
{
	md = object@dist$dist
	ans = switch(md,
			manig = .nmoments.fft(object, subdivisions = subdivisions, 
					rel.tol = rel.tol, trace = trace, ...),
			magh = .nmoments.fft(object, subdivisions = subdivisions, 
					rel.tol = rel.tol, trace = trace, ...),
			mvnorm = .nmoments.lin(object))
	return(ans)
}

setMethod("nportmoments", signature(object = "goGARCHfft"), .nportmoments)

.nmoments.fft = function(object, subdivisions = 200, rel.tol = .Machine$double.eps^0.5, trace = 0, ...)
{
	
	if( object@dist$support.method == "user" ){
		xpdf = seq(object@dist$fft.support[1], object@dist$fft.support[2], by = object@dist$fft.by)
		yd = object@dist$y
		N = dim(yd)[2]
		nmom = matrix(NA, ncol = 4, nrow = N)
		
		for(i in 1:N){
			dmad = .makeDNew(x = xpdf, dx = yd[,i], h = object@dist$fft.by, standM = "sum")
			fm1 = function(x) x * dmad(x)
			nmom[i,1] = integrate(fm1, -Inf, Inf, subdivisions = 150, rel.tol = .Machine$double.eps^0.25)$value
			m1 = nmom[i,1]
			fm2 = function(x) ((x-m1)^2) * dmad(x)
			nmom[i,2] = sqrt(integrate(fm2, -Inf, Inf, subdivisions = subdivisions, rel.tol = rel.tol, stop.on.error = FALSE)$value)
			fm3 = function(x) ((x-m1)^3) * dmad(x)
			nmom[i,3] = (integrate(fm3, -Inf, Inf, subdivisions = subdivisions, rel.tol = rel.tol, stop.on.error = FALSE)$value)
			nmom[i,3] = nmom[i,3]/nmom[i,2]^3
			fm4 = function(x) ((x-m1)^4) * dmad(x)
			nmom[i,4] = (integrate(fm4, -Inf, Inf, subdivisions = subdivisions, rel.tol = rel.tol, stop.on.error = FALSE)$value)
			nmom[i,4] = (nmom[i,4]/nmom[i,2]^4)-3
			if(trace) print(i)
		}
		colnames(nmom) = c("mean", "sigma", "m3", "m4")
	} else{
		yd = object@dist$y
		N = length(yd)
		nmom = matrix(NA, ncol = 4, nrow = N)
		
		for(i in 1:N){
			xpdf = seq(object@dist$support.user[i, 1], object@dist$support.user[i, 2], by = object@dist$fft.by)
			dmad = .makeDNew(x = xpdf, dx = yd[[i]], h = object@dist$fft.by, standM = "sum")
			fm1 = function(x) x * dmad(x)
			nmom[i,1] = integrate(fm1, -Inf, Inf, subdivisions = 150, rel.tol = .Machine$double.eps^0.25)$value
			m1 = nmom[i,1]
			fm2 = function(x) ((x-m1)^2) * dmad(x)
			nmom[i,2] = sqrt(integrate(fm2, -Inf, Inf, subdivisions = subdivisions, rel.tol = rel.tol, stop.on.error = FALSE)$value)
			fm3 = function(x) ((x-m1)^3) * dmad(x)
			nmom[i,3] = (integrate(fm3, -Inf, Inf, subdivisions = subdivisions, rel.tol = rel.tol, stop.on.error = FALSE)$value)
			nmom[i,3] = nmom[i,3]/nmom[i,2]^3
			fm4 = function(x) ((x-m1)^4) * dmad(x)
			nmom[i,4] = (integrate(fm4, -Inf, Inf, subdivisions = subdivisions, rel.tol = rel.tol, stop.on.error = FALSE)$value)
			nmom[i,4] = (nmom[i,4]/nmom[i,2]^4)-3
			if(trace) print(i)
		}
		colnames(nmom) = c("mean", "sigma", "m3", "m4")
	}
	return(nmom)
}

.nmoments.lin = function(object)
{
	return( object@dist$y )
}

#------------------------------------------------------------------------------
###############################################################################


###############################################################################
#------------------------------------------------------------------------------
# Convolution
convolution = function(object, ...)
{
	UseMethod("convolution")
}


.convolution.gogarchfit = function(object, weights, fft.step = 0.01, fft.by = 0.0001, fft.support = c(-1, 1), 
		support.method = c("user", "adaptive"), use.ff = TRUE, parallel = FALSE, 
		parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2), trace = 0, ...)
{
	dist = object@mfit$mdist
	pars = object@mfit$distr$z.D
	Sigma = rcov(object, type = 2)
	Mu = object@mfit$x.M
	A = object@mfit$A
	support.method = support.method[1]
	
	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])
		}
	}		
	debug = as.logical(trace)
	ans = switch(dist,
			manig = .convolve.nig(nigpars = pars, Sigma = Sigma, Mu = Mu, A = A, 
					weights = weights, fft.step = fft.step, fft.by = fft.by, fft.support = fft.support, 
					use.ff = use.ff, support.method = support.method, parallel = parallel, 
					parallel.control = parallel.control, debug = debug),
			magh = .convolve.gh(ghpars = pars, Sigma = Sigma, Mu = Mu, A = A, 
					weights = weights, fft.step = fft.step, fft.by = fft.by, fft.support = fft.support, 
					use.ff = use.ff, support.method = support.method, parallel = parallel, 
					parallel.control = parallel.control, debug = debug),
			mvnorm = .convolve.norm(Sigma  = Sigma, Mu = Mu, A = A, weights = weights))
	return(ans)
}

setMethod("convolution", signature(object = "goGARCHfit"), .convolution.gogarchfit)

.convolution.gogarchfilter = function(object, weights, fft.step = 0.001, fft.by = 0.0001, fft.support = c(-1, 1), 
		support.method = c("user", "adaptive"), use.ff = TRUE, parallel = FALSE, 
		parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2), debug = 0, ...)
{
	dist = object@mfilter$mdist
	pars = object@mfilter$distr$z.D
	Sigma = rcov(object, type = 2)
	Mu	= object@mfilter$x.M
	A = object@mfilter$A
	support.method = support.method[1]
	
	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])
		}
	}
	debug = FALSE
	ans = switch(dist,
			manig = .convolve.nig(nigpars = pars, Sigma = Sigma, Mu = Mu, A = A, 
					weights = weights, fft.step = fft.step, fft.by = fft.by, fft.support = fft.support, 
					use.ff = use.ff, support.method = support.method, parallel = parallel, 
					parallel.control = parallel.control, debug = debug),
			magh = .convolve.gh(ghpars = pars, Sigma = Sigma, Mu = Mu, A = A, 
					weights = weights, fft.step = fft.step, fft.by = fft.by, fft.support = fft.support, 
					use.ff = use.ff, support.method = support.method, parallel = parallel, 
					parallel.control = parallel.control, debug = debug),
			mvnorm = .convolve.norm(Sigma  = Sigma, Mu = Mu, A = A, weights = weights))
	return(ans)
}

setMethod("convolution", signature(object = "goGARCHfilter"), .convolution.gogarchfilter)

.convolution.gogarchforecast = function(object, weights, fft.step = 0.001, fft.by = 0.0001, fft.support = c(-1, 1), 
		support.method = c("user", "adaptive"), use.ff = TRUE, parallel = FALSE, 
		parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2), debug = 0, ...)
{
	dist = object@mforecast$mdist
	pars = object@mforecast$distr$z.D
	Sigma = rcov(object, type = 2)
	Mu = object@mforecast$x.M
	A = object@mforecast$A
	support.method = support.method[1]
	
	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])
		}
	}
	debug = FALSE
	ans = switch(dist,
			manig = .convolve.nig(nigpars = pars, Sigma = Sigma, Mu = Mu, A = A, 
					weights = weights, fft.step = fft.step, fft.by = fft.by, fft.support = fft.support, 
					use.ff = use.ff, support.method = support.method, parallel = parallel, 
					parallel.control = parallel.control, debug = debug),
			magh = .convolve.gh(ghpars = pars, Sigma = Sigma, Mu = Mu, A = A, 
					weights = weights, fft.step = fft.step, fft.by = fft.by, fft.support = fft.support, 
					use.ff = use.ff, support.method = support.method, parallel = parallel, 
					parallel.control = parallel.control, debug = debug),
			mvnorm = .convolve.norm(Sigma  = Sigma, Mu = Mu, A = A, weights = weights))
	return(ans)
}

setMethod("convolution", signature(object = "goGARCHforecast"), .convolution.gogarchforecast)


.convolution.gogarchsim = function(object, weights, fft.step = 0.001, fft.by = 0.0001, fft.support = c(-1, 1), 
		support.method = c("user", "adaptive"), use.ff = TRUE, sim = 1, parallel = FALSE, 
		parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2), debug = 0, ...)
{
	dist = object@msim$mdist
	pars = object@msim$z.D
	Sigma = rcov(object, type = 2, sim = sim)
	Mu	= object@msim$seriesSim[[sim]]
	A = object@msim$A
	support.method = support.method[1]
	
	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])
		}
	}
	
	debug = FALSE
	ans = switch(dist,
			manig = .convolve.nig(nigpars = pars, Sigma = Sigma, Mu = Mu, A = A, 
					weights = weights, fft.step = fft.step, fft.by = fft.by, fft.support = fft.support, 
					use.ff = use.ff, support.method = support.method, parallel = parallel, 
					parallel.control = parallel.control, debug = debug),
			magh = .convolve.gh(ghpars = pars, Sigma = Sigma, Mu = Mu, A = A, 
					weights = weights, fft.step = fft.step, fft.by = fft.by, fft.support = fft.support, 
					use.ff = use.ff, support.method = support.method, parallel = parallel, 
					parallel.control = parallel.control, debug = debug),
			mvnorm = .convolve.norm(Sigma  = Sigma, Mu = Mu, A = A, weights = weights))
	return(ans)
}

setMethod("convolution", signature(object = "goGARCHsim"), .convolution.gogarchsim)

.convolution.gogarchroll = function(object, weights, fft.step = 0.001, fft.by = 0.0001, fft.support = c(-1, 1), 
		support.method = c("user", "adaptive"), use.ff = TRUE, parallel = FALSE, 
		parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2), debug = 0, ...)
{
	dist = object@spec$gospec@mspec$dist
	support.method = support.method[1]
	
	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])
		}
	}
	ans = switch(dist,
			mvnorm = .convolve.roll.mn(object, weights, ...),
			manig = .convolve.roll.nig(object, weights, fft.step, fft.by, fft.support, support.method, use.ff, 
					parallel, parallel.control, debug, ...),
			magh = .convolve.roll.nig(object, weights, fft.step, fft.by, fft.support, support.method, use.ff, 
					parallel, parallel.control, debug, ...))
	return(ans)
}

.convolve.roll.mn = function(object, weights, ...)
{
	forecast.length = object@spec$forecast.length
	refit.every = object@spec$refit.every
	forecast.length = floor(forecast.length/refit.every) * refit.every
	m = dim(object@forecast[[1]]@mforecast$x.M)[2]
	
	if(is.matrix(weights)){
		if(dim(weights)[2] != m) stop("\nconvolution-->error: Weights column dimension not equal to data column dimension!")
		if(dim(weights)[1] != forecast.length) stop("\nconvolution-->error: Weights row dimension must equal floor(forecast.length/refit.every) * refit.every")
	} else{
		weights = matrix(weights, ncol = m, nrow = forecast.length, byrow = TRUE)
	}
	nf = length(object@forecast)
	idx = matrix(NA, ncol = 2, nrow = nf)
	idx[,1] = seq(1, forecast.length, by = refit.every)
	idx[,2] = seq(refit.every, forecast.length, by = refit.every)
	mi = 0
	rpdf = matrix(NA, ncol = 2, nrow = forecast.length)
	for(i in 1:nf){
		tmp = .convolution.gogarchforecast(object@forecast[[i]], 
				weights = weights[( (mi*refit.every) +1):(i*refit.every),,drop = FALSE])
		rpdf[idx[i,1]:idx[i,2], ] = tmp@dist$y
		mi = mi + 1
	}
	sol = list()
	sol$dist = "mvnorm"
	sol$y = rpdf
	sol$fft.support = NULL
	sol$fft.step = NULL
	sol$fft.by = NULL
	ans = new("goGARCHfft",
			dist = sol)
	return(ans)
}

.convolve.roll.nig = function(object, weights, fft.step, fft.by, fft.support, support.method = c("user", "adaptive"), 
		use.ff = TRUE, parallel = FALSE, parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2), debug = 0, ...)
{
	dist = object@spec$gospec@mspec$dist
	forecast.length = object@spec$forecast.length
	refit.every = object@spec$refit.every
	forecast.length = floor(forecast.length/refit.every) * refit.every
	m = dim(object@forecast[[1]]@mforecast$x.M)[2]
	
	if(is.matrix(weights)){
		if(dim(weights)[2] != m) stop("\nconvolution-->error: Weights column dimension not equal to data column dimension!")
		if(dim(weights)[1] != forecast.length) stop("\nconvolution-->error: Weights row dimension must equal floor(forecast.length/refit.every) * refit.every")
	} else{
		weights = matrix(weights, ncol = m, nrow = forecast.length, byrow = TRUE)
	}
	nf = length(object@forecast)
	
	idx = matrix(NA, ncol = 2, nrow = nf)
	idx[,1] = seq(1, forecast.length, by = refit.every)
	idx[,2] = seq(refit.every, forecast.length, by = refit.every)
	
	
	if( support.method == "user" ){
		zz = seq(fft.support[1], fft.support[2], by = fft.by)
		if(use.ff) rpdf = ff(vmode="double", dim=c(length(zz), forecast.length)) else rpdf = matrix(NA, nrow = length(zz), ncol = forecast.length)
		mi = 0
		for(i in 1:nf){
			tmp = .convolution.gogarchforecast(object@forecast[[i]], 
					weights = weights[( (mi*refit.every) +1):(i*refit.every),,drop = FALSE], 
					fft.step = fft.step, fft.by = fft.by, fft.support = fft.support, support.method = support.method, 
					use.ff = use.ff, parallel = parallel, parallel.control = parallel.control, debug = debug)
			
			rpdf[1:length(zz), idx[i,1]:idx[i,2]] = tmp@dist$y[1:length(zz), 1:refit.every]
			mi = mi + 1
		}
	} else{
		rpdf = NULL
		mi = 0
		for(i in 1:nf){
			tmp = .convolution.gogarchforecast(object@forecast[[i]], 
					weights = weights[( (mi*refit.every) +1):(i*refit.every),,drop = FALSE], 
					fft.step = fft.step, fft.by = fft.by, fft.support = fft.support, support.method = support.method, 
					use.ff = use.ff, parallel = parallel, parallel.control = parallel.control, debug = debug)
			
			rpdf = c(rpdf, tmp@dist$y)
			mi = mi + 1
		}
	}
	sol = list()
	sol$dist = dist
	sol$y = rpdf
	sol$support.method = support.method
	sol$fft.support = fft.support
	sol$fft.step = fft.step
	sol$fft.by = fft.by
	ans = new("goGARCHfft",
			dist = sol)
	return(ans)
}

setMethod("convolution", signature(object = "goGARCHroll"), .convolution.gogarchroll)


# manig convolution of independent densities
.convolve.nig = function(nigpars, Sigma, Mu, A, weights, fft.step = 0.001, fft.by = 0.0001, 
		fft.support = c(-1, 1), support.method = c("user", "adaptive"), use.ff = TRUE, parallel = FALSE, 
		parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2), debug = 0)
{
	sol = list()
	n = dim(Sigma)[3]
	m = dim(Sigma)[2]
	if(!is.matrix(weights)) weights = matrix(weights, ncol = m, nrow = n, byrow = TRUE)
	if(dim(weights)[2] != m) stop("\nconvolution-->error: weights matrix dimension [2] not of same length as Sigma dimension [2].", call. = FALSE)
	if(dim(weights)[1] != n) stop("\nconvolution-->error: weights matrix dimension [1] not of same length as Sigma dimension [1].", call. = FALSE)
	
	# The weighting vector for the distribution
	w.hat = matrix(NA, ncol = m, nrow = n)
	for(i in 1:n){
		dS = sqrt(Sigma[,,i])
		w.hat[i,] = weights[i,] %*% (t(A) %*% dS)
	}
	# weights[i,] * (t(A) %*% diag(dS))
	# port = rep(0,n)
	# for(i in 1:n)  port[i] = weights[i,]%*%Mu[i,] + weights[i,]%*%(t(A) %*% sqrt(Sigma[,,i]) %*% (zres[i,]))
	# Mu[i,] + (t(A) %*% dS %*% (zres[i,]))
	wpars = array(data = NA, dim = c(m, 4, n))
	
	# Scaling of parameters (Blaesild proof)
	for(i in 1:n){

		for(j in 1:m){
			wpars[j,1:4,i] 	= nigpars[i,j,1:4]*c(1/abs(w.hat[i,j]),1/w.hat[i,j],abs(w.hat[i,j]),w.hat[i,j])
			wpars[j,4,i] 	= wpars[j,4,i] + Mu[i,j] * weights[i,j]
		}
	}
	
	if( support.method == "user" ){
		support.user = NULL
		zz = seq(fft.support[1], fft.support[2], by = fft.by)
		zzLength = length(zz)
		if(use.ff){
			nigpdf = ff(vmode="double", dim = c(zzLength, n))
		} else{
			nigpdf = matrix(NA, nrow = zzLength, ncol = n)
		}
		if( parallel ){
			if( parallel.control$pkg == "multicore" ){
				if(!exists("mclapply")){
					require('multicore')
				}
				xpdf = mclapply(1:n, FUN = function(i) cfinv(z = zz, f = nigmvcf, step = fft.step, 
									alpha = wpars[,1,i], beta = wpars[,2,i], delta = wpars[,3,i], mu = wpars[,4,i]), 
						mc.cores = parallel.control$cores)
				for(i in 1:n) nigpdf[,i] = xpdf[[i]]
				rm(xpdf)
			} else{
				if(!exists("sfLapply")){
					require('snowfall')
				}
				sfInit(parallel=TRUE, cpus = parallel.control$cores)
				sfExport("zz", "rgarch:::nigmvcf", "fft.step", "wpars",local = TRUE)
				xpdf = sfLapply(1:n, fun = function(i) rgarch:::cfinv(z = zz, f = nigmvcf, step = fft.step, 
									alpha = wpars[,1,i], beta = wpars[,2,i], delta = wpars[,3,i], mu = wpars[,4,i]))
				for(i in 1:n) nigpdf[,i] = xpdf[[i]]
				rm(xpdf)
				sfStop()
			}
		} else{	
			for(i in 1:n){
				nigpdf[,i] = cfinv(z = zz, f = nigmvcf, step = fft.step, alpha = wpars[,1,i], beta = wpars[,2,i], 
						delta = wpars[,3,i], mu = wpars[,4,i])
				if(debug) print(i)
			}
		}
		
	} else{
		nigpdf = vector(mode = "list", length = n)
		support.user = matrix(NA, ncol = 2, nrow = n)
		if( parallel ){
			if( parallel.control$pkg == "multicore" ){
				if(!exists("mclapply")){
					require('multicore')
				}
				xsol = mclapply(1:n, FUN = function(i){
							xmin = min(apply(cbind(wpars[,1,i], wpars[,2,i], wpars[,3,i], wpars[,4,i]), 1, FUN = function(x) qnig(0.0000001, alpha = x[1], beta = x[2], delta = x[3], mu = x[4])))
							xmax = max(apply(cbind(wpars[,1,i], wpars[,2,i], wpars[,3,i], wpars[,4,i]), 1, FUN = function(x) qnig(1-0.0000001, alpha = x[1], beta = x[2], delta = x[3], mu = x[4])))
							zz = seq(xmin, xmax, by = fft.by)
							ans = cfinv(z = zz, f = nigmvcf, step = fft.step, 
									alpha = wpars[,1,i], beta = wpars[,2,i], delta = wpars[,3,i], mu = wpars[,4,i])
							sol = list(d = ans, support = c(xmin, xmax))
							if(debug) print(i)
							return(sol)
						}, mc.cores = parallel.control$cores)
				for(i in 1:n){
					nigpdf[[i]] = xsol[[i]]$d
					support.user[i, ] = xsol[[i]]$support
				}
			} else{
				if(!exists("sfLapply")){
					require('snowfall')
				}
				sfInit(parallel=TRUE, cpus = parallel.control$cores)
				sfExport("zz", "rgarch:::nigmvcf", "fft.step", "fft.by", "wpars", "rgarch:::qnig", local = TRUE)
				xsol = sfLapply(1:n, fun = function(i){
							xmin = min(apply(cbind(wpars[,1,i], wpars[,2,i], wpars[,3,i], wpars[,4,i]), 1, FUN = function(x) qnig(0.0000001, alpha = x[1], beta = x[2], delta = x[3], mu = x[4])))
							xmax = max(apply(cbind(wpars[,1,i], wpars[,2,i], wpars[,3,i], wpars[,4,i]), 1, FUN = function(x) qnig(1-0.0000001, alpha = x[1], beta = x[2], delta = x[3], mu = x[4])))
							zz = seq(xmin, xmax, by = fft.by)
							ans = rgarch:::cfinv(z = zz, f = nigmvcf, step = fft.step, alpha = wpars[,1,i], beta = wpars[,2,i], delta = wpars[,3,i], mu = wpars[,4,i])
							sol = list(d = ans, support = c(xmin, xmax))
							if(debug) print(i)
							return(sol)
						})
				sfStop()
				for(i in 1:n){
					nigpdf[[i]] = xsol[[i]]$d
					support.user[i, ] = xsol[[i]]$support
				}
			}
		} else{
			for(i in 1:n){
				xmin = min(apply(cbind(wpars[,1,i], wpars[,2,i], wpars[,3,i], wpars[,4,i]), 1, FUN = function(x) qnig(0.0000001, alpha = x[1], beta = x[2], delta = x[3], mu = x[4])))
				xmax = max(apply(cbind(wpars[,1,i], wpars[,2,i], wpars[,3,i], wpars[,4,i]), 1, FUN = function(x) qnig(1-0.0000001, alpha = x[1], beta = x[2], delta = x[3], mu = x[4])))
				zz = seq(xmin, xmax, by = fft.by)
				nigpdf[[i]] = cfinv(z = zz, f = nigmvcf, step = fft.step, 
						alpha = wpars[,1,i], beta = wpars[,2,i], delta = wpars[,3,i], mu = wpars[,4,i])
				support.user[i, ] = c(xmin, xmax)
				if(debug) print(i)
			}
		}
	}
	# i.e. for each data point(n points) there is a density of size zzLength
	#nigpdf = big.matrix(nrow = zzLength, ncol = n, type = "double", init = NA, dimnames = NULL)


	sol$dist = "manig"
	sol$y = nigpdf
	sol$support.method = support.method
	sol$support.user = support.user
	sol$fft.support = fft.support
	sol$fft.step = fft.step
	sol$fft.by = fft.by
	ans = new("goGARCHfft",
			dist = sol)
	return(ans)
}


.convolve.gh = function(ghpars, Sigma, Mu, A, weights, fft.step = 0.001, fft.by = 0.0001, 
		fft.support = c(-1, 1), support.method = c("user", "adaptive"), use.ff = TRUE, parallel = FALSE, 
		parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2), debug = 0)
{
	if(!exists("BesselK")) {
        	require('Bessel')
    }
	sol = list()
	n = dim(Sigma)[3]
	m = dim(Sigma)[2]
	if(!is.matrix(weights)) weights = matrix(weights, ncol = m, nrow = n, byrow = TRUE)
	if(dim(weights)[2] != m) stop("\nconvolution-->error: weights matrix dimension [2] not of same length as Sigma dimension [2].", call. = FALSE)
	if(dim(weights)[1] != n) stop("\nconvolution-->error: weights matrix dimension [1] not of same length as Sigma dimension [1].", call. = FALSE)
	
	# The weighting vector for the distribution
	w.hat = matrix(NA, ncol = m, nrow = n)
	for(i in 1:n){
		dS = sqrt(Sigma[,,i])
		w.hat[i,] = weights[i,] %*% (t(A) %*% dS)
	}
	wpars = array(data = NA, dim = c(m, 5, n))
	
	# Scaling of parameters (Blaesild proof)
	for(i in 1:n){
		for(j in 1:m){
			wpars[j,1:4,i] 	= ghpars[i,j,1:4]*c(1/abs(w.hat[i,j]),1/w.hat[i,j],abs(w.hat[i,j]),w.hat[i,j])
			wpars[j,4,i] 	= wpars[j,4,i] + Mu[i,j] * weights[i,j]
			wpars[j,5,i]	= ghpars[i, j, 5]
		}
	}
	
	
	if( support.method == "user" ){
		support.user = NULL
		zz = seq(fft.support[1], fft.support[2], by = fft.by)
		zzLength = length(zz)
		if(use.ff){
			ghpdf = ff(vmode="double", dim = c(zzLength, n))
		} else{
			ghpdf = matrix(NA, nrow = zzLength, ncol = n)
		}
		if( parallel ){
			if( parallel.control$pkg == "multicore" ){
				if(!exists("mclapply")){
					require('multicore')
				}
				xpdf = mclapply(1:n, FUN = function(i) cfinv(z = zz, f = ghypmvcf, step = fft.step, lambda = wpars[,5,i], 
									alpha = wpars[,1,i], beta = wpars[,2,i], delta = wpars[,3,i], mu = wpars[,4,i]), 
						mc.cores = parallel.control$cores)
				for(i in 1:n) ghpdf[,i] = xpdf[[i]]
				rm(xpdf)
			} else{
				if(!exists("sfLapply")){
					require('snowfall')
				}
				sfInit(parallel=TRUE, cpus = parallel.control$cores)
				sfExport("zz", "rgarch:::ghypmvcf", "fft.step", "wpars",local = TRUE)
				xpdf = sfLapply(1:n, fun = function(i) rgarch:::cfinv(z = zz, f = ghypmvcf, step = fft.step, lambda = wpars[,5,i], 
									alpha = wpars[,1,i], beta = wpars[,2,i], delta = wpars[,3,i], mu = wpars[,4,i]))
				for(i in 1:n) ghpdf[,i] = xpdf[[i]]
				rm(xpdf)
				sfStop()
			}
		} else{	
			for(i in 1:n){
				ghpdf[,i] = cfinv(z = zz, f = ghypmvcf, step = fft.step, lambda = wpars[,5,i], 
						alpha = wpars[,1,i], beta = wpars[,2,i], delta = wpars[,3,i], mu = wpars[,4,i])
				if(debug) print(i)
			}
		}
		
	} else{
		ghpdf = vector(mode = "list", length = n)
		support.user = matrix(NA, ncol = 2, nrow = n)
		if( parallel ){
			if( parallel.control$pkg == "multicore" ){
				if(!exists("mclapply")){
					require('multicore')
				}
				xsol = mclapply(1:n, FUN = function(i){
							xmin = min(apply(cbind(wpars[,1,i], wpars[,2,i], wpars[,3,i], wpars[,4,i], wpars[,5,i]), 1, FUN = function(x) qgh(0.0000001, alpha = x[1], beta = x[2], delta = x[3], mu = x[4], lambda = x[5])))
							xmax = max(apply(cbind(wpars[,1,i], wpars[,2,i], wpars[,3,i], wpars[,4,i], wpars[,5,i]), 1, FUN = function(x) qgh(1-0.0000001, alpha = x[1], beta = x[2], delta = x[3], mu = x[4], lambda = x[5])))
							zz = seq(xmin, xmax, by = fft.by)
							ans = cfinv(z = zz, f = ghypmvcf, step = fft.step, lambda = wpars[,5,i], 
									alpha = wpars[,1,i], beta = wpars[,2,i], delta = wpars[,3,i], mu = wpars[,4,i])
							sol = list(d = ans, support = c(xmin, xmax))
							if(debug) print(i)
							return(sol)
						}, mc.cores = parallel.control$cores)
				for(i in 1:n){
					ghpdf[[i]] = xsol[[i]]$d
					support.user[i, ] = xsol[[i]]$support
				}
			} else{
				if(!exists("sfLapply")){
					require('snowfall')
				}
				sfInit(parallel=TRUE, cpus = parallel.control$cores)
				sfExport("zz", "rgarch:::ghypmvcf", "fft.step", "fft.by", "wpars", "rgarch:::qgh", local = TRUE)
				xsol = sfLapply(1:n, fun = function(i){
							xmin = min(apply(cbind(wpars[,1,i], wpars[,2,i], wpars[,3,i], wpars[,4,i], wpars[,5,i]), 1, FUN = function(x) qgh(0.0000001, alpha = x[1], beta = x[2], delta = x[3], mu = x[4], lambda = x[5])))
							xmax = max(apply(cbind(wpars[,1,i], wpars[,2,i], wpars[,3,i], wpars[,4,i], wpars[,5,i]), 1, FUN = function(x) qgh(1-0.0000001, alpha = x[1], beta = x[2], delta = x[3], mu = x[4], lambda = x[5])))
							zz = seq(xmin, xmax, by = fft.by)
							ans = rgarch:::cfinv(z = zz, f = ghypmvcf, step = fft.step, lambda = wpars[,5,i], 
									alpha = wpars[,1,i], beta = wpars[,2,i], delta = wpars[,3,i], mu = wpars[,4,i])
							sol = list(d = ans, support = c(xmin, xmax))
							if(debug) print(i)
							return(sol)
						})
				sfStop()
				for(i in 1:n){
					ghpdf[[i]] = xsol[[i]]$d
					support.user[i, ] = xsol[[i]]$support
				}
			}
		} else{
			for(i in 1:n){
				xmin = min(apply(cbind(wpars[,1,i], wpars[,2,i], wpars[,3,i], wpars[,4,i], wpars[,5,i]), 1, FUN = function(x) qgh(0.0000001, alpha = x[1], beta = x[2], delta = x[3], mu = x[4], lambda = x[5])))
				xmax = max(apply(cbind(wpars[,1,i], wpars[,2,i], wpars[,3,i], wpars[,4,i], wpars[,5,i]), 1, FUN = function(x) qgh(1-0.0000001, alpha = x[1], beta = x[2], delta = x[3], mu = x[4], lambda = x[5])))
				zz = seq(xmin, xmax, by = fft.by)
				ghpdf[[i]] = cfinv(z = zz, f = ghypmvcf, step = fft.step, lambda = wpars[,5,i], 
						alpha = wpars[,1,i], beta = wpars[,2,i], delta = wpars[,3,i], mu = wpars[,4,i])
				support.user[i, ] = c(xmin, xmax)
				if(debug) print(i)
			}
		}
	}
	sol$dist = "magh"
	sol$y = ghpdf
	sol$fft.support = fft.support
	sol$support.method = support.method
	sol$support.user = support.user
	sol$fft.step = fft.step
	sol$fft.by = fft.by
	ans = new("goGARCHfft",
			dist = sol)
	return(ans)
}

.convolve.norm = function(Sigma, Mu, A, weights, debug = FALSE)
{
	sol = list()
	n = dim(Sigma)[3]
	m = dim(Sigma)[2]
	if(!is.matrix(weights)) weights = matrix(weights, ncol = m, nrow = n, byrow = TRUE)
	if(dim(weights)[2] != m) stop("\nconvolution-->error: weights matrix dimension [2] not of same length as Sigma dimension [2].", call. = FALSE)
	if(dim(weights)[1] != n) stop("\nconvolution-->error: weights matrix dimension [1] not of same length as Sigma dimension [1].", call. = FALSE)
	
	# A matrix with the weighte portfolio sigma and mu
	normpdf = matrix(NA, ncol = 2, nrow = n)
	w = weights
	for(i in 1:n){
		normpdf[i, 1] = w[i,]%*%Mu[i,]
		normpdf[i, 2] = t(w[i,])%*%(t(A)%*%Sigma[,,i]%*%A)%*%w[i,]
		if(debug) print(i)
	}
	sol$dist = "mvnorm"
	sol$y = normpdf
	sol$fft.support = NULL
	sol$fft.step = NULL
	sol$fft.by = NULL
	ans = new("goGARCHfft",
			dist = sol)
	return(ans)
}
#------------------------------------------------------------------------------
# interpolation of d, p, q portfolio density
dfft = function(object, index = 1)
{
	UseMethod("dfft")
}

setMethod("dfft", signature(object = "goGARCHfft"), .dfft)

pfft = function(object, index = 1)
{
	UseMethod("pfft")
}

setMethod("pfft", signature(object = "goGARCHfft"), .pfft)


qfft = function(object, index = 1)
{
	UseMethod("qfft")
}

setMethod("qfft", signature(object = "goGARCHfft"), .qfft)

#------------------------------------------------------------------------------
###############################################################################

nisurface = function(object, ...)
{
	UseMethod("nisurface")
}

.newsimpact.gogarch = function(object, type = "cov", pair = c(1,2), factor = c(1,2), plot = TRUE)
{
	ans = switch(type,
			cov = .newsimpact.gogarch.cov(object = object, pair = pair, plot = plot),
			cor = .newsimpact.gogarch.cor(object = object, pair = pair, plot = plot)
			# coskew = .newsimpact.gogarch.coskew(object = object, pair = pair, plot = plot),
			# cokurt = .newsimpact.gogarch.cokurt(object = object, pair = pair, plot = plot)
			)
	return(ans)
}

setMethod("nisurface", signature(object = "goGARCHfit"), .newsimpact.gogarch)

.newsimpact.gogarch.cov = function(object, pair = c(1,2), factor = c(1,2), plot = TRUE)
{
	A = object@mfit$A
	m = dim(A)[1]
	nilist = vector(mode = "list", length = m)
	nilist[[1]] = newsimpact(z = NULL, object = object@ufit@fit[[1]])
	# we want a common epsilon
	zeps = nilist[[1]]$zx
	for(i in 2:m) nilist[[i]]  = newsimpact(z = zeps, object = object@ufit@fit[[i]])
	sigmas = sapply(nilist, FUN = function(x) x$zy)
	N = dim(sigmas)[1]
	H = array(NA, dim = c(dim(A)[1], dim(A)[1], N))
	nicurve = matrix(NA, ncol = N, nrow = N)
	zr = as.numeric(which(nilist[[1]]$zx == 0))
	# to get the off diagonals we need to keep 1 fixed and
	# revolve the epsilon of the other
	s2 = sapply(nilist, FUN = function(x) x$zy)
	rownames(s2) = round(as.numeric(zeps),4)
	np = 1:m
	for(j in 1:N){
		for(i in 1:N){	
			sigmas = s2[zr,]
			sigmas1 = as.numeric(s2[j, factor[1]])
			sigmas2 = as.numeric(s2[i, factor[2]])
			sigmas[factor[1]] = sigmas1
			sigmas[factor[2]] = sigmas2
			H = t(A)%*%diag(sigmas)%*%A
			nicurve[j,i] = H[pair[1], pair[2]]
		}
	}
	if(plot){
		tcol <- terrain.colors(12)
		par(mfrow = c(1,2))
		persp(  x = zeps,
			y = zeps,
			z = nicurve,  col = "lightblue", theta = 30, phi = 30,
			ltheta = 120, shade = 0.75, ticktype = "detailed" )	
		contour(  x = zeps,
			y = zeps,
			z = nicurve,  col = tcol[2], lty = "solid")
	}
	return(list(nicurve = nicurve, axis = zeps))
}


.newsimpact.gogarch.cor = function(object, pair = c(1,2), factor = c(1,2), plot = TRUE)
{
	A = object@mfit$A
	m = dim(A)[1]
	nilist = vector(mode = "list", length = m)
	nilist[[1]] = newsimpact(z = NULL, object = object@ufit@fit[[1]])
	zeps = nilist[[1]]$zx
	for(i in 2:m) nilist[[i]]  = newsimpact(z = zeps, object = object@ufit@fit[[i]])
	sigmas = sapply(nilist, FUN = function(x) x$zy)
	N = dim(sigmas)[1]
	H = array(NA, dim = c(dim(A)[1], dim(A)[1], N))
	nicurve = matrix(NA, ncol = N, nrow = N)
	zr = as.numeric(which(nilist[[1]]$zx == 0))
	s2 = sapply(nilist, FUN = function(x) x$zy)
	rownames(s2) = round(as.numeric(zeps),4)
	np = 1:m
	for(j in 1:N){
		for(i in 1:N){	
			sigmas = s2[zr,]
			sigmas1 = as.numeric(s2[j, factor[1]])
			sigmas2 = as.numeric(s2[i, factor[2]])
			sigmas[factor[1]] = sigmas1
			sigmas[factor[2]] = sigmas2
			H = cov2cor(t(A)%*%diag(sigmas)%*%A)
			nicurve[j,i] = H[pair[1], pair[2]]
		}
	}
	if(plot){
		tcol <- terrain.colors(12)
		par(mfrow = c(1,2))
		persp(  x = zeps,
				y = zeps,
				z = nicurve,  col = "lightblue", theta = 30, phi = 30,
				ltheta = 120, shade = 0.75, ticktype = "detailed" )	
		contour(  x = zeps,
				y = zeps,
				z = nicurve,  col = tcol[2], lty = "solid")
	}
	return(list(nicurve = nicurve, axis = zeps))
}

#-----------------------------------------------------------------------------

Try the rgarch package in your browser

Any scripts or data that you put into this service are public.

rgarch documentation built on May 2, 2019, 5:22 p.m.