Nothing
#################################################################################
##
## 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))
}
#-----------------------------------------------------------------------------
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.