Nothing
R/copula-garch.R
In rgarch: Flexible GARCH modelling in R
#################################################################################
##
## R package rgarch by Alexios Ghalanos Copyright (C) 2008, 2009, 2010
## 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.
##
#################################################################################
#################################################################################
# Fit, Filter, Forecast and Simulation
#-------------------------------------
.cgarchspec = function(uspec, VAR = FALSE, VAR.opt = list(robust = FALSE, lag = 1, lag.max = NULL,
lag.criterion = c("AIC", "HQ", "SC", "FPE"), external.regressors = NULL,
robust.control = list("gamma" = 0.25, "delta" = 0.01, "nc" = 10, "ns" = 500)),
dccOrder = c(1,1), distribution.model = list(copula = c("mvnorm", "mvt"), method = c("Kendall", "ML"),
time.varying = FALSE, transformation = c("parametric", "empirical", "spd")),
start.pars = list(), fixed.pars = list())
{
m = length(uspec@spec)
if(is.null(distribution.model$copula)) distribution = "mvnorm" else distribution = tolower(distribution.model$copula)
distribution = distribution[1]
valid.distributions = c("mvnorm", "mvt")
if(!any(distribution == valid.distributions)) stop("\nInvalid Copula Distribution Choice\n", call. = FALSE)
if(is.null(distribution.model$method)) method = "kendall" else method = tolower(distribution.model$method)
method = method[1]
valid.methods = c("kendall", "ml")
if(!any(method == valid.methods)){
if(distribution == "mvt") warning("\nInvalid Rho Method Estimation Choice\n", call. = FALSE)
method = "kendall"
}
if(is.null(distribution.model$time.varying)) timecopula = FALSE else timecopula = as.logical(distribution.model$time.varying)
timecopula = timecopula[1]
if(is.null(distribution.model$transformation)) transformation = "parametric" else transformation = tolower(distribution.model$transformation)
transformation = transformation[1]
valid.transformations = c("parametric", "empirical", "spd")
if(!any(transformation == valid.transformations)) stop("\nInvalid Copula Transformation Choice\n", call. = FALSE)
if(is.null(dccOrder) && timecopula) dccOrder = c(1,1)
if(!timecopula) dccOrder = c(0,0)
if( length(dccOrder) != 2) stop("\nInvalid dccOrder\n", call. = FALSE)
dccOrder = as.integer(dccOrder)
maxdccOrder = max(dccOrder)
maxgarchOrder = max( sapply(uspec@spec, FUN = function(x) max(max(x@mean.model$armaOrder), max(x@variance.model$garchOrder) ) ) )
maxOrder = max( maxdccOrder, maxgarchOrder )
shape.distributions = c("mvt")
skew.distributions = NULL
include.skew = include.shape = FALSE
# we will have to change this eventually to accomodate vector
# skewness and shape parameters for more interesting distributions
if(any(distribution == shape.distributions)) include.shape = TRUE
if(any(distribution == skew.distributions)) include.skew = TRUE
mspec = list()
mspec$timecopula = timecopula
mspec$method = method
mspec$transformation = transformation
mspec$dccOrder = dccOrder
mspec$start.pars = start.pars
mspec$fixed.pars = fixed.pars
mspec$distribution = distribution
mspec$include.shape = include.shape
mspec$include.skew = include.skew
modelnames = c(
if(dccOrder[1] > 0 && timecopula) paste("dcca", 1:dccOrder[1], sep = ""),
if(dccOrder[2] > 0 && timecopula) paste("dccb", 1:dccOrder[2], sep = ""),
if(include.skew) paste("skew"),
if(include.shape) paste("shape"),
if(!timecopula && method == "ml") paste("Cov", 1:((m*m - m)/2), sep = ""))
pos = 1
pos.matrix = matrix(0, ncol = 4, nrow = 5)
colnames(pos.matrix) = c("start", "stop", "include", "npars")
rownames(pos.matrix) = c("dcca","dccb","skew","shape", "Cov")
if(dccOrder[1] > 0 && timecopula){
pos.matrix[1,1:3] = c(pos, pos+dccOrder[1]-1, 1)
pos = max(pos.matrix[1:1, 2]) + 1
}
if(dccOrder[2] > 0 && timecopula){
pos.matrix[2,1:3] = c(pos, pos+dccOrder[2]-1, 1)
pos = max(pos.matrix[1:2, 2]) + 1
}
if(include.skew)
{
pos.matrix[3,1:3] = c(pos, pos, 1)
pos = max(pos.matrix[1:3, 2]) + 1
}
if(include.shape)
{
pos.matrix[4,1:3] = c(pos, pos, 1)
pos = max(pos.matrix[1:4, 2]) + 1
}
if(!timecopula && method == "ml"){
pos.matrix[5,1:3] = c(pos, pos + ((m*m - m)/2) - 1, 1)
}
nn = length(modelnames)
modelmatrix = matrix(0, ncol = 3, nrow = nn)
rownames(modelmatrix) = modelnames
colnames(modelmatrix) = c("opt", "fixed", "start")
fixed.names = names(fixed.pars)
fp = charmatch(fixed.names, modelnames)
if(!is.null(fixed.names) && any(!is.na(fp))){
fixed = fp[!is.na(fp)]
modelmatrix[fixed,2] = 1
fz = charmatch(modelnames, fixed.names)
fz = fz[!is.na(fz)]
fixed.pars = fixed.pars[fz]
names(fixed.pars) = fixed.names[fz]
} else{
fixed.pars = NULL
}
modelmatrix[,1] = 1 - modelmatrix[,2]
start.names = names(start.pars)
sp = charmatch(start.names, modelnames)
if(!is.null(start.names) && any(!is.na(sp))){
start = sp[!is.na(sp)]
modelmatrix[start,3] = 1
sz = charmatch(modelnames, start.names)
sz = sz[!is.na(sz)]
start.pars = start.pars[sz]
} else{
start.pars = NULL
}
opt.pars = vector(mode = "numeric", length = length(modelnames))
names(opt.pars) = modelnames
# optimization model
omodel = list(modelnames = modelnames, opt.pars = opt.pars, start.pars = start.pars,
fixed.pars = fixed.pars, modelmatrix = modelmatrix, maxOrder = maxOrder,
maxgarchOrder = maxgarchOrder, maxdccOrder = maxdccOrder,
pos.matrix = pos.matrix)
mspec$optimization.model = omodel
if( !is(uspec, "uGARCHmultispec") ) stop("\ncgarchspec-->error: uspec must be a uGARCHmultispec object")
mspec$varmodel = list()
if( VAR ){
mspec$varmodel$VAR = TRUE
mspec$varmodel$robust = VAR.opt$robust
mspec$varmodel$lag = VAR.opt$lag
mspec$varmodel$lag.max = VAR.opt$lag.max
mspec$varmodel$lag.criterion = VAR.opt$lag.criterion
mspec$varmodel$external.regressors = VAR.opt$external.regressors
mspec$varmodel$robust.control = VAR.opt$robust.control
if( !is.null(mspec$varmodel$external.regressors) ){
mspec$varmodel$mxn = dim( as.matrix(VAR.opt$external.regressors) )[2]
} else{
mspec$varmodel$mxn = 0
}
# loop and reconstruct the spec to exclude a mean specification
nl = length(uspec@spec)
tmpspec = vector(mode ="list", length = nl)
for(i in 1:nl){
tmpspec[[i]] = ugarchspec(mean.model = list(include.mean = FALSE, armaOrder = c(0,0),
external.regressors = NULL, arfima = FALSE, garchInMean = FALSE),
variance.model = list(model = uspec@spec[[i]]@variance.model$model,
submodel = uspec@spec[[i]]@variance.model$submodel,
garchOrder = uspec@spec[[i]]@variance.model$garchOrder,
external.regressors = uspec@spec[[i]]@variance.model$vexdata,
variance.targeting = !uspec@spec[[i]]@variance.model$include.omega),
distribution.model = uspec@spec[[i]]@distribution.model$distribution,
start.pars = uspec@spec[[i]]@optimization.model$start.pars,
fixed.pars = uspec@spec[[i]]@optimization.model$fixed.pars)
}
uspec = multispec( tmpspec )
} else{
mspec$varmodel$VAR = FALSE
mspec$varmodel$lag = 1
mspec$varmodel$lag.max = 1
mspec$varmodel$lag.criterion = "HQ"
mspec$varmodel$external.regressors = NULL
mspec$varmodel$mxn = 0
}
ans = new("cGARCHspec",
mspec = mspec,
uspec = uspec)
return(ans)
}
setMethod("cgarchspec", signature(uspec = "uGARCHmultispec"), .cgarchspec)
.cgarchfit = function(spec, data, spd.control = list(lower = 0.1, upper = 0.9, type = "pwm", kernel = "epanech"),
fit.control = list(eval.se = TRUE, trace = TRUE, stationarity = TRUE), solver = "solnp", solver.control = list(), out.sample = 0, parallel = FALSE,
parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2), fit = NULL, VAR.fit = NULL, ...)
{
dist = spec@mspec$distribution
tv = ifelse(spec@mspec$timecopula, 1, 0)
dx = paste(dist, tv, sep = "")
ans = switch(dx,
mvt0 = .cgarchfit.student(spec = spec, data = data, spd.control = spd.control, fit.control = fit.control, solver = solver,
solver.control = solver.control, out.sample = out.sample, parallel = parallel,
parallel.control = parallel.control, fit = fit, VAR.fit = VAR.fit, ...),
mvt1 = .cgarchfit.tvstudent(spec = spec, data = data, spd.control = spd.control, fit.control = fit.control, solver = solver,
solver.control = solver.control, out.sample = out.sample, parallel = parallel,
parallel.control = parallel.control, fit = fit, VAR.fit = VAR.fit, ...),
mvnorm0 = .cgarchfit.norm(spec = spec, data = data, spd.control = spd.control, fit.control = fit.control, solver = solver,
solver.control = solver.control, out.sample = out.sample, parallel = parallel,
parallel.control = parallel.control, fit = fit, VAR.fit = VAR.fit, ...),
mvnorm1 = .cgarchfit.tvnorm(spec = spec, data = data, spd.control = spd.control, fit.control = fit.control, solver = solver,
solver.control = solver.control, out.sample = out.sample, parallel = parallel,
parallel.control = parallel.control, fit = fit, VAR.fit = VAR.fit, ...))
return( ans )
}
setMethod("cgarchfit", signature(spec = "cGARCHspec"), .cgarchfit)
.cgarchfit.tvnorm = function(spec, data, spd.control = list(lower = 0.1, upper = 0.9, type = "pwm", kernel = "epanech"),
fit.control = list(eval.se = TRUE, trace = TRUE, stationarity = TRUE), solver = "solnp", solver.control = list(), out.sample = 0, parallel = FALSE,
parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2), fit = NULL, VAR.fit = NULL, ...)
{
tic = Sys.time()
ufit.control = list()
if(is.null(fit.control$stationarity)){
ufit.control$stationarity = TRUE
fit.control$stationarity = TRUE
} else {
ufit.control$stationarity = fit.control$stationarity
fit.control$stationarity = fit.control$stationarity
}
if(is.null(fit.control$scale)){
ufit.control$scale = TRUE
} else{
ufit.control$scale = fit.control$scale
fit.control$scale = NULL
}
if(is.null(fit.control$eval.se)){
eval.se = TRUE
} else{
eval.se = as.logical(fit.control$eval.se)
}
# Parallel Execution Prelim Checks (we need them because of the hessian estimation)
if( parallel ){
os = .Platform$OS.type
if(is.null(parallel.control$pkg)){
if( os == "windows" ) parallel.control$pkg = "snowfall" else parallel.control$pkg = "multicore"
if( is.null(parallel.control$cores) ) parallel.control$cores = 2
} else{
mtype = match(tolower(parallel.control$pkg[1]), c("multicore", "snowfall"))
if(is.na(mtype)) stop("\nParallel Package type not recognized in parallel.control\n")
parallel.control$pkg = tolower(parallel.control$pkg[1])
if( os == "windows" && parallel.control$pkg == "multicore" ) stop("\nmulticore not supported on windows O/S\n")
if( is.null(parallel.control$cores) ) parallel.control$cores = 2 else parallel.control$cores = as.integer(parallel.control$cores[1])
}
}
if( is.null( colnames(data) ) ) cnames = paste("Asset_", 1:k, sep = "") else cnames = colnames(data)
xdata = .extractmdata(data)
if(!is.numeric(out.sample)) stop("\ncgarchfit-->error: out.sample must be numeric\n")
if(as.numeric(out.sample) < 0) stop("\ncgarchfit-->error: out.sample must be positive\n")
n.start = round(out.sample, 0)
n = dim(xdata$data)[1]
k = dim(xdata$data)[2]
if( (n-n.start) < 100) stop("\ncgarchfit-->error: function requires at least 100 data\n points to run\n")
data = xdata$data[1:(n-n.start), , drop = FALSE]
dates = xdata$pos[1:(n-n.start), drop = FALSE]
origdata = xdata$data
origdates = xdata$pos
dformat = xdata$dformat
vrmodel = NULL
if( spec@mspec$varmodel$VAR ){
if(spec@mspec$varmodel$mxn>0){
ex = spec@mspec$varmodel$external.regressors
if( dim(ex)[1] < dim(data)[1] ) stop("\ncgarchfit-->error: external regressor matrix has less points than data!...", call. = FALSE)
orex = ex
ex = ex[1:(n - n.start), , drop = FALSE]
mxn = dim(orex)[2]
} else{
orex = NULL
ex = NULL
mxn = 0
}
if( !is.null(VAR.fit) ){
zdata = VAR.fit$xresiduals
# should be N - p
p = mp = VAR.fit$lag
N = dim(zdata)[1]
K = dim(zdata)[2]
mu = VAR.fit$xfitted
vrmodel = VAR.fit
#spec@mspec$optimization.model$maxOrder = max( mp, spec@mspec$optimization.model$maxdccOrder )
#spec@mspec$optimization.model$maxgarchOrder = mp
} else{
cat("\nestimating VAR model...")
ic = spec@mspec$varmodel$lag.criterion[1]
if( !is.null(spec@mspec$varmodel$lag.max) ){
mp = .varxselect(y = data, lag.max = spec@mspec$varmodel$lag.max, exogen = ex)$selection
mp = as.numeric( mp[which(substr(names(mp), 1, 2) == substr(ic, 1, 2))] )
} else {
mp = spec@mspec$varmodel$lag
}
vrmodel = varxfilter(X = data, p = mp, exogen = ex, Bcoef = NULL, robust = spec@mspec$varmodel$robust,
gamma = spec@mspec$varmodel$robust.control$gamma, delta = spec@mspec$varmodel$robust.control$delta,
nc = spec@mspec$varmodel$robust.control$nc, ns = spec@mspec$varmodel$robust.control$ns)
cat("done!\n")
zdata = vrmodel$xresiduals
# should be N - p
N = dim(zdata)[1]
K = dim(zdata)[2]
mu = vrmodel$xfitted
vrmodel$lag = mp
vrmodel$mxn = mxn
#spec@mspec$optimization.model$maxOrder = max( mp, spec@mspec$optimization.model$maxdccOrder )
#spec@mspec$optimization.model$maxgarchOrder = mp
}
# if we have out of sample data we need to filter
if(out.sample>0){
Bcoef = vrmodel$Bcoef
vrmodel2 = varxfilter(X = origdata, p = mp, exogen = orex, Bcoef = Bcoef)
zorigdata = vrmodel2$xresiduals
} else{
zorigdata = zdata
}
} else{
zdata = data
zorigdata = origdata
orex = NULL
ex = NULL
}
if( is(spec@uspec, "uGARCHmultispec") ){
speclist = spec@uspec
} else{
stop("\ncgarchfit-->error: cGARCHspec must contain multispec object with fixed parameters")
}
maxgarchOrder = spec@mspec$optimization.model$maxgarchOrder
tmpidx = .dccparidx( speclist )
paridx = tmpidx$idx
paridxi = tmpidx$idxi
paridxa = tmpidx$idxa
if( !is.null(fit) && is(fit, "uGARCHmultifit") ){
fitlist = fit
} else{
fitlist = multifit(multispec = speclist, data = zorigdata, out.sample = n.start, solver = solver,
solver.control = solver.control, fit.control = ufit.control, parallel = parallel,
parallel.control = parallel.control)
assign(".fitlist", fitlist, envir = .GlobalEnv)
}
garchpars = as.numeric( unlist(coef(fitlist) ) )
converge = sapply(fitlist@fit, FUN = function(x) x@fit$convergence)
if( any( converge == 1 ) ){
pr = which(converge != 1)
cat("\nNon-Converged:\n")
print(pr)
cat("\ncgarchfit-->error: convergence problem in univariate fit...")
cat("\n...returning uGARCHmultifit object instead...check and resubmit...")
return( fitlist )
}
# create a temporary environment to store values (deleted at end of function)
tempenvir = new.env(hash = TRUE, parent = .GlobalEnv)
if(is.null(fit.control$trace)) trace = FALSE else trace = as.logical(fit.control$trace)
assign("trace", trace, envir = tempenvir)
assign("eval.se", eval.se, envir = tempenvir)
assign("cnames", cnames, envir = tempenvir)
assign("spec", spec, envir = tempenvir)
assign("parallel", parallel, envir = tempenvir)
assign("parallel.control", parallel.control, envir = tempenvir)
assign("paridx", paridx, envir = tempenvir)
assign("paridxi", paridxi, envir = tempenvir)
assign("paridxa", paridxa, envir = tempenvir)
assign("speclist", speclist, envir = tempenvir)
assign("fitlist", fitlist, envir = tempenvir)
omodel = spec@mspec$optimization.model
distribution = spec@mspec$distribution
transformation = spec@mspec$transformation
sig = sigma(fitlist)
res = residuals(fitlist)
#stdresid = res/sig
zres = sapply(fitlist@fit, FUN = function(x) x@fit$z)
ures = matrix(NA, ncol = 5, nrow = dim(zres)[1])
if( maxgarchOrder > 0 ) zres = zres[-(1:maxgarchOrder), , drop = FALSE]
# Transformation of z to [0,1] domain
ans = switch(transformation,
parametric = .pparametric(fitlist, zres),
empirical = .pempirical(zres),
spd = .pspd(zres, spd.control))
if(transformation == "spd"){
ures = ans$ures
sfit = ans$sfit
} else{
ures = ans
sfit = NULL
}
assign("spd.control", spd.control, envir = tempenvir)
# make tail adjustments in order to avoid problem with the quantile functions in
# optimization
if(any(ures > 0.99999)){
xn = which(ures > 0.99999)
ures[xn] = 0.99999
}
if(any(ures < eps)){
xn = which(ures < (1.5*eps))
ures[xn] = eps
}
assign("omodel", omodel, envir = tempenvir)
assign("mspec", spec@mspec, envir = tempenvir)
assign("ures", ures, envir = tempenvir)
assign("solver", solver, envir = tempenvir)
assign("fit.control", fit.control, envir = tempenvir)
Ifn = .dcccon
ILB = 0.0001
IUB = 0.9998
if( solver == "solnp" ) fit.control$stationarity = FALSE else fit.control$stationarity = TRUE
assign("fit.control", fit.control, envir = tempenvir)
ipars = .tvcopulastart(data = ures, garchenv = tempenvir)
pars = ipars$pars
LB = ipars$modelLB
UB = ipars$modelUB
assign("npar", length(pars), envir = tempenvir)
if( any( LB == UB ) )
{
fixid = which(LB == UB)
if( length(fixid) == length(pars) )
{
if( !eval.se )
{
# if all parameters are fixed an no standard erros are to
# be calculated then we return a ugarchfilter object
cat("\ncgarchfit-->warning: all parameters fixed...returning cgarchfilter object instead\n")
return(cgarchfilter(spec = spec, data = ures, out.sample = out.sample,
parallel = parallel, parallel.control = parallel.control, VAR.fit = VAR.fit,
fit = fitlist))
} else{
# if all parameters are fixed but we require standard errors, we
# skip the solver
use.solver = FALSE
fixed = pars[fixid]
pars = pars[-fixid]
}
} else{
# with some parameters fixed we extract them (to be rejoined at end)
# so that they do not enter the solver
use.solver = TRUE
fixed = pars[fixid]
pars = pars[-fixid]
LB = LB[-fixid]
UB = UB[-fixid]
}
} else{
# no fixed parameters, all normal
fixed = NULL
fixid = NULL
use.solver = TRUE
}
assign("fixed", fixed, envir = tempenvir)
assign("fixid", fixid, envir = tempenvir)
assign(".llh", 1, envir = tempenvir)
# get
if( use.solver )
{
solution = .copulasolver(solver, pars, fun = copula.tvnormalLLH, Ifn, ILB,
IUB, gr = NULL, hessian = NULL, control = solver.control,
LB = LB, UB = UB, data = ures, returnType = "llh", garchenv = tempenvir)
sol = solution$sol
hess = solution$hess
opt.pars = sol$par
names(opt.pars) = names(pars)
convergence = sol$convergence
} else{
hess = NULL
timer = Sys.time()-tic
opt.pars = NULL
convergence = 0
sol = list()
sol$message = "all parameters fixed"
}
mfit = list()
mfit$dates = dates
mfit$date.format = dformat
mfit$origdata = origdata
mfit$origdates = origdates
mfit$vrmodel = vrmodel
# convert back from uniform to standard variates
model = list()
model$copula.distribution = distribution
model$out.sample = out.sample
model$transformation = transformation
model$spd.control = spd.control
model$sfit = sfit
model$time.varying = TRUE
# add some tail data of the univariate fits to be used with simulation (we assume that noone is
# going to add more than 25 lags in a univariate fit!
mfit$tailusigma = tail(sigma(fitlist), 25)
mfit$tailuresids = tail(residuals(fitlist), 25)
mfit$tailuret = tail(data, 25)
cfit = copula.tvnormalLLH(opt.pars, data = ures, returnType = "ALL", garchenv = tempenvir)
if(eval.se){
nderiv = .cgarchcvar(f = copula.tvnormalLLH, pars = opt.pars, garchfit = fitlist, model = model,
env = tempenvir, fname = "copula.tvnormalLLH")
mfit$pars = nderiv$pars
mfit$garchpars = coef(fitlist)
# stpars = second stage parameters
mfit$stpars = opt.pars
mfit$scores = nderiv$scores
mfit$A = nderiv$A
mfit$cvar = nderiv$cvar
mfit$matcoef = nderiv$matcoef
} else{
garchNames = vector(mode = "character")
assetNames = cnames
dccNames = spec@mspec$optimization.model$modelnames
m = length(spec@uspec@spec)
for(i in 1:m){
cf = coef(fitlist@fit[[i]])
cfnames = names(cf)
garchNames = c(garchNames, paste(assetNames[i], cfnames, sep = ".") )
}
allnames = c(garchNames, dccNames)
mfit$pars = c(as.numeric(unlist(coef(fitlist))), opt.pars)
mfit$garchpars = coef(fitlist)
mfit$stpars = opt.pars
names(mfit$pars) = allnames
mfit$scores = NULL
mfit$A = NULL
mfit$cvar = NULL
matcoef = matrix(NA, ncol = 4, nrow = length(mfit$pars))
matcoef[,1] = mfit$pars
dimnames(matcoef) = list(allnames, c(" Estimate",
" Std. Error", " t value", "Pr(>|t|)"))
mfit$matcoef = matcoef
}
mfit$paridx = tmpidx
garchllh = sapply(fitlist@fit, FUN = function(x) -x@fit$log.likelihoods)
copllh = cfit$lik
if( maxgarchOrder > 0 ) copllh = c(rep(0, maxgarchOrder), copllh)
mfit$llh = sum( apply(cbind(copllh, garchllh), 1, "sum") )
mfit$timer = Sys.time() - tic
mfit$zres = zres
ans = new("cGARCHfit",
mfit = mfit,
mspec = spec@mspec,
uspec = spec@uspec,
copula = cfit,
model = model)
return(ans)
}
.cgarchfit.norm = function(spec, data, spd.control = list(lower = 0.1, upper = 0.9, type = "pwm", kernel = "epanech"),
fit.control = list(eval.se = TRUE, trace = TRUE, stationarity = TRUE), solver = "solnp", solver.control = list(), out.sample = 0, parallel = FALSE,
parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2), fit = NULL, VAR.fit = NULL, ...)
{
tic = Sys.time()
ufit.control = list()
if(is.null(fit.control$stationarity)){
ufit.control$stationarity = TRUE
} else {
ufit.control$stationarity = fit.control$stationarity
fit.control$stationarity = NULL
}
if(is.null(fit.control$scale)){
ufit.control$scale = TRUE
} else{
ufit.control$scale = fit.control$scale
fit.control$scale = NULL
}
if(is.null(fit.control$eval.se)){
eval.se = TRUE
} else{
eval.se = as.logical(fit.control$eval.se)
}
# Parallel Execution Prelim Checks (we need them because of the hessian estimation)
if( parallel ){
os = .Platform$OS.type
if(is.null(parallel.control$pkg)){
if( os == "windows" ) parallel.control$pkg = "snowfall" else parallel.control$pkg = "multicore"
if( is.null(parallel.control$cores) ) parallel.control$cores = 2
} else{
mtype = match(tolower(parallel.control$pkg[1]), c("multicore", "snowfall"))
if(is.na(mtype)) stop("\nParallel Package type not recognized in parallel.control\n")
parallel.control$pkg = tolower(parallel.control$pkg[1])
if( os == "windows" && parallel.control$pkg == "multicore" ) stop("\nmulticore not supported on windows O/S\n")
if( is.null(parallel.control$cores) ) parallel.control$cores = 2 else parallel.control$cores = as.integer(parallel.control$cores[1])
}
}
if( is.null( colnames(data) ) ) cnames = paste("Asset_", 1:k, sep = "") else cnames = colnames(data)
xdata = .extractmdata(data)
if(!is.numeric(out.sample)) stop("\ncgarchfit-->error: out.sample must be numeric\n")
if(as.numeric(out.sample) < 0) stop("\ncgarchfit-->error: out.sample must be positive\n")
n.start = round(out.sample, 0)
n = dim(xdata$data)[1]
k = dim(xdata$data)[2]
if( (n-n.start) < 100) stop("\ncgarchfit-->error: function requires at least 100 data\n points to run\n")
data = xdata$data[1:(n-n.start), , drop = FALSE]
dates = xdata$pos[1:(n-n.start), drop = FALSE]
origdata = xdata$data
origdates = xdata$pos
dformat = xdata$dformat
vrmodel = NULL
if( spec@mspec$varmodel$VAR ){
if(spec@mspec$varmodel$mxn>0){
ex = spec@mspec$varmodel$external.regressors
if( dim(ex)[1] < dim(data)[1] ) stop("\ncgarchfit-->error: external regressor matrix has less points than data!...", call. = FALSE)
orex = ex
ex = ex[1:(n - n.start), , drop = FALSE]
mxn = dim(orex)[2]
} else{
orex = NULL
ex = NULL
mxn = 0
}
if( !is.null(VAR.fit) ){
zdata = VAR.fit$xresiduals
# should be N - p
p = mp = VAR.fit$lag
N = dim(zdata)[1]
K = dim(zdata)[2]
mu = VAR.fit$xfitted
vrmodel = VAR.fit
#spec@mspec$optimization.model$maxOrder = max( mp, spec@mspec$optimization.model$maxdccOrder )
#spec@mspec$optimization.model$maxgarchOrder = mp
} else{
cat("\nestimating VAR model...")
ic = spec@mspec$varmodel$lag.criterion[1]
if( !is.null(spec@mspec$varmodel$lag.max) ){
mp = .varxselect(y = data, lag.max = spec@mspec$varmodel$lag.max, exogen = ex)$selection
mp = as.numeric( mp[which(substr(names(mp), 1, 2) == substr(ic, 1, 2))] )
} else {
mp = spec@mspec$varmodel$lag
}
vrmodel = varxfilter(X = data, p = mp, exogen = ex, Bcoef = NULL, robust = spec@mspec$varmodel$robust,
gamma = spec@mspec$varmodel$robust.control$gamma, delta = spec@mspec$varmodel$robust.control$delta,
nc = spec@mspec$varmodel$robust.control$nc, ns = spec@mspec$varmodel$robust.control$ns)
cat("done!\n")
zdata = vrmodel$xresiduals
# should be N - p
N = dim(zdata)[1]
K = dim(zdata)[2]
mu = vrmodel$xfitted
vrmodel$lag = mp
vrmodel$mxn = mxn
#spec@mspec$optimization.model$maxOrder = max( mp, spec@mspec$optimization.model$maxdccOrder )
#spec@mspec$optimization.model$maxgarchOrder = mp
}
# if we have out of sample data we need to filter
if(out.sample>0){
Bcoef = vrmodel$Bcoef
vrmodel2 = varxfilter(X = origdata, p = mp, exogen = orex, Bcoef = Bcoef)
zorigdata = vrmodel2$xresiduals
} else{
zorigdata = zdata
}
} else{
zdata = data
zorigdata = origdata
orex = NULL
ex = NULL
}
if( is(spec@uspec, "uGARCHmultispec") ){
speclist = spec@uspec
} else{
stop("\ncgarchfit-->error: cGARCHspec must contain multispec object with fixed parameters")
}
maxgarchOrder = spec@mspec$optimization.model$maxgarchOrder
tmpidx = .dccparidx( speclist )
paridx = tmpidx$idx
paridxi = tmpidx$idxi
paridxa = tmpidx$idxa
if( !is.null(fit) && is(fit, "uGARCHmultifit") ){
fitlist = fit
} else{
fitlist = multifit(multispec = speclist, data = zorigdata, out.sample = n.start, solver = solver,
solver.control = solver.control, fit.control = ufit.control, parallel = parallel,
parallel.control = parallel.control)
assign(".fitlist", fitlist, envir = .GlobalEnv)
}
garchpars = as.numeric( unlist(coef(fitlist) ) )
converge = sapply(fitlist@fit, FUN = function(x) x@fit$convergence)
if( any( converge == 1 ) ){
pr = which(converge != 1)
cat("\nNon-Converged:\n")
print(pr)
cat("\ncgarchfit-->error: convergence problem in univariate fit...")
cat("\n...returning uGARCHmultifit object instead...check and resubmit...")
return( fitlist )
}
# create a temporary environment to store values (deleted at end of function)
tempenvir = new.env(hash = TRUE, parent = .GlobalEnv)
if(is.null(fit.control$trace)) trace = FALSE else trace = as.logical(fit.control$trace)
assign("trace", trace, envir = tempenvir)
assign("method", spec@mspec$method, envir = tempenvir)
assign("eval.se", eval.se, envir = tempenvir)
assign("cnames", cnames, envir = tempenvir)
assign("spec", spec, envir = tempenvir)
assign("parallel", parallel, envir = tempenvir)
assign("parallel.control", parallel.control, envir = tempenvir)
assign("paridx", paridx, envir = tempenvir)
assign("paridxi", paridxi, envir = tempenvir)
assign("paridxa", paridxa, envir = tempenvir)
assign("speclist", speclist, envir = tempenvir)
assign("fitlist", fitlist, envir = tempenvir)
omodel = spec@mspec$optimization.model
distribution = spec@mspec$distribution
transformation = spec@mspec$transformation
sig = sigma(fitlist)
res = residuals(fitlist)
#stdresid = res/sig
assign("transformation", transformation, envir = tempenvir)
zres = sapply(fitlist@fit, FUN = function(x) x@fit$z)
ures = matrix(NA, ncol = 5, nrow = dim(zres)[1])
if( maxgarchOrder > 0 ) zres = zres[-(1:maxgarchOrder), , drop = FALSE]
# Transformation of z to [0,1] domain
ans = switch(transformation,
parametric = .pparametric(fitlist, zres),
empirical = .pempirical(zres),
spd = .pspd(zres, spd.control))
if(transformation == "spd"){
ures = ans$ures
sfit = ans$sfit
} else{
ures = ans
sfit = NULL
}
assign("spd.control", spd.control, envir = tempenvir)
# make tail adjustments in order to avoid problem with the quantile functions in
# optimization
if(any(ures > 0.99999)){
xn = which(ures > 0.99999)
ures[xn] = 0.99999
}
if(any(ures < eps)){
xn = which(ures < (1.5*eps))
ures[xn] = eps
}
assign("omodel", omodel, envir = tempenvir)
assign("mspec", spec@mspec, envir = tempenvir)
assign("ures", ures, envir = tempenvir)
assign("solver", solver, envir = tempenvir)
assign("fit.control", fit.control, envir = tempenvir)
Ifn = NULL
ILB = NULL
IUB = NULL
if( solver == "solnp" ) fit.control$stationarity = FALSE else fit.control$stationarity = TRUE
assign("fit.control", fit.control, envir = tempenvir)
if( spec@mspec$method == "ml" ){
ipars = .copulastart(data = ures, garchenv = tempenvir)
pars = ipars$pars
LB = ipars$modelLB
UB = ipars$modelUB
assign("npar", length(pars), envir = tempenvir)
if( any( LB == UB ) )
{
fixid = which(LB == UB)
if( length(fixid) == length(pars) )
{
if( !eval.se )
{
cat("\ncgarchfit-->warning: all parameters fixed...returning cgarchfilter object instead\n")
spec2 = spec
cf = coef(fitlist)
for(i in 1:k){
if(spec@uspec@type == "equal"){
spec2@uspec@spec[[i]]@optimization.model$fixed.pars = as.list(cf[,i])
} else{
spec2@uspec@spec[[i]]@optimization.model$fixed.pars = as.list(cf[[i]])
}
}
return(cgarchfilter(spec = spec2, filter.control = list(n.old = NULL), spd.control = spd.control,
out.sample = out.sample, parallel = parallel, parallel.control = parallel.control,
VAR.fit = vrmodel))
} else{
# if all parameters are fixed but we require standard errors, we
# skip the solver
use.solver = FALSE
fixed = pars[fixid]
pars = pars[-fixid]
}
} else{
# with some parameters fixed we extract them (to be rejoined at end)
# so that they do not enter the solver
use.solver = TRUE
fixed = pars[fixid]
pars = pars[-fixid]
LB = LB[-fixid]
UB = UB[-fixid]
}
} else{
# no fixed parameters, all normal
fixed = NULL
fixid = NULL
use.solver = TRUE
}
assign("fixed", fixed, envir = tempenvir)
assign("fixid", fixid, envir = tempenvir)
assign(".llh", 1, envir = tempenvir)
# get
if( use.solver )
{
# overide solver control for static R
solution = .copulasolver(solver, pars, fun = copula.normalLLH, Ifn, ILB,
IUB, gr = NULL, hessian = NULL, control = solver.control,
LB = LB, UB = UB, data = ures, returnType = "llh", garchenv = tempenvir)
sol = solution$sol
hess = solution$hess
opt.pars = sol$par
names(opt.pars) = names(pars)
convergence = sol$convergence
} else{
hess = NULL
timer = Sys.time()-tic
opt.pars = NULL
convergence = 0
sol = list()
sol$message = "all parameters fixed"
}
} else{
assign("npar", 0, envir = tempenvir)
fixed = NULL
fixid = NULL
use.solver = FALSE
assign("fixed", fixed, envir = tempenvir)
assign("fixid", fixid, envir = tempenvir)
assign(".llh", 1, envir = tempenvir)
hess = NULL
timer = Sys.time()-tic
opt.pars = NULL
convergence = 0
sol = list()
sol$message = NULL
}
mfit = list()
mfit$dates = dates
mfit$date.format = dformat
mfit$origdata = origdata
mfit$origdates = origdates
mfit$vrmodel = vrmodel
# convert back from uniform to standard variates
model = list()
model$copula.distribution = distribution
model$out.sample = out.sample
model$transformation = transformation
model$spd.control = spd.control
model$sfit = sfit
model$time.varying = FALSE
# add some tail data of the univariate fits to be used with simulation
mfit$tailusigma = tail(sigma(fitlist), 25)
mfit$tailuresids = tail(residuals(fitlist), 25)
mfit$tailuret = tail(data, 25)
cfit = copula.normalLLH(opt.pars, data = ures, returnType = "ALL", garchenv = tempenvir)
if(eval.se && spec@mspec$method == "ml"){
nderiv = .cgarchcvar(f = copula.normalLLH, pars = opt.pars, garchfit = fitlist, model = model,
env = tempenvir, fname = "copula.normalLLH")
mfit$pars = nderiv$pars
mfit$garchpars = coef(fitlist)
# stpars = second stage parameters
mfit$stpars = opt.pars
mfit$scores = nderiv$scores
mfit$A = nderiv$A
mfit$cvar = nderiv$cvar
mfit$matcoef = nderiv$matcoef
} else{
if( spec@mspec$method == "ml" ){
garchNames = vector(mode = "character")
assetNames = cnames
dccNames = spec@mspec$optimization.model$modelnames
m = length(spec@uspec@spec)
for(i in 1:m){
cf = coef(fitlist@fit[[i]])
cfnames = names(cf)
garchNames = c(garchNames, paste(assetNames[i], cfnames, sep = ".") )
}
allnames = c(garchNames, dccNames)
mfit$pars = c(as.numeric(unlist(coef(fitlist))), opt.pars)
mfit$garchpars = coef(fitlist)
# stpars = second stage parameters
mfit$stpars = opt.pars
names(mfit$pars) = allnames
mfit$scores = NULL
mfit$A = NULL
mfit$cvar = NULL
matcoef = matrix(NA, ncol = 4, nrow = length(mfit$pars))
matcoef[,1] = mfit$pars
dimnames(matcoef) = list(allnames, c(" Estimate",
" Std. Error", " t value", "Pr(>|t|)"))
mfit$matcoef = matcoef
} else{
garchNames = vector(mode = "character")
assetNames = cnames
dccNames = NULL
m = length(spec@uspec@spec)
for(i in 1:m){
cf = coef(fitlist@fit[[i]])
cfnames = names(cf)
garchNames = c(garchNames, paste(assetNames[i], cfnames, sep = ".") )
}
allnames = c(garchNames, dccNames)
mfit$pars = c(as.numeric(unlist(coef(fitlist))), opt.pars)
mfit$garchpars = coef(fitlist)
# stpars = second stage parameters
mfit$stpars = NULL
names(mfit$pars) = allnames
matcoef = NULL
mfit$scores = NULL
mfit$A = NULL
mfit$cvar = NULL
for(i in 1:m){
matcoef = rbind(matcoef, fitlist@fit[[i]]@fit$robust.matcoef)
}
dimnames(matcoef) = list(allnames, c(" Estimate",
" Std. Error", " t value", "Pr(>|t|)"))
mfit$matcoef = matcoef
}
}
mfit$paridx = tmpidx
garchllh = sapply(fitlist@fit, FUN = function(x) -x@fit$log.likelihoods)
copllh = cfit$lik
if( maxgarchOrder > 0 ) copllh = c(rep(0, maxgarchOrder), copllh)
mfit$llh = sum( apply(cbind(copllh, garchllh), 1, "sum") )
mfit$zres = zres
mfit$timer = Sys.time() - tic
ans = new("cGARCHfit",
mfit = mfit,
mspec = spec@mspec,
uspec = spec@uspec,
copula = cfit,
model = model)
return(ans)
}
.cgarchfit.tvstudent = function(spec, data, spd.control = list(lower = 0.1, upper = 0.9, type = "pwm", kernel = "epanech"),
fit.control = list(eval.se = TRUE, trace = TRUE, stationarity = TRUE), solver = "solnp", solver.control = list(), out.sample = 0, parallel = FALSE,
parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2), fit = NULL, VAR.fit = NULL, ...)
{
tic = Sys.time()
ufit.control = list()
if(is.null(fit.control$stationarity)){
ufit.control$stationarity = TRUE
fit.control$stationarity = TRUE
} else {
ufit.control$stationarity = fit.control$stationarity
fit.control$stationarity = fit.control$stationarity
}
if(is.null(fit.control$scale)){
ufit.control$scale = TRUE
} else{
ufit.control$scale = fit.control$scale
fit.control$scale = NULL
}
if(is.null(fit.control$eval.se)){
eval.se = TRUE
} else{
eval.se = as.logical(fit.control$eval.se)
}
# Parallel Execution Prelim Checks (we need them because of the hessian estimation)
if( parallel ){
os = .Platform$OS.type
if(is.null(parallel.control$pkg)){
if( os == "windows" ) parallel.control$pkg = "snowfall" else parallel.control$pkg = "multicore"
if( is.null(parallel.control$cores) ) parallel.control$cores = 2
} else{
mtype = match(tolower(parallel.control$pkg[1]), c("multicore", "snowfall"))
if(is.na(mtype)) stop("\nParallel Package type not recognized in parallel.control\n")
parallel.control$pkg = tolower(parallel.control$pkg[1])
if( os == "windows" && parallel.control$pkg == "multicore" ) stop("\nmulticore not supported on windows O/S\n")
if( is.null(parallel.control$cores) ) parallel.control$cores = 2 else parallel.control$cores = as.integer(parallel.control$cores[1])
}
}
if( is.null( colnames(data) ) ) cnames = paste("Asset_", 1:k, sep = "") else cnames = colnames(data)
xdata = .extractmdata(data)
if(!is.numeric(out.sample)) stop("\ncgarchfit-->error: out.sample must be numeric\n")
if(as.numeric(out.sample) < 0) stop("\ncgarchfit-->error: out.sample must be positive\n")
n.start = round(out.sample, 0)
n = dim(xdata$data)[1]
k = dim(xdata$data)[2]
if( (n-n.start) < 100) stop("\ncgarchfit-->error: function requires at least 100 data\n points to run\n")
data = xdata$data[1:(n-n.start), , drop = FALSE]
dates = xdata$pos[1:(n-n.start), drop = FALSE]
origdata = xdata$data
origdates = xdata$pos
dformat = xdata$dformat
vrmodel = NULL
if( spec@mspec$varmodel$VAR ){
if(spec@mspec$varmodel$mxn>0){
ex = spec@mspec$varmodel$external.regressors
if( dim(ex)[1] < dim(data)[1] ) stop("\ncgarchfit-->error: external regressor matrix has less points than data!...", call. = FALSE)
orex = ex
ex = ex[1:(n - n.start), , drop = FALSE]
mxn = dim(orex)[2]
} else{
orex = NULL
ex = NULL
mxn = 0
}
if( !is.null(VAR.fit) ){
zdata = VAR.fit$xresiduals
# should be N - p
p = mp = VAR.fit$lag
N = dim(zdata)[1]
K = dim(zdata)[2]
mu = VAR.fit$xfitted
vrmodel = VAR.fit
#spec@mspec$optimization.model$maxOrder = max( mp, spec@mspec$optimization.model$maxdccOrder )
#spec@mspec$optimization.model$maxgarchOrder = mp
} else{
cat("\nestimating VAR model...")
ic = spec@mspec$varmodel$lag.criterion[1]
if( !is.null(spec@mspec$varmodel$lag.max) ){
mp = .varxselect(y = data, lag.max = spec@mspec$varmodel$lag.max, exogen = ex)$selection
mp = as.numeric( mp[which(substr(names(mp), 1, 2) == substr(ic, 1, 2))] )
} else {
mp = spec@mspec$varmodel$lag
}
vrmodel = varxfilter(X = data, p = mp, exogen = ex, Bcoef = NULL, robust = spec@mspec$varmodel$robust,
gamma = spec@mspec$varmodel$robust.control$gamma, delta = spec@mspec$varmodel$robust.control$delta,
nc = spec@mspec$varmodel$robust.control$nc, ns = spec@mspec$varmodel$robust.control$ns)
cat("done!\n")
zdata = vrmodel$xresiduals
# should be N - p
N = dim(zdata)[1]
K = dim(zdata)[2]
mu = vrmodel$xfitted
vrmodel$lag = mp
vrmodel$mxn = mxn
#spec@mspec$optimization.model$maxOrder = max( mp, spec@mspec$optimization.model$maxdccOrder )
#spec@mspec$optimization.model$maxgarchOrder = mp
}
# if we have out of sample data we need to filter
if(out.sample>0){
Bcoef = vrmodel$Bcoef
vrmodel2 = varxfilter(X = origdata, p = mp, exogen = orex, Bcoef = Bcoef)
zorigdata = vrmodel2$xresiduals
} else{
zorigdata = zdata
}
} else{
zdata = data
zorigdata = origdata
orex = NULL
ex = NULL
}
if( is(spec@uspec, "uGARCHmultispec") ){
speclist = spec@uspec
} else{
stop("\ncgarchfit-->error: cGARCHspec must contain multispec object with fixed parameters")
}
maxgarchOrder = spec@mspec$optimization.model$maxgarchOrder
tmpidx = .dccparidx( speclist )
paridx = tmpidx$idx
paridxi = tmpidx$idxi
paridxa = tmpidx$idxa
if( !is.null(fit) && is(fit, "uGARCHmultifit") ){
fitlist = fit
} else{
fitlist = multifit(multispec = speclist, data = zorigdata, out.sample = n.start, solver = solver,
solver.control = solver.control, fit.control = ufit.control, parallel = parallel,
parallel.control = parallel.control)
assign(".fitlist", fitlist, envir = .GlobalEnv)
}
garchpars = as.numeric( unlist(coef(fitlist) ) )
converge = sapply(fitlist@fit, FUN = function(x) x@fit$convergence)
if( any( converge == 1 ) ){
pr = which(converge != 1)
cat("\nNon-Converged:\n")
print(pr)
cat("\ncgarchfit-->error: convergence problem in univariate fit...")
cat("\n...returning uGARCHmultifit object instead...check and resubmit...")
return( fitlist )
}
# create a temporary environment to store values (deleted at end of function)
tempenvir = new.env(hash = TRUE, parent = .GlobalEnv)
if(is.null(fit.control$trace)) trace = FALSE else trace = as.logical(fit.control$trace)
assign("trace", trace, envir = tempenvir)
assign("eval.se", eval.se, envir = tempenvir)
assign("cnames", cnames, envir = tempenvir)
assign("spec", spec, envir = tempenvir)
assign("parallel", parallel, envir = tempenvir)
assign("parallel.control", parallel.control, envir = tempenvir)
assign("paridx", paridx, envir = tempenvir)
assign("paridxi", paridxi, envir = tempenvir)
assign("paridxa", paridxa, envir = tempenvir)
assign("speclist", speclist, envir = tempenvir)
assign("fitlist", fitlist, envir = tempenvir)
omodel = spec@mspec$optimization.model
distribution = spec@mspec$distribution
transformation = spec@mspec$transformation
sig = sigma(fitlist)
res = residuals(fitlist)
#stdresid = res/sig
assign("transformation", transformation, envir = tempenvir)
zres = sapply(fitlist@fit, FUN = function(x) x@fit$z)
ures = matrix(NA, ncol = 5, nrow = dim(zres)[1])
if( maxgarchOrder > 0 ) zres = zres[-(1:maxgarchOrder), , drop = FALSE]
# Transformation of z to [0,1] domain
ans = switch(transformation,
parametric = .pparametric(fitlist, zres),
empirical = .pempirical(zres),
spd = .pspd(zres, spd.control))
if(transformation == "spd"){
ures = ans$ures
sfit = ans$sfit
} else{
ures = ans
sfit = NULL
}
assign("spd.control", spd.control, envir = tempenvir)
# make tail adjustments in order to avoid problem with the quantile functions in
# optimization
if(any(ures > 0.99999)){
xn = which(ures > 0.99999)
ures[xn] = 0.99999
}
if(any(ures < eps)){
xn = which(ures < (1.5*eps))
ures[xn] = eps
}
assign("omodel", omodel, envir = tempenvir)
assign("mspec", spec@mspec, envir = tempenvir)
assign("ures", ures, envir = tempenvir)
assign("solver", solver, envir = tempenvir)
assign("fit.control", fit.control, envir = tempenvir)
Ifn = .dcccon
ILB = 0.0001
IUB = 0.9998
if( solver == "solnp" ) fit.control$stationarity = FALSE else fit.control$stationarity = TRUE
assign("fit.control", fit.control, envir = tempenvir)
ipars = .tvcopulastart(data = ures, garchenv = tempenvir)
pars = ipars$pars
LB = ipars$modelLB
UB = ipars$modelUB
assign("npar", length(pars), envir = tempenvir)
if( any( LB == UB ) )
{
fixid = which(LB == UB)
if( length(fixid) == length(pars) )
{
if( !eval.se )
{
cat("\ncgarchfit-->warning: all parameters fixed...returning cgarchfilter object instead\n")
spec2 = spec
cf = coef(fitlist)
for(i in 1:k){
if(spec@uspec@type == "equal"){
spec2@uspec@spec[[i]]@optimization.model$fixed.pars = as.list(cf[,i])
} else{
spec2@uspec@spec[[i]]@optimization.model$fixed.pars = as.list(cf[[i]])
}
}
return(cgarchfilter(spec = spec2, filter.control = list(n.old = NULL), spd.control = spd.control,
out.sample = out.sample, parallel = parallel, parallel.control = parallel.control,
VAR.fit = vrmodel))
} else{
# if all parameters are fixed but we require standard errors, we
# skip the solver
use.solver = FALSE
fixed = pars[fixid]
pars = pars[-fixid]
}
} else{
# with some parameters fixed we extract them (to be rejoined at end)
# so that they do not enter the solver
use.solver = TRUE
fixed = pars[fixid]
pars = pars[-fixid]
LB = LB[-fixid]
UB = UB[-fixid]
}
} else{
# no fixed parameters, all normal
fixed = NULL
fixid = NULL
use.solver = TRUE
}
assign("fixed", fixed, envir = tempenvir)
assign("fixid", fixid, envir = tempenvir)
assign(".llh", 1, envir = tempenvir)
# get
if( use.solver )
{
solution = .copulasolver(solver, pars, fun = copula.tvstudentLLH, Ifn, ILB,
IUB, gr = NULL, hessian = NULL, control = solver.control,
LB = LB, UB = UB, data = ures, returnType = "llh", garchenv = tempenvir)
sol = solution$sol
hess = solution$hess
opt.pars = sol$par
names(opt.pars) = names(pars)
convergence = sol$convergence
} else{
hess = NULL
timer = Sys.time()-tic
opt.pars = NULL
convergence = 0
sol = list()
sol$message = "all parameters fixed"
}
mfit = list()
mfit$dates = dates
mfit$date.format = dformat
mfit$origdata = origdata
mfit$origdates = origdates
mfit$vrmodel = vrmodel
# convert back from uniform to standard variates
model = list()
model$copula.distribution = distribution
model$out.sample = out.sample
model$transformation = transformation
model$spd.control = spd.control
model$sfit = sfit
model$time.varying = TRUE
# add some tail data of the univariate fits to be used with simulation
mfit$tailusigma = tail(sigma(fitlist), 25)
mfit$tailuresids = tail(residuals(fitlist), 25)
mfit$tailuret = tail(data, 25)
cfit = copula.tvstudentLLH(opt.pars, data = ures, returnType = "ALL", garchenv = tempenvir)
if(eval.se){
nderiv = .cgarchcvar(f = copula.tvstudentLLH, pars = opt.pars, garchfit = fitlist, model = model,
env = tempenvir, fname = "copula.tvstudentLLH")
mfit$pars = nderiv$pars
mfit$garchpars = coef(fitlist)
# stpars = second stage parameters
mfit$stpars = opt.pars
mfit$scores = nderiv$scores
mfit$A = nderiv$A
mfit$cvar = nderiv$cvar
mfit$matcoef = nderiv$matcoef
} else{
garchNames = vector(mode = "character")
assetNames = cnames
dccNames = spec@mspec$optimization.model$modelnames
m = length(spec@uspec@spec)
for(i in 1:m){
cf = coef(fitlist@fit[[i]])
cfnames = names(cf)
garchNames = c(garchNames, paste(assetNames[i], cfnames, sep = ".") )
}
allnames = c(garchNames, dccNames)
mfit$pars = c(as.numeric(unlist(coef(fitlist))), opt.pars)
mfit$garchpars = coef(fitlist)
# stpars = second stage parameters
mfit$stpars = opt.pars
names(mfit$pars) = allnames
mfit$scores = NULL
mfit$A = NULL
mfit$cvar = NULL
matcoef = matrix(NA, ncol = 4, nrow = length(mfit$pars))
matcoef[,1] = mfit$pars
dimnames(matcoef) = list(allnames, c(" Estimate",
" Std. Error", " t value", "Pr(>|t|)"))
mfit$matcoef = matcoef
}
mfit$paridx = tmpidx
garchllh = sapply(fitlist@fit, FUN = function(x) -x@fit$log.likelihoods)
copllh = cfit$lik
if( maxgarchOrder > 0 ) copllh = c(rep(0, maxgarchOrder), copllh)
mfit$llh = sum( apply(cbind(copllh, garchllh), 1, "sum") )
mfit$zres = zres
mfit$timer = Sys.time() - tic
ans = new("cGARCHfit",
mfit = mfit,
mspec = spec@mspec,
uspec = spec@uspec,
copula = cfit,
model = model)
return(ans)
}
.cgarchfit.student = function(spec, data, spd.control = list(lower = 0.1, upper = 0.9, type = "pwm", kernel = "epanech"),
fit.control = list(eval.se = TRUE, trace = TRUE, stationarity = TRUE), solver = "solnp", solver.control = list(), out.sample = 0, parallel = FALSE,
parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2), fit = NULL, VAR.fit = NULL, ...)
{
tic = Sys.time()
ufit.control = list()
if(is.null(fit.control$stationarity)){
ufit.control$stationarity = TRUE
} else {
ufit.control$stationarity = fit.control$stationarity
fit.control$stationarity = NULL
}
if(is.null(fit.control$scale)){
ufit.control$scale = TRUE
} else{
ufit.control$scale = fit.control$scale
fit.control$scale = NULL
}
if(is.null(fit.control$eval.se)){
eval.se = TRUE
} else{
eval.se = as.logical(fit.control$eval.se)
}
# Parallel Execution Prelim Checks (we need them because of the hessian estimation)
if( parallel ){
os = .Platform$OS.type
if(is.null(parallel.control$pkg)){
if( os == "windows" ) parallel.control$pkg = "snowfall" else parallel.control$pkg = "multicore"
if( is.null(parallel.control$cores) ) parallel.control$cores = 2
} else{
mtype = match(tolower(parallel.control$pkg[1]), c("multicore", "snowfall"))
if(is.na(mtype)) stop("\nParallel Package type not recognized in parallel.control\n")
parallel.control$pkg = tolower(parallel.control$pkg[1])
if( os == "windows" && parallel.control$pkg == "multicore" ) stop("\nmulticore not supported on windows O/S\n")
if( is.null(parallel.control$cores) ) parallel.control$cores = 2 else parallel.control$cores = as.integer(parallel.control$cores[1])
}
}
if( is.null( colnames(data) ) ) cnames = paste("Asset_", 1:k, sep = "") else cnames = colnames(data)
xdata = .extractmdata(data)
if(!is.numeric(out.sample)) stop("\ncgarchfit-->error: out.sample must be numeric\n")
if(as.numeric(out.sample) < 0) stop("\ncgarchfit-->error: out.sample must be positive\n")
n.start = round(out.sample, 0)
n = dim(xdata$data)[1]
k = dim(xdata$data)[2]
if( (n-n.start) < 100) stop("\ncgarchfit-->error: function requires at least 100 data\n points to run\n")
data = xdata$data[1:(n-n.start), , drop = FALSE]
dates = xdata$pos[1:(n-n.start), drop = FALSE]
origdata = xdata$data
origdates = xdata$pos
dformat = xdata$dformat
vrmodel = NULL
if( spec@mspec$varmodel$VAR ){
if(spec@mspec$varmodel$mxn>0){
ex = spec@mspec$varmodel$external.regressors
if( dim(ex)[1] < dim(data)[1] ) stop("\ncgarchfit-->error: external regressor matrix has less points than data!...", call. = FALSE)
orex = ex
ex = ex[1:(n - n.start), , drop = FALSE]
mxn = dim(orex)[2]
} else{
orex = NULL
ex = NULL
mxn = 0
}
if( !is.null(VAR.fit) ){
zdata = VAR.fit$xresiduals
# should be N - p
p = mp = VAR.fit$lag
N = dim(zdata)[1]
K = dim(zdata)[2]
mu = VAR.fit$xfitted
vrmodel = VAR.fit
#spec@mspec$optimization.model$maxOrder = max( mp, spec@mspec$optimization.model$maxdccOrder )
#spec@mspec$optimization.model$maxgarchOrder = mp
} else{
cat("\nestimating VAR model...")
ic = spec@mspec$varmodel$lag.criterion[1]
if( !is.null(spec@mspec$varmodel$lag.max) ){
mp = .varxselect(y = data, lag.max = spec@mspec$varmodel$lag.max, exogen = ex)$selection
mp = as.numeric( mp[which(substr(names(mp), 1, 2) == substr(ic, 1, 2))] )
} else {
mp = spec@mspec$varmodel$lag
}
vrmodel = varxfilter(X = data, p = mp, exogen = ex, Bcoef = NULL, robust = spec@mspec$varmodel$robust,
gamma = spec@mspec$varmodel$robust.control$gamma, delta = spec@mspec$varmodel$robust.control$delta,
nc = spec@mspec$varmodel$robust.control$nc, ns = spec@mspec$varmodel$robust.control$ns)
cat("done!\n")
zdata = vrmodel$xresiduals
# should be N - p
N = dim(zdata)[1]
K = dim(zdata)[2]
mu = vrmodel$xfitted
vrmodel$lag = mp
vrmodel$mxn = mxn
#spec@mspec$optimization.model$maxOrder = max( mp, spec@mspec$optimization.model$maxdccOrder )
#spec@mspec$optimization.model$maxgarchOrder = mp
}
# if we have out of sample data we need to filter
if(out.sample>0){
Bcoef = vrmodel$Bcoef
vrmodel2 = varxfilter(X = origdata, p = mp, exogen = orex, Bcoef = Bcoef)
zorigdata = vrmodel2$xresiduals
} else{
zorigdata = zdata
}
} else{
zdata = data
zorigdata = origdata
orex = NULL
ex = NULL
}
if( is(spec@uspec, "uGARCHmultispec") ){
speclist = spec@uspec
} else{
stop("\ncgarchfit-->error: cGARCHspec must contain multispec object with fixed parameters")
}
maxgarchOrder = spec@mspec$optimization.model$maxgarchOrder
tmpidx = .dccparidx( speclist )
paridx = tmpidx$idx
paridxi = tmpidx$idxi
paridxa = tmpidx$idxa
if( !is.null(fit) && is(fit, "uGARCHmultifit") ){
fitlist = fit
} else{
fitlist = multifit(multispec = speclist, data = zorigdata, out.sample = n.start, solver = solver,
solver.control = solver.control, fit.control = ufit.control, parallel = parallel,
parallel.control = parallel.control)
assign(".fitlist", fitlist, envir = .GlobalEnv)
}
garchpars = as.numeric( unlist(coef(fitlist) ) )
converge = sapply(fitlist@fit, FUN = function(x) x@fit$convergence)
if( any( converge == 1 ) ){
pr = which(converge != 1)
cat("\nNon-Converged:\n")
print(pr)
cat("\ncgarchfit-->error: convergence problem in univariate fit...")
cat("\n...returning uGARCHmultifit object instead...check and resubmit...")
return( fitlist )
}
# create a temporary environment to store values (deleted at end of function)
tempenvir = new.env(hash = TRUE, parent = .GlobalEnv)
if(is.null(fit.control$trace)) trace = FALSE else trace = as.logical(fit.control$trace)
assign("trace", trace, envir = tempenvir)
assign("method", spec@mspec$method, envir = tempenvir)
assign("eval.se", eval.se, envir = tempenvir)
assign("cnames", cnames, envir = tempenvir)
assign("spec", spec, envir = tempenvir)
assign("parallel", parallel, envir = tempenvir)
assign("parallel.control", parallel.control, envir = tempenvir)
assign("paridx", paridx, envir = tempenvir)
assign("paridxi", paridxi, envir = tempenvir)
assign("paridxa", paridxa, envir = tempenvir)
assign("speclist", speclist, envir = tempenvir)
assign("fitlist", fitlist, envir = tempenvir)
omodel = spec@mspec$optimization.model
distribution = spec@mspec$distribution
transformation = spec@mspec$transformation
sig = sigma(fitlist)
res = residuals(fitlist)
#stdresid = res/sig
assign("transformation", transformation, envir = tempenvir)
zres = sapply(fitlist@fit, FUN = function(x) x@fit$z)
ures = matrix(NA, ncol = 5, nrow = dim(zres)[1])
if( maxgarchOrder > 0 ) zres = zres[-(1:maxgarchOrder), , drop = FALSE]
# Transformation of z to [0,1] domain
ans = switch(transformation,
parametric = .pparametric(fitlist, zres),
empirical = .pempirical(zres),
spd = .pspd(zres, spd.control))
if(transformation == "spd"){
ures = ans$ures
sfit = ans$sfit
} else{
ures = ans
sfit = NULL
}
assign("spd.control", spd.control, envir = tempenvir)
# make tail adjustments in order to avoid problem with the quantile functions in
# optimization
if(any(ures > 0.99999)){
xn = which(ures > 0.99999)
ures[xn] = 0.99999
}
if(any(ures < eps)){
xn = which(ures < (1.5*eps))
ures[xn] = eps
}
assign("omodel", omodel, envir = tempenvir)
assign("mspec", spec@mspec, envir = tempenvir)
assign("ures", ures, envir = tempenvir)
assign("solver", solver, envir = tempenvir)
assign("fit.control", fit.control, envir = tempenvir)
Ifn = NULL
ILB = NULL
IUB = NULL
if( solver == "solnp" ) fit.control$stationarity = FALSE else fit.control$stationarity = TRUE
assign("fit.control", fit.control, envir = tempenvir)
ipars = .copulastart(data = ures, garchenv = tempenvir)
pars = ipars$pars
LB = ipars$modelLB
UB = ipars$modelUB
assign("npar", length(pars), envir = tempenvir)
if( any( LB == UB ) )
{
fixid = which(LB == UB)
if( length(fixid) == length(pars) )
{
if( !eval.se )
{
# if all parameters are fixed an no standard erros are to
# be calculated then we return a ugarchfilter object
cat("\ncgarchfit-->warning: all parameters fixed...returning cgarchfilter object instead\n")
spec2 = spec
cf = coef(fitlist)
for(i in 1:k){
if(spec@uspec@type == "equal"){
spec2@uspec@spec[[i]]@optimization.model$fixed.pars = as.list(cf[,i])
} else{
spec2@uspec@spec[[i]]@optimization.model$fixed.pars = as.list(cf[[i]])
}
}
return(cgarchfilter(spec = spec2, filter.control = list(n.old = NULL), spd.control = spd.control,
out.sample = out.sample, parallel = parallel, parallel.control = parallel.control,
VAR.fit = vrmodel))
} else{
# if all parameters are fixed but we require standard errors, we
# skip the solver
use.solver = FALSE
fixed = pars[fixid]
pars = pars[-fixid]
}
} else{
# with some parameters fixed we extract them (to be rejoined at end)
# so that they do not enter the solver
use.solver = TRUE
fixed = pars[fixid]
pars = pars[-fixid]
LB = LB[-fixid]
UB = UB[-fixid]
}
} else{
# no fixed parameters, all normal
fixed = NULL
fixid = NULL
use.solver = TRUE
}
assign("fixed", fixed, envir = tempenvir)
assign("fixid", fixid, envir = tempenvir)
assign(".llh", 1, envir = tempenvir)
# get
if( use.solver )
{
# overide solver control for static R and use plain optim (L-BFGS-B)
if(spec@mspec$method == "kendall"){
solution = .copulasolver(solver = "lbfgs", pars, fun = copula.studentLLH, Ifn, ILB,
IUB, gr = NULL, hessian = NULL, control = list(),
LB = LB, UB = UB, data = ures, returnType = "llh", garchenv = tempenvir)
sol = solution$sol
hess = solution$hess
opt.pars = sol$par
names(opt.pars) = names(pars)
convergence = sol$convergence
} else{
solution = .copulasolver(solver, pars, fun = copula.studentLLH, Ifn, ILB,
IUB, gr = NULL, hessian = NULL, control = solver.control,
LB = LB, UB = UB, data = ures, returnType = "llh", garchenv = tempenvir)
sol = solution$sol
hess = solution$hess
opt.pars = sol$par
names(opt.pars) = names(pars)
convergence = sol$convergence
}
} else{
hess = NULL
timer = Sys.time()-tic
opt.pars = NULL
convergence = 0
sol = list()
sol$message = "all parameters fixed"
}
mfit = list()
mfit$dates = dates
mfit$date.format = dformat
mfit$origdata = origdata
mfit$origdates = origdates
mfit$vrmodel = vrmodel
# convert back from uniform to standard variates
model = list()
model$copula.distribution = distribution
model$out.sample = out.sample
model$transformation = transformation
model$spd.control = spd.control
model$sfit = sfit
model$time.varying = FALSE
# add some tail data of the univariate fits to be used with simulation
mfit$tailusigma = tail(sigma(fitlist), 25)
mfit$tailuresids = tail(residuals(fitlist), 25)
mfit$tailuret = tail(data, 25)
cfit = copula.studentLLH(opt.pars, data = ures, returnType = "ALL", garchenv = tempenvir)
if(eval.se){
nderiv = .cgarchcvar(f = copula.studentLLH, pars = opt.pars, garchfit = fitlist, model = model,
env = tempenvir, fname = "copula.studentLLH")
mfit$pars = nderiv$pars
mfit$garchpars = coef(fitlist)
mfit$stpars = opt.pars
mfit$scores = nderiv$scores
mfit$A = nderiv$A
mfit$cvar = nderiv$cvar
mfit$matcoef = nderiv$matcoef
} else{
garchNames = vector(mode = "character")
assetNames = cnames
dccNames = spec@mspec$optimization.model$modelnames
m = length(spec@uspec@spec)
for(i in 1:m){
cf = coef(fitlist@fit[[i]])
cfnames = names(cf)
garchNames = c(garchNames, paste(assetNames[i], cfnames, sep = ".") )
}
allnames = c(garchNames, dccNames)
mfit$pars = c(as.numeric(unlist(coef(fitlist))), opt.pars)
names(mfit$pars) = allnames
mfit$garchpars = coef(fitlist)
mfit$stpars = opt.pars
mfit$scores = NULL
mfit$A = NULL
mfit$cvar = NULL
matcoef = matrix(NA, ncol = 4, nrow = length(mfit$pars))
matcoef[,1] = mfit$pars
dimnames(matcoef) = list(allnames, c(" Estimate",
" Std. Error", " t value", "Pr(>|t|)"))
mfit$matcoef = matcoef
}
mfit$paridx = tmpidx
garchllh = sapply(fitlist@fit, FUN = function(x) -x@fit$log.likelihoods)
copllh = cfit$lik
if( maxgarchOrder > 0 ) copllh = c(rep(0, maxgarchOrder), copllh)
mfit$llh = sum( apply(cbind(copllh, garchllh), 1, "sum") )
mfit$zres = zres
mfit$timer = Sys.time() - tic
ans = new("cGARCHfit",
mfit = mfit,
mspec = spec@mspec,
uspec = spec@uspec,
copula = cfit,
model = model)
return(ans)
}
.cgarchfilter = function(spec, data, out.sample = 0, filter.control = list(n.old = NULL), spd.control = list(lower = 0.1, upper = 0.9, type = "pwm", kernel = "epanech"),
parallel = FALSE, parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2), VAR.fit = NULL, ...)
{
dist = spec@mspec$distribution
tv = ifelse(spec@mspec$timecopula, 1, 0)
dx = paste(dist, tv, sep = "")
ans = switch(dx,
mvt0 = .cgarchfilter.student(spec = spec, data = data, spd.control = spd.control, out.sample = out.sample, parallel = parallel,
parallel.control = parallel.control, VAR.fit = VAR.fit, ...),
mvt1 = .cgarchfilter.tvstudent(spec = spec, data = data, spd.control = spd.control, out.sample = out.sample, parallel = parallel,
parallel.control = parallel.control, VAR.fit = VAR.fit, ...),
mvnorm0 = .cgarchfilter.norm(spec = spec, data = data, spd.control = spd.control, out.sample = out.sample, parallel = parallel,
parallel.control = parallel.control, VAR.fit = VAR.fit, ...),
mvnorm1 = .cgarchfilter.tvnorm(spec = spec, data = data, spd.control = spd.control, out.sample = out.sample, parallel = parallel,
parallel.control = parallel.control, VAR.fit = VAR.fit, ...))
return( ans )
}
setMethod("cgarchfilter", signature(spec = "cGARCHspec"), .cgarchfilter)
.cgarchfilter.tvstudent = function(spec, data, filter.control = list(n.old = NULL), spd.control = list(lower = 0.1, upper = 0.9, type = "pwm", kernel = "epanech"),
out.sample = 0, parallel = FALSE, parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2),
VAR.fit = NULL, ...)
{
tic = Sys.time()
if( parallel ){
os = .Platform$OS.type
if(is.null(parallel.control$pkg)){
if( os == "windows" ) parallel.control$pkg = "snowfall" else parallel.control$pkg = "multicore"
if( is.null(parallel.control$cores) ) parallel.control$cores = 2
} else{
mtype = match(tolower(parallel.control$pkg[1]), c("multicore", "snowfall"))
if(is.na(mtype)) stop("\nParallel Package type not recognized in parallel.control\n")
parallel.control$pkg = tolower(parallel.control$pkg[1])
if( os == "windows" && parallel.control$pkg == "multicore" ) stop("\nmulticore not supported on windows O/S\n")
if( is.null(parallel.control$cores) ) parallel.control$cores = 2 else parallel.control$cores = as.integer(parallel.control$cores[1])
}
}
n.old = filter.control$n.old
if( is.null( colnames(data) ) ) cnames = paste("Asset_", 1:k, sep = "") else cnames = colnames(data)
xdata = .extractmdata(data)
if(!is.numeric(out.sample)) stop("\ncgarchfilter-->error: out.sample must be numeric\n")
if(as.numeric(out.sample) < 0) stop("\ncgarchfilter-->error: out.sample must be positive\n")
n.start = round(out.sample, 0)
n = dim(xdata$data)[1]
k = dim(xdata$data)[2]
if( (n - n.start) < 100) stop("\ncgarchfilter-->error: function requires at least 100 data\n points to run\n")
data = xdata$data[1:(n - n.start), , drop = FALSE]
dates = xdata$pos[1:(n - n.start)]
origdata = xdata$data
origdates = xdata$pos
dformat = xdata$dformat
vrmodel = NULL
if( spec@mspec$varmodel$VAR ){
if(spec@mspec$varmodel$mxn>0){
ex = spec@mspec$varmodel$external.regressors
if( dim(ex)[1] < dim(data)[1] ) stop("\ndccfilter-->error: external regressor matrix has less points than data!...", call. = FALSE)
orex = ex
ex = ex[1:(n - n.start), , drop = FALSE]
mxn = dim(orex)[2]
} else{
orex = NULL
ex = NULL
mxn = 0
}
if( !is.null(VAR.fit) ){
zdata = VAR.fit$xresiduals
# should be N - p
p = mp = VAR.fit$lag
N = dim(zdata)[1]
K = dim(zdata)[2]
mu = VAR.fit$xfitted
vrmodel = VAR.fit
if(out.sample>0){
Bcoef = VAR.fit$Bcoef
vrmodel2 = varxfilter(X = origdata, p = mp, exogen = orex, Bcoef = Bcoef)
zorigdata = vrmodel2$xresiduals
} else{
zorigdata = zdata
}
} else{
cat("\nestimating VAR model...")
ic = spec@mspec$varmodel$lag.criterion[1]
if( !is.null(spec@mspec$varmodel$lag.max) ){
mp = .varxselect(y = data, lag.max = spec@mspec$varmodel$lag.max, exogen = ex)$selection
mp = as.numeric( mp[which(substr(names(mp), 1, 2) == substr(ic, 1, 2))] )
} else {
mp = spec@mspec$varmodel$lag
}
vrmodel = varxfilter(X = data, p = mp, exogen = ex, Bcoef = NULL, robust = spec@mspec$varmodel$robust,
gamma = spec@mspec$varmodel$robust.control$gamma, delta = spec@mspec$varmodel$robust.control$delta,
nc = spec@mspec$varmodel$robust.control$nc, ns = spec@mspec$varmodel$robust.control$ns)
cat("done!\n")
zdata = vrmodel$xresiduals
# should be N - p
N = dim(zdata)[1]
K = dim(zdata)[2]
mu = vrmodel$xfitted
vrmodel$lag = mp
vrmodel$mxn = mxn
# if we have out of sample data we need to filter
if(out.sample>0){
Bcoef = vrmodel$Bcoef
vrmodel2 = varxfilter(X = origdata, p = mp, exogen = orex, Bcoef = Bcoef)
zorigdata = vrmodel2$xresiduals
} else{
zorigdata = zdata
}
}
} else{
zdata = data
zorigdata = origdata
orex = NULL
ex = NULL
}
if( is(spec@uspec, "uGARCHmultispec") ){
speclist = spec@uspec
} else{
stop("\ncgarchfilter-->error: cGARCHspec must contain multispec object with fixed parameters")
}
maxgarchOrder = spec@mspec$optimization.model$maxgarchOrder
filterlist = multifilter(multifitORspec = spec@uspec, data = as.data.frame(zorigdata), out.sample = out.sample,
parallel = parallel, parallel.control = parallel.control, n.old = n.old)
tmpidx = .dccparidx( speclist )
paridx = tmpidx$idx
paridxi = tmpidx$idxi
paridxa = tmpidx$idxa
garchpars = as.numeric( unlist(coef(filterlist) ) )
# create a temporary environment to store values (deleted at end of function)
tempenvir = new.env(hash = TRUE, parent = .GlobalEnv)
assign("trace", FALSE, envir = tempenvir)
assign("cnames", cnames, envir = tempenvir)
assign("spec", spec, envir = tempenvir)
assign("parallel", parallel, envir = tempenvir)
assign("parallel.control", parallel.control, envir = tempenvir)
assign("paridx", paridx, envir = tempenvir)
assign("paridxi", paridxi, envir = tempenvir)
assign("paridxa", paridxa, envir = tempenvir)
assign("speclist", speclist, envir = tempenvir)
assign("filterlist", filterlist, envir = tempenvir)
omodel = spec@mspec$optimization.model
distribution = spec@mspec$distribution
transformation = spec@mspec$transformation
sig = sigma(filterlist)
res = residuals(filterlist)
#stdresid = res/sig
zres = sapply(filterlist@filter, FUN = function(x) x@filter$z)
ures = matrix(NA, ncol = 5, nrow = dim(zres)[1])
if( maxgarchOrder > 0 ) zres = zres[-(1:maxgarchOrder), , drop = FALSE]
assign("spd.control", spd.control, envir = tempenvir)
# Transformation of z to [0,1] domain
ans = switch(transformation,
parametric = .pparametric.filter(filterlist@filter, zres),
empirical = .pempirical(zres),
spd = .pspd(zres, spd.control))
if(transformation == "spd"){
ures = ans$ures
sfilter = ans$sfit
} else{
ures = ans
sfilter = NULL
}
# make tail adjustments in order to avoid problem with the quantile functions in
# optimization
if(any(ures > 0.99999)){
xn = which(ures > 0.99999)
ures[xn] = 0.99999
}
if(any(ures < eps)){
xn = which(ures < (1.5*eps))
ures[xn] = eps
}
assign("omodel", omodel, envir = tempenvir)
assign("mspec", spec@mspec, envir = tempenvir)
assign("ures", ures, envir = tempenvir)
assign("solver", "solnp", envir = tempenvir)
ipars = .tvcopulastart(data = ures, garchenv = tempenvir)
pars = ipars$pars
assign("npar", length(pars), envir = tempenvir)
fixed = NULL
fixid = NULL
assign("fixed", fixed, envir = tempenvir)
assign("fixid", fixid, envir = tempenvir)
assign(".llh", 1, envir = tempenvir)
mfilter = list()
mfilter$dates = dates
mfilter$date.format = dformat
mfilter$origdata = origdata
mfilter$origdates = origdates
mfilter$vrmodel = vrmodel
# convert back from uniform to standard variates
model = list()
model$copula.distribution = distribution
model$out.sample = out.sample
model$transformation = transformation
model$spd.control = spd.control
model$sfilter = sfilter
model$time.varying = TRUE
# add some tail data of the univariate fits to be used with simulation
mfilter$tailusigma = tail(sigma(filterlist), 25)
mfilter$tailuresids = tail(residuals(filterlist), 25)
mfilter$tailuret = tail(data, 25)
cfilter = copula.tvstudentLLH(pars, data = ures, returnType = "ALL", garchenv = tempenvir)
garchNames = vector(mode = "character")
assetNames = cnames
dccNames = spec@mspec$optimization.model$modelnames
m = length(spec@uspec@spec)
for(i in 1:m){
cf = coef(filterlist@filter[[i]])
cfnames = names(cf)
garchNames = c(garchNames, paste(assetNames[i], cfnames, sep = ".") )
}
allnames = c(garchNames, dccNames)
mfilter$pars = c(as.numeric(coef(filterlist)), pars)
mfilter$garchpars = coef(filterlist)
mfilter$stpars = pars
names(mfilter$pars) = allnames
mfilter$scores = NULL
mfilter$A = NULL
mfilter$cvar = NULL
matcoef = matrix(NA, ncol = 4, nrow = length(mfilter$pars))
matcoef[,1] = mfilter$pars
dimnames(matcoef) = list(allnames, c(" Estimate",
" Std. Error", " t value", "Pr(>|t|)"))
mfilter$matcoef = matcoef
mfilter$paridx = tmpidx
garchllh = sapply(filterlist@filter, FUN = function(x) x@filter$log.likelihoods)
copllh = cfilter$lik
if( maxgarchOrder > 0 ) copllh = c(rep(0, maxgarchOrder), copllh)
mfilter$llh = sum( apply(cbind(copllh, garchllh), 1, "sum") )
mfilter$zres = zres
mfilter$timer = Sys.time() - tic
ans = new("cGARCHfilter",
mfilter = mfilter,
mspec = spec@mspec,
uspec = spec@uspec,
copula = cfilter,
model = model)
return(ans)
}
.cgarchfilter.student = function(spec, data, filter.control = list(n.old = NULL), spd.control = list(lower = 0.1, upper = 0.9, type = "pwm", kernel = "epanech"),
out.sample = 0, parallel = FALSE, parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2),
VAR.fit = NULL, ...)
{
tic = Sys.time()
if( parallel ){
os = .Platform$OS.type
if(is.null(parallel.control$pkg)){
if( os == "windows" ) parallel.control$pkg = "snowfall" else parallel.control$pkg = "multicore"
if( is.null(parallel.control$cores) ) parallel.control$cores = 2
} else{
mtype = match(tolower(parallel.control$pkg[1]), c("multicore", "snowfall"))
if(is.na(mtype)) stop("\nParallel Package type not recognized in parallel.control\n")
parallel.control$pkg = tolower(parallel.control$pkg[1])
if( os == "windows" && parallel.control$pkg == "multicore" ) stop("\nmulticore not supported on windows O/S\n")
if( is.null(parallel.control$cores) ) parallel.control$cores = 2 else parallel.control$cores = as.integer(parallel.control$cores[1])
}
}
n.old = filter.control$n.old
if( is.null( colnames(data) ) ) cnames = paste("Asset_", 1:k, sep = "") else cnames = colnames(data)
xdata = .extractmdata(data)
if(!is.numeric(out.sample)) stop("\ncgarchfilter-->error: out.sample must be numeric\n")
if(as.numeric(out.sample) < 0) stop("\ncgarchfilter-->error: out.sample must be positive\n")
n.start = round(out.sample, 0)
n = dim(xdata$data)[1]
k = dim(xdata$data)[2]
if( (n - n.start) < 100) stop("\ncgarchfilter-->error: function requires at least 100 data\n points to run\n")
data = xdata$data[1:(n - n.start), , drop = FALSE]
dates = xdata$pos[1:(n - n.start)]
origdata = xdata$data
origdates = xdata$pos
dformat = xdata$dformat
vrmodel = NULL
if( spec@mspec$varmodel$VAR ){
if(spec@mspec$varmodel$mxn>0){
ex = spec@mspec$varmodel$external.regressors
if( dim(ex)[1] < dim(data)[1] ) stop("\ndccfilter-->error: external regressor matrix has less points than data!...", call. = FALSE)
orex = ex
ex = ex[1:(n - n.start), , drop = FALSE]
mxn = dim(orex)[2]
} else{
orex = NULL
ex = NULL
mxn = 0
}
if( !is.null(VAR.fit) ){
zdata = VAR.fit$xresiduals
# should be N - p
p = mp = VAR.fit$lag
N = dim(zdata)[1]
K = dim(zdata)[2]
mu = VAR.fit$xfitted
vrmodel = VAR.fit
if(out.sample>0){
Bcoef = VAR.fit$Bcoef
vrmodel2 = varxfilter(X = origdata, p = mp, exogen = orex, Bcoef = Bcoef)
zorigdata = vrmodel2$xresiduals
} else{
zorigdata = zdata
}
} else{
cat("\nestimating VAR model...")
ic = spec@mspec$varmodel$lag.criterion[1]
if( !is.null(spec@mspec$varmodel$lag.max) ){
mp = .varxselect(y = data, lag.max = spec@mspec$varmodel$lag.max, exogen = ex)$selection
mp = as.numeric( mp[which(substr(names(mp), 1, 2) == substr(ic, 1, 2))] )
} else {
mp = spec@mspec$varmodel$lag
}
vrmodel = varxfilter(X = data, p = mp, exogen = ex, Bcoef = NULL, robust = spec@mspec$varmodel$robust,
gamma = spec@mspec$varmodel$robust.control$gamma, delta = spec@mspec$varmodel$robust.control$delta,
nc = spec@mspec$varmodel$robust.control$nc, ns = spec@mspec$varmodel$robust.control$ns)
cat("done!\n")
zdata = vrmodel$xresiduals
# should be N - p
N = dim(zdata)[1]
K = dim(zdata)[2]
mu = vrmodel$xfitted
vrmodel$lag = mp
vrmodel$mxn = mxn
# if we have out of sample data we need to filter
if(out.sample>0){
Bcoef = vrmodel$Bcoef
vrmodel2 = varxfilter(X = origdata, p = mp, exogen = orex, Bcoef = Bcoef)
zorigdata = vrmodel2$xresiduals
} else{
zorigdata = zdata
}
}
} else{
zdata = data
zorigdata = origdata
orex = NULL
ex = NULL
}
if( is(spec@uspec, "uGARCHmultispec") ){
speclist = spec@uspec
} else{
stop("\ncgarchfilter-->error: cGARCHspec must contain multispec object with fixed parameters")
}
maxgarchOrder = spec@mspec$optimization.model$maxgarchOrder
filterlist = multifilter(multifitORspec = spec@uspec, data = as.data.frame(zorigdata), out.sample = out.sample,
parallel = parallel, parallel.control = parallel.control, n.old = n.old)
assign(".filterlist", filterlist, envir = .GlobalEnv)
garchpars = as.numeric( unlist(coef(filterlist) ) )
tempenvir = new.env(hash = TRUE, parent = .GlobalEnv)
assign("trace", FALSE, envir = tempenvir)
assign("method", spec@mspec$method, envir = tempenvir)
assign("cnames", cnames, envir = tempenvir)
assign("spec", spec, envir = tempenvir)
assign("parallel", parallel, envir = tempenvir)
assign("parallel.control", parallel.control, envir = tempenvir)
#assign("paridx", paridx, envir = tempenvir)
#assign("paridxi", paridxi, envir = tempenvir)
#assign("paridxa", paridxa, envir = tempenvir)
assign("speclist", speclist, envir = tempenvir)
assign("filterlist", filterlist, envir = tempenvir)
omodel = spec@mspec$optimization.model
distribution = spec@mspec$distribution
transformation = spec@mspec$transformation
sig = sigma(filterlist)
res = residuals(filterlist)
#stdresid = res/sig
assign("transformation", transformation, envir = tempenvir)
zres = sapply(filterlist@filter, FUN = function(x) x@filter$z)
ures = matrix(NA, ncol = 5, nrow = dim(zres)[1])
if( maxgarchOrder > 0 ) zres = zres[-(1:maxgarchOrder), , drop = FALSE]
assign("spd.control", spd.control, envir = tempenvir)
# Transformation of z to [0,1] domain
ans = switch(transformation,
parametric = .pparametric.filter(filterlist@filter, zres),
empirical = .pempirical(zres),
spd = .pspd(zres, spd.control))
if(transformation == "spd"){
ures = ans$ures
sfilter = ans$sfit
} else{
ures = ans
sfilter = NULL
}
# make tail adjustments in order to avoid problem with the quantile functions in
# optimization
if(any(ures > 0.99999)){
xn = which(ures > 0.99999)
ures[xn] = 0.99999
}
if(any(ures < eps)){
xn = which(ures < (1.5*eps))
ures[xn] = eps
}
assign("omodel", omodel, envir = tempenvir)
assign("mspec", spec@mspec, envir = tempenvir)
assign("ures", ures, envir = tempenvir)
ipars = .copulastart(data = ures, garchenv = tempenvir)
pars = ipars$pars
assign("npar", length(pars), envir = tempenvir)
# no fixed parameters, all normal
fixed = NULL
fixid = NULL
use.solver = TRUE
assign("fixed", fixed, envir = tempenvir)
assign("fixid", fixid, envir = tempenvir)
assign(".llh", 1, envir = tempenvir)
mfilter = list()
mfilter$dates = dates
mfilter$date.format = dformat
mfilter$origdata = origdata
mfilter$origdates = origdates
mfilter$vrmodel = vrmodel
# convert back from uniform to standard variates
model = list()
model$copula.distribution = distribution
model$out.sample = out.sample
model$transformation = transformation
model$spd.control = spd.control
model$sfit = sfilter
model$time.varying = FALSE
# add some tail data of the univariate fits to be used with simulation
mfilter$tailusigma = tail(sigma(filterlist), 25)
mfilter$tailuresids = tail(residuals(filterlist), 25)
mfilter$tailuret = tail(data, 25)
cfilter = copula.studentLLH(pars, data = ures, returnType = "ALL", garchenv = tempenvir)
garchNames = vector(mode = "character")
assetNames = cnames
dccNames = spec@mspec$optimization.model$modelnames
m = length(spec@uspec@spec)
for(i in 1:m){
cf = coef(filterlist@filter[[i]])
cfnames = names(cf)
garchNames = c(garchNames, paste(assetNames[i], cfnames, sep = ".") )
}
allnames = c(garchNames, dccNames)
mfilter$pars = c(as.numeric(coef(filterlist)), pars)
names(mfilter$pars) = allnames
mfilter$garchpars = coef(filterlist)
mfilter$stpars = pars
mfilter$scores = NULL
mfilter$A = NULL
mfilter$cvar = NULL
matcoef = matrix(NA, ncol = 4, nrow = length(mfilter$pars))
matcoef[,1] = mfilter$pars
dimnames(matcoef) = list(allnames, c(" Estimate",
" Std. Error", " t value", "Pr(>|t|)"))
mfilter$matcoef = matcoef
# remember that the filter returns the correct sign for the log.likelihoods (unlike fit which returns
# the negative of.
garchllh = sapply(filterlist@filter, FUN = function(x) x@filter$log.likelihoods)
copllh = cfilter$lik
if( maxgarchOrder > 0 ) copllh = c(rep(0, maxgarchOrder), copllh)
mfilter$llh = sum( apply(cbind(copllh, garchllh), 1, "sum") )
mfilter$zres = zres
mfilter$timer = Sys.time() - tic
ans = new("cGARCHfilter",
mfilter = mfilter,
mspec = spec@mspec,
uspec = spec@uspec,
copula = cfilter,
model = model)
return(ans)
}
.cgarchfilter.tvnorm = function(spec, data, filter.control = list(n.old = NULL), spd.control = list(lower = 0.1, upper = 0.9, type = "pwm", kernel = "epanech"),
out.sample = 0, parallel = FALSE, parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2),
VAR.fit = NULL, ...)
{
tic = Sys.time()
if( parallel ){
os = .Platform$OS.type
if(is.null(parallel.control$pkg)){
if( os == "windows" ) parallel.control$pkg = "snowfall" else parallel.control$pkg = "multicore"
if( is.null(parallel.control$cores) ) parallel.control$cores = 2
} else{
mtype = match(tolower(parallel.control$pkg[1]), c("multicore", "snowfall"))
if(is.na(mtype)) stop("\nParallel Package type not recognized in parallel.control\n")
parallel.control$pkg = tolower(parallel.control$pkg[1])
if( os == "windows" && parallel.control$pkg == "multicore" ) stop("\nmulticore not supported on windows O/S\n")
if( is.null(parallel.control$cores) ) parallel.control$cores = 2 else parallel.control$cores = as.integer(parallel.control$cores[1])
}
}
n.old = filter.control$n.old
if( is.null( colnames(data) ) ) cnames = paste("Asset_", 1:k, sep = "") else cnames = colnames(data)
xdata = .extractmdata(data)
if(!is.numeric(out.sample)) stop("\ncgarchfilter-->error: out.sample must be numeric\n")
if(as.numeric(out.sample) < 0) stop("\ncgarchfilter-->error: out.sample must be positive\n")
n.start = round(out.sample, 0)
n = dim(xdata$data)[1]
k = dim(xdata$data)[2]
if( (n - n.start) < 100) stop("\ncgarchfilter-->error: function requires at least 100 data\n points to run\n")
data = xdata$data[1:(n - n.start), , drop = FALSE]
dates = xdata$pos[1:(n - n.start)]
origdata = xdata$data
origdates = xdata$pos
dformat = xdata$dformat
vrmodel = NULL
if( spec@mspec$varmodel$VAR ){
if(spec@mspec$varmodel$mxn>0){
ex = spec@mspec$varmodel$external.regressors
if( dim(ex)[1] < dim(data)[1] ) stop("\ndccfilter-->error: external regressor matrix has less points than data!...", call. = FALSE)
orex = ex
ex = ex[1:(n - n.start), , drop = FALSE]
mxn = dim(orex)[2]
} else{
orex = NULL
ex = NULL
mxn = 0
}
if( !is.null(VAR.fit) ){
zdata = VAR.fit$xresiduals
# should be N - p
p = mp = VAR.fit$lag
N = dim(zdata)[1]
K = dim(zdata)[2]
mu = VAR.fit$xfitted
vrmodel = VAR.fit
if(out.sample>0){
Bcoef = VAR.fit$Bcoef
vrmodel2 = varxfilter(X = origdata, p = mp, exogen = orex, Bcoef = Bcoef)
zorigdata = vrmodel2$xresiduals
} else{
zorigdata = zdata
}
} else{
cat("\nestimating VAR model...")
ic = spec@mspec$varmodel$lag.criterion[1]
if( !is.null(spec@mspec$varmodel$lag.max) ){
mp = .varxselect(y = data, lag.max = spec@mspec$varmodel$lag.max, exogen = ex)$selection
mp = as.numeric( mp[which(substr(names(mp), 1, 2) == substr(ic, 1, 2))] )
} else {
mp = spec@mspec$varmodel$lag
}
vrmodel = varxfilter(X = data, p = mp, exogen = ex, Bcoef = NULL, robust = spec@mspec$varmodel$robust,
gamma = spec@mspec$varmodel$robust.control$gamma, delta = spec@mspec$varmodel$robust.control$delta,
nc = spec@mspec$varmodel$robust.control$nc, ns = spec@mspec$varmodel$robust.control$ns)
cat("done!\n")
zdata = vrmodel$xresiduals
# should be N - p
N = dim(zdata)[1]
K = dim(zdata)[2]
mu = vrmodel$xfitted
vrmodel$lag = mp
vrmodel$mxn = mxn
# if we have out of sample data we need to filter
if(out.sample>0){
Bcoef = vrmodel$Bcoef
vrmodel2 = varxfilter(X = origdata, p = mp, exogen = orex, Bcoef = Bcoef)
zorigdata = vrmodel2$xresiduals
} else{
zorigdata = zdata
}
}
} else{
zdata = data
zorigdata = origdata
orex = NULL
ex = NULL
}
if( is(spec@uspec, "uGARCHmultispec") ){
speclist = spec@uspec
} else{
stop("\ncgarchfilter-->error: cGARCHspec must contain multispec object with fixed parameters")
}
maxgarchOrder = spec@mspec$optimization.model$maxgarchOrder
filterlist = multifilter(multifitORspec = spec@uspec, data = as.data.frame(zorigdata), out.sample = out.sample,
parallel = parallel, parallel.control = parallel.control, n.old = n.old)
tmpidx = .dccparidx( speclist )
paridx = tmpidx$idx
paridxi = tmpidx$idxi
paridxa = tmpidx$idxa
garchpars = as.numeric( unlist(coef(filterlist) ) )
# create a temporary environment to store values (deleted at end of function)
tempenvir = new.env(hash = TRUE, parent = .GlobalEnv)
assign("trace", FALSE, envir = tempenvir)
assign("cnames", cnames, envir = tempenvir)
assign("spec", spec, envir = tempenvir)
assign("parallel", parallel, envir = tempenvir)
assign("parallel.control", parallel.control, envir = tempenvir)
assign("paridx", paridx, envir = tempenvir)
assign("paridxi", paridxi, envir = tempenvir)
assign("paridxa", paridxa, envir = tempenvir)
assign("speclist", speclist, envir = tempenvir)
assign("filterlist", filterlist, envir = tempenvir)
omodel = spec@mspec$optimization.model
distribution = spec@mspec$distribution
transformation = spec@mspec$transformation
sig = sigma(filterlist)
res = residuals(filterlist)
#stdresid = res/sig
zres = sapply(filterlist@filter, FUN = function(x) x@filter$z)
ures = matrix(NA, ncol = 5, nrow = dim(zres)[1])
if( maxgarchOrder > 0 ) zres = zres[-(1:maxgarchOrder), , drop = FALSE]
assign("spd.control", spd.control, envir = tempenvir)
# Transformation of z to [0,1] domain
ans = switch(transformation,
parametric = .pparametric.filter(filterlist@filter, zres),
empirical = .pempirical(zres),
spd = .pspd(zres, spd.control))
if(transformation == "spd"){
ures = ans$ures
sfilter = ans$sfit
} else{
ures = ans
sfilter = NULL
}
# make tail adjustments in order to avoid problem with the quantile functions in
# optimization
if(any(ures > 0.99999)){
xn = which(ures > 0.99999)
ures[xn] = 0.99999
}
if(any(ures < eps)){
xn = which(ures < (1.5*eps))
ures[xn] = eps
}
assign("omodel", omodel, envir = tempenvir)
assign("mspec", spec@mspec, envir = tempenvir)
assign("ures", ures, envir = tempenvir)
assign("solver", "solnp", envir = tempenvir)
ipars = .tvcopulastart(data = ures, garchenv = tempenvir)
pars = ipars$pars
assign("npar", length(pars), envir = tempenvir)
fixed = NULL
fixid = NULL
assign("fixed", fixed, envir = tempenvir)
assign("fixid", fixid, envir = tempenvir)
assign(".llh", 1, envir = tempenvir)
mfilter = list()
mfilter$dates = dates
mfilter$date.format = dformat
mfilter$origdata = origdata
mfilter$origdates = origdates
mfilter$vrmodel = vrmodel
# convert back from uniform to standard variates
model = list()
model$copula.distribution = distribution
model$out.sample = out.sample
model$transformation = transformation
model$spd.control = spd.control
model$sfilter = sfilter
model$time.varying = TRUE
# add some tail data of the univariate fits to be used with simulation
mfilter$tailusigma = tail(sigma(filterlist), 25)
mfilter$tailuresids = tail(residuals(filterlist), 25)
mfilter$tailuret = tail(data, 25)
cfilter = copula.tvnormalLLH(pars, data = ures, returnType = "ALL", garchenv = tempenvir)
garchNames = vector(mode = "character")
assetNames = cnames
dccNames = spec@mspec$optimization.model$modelnames
m = length(spec@uspec@spec)
for(i in 1:m){
cf = coef(filterlist@filter[[i]])
cfnames = names(cf)
garchNames = c(garchNames, paste(assetNames[i], cfnames, sep = ".") )
}
allnames = c(garchNames, dccNames)
mfilter$pars = c(as.numeric(coef(filterlist)), pars)
mfilter$garchpars = coef(filterlist)
mfilter$stpars = pars
names(mfilter$pars) = allnames
mfilter$scores = NULL
mfilter$A = NULL
mfilter$cvar = NULL
matcoef = matrix(NA, ncol = 4, nrow = length(mfilter$pars))
matcoef[,1] = mfilter$pars
dimnames(matcoef) = list(allnames, c(" Estimate",
" Std. Error", " t value", "Pr(>|t|)"))
mfilter$matcoef = matcoef
mfilter$paridx = tmpidx
garchllh = sapply(filterlist@filter, FUN = function(x) x@filter$log.likelihoods)
copllh = cfilter$lik
if( maxgarchOrder > 0 ) copllh = c(rep(0, maxgarchOrder), copllh)
mfilter$llh = sum( apply(cbind(copllh, garchllh), 1, "sum") )
mfilter$zres = zres
mfilter$timer = Sys.time() - tic
ans = new("cGARCHfilter",
mfilter = mfilter,
mspec = spec@mspec,
uspec = spec@uspec,
copula = cfilter,
model = model)
return(ans)
}
.cgarchfilter.norm = function(spec, data, filter.control = list(n.old = NULL), spd.control = list(lower = 0.1, upper = 0.9, type = "pwm", kernel = "epanech"),
out.sample = 0, parallel = FALSE, parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2),
VAR.fit = NULL, ...)
{
tic = Sys.time()
if( parallel ){
os = .Platform$OS.type
if(is.null(parallel.control$pkg)){
if( os == "windows" ) parallel.control$pkg = "snowfall" else parallel.control$pkg = "multicore"
if( is.null(parallel.control$cores) ) parallel.control$cores = 2
} else{
mtype = match(tolower(parallel.control$pkg[1]), c("multicore", "snowfall"))
if(is.na(mtype)) stop("\nParallel Package type not recognized in parallel.control\n")
parallel.control$pkg = tolower(parallel.control$pkg[1])
if( os == "windows" && parallel.control$pkg == "multicore" ) stop("\nmulticore not supported on windows O/S\n")
if( is.null(parallel.control$cores) ) parallel.control$cores = 2 else parallel.control$cores = as.integer(parallel.control$cores[1])
}
}
n.old = filter.control$n.old
if( is.null( colnames(data) ) ) cnames = paste("Asset_", 1:k, sep = "") else cnames = colnames(data)
xdata = .extractmdata(data)
if(!is.numeric(out.sample)) stop("\ncgarchfilter-->error: out.sample must be numeric\n")
if(as.numeric(out.sample) < 0) stop("\ncgarchfilter-->error: out.sample must be positive\n")
n.start = round(out.sample, 0)
n = dim(xdata$data)[1]
k = dim(xdata$data)[2]
if( (n - n.start) < 100) stop("\ncgarchfilter-->error: function requires at least 100 data\n points to run\n")
data = xdata$data[1:(n - n.start), , drop = FALSE]
dates = xdata$pos[1:(n - n.start)]
origdata = xdata$data
origdates = xdata$pos
dformat = xdata$dformat
vrmodel = NULL
if( spec@mspec$varmodel$VAR ){
if(spec@mspec$varmodel$mxn>0){
ex = spec@mspec$varmodel$external.regressors
if( dim(ex)[1] < dim(data)[1] ) stop("\ndccfilter-->error: external regressor matrix has less points than data!...", call. = FALSE)
orex = ex
ex = ex[1:(n - n.start), , drop = FALSE]
mxn = dim(orex)[2]
} else{
orex = NULL
ex = NULL
mxn = 0
}
if( !is.null(VAR.fit) ){
zdata = VAR.fit$xresiduals
# should be N - p
p = mp = VAR.fit$lag
N = dim(zdata)[1]
K = dim(zdata)[2]
mu = VAR.fit$xfitted
vrmodel = VAR.fit
if(out.sample>0){
Bcoef = VAR.fit$Bcoef
vrmodel2 = varxfilter(X = origdata, p = mp, exogen = orex, Bcoef = Bcoef)
zorigdata = vrmodel2$xresiduals
} else{
zorigdata = zdata
}
} else{
cat("\nestimating VAR model...")
ic = spec@mspec$varmodel$lag.criterion[1]
if( !is.null(spec@mspec$varmodel$lag.max) ){
mp = .varxselect(y = data, lag.max = spec@mspec$varmodel$lag.max, exogen = ex)$selection
mp = as.numeric( mp[which(substr(names(mp), 1, 2) == substr(ic, 1, 2))] )
} else {
mp = spec@mspec$varmodel$lag
}
vrmodel = varxfilter(X = data, p = mp, exogen = ex, Bcoef = NULL, robust = spec@mspec$varmodel$robust,
gamma = spec@mspec$varmodel$robust.control$gamma, delta = spec@mspec$varmodel$robust.control$delta,
nc = spec@mspec$varmodel$robust.control$nc, ns = spec@mspec$varmodel$robust.control$ns)
cat("done!\n")
zdata = vrmodel$xresiduals
# should be N - p
N = dim(zdata)[1]
K = dim(zdata)[2]
mu = vrmodel$xfitted
vrmodel$lag = mp
vrmodel$mxn = mxn
# if we have out of sample data we need to filter
if(out.sample>0){
Bcoef = vrmodel$Bcoef
vrmodel2 = varxfilter(X = origdata, p = mp, exogen = orex, Bcoef = Bcoef)
zorigdata = vrmodel2$xresiduals
} else{
zorigdata = zdata
}
}
} else{
zdata = data
zorigdata = origdata
orex = NULL
ex = NULL
}
if( is(spec@uspec, "uGARCHmultispec") ){
speclist = spec@uspec
} else{
stop("\ncgarchfilter-->error: cGARCHspec must contain multispec object with fixed parameters")
}
maxgarchOrder = spec@mspec$optimization.model$maxgarchOrder
filterlist = multifilter(multifitORspec = spec@uspec, data = as.data.frame(zorigdata), out.sample = out.sample,
parallel = parallel, parallel.control = parallel.control, n.old = n.old)
assign(".filterlist", filterlist, envir = .GlobalEnv)
garchpars = as.numeric( unlist(coef(filterlist) ) )
tempenvir = new.env(hash = TRUE, parent = .GlobalEnv)
assign("trace", FALSE, envir = tempenvir)
assign("method", spec@mspec$method, envir = tempenvir)
assign("cnames", cnames, envir = tempenvir)
assign("spec", spec, envir = tempenvir)
assign("parallel", parallel, envir = tempenvir)
assign("parallel.control", parallel.control, envir = tempenvir)
#assign("paridx", paridx, envir = tempenvir)
#assign("paridxi", paridxi, envir = tempenvir)
#assign("paridxa", paridxa, envir = tempenvir)
assign("speclist", speclist, envir = tempenvir)
assign("filterlist", filterlist, envir = tempenvir)
omodel = spec@mspec$optimization.model
distribution = spec@mspec$distribution
transformation = spec@mspec$transformation
sig = sigma(filterlist)
res = residuals(filterlist)
#stdresid = res/sig
assign("transformation", transformation, envir = tempenvir)
zres = sapply(filterlist@filter, FUN = function(x) x@filter$z)
ures = matrix(NA, ncol = 5, nrow = dim(zres)[1])
if( maxgarchOrder > 0 ) zres = zres[-(1:maxgarchOrder), , drop = FALSE]
assign("spd.control", spd.control, envir = tempenvir)
# Transformation of z to [0,1] domain
ans = switch(transformation,
parametric = .pparametric.filter(filterlist@filter, zres),
empirical = .pempirical(zres),
spd = .pspd(zres, spd.control))
if(transformation == "spd"){
ures = ans$ures
sfilter = ans$sfit
} else{
ures = ans
sfilter = NULL
}
# make tail adjustments in order to avoid problem with the quantile functions in
# optimization
if(any(ures > 0.99999)){
xn = which(ures > 0.99999)
ures[xn] = 0.99999
}
if(any(ures < eps)){
xn = which(ures < (1.5*eps))
ures[xn] = eps
}
assign("omodel", omodel, envir = tempenvir)
assign("mspec", spec@mspec, envir = tempenvir)
assign("ures", ures, envir = tempenvir)
ipars = .copulastart(data = ures, garchenv = tempenvir)
pars = ipars$pars
assign("npar", length(pars), envir = tempenvir)
# no fixed parameters, all normal
fixed = NULL
fixid = NULL
use.solver = TRUE
assign("fixed", fixed, envir = tempenvir)
assign("fixid", fixid, envir = tempenvir)
assign(".llh", 1, envir = tempenvir)
mfilter = list()
mfilter$dates = dates
mfilter$date.format = dformat
mfilter$origdata = origdata
mfilter$origdates = origdates
mfilter$vrmodel = vrmodel
# convert back from uniform to standard variates
model = list()
model$copula.distribution = distribution
model$out.sample = out.sample
model$transformation = transformation
model$spd.control = spd.control
model$sfit = sfilter
model$time.varying = FALSE
# add some tail data of the univariate fits to be used with simulation
mfilter$tailusigma = tail(sigma(filterlist), 25)
mfilter$tailuresids = tail(residuals(filterlist), 25)
mfilter$tailuret = tail(data, 25)
cfilter = copula.normalLLH(pars, data = ures, returnType = "ALL", garchenv = tempenvir)
garchNames = vector(mode = "character")
assetNames = cnames
dccNames = spec@mspec$optimization.model$modelnames
m = length(spec@uspec@spec)
for(i in 1:m){
cf = coef(filterlist@filter[[i]])
cfnames = names(cf)
garchNames = c(garchNames, paste(assetNames[i], cfnames, sep = ".") )
}
allnames = c(garchNames, dccNames)
mfilter$pars = c(as.numeric(coef(filterlist)), pars)
names(mfilter$pars) = allnames
mfilter$garchpars = coef(filterlist)
mfilter$stpars = pars
mfilter$scores = NULL
mfilter$A = NULL
mfilter$cvar = NULL
matcoef = matrix(NA, ncol = 4, nrow = length(mfilter$pars))
matcoef[,1] = mfilter$pars
dimnames(matcoef) = list(allnames, c(" Estimate",
" Std. Error", " t value", "Pr(>|t|)"))
mfilter$matcoef = matcoef
# remember that the filter returns the correct sign for the log.likelihoods (unlike fit which returns
# the negative of.
garchllh = sapply(filterlist@filter, FUN = function(x) x@filter$log.likelihoods)
copllh = cfilter$lik
if( maxgarchOrder > 0 ) copllh = c(rep(0, maxgarchOrder), copllh)
mfilter$llh = sum( apply(cbind(copllh, garchllh), 1, "sum") )
mfilter$zres = zres
mfilter$timer = Sys.time() - tic
ans = new("cGARCHfilter",
mfilter = mfilter,
mspec = spec@mspec,
uspec = spec@uspec,
copula = cfilter,
model = model)
return(ans)
}
.cgarchsim = function(fit, n.sim = 1000, n.start = 0, m.sim = 1, startMethod = c("unconditional", "sample"),
presigma = NULL, preresiduals = NULL, prereturns = NULL, preR = NULL, rseed = NULL, mexsimdata = NULL,
vexsimdata = NULL, parallel = FALSE, parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2),
VAR.fit = NULL, ...)
{
dist = fit@model$copula.distribution
tv = ifelse(fit@model$time.varying, 1, 0)
dx = paste(dist, tv, sep = "")
ans = switch(dx,
mvt0 = .cgarchsim1(fit = fit, n.sim = n.sim, n.start = n.start, m.sim = m.sim, startMethod = startMethod,
presigma = presigma, preresiduals = preresiduals, prereturns = prereturns, preR = preR, rseed = rseed,
mexsimdata = mexsimdata, vexsimdata = vexsimdata, parallel = parallel, parallel.control = parallel.control,
VAR.fit = VAR.fit, ...),
mvt1 = .cgarchsim2(fit = fit, n.sim = n.sim, n.start = n.start, m.sim = m.sim, startMethod = startMethod,
presigma = presigma, preresiduals = preresiduals, prereturns = prereturns, preR = preR, rseed = rseed,
mexsimdata = mexsimdata, vexsimdata = vexsimdata, parallel = parallel, parallel.control = parallel.control,
VAR.fit = VAR.fit, ...),
mvnorm0 = .cgarchsim1(fit = fit, n.sim = n.sim, n.start = n.start, m.sim = m.sim, startMethod = startMethod,
presigma = presigma, preresiduals = preresiduals, prereturns = prereturns, preR = preR, rseed = rseed,
mexsimdata = mexsimdata, vexsimdata = vexsimdata, parallel = parallel, parallel.control = parallel.control,
VAR.fit = VAR.fit, ...),
mvnorm1 = .cgarchsim2(fit = fit, n.sim = n.sim, n.start = n.start, m.sim = m.sim, startMethod = startMethod,
presigma = presigma, preresiduals = preresiduals, prereturns = prereturns, preR = preR, rseed = rseed,
mexsimdata = mexsimdata, vexsimdata = vexsimdata, parallel = parallel, parallel.control = parallel.control,
VAR.fit = VAR.fit, ...))
return( ans )
}
setMethod("cgarchsim", signature(fit = "cGARCHfit"), .cgarchsim)
# non time varying copula simulation
.cgarchsim1 = function(fit, n.sim = 1000, n.start = 0, m.sim = 1, startMethod = c("unconditional", "sample"),
presigma = NULL, preresiduals = NULL, prereturns = NULL, preR = NULL, rseed = NULL, mexsimdata = NULL,
vexsimdata = NULL, parallel = FALSE, parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2),
VAR.fit = NULL, ...)
{
# first generate the copula random uniform numbers (static copula)
# ures --> transform --> zres
# pass zres to GARCH simulation
# pass if needed to VAR model
if( is.null(rseed) ){
rseed = as.integer(runif(m.sim, 1, Sys.time()))
} else {
if(length(rseed) == 1) c(rseed, as.integer(runif(m.sim-1, 1, Sys.time()))) else rseed = as.integer( rseed[1:m.sim] )
}
Data = head(fit@mfit$origdata, dim(fit@mfit$origdata)[1] - fit@model$out.sample)
m = dim(Data)[2]
n = dim(Data)[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])
}
}
nsim = n.sim + n.start
spec = fit@uspec
cf = coef(fit, "garch")
cfl = vector(mode = "list", length = m)
for(i in 1:m){
if(spec@type == "equal"){
spec@spec[[i]]@optimization.model$fixed.pars = as.list(cf[,i])
cfl[[i]] = cf[,i]
} else{
spec@spec[[i]]@optimization.model$fixed.pars = as.list(cf[[i]])
cfl[[i]] = cf[[i]]
}
}
if(is.null(preR)){
Rbar = fit@copula$Rt
} else{
if(dim(preR)[1] != m) stop("\ncgarchsim-->error: wrong dimension for preR\n")
if(dim(preR)[2] != m) stop("\ncgarchsim-->error: wrong dimension for preR\n")
if(any(diag(preR) != 1)) stop("\ncgarchsim-->error: preR diagonals must be 1.\n")
Rbar = preR
}
dist = fit@model$copula.distribution
dx = paste(dist, 0, sep = "")
ures = .sample.copula(fit, Rbar, dist = dx, n.sim = n.sim, n.start = n.start, m.sim = m.sim, rseed = rseed)
# ures = [0,1] copula random numbers
# now transform back into margins for use in garch
transformation = fit@model$transformation
include.skew = sapply(fit@uspec@spec, FUN = function(x) x@distribution.model$include.skew)
include.shape = sapply(fit@uspec@spec, FUN = function(x) x@distribution.model$include.shape)
distx = sapply(fit@uspec@spec, FUN = function(x) x@distribution.model$distribution)
zres = array(NA, dim = c(nsim, m, m.sim))
# parallel:
if( parallel ){
if( parallel.control$pkg == "multicore" ){
mtmp = mclapply(as.list(1:m.sim), FUN = function(i) switch(transformation,
parametric = .qparametric(ures[,,i], ucoef = cfl, include.skew, include.shape, dist = distx),
empirical = .qempirical(matrix(ures[,,i], ncol = m), fit@mfit$zres),
spd = .qspd(ures[,,i], sfit = fit@model$sfit)), mc.cores = parallel.control$cores)
for(i in 1:m.sim) zres[,,i] = mtmp[[i]]
} else{
ssfit = fit@model$sfit
sfInit(parallel = TRUE, cpus = parallel.control$cores)
sfExport("ures", "cfl", "include.skew", "include.shape", "distx", "ssfit", "transformation", local = TRUE)
mtmp = sfLapply(as.list(1:m.sim), fun = function(i)
switch(transformation,
parametric = rgarch:::.qparametric(ures[,,i], ucoef = cfl, include.skew, include.shape, dist = distx),
empirical = rgarch:::.qempirical(matrix(ures[,,i], ncol = m), fit@mfit$zres),
spd = rgarch:::.qspd(matrix(ures[,,i], ncol = m), sfit = ssfit)))
for(i in 1:m.sim) zres[,,i] = mtmp[[i]]
#sfStop()
}
} else{
mtmp = lapply(as.list(1:m.sim), FUN = function(i) switch(transformation,
parametric = .qparametric(matrix(ures[,,i], ncol = m), ucoef = cfl, include.skew, include.shape, dist = distx),
empirical = .qempirical(matrix(ures[,,i], ncol = m), fit@mfit$zres),
spd = .qspd(matrix(ures[,,i], ncol = m), sfit = fit@model$sfit)))
for(i in 1:m.sim) zres[,,i] = mtmp[[i]]
}
# Now we have the innovations which we feed back into the garch routine to simulate
simRes = simX = simR = simQ = simH = simSeries = vector(mode = "list", length = m.sim)
if( is.null(fit@mfit$vrmodel) ){
if(is.null(prereturns)) prereturns = fit@mfit$tailuret
}
if( parallel ){
if( parallel.control$pkg == "multicore" ){
simlist = mclapply(as.list(1:m), FUN = function(i){
maxx = spec@spec[[i]]@optimization.model$maxOrder2;
htmp = ugarchpath(spec@spec[[i]], n.sim = n.sim + n.start, n.start = 0, m.sim = m.sim,
custom.dist = list(name = "sample", distfit = matrix(zres[,i,1:m.sim], ncol = m.sim)),
presigma = if( is.null(presigma) ) tail(fit@mfit$tailusigma[,i], maxx) else tail(presigma[,i], maxx),
preresiduals = if( is.null(preresiduals) ) tail(fit@mfit$tailuresids[,i], maxx) else tail(preresiduals[,i], maxx),
prereturns = if(is.null(fit@mfit$vrmodel)) tail(prereturns[,i], maxx) else NA,
mexsimdata = if( is.null(fit@mfit$vrmodel) ) mexsimdata[[i]] else NULL, vexsimdata = vexsimdata[[i]] )
h = matrix(tail(htmp@path$sigmaSim^2, n.sim), nrow = n.sim)
x = matrix(htmp@path$seriesSim, nrow = n.sim + n.start)
return(list(h = h, x = x))}, mc.cores = parallel.control$cores)
H = array(NA, dim = c(n.sim, m, m.sim))
tmpH = array(NA, dim = c(m, m, n.sim))
for(i in 1:n.sim) H[i,,] = t(sapply(simlist, FUN = function(x) as.numeric(x$h[i,])))
for(i in 1:m.sim){
for(j in 1:n.sim){
tmpH[ , , j] = diag(sqrt( H[j, , i]) ) %*% fit@copula$Rt %*% diag(sqrt( H[j, , i] ) )
}
simH[[i]] = tmpH
}
simxX = array(NA, dim = c(n.sim, m, m.sim))
for(i in 1:n.sim) simxX[i,,] = t(sapply(simlist, FUN = function(x) as.numeric(x$x[i,])))
simX = vector(mode = "list", length = m.sim)
for(i in 1:m.sim) simX[[i]] = matrix(simxX[,,i], nrow = n.sim)
} else{
if(is.null(fit@mfit$vrmodel)) include.vrmodel = FALSE else include.vrmodel = TRUE
sfExport("spec", "n.sim", "n.start", "m.sim", "startMethod", "zres", "presigma",
"preresiduals", "prereturns", "fit", "mexsimdata", "vexsimdata", "include.vrmodel", local = TRUE)
simlist = sfLapply(as.list(1:m), fun = function(i){
maxx = spec@spec[[i]]@optimization.model$maxOrder2;
htmp = rgarch::ugarchpath(spec@spec[[i]], n.sim = n.sim + n.start, n.start = 0, m.sim = m.sim,
custom.dist = list(name = "sample", distfit = matrix(zres[,i,1:m.sim], ncol = m.sim)),
presigma = if( is.null(presigma) ) tail(fit@mfit$tailusigma[,i], maxx) else tail(presigma[,i], maxx),
preresiduals = if( is.null(preresiduals) ) tail(fit@mfit$tailuresids[,i], maxx) else tail(preresiduals[,i], maxx),
prereturns = if(is.null(fit@mfit$vrmodel)) tail(prereturns[,i], maxx) else NA,
mexsimdata = if( is.null(fit@mfit$vrmodel) ) mexsimdata[[i]] else NULL, vexsimdata = vexsimdata[[i]] )
h = matrix(tail(htmp@path$sigmaSim^2, n.sim), nrow = n.sim)
x = matrix(htmp@path$seriesSim, nrow = n.sim + n.start)
return(list(h = h, x = x))})
if( !include.vrmodel ) sfStop()
H = array(NA, dim = c(n.sim, m, m.sim))
tmpH = array(NA, dim = c(m, m, n.sim))
for(i in 1:n.sim) H[i,,] = t(sapply(simlist, FUN = function(x) as.numeric(x$h[i,])))
for(i in 1:m.sim){
for(j in 1:n.sim){
tmpH[ , , j] = diag(sqrt( H[j, , i]) ) %*% fit@copula$Rt %*% diag(sqrt( H[j, , i] ) )
}
simH[[i]] = tmpH
}
simxX = array(NA, dim = c(n.sim, m, m.sim))
for(i in 1:n.sim) simxX[i,,] = t(sapply(simlist, FUN = function(x) as.numeric(x$x[i,])))
simX = vector(mode = "list", length = m.sim)
for(i in 1:m.sim) simX[[i]] = matrix(simxX[,,i], nrow = n.sim)
}
} else{
simlist = lapply(as.list(1:m), FUN = function(i){
maxx = spec@spec[[i]]@optimization.model$maxOrder2;
htmp = ugarchpath(spec@spec[[i]], n.sim = n.sim + n.start, n.start = 0, m.sim = m.sim,
custom.dist = list(name = "sample", distfit = matrix(zres[,i,1:m.sim], ncol = m.sim)),
presigma = if( is.null(presigma) ) tail(fit@mfit$tailusigma[,i], maxx) else tail(presigma[,i], maxx),
preresiduals = if( is.null(preresiduals) ) tail(fit@mfit$tailuresids[,i], maxx) else tail(preresiduals[,i], maxx),
prereturns = if(is.null(fit@mfit$vrmodel)) tail(prereturns[,i], maxx) else NA,
mexsimdata = if( is.null(fit@mfit$vrmodel) ) mexsimdata[[i]] else NULL, vexsimdata = vexsimdata[[i]] )
h = matrix(tail(htmp@path$sigmaSim^2, n.sim), nrow = n.sim)
x = matrix(htmp@path$seriesSim, nrow = n.sim + n.start)
return(list(h = h, x = x))})
H = array(NA, dim = c(n.sim, m, m.sim))
tmpH = array(NA, dim = c(m, m, n.sim))
for(i in 1:n.sim) H[i,,] = t(sapply(simlist, FUN = function(x) as.numeric(x$h[i,])))
for(i in 1:m.sim){
for(j in 1:n.sim){
tmpH[ , , j] = diag(sqrt( H[j, , i]) ) %*% fit@copula$Rt %*% diag(sqrt( H[j, , i] ) )
}
simH[[i]] = tmpH
}
simxX = array(NA, dim = c(n.sim, m, m.sim))
for(i in 1:n.sim) simxX[i,,] = t(sapply(simlist, FUN = function(x) as.numeric(x$x[i,])))
simX = vector(mode = "list", length = m.sim)
for(i in 1:m.sim) simX[[i]] = matrix(simxX[,,i], nrow = n.sim)
}
if( !is.null(fit@mfit$vrmodel) ){
# mexsimdata check
mxn = fit@mfit$vrmodel$mxn
if(mxn > 0 ){
if( !is.null(mexsimdata) ){
if( !is.list(mexsimdata) | length(mexsimdata) != m.sim ) stop("\ncgarchsim-->error: mexsimdata must be a list of length equal to m.sim...", call. = FALSE)
for(j in 1:m.sim){
if( !is.matrix(mexsimdata[[j]]) ) stop("\ncgarchsim-->error: mexsimdata list submatrix must be a matrix...", call. = FALSE)
if( dim(mexsimdata[[j]])[2] != mxn ) stop("\ncgarchsim-->error: wrong number of external regressors in list submatrix!...", call. = FALSE)
if( dim(mexsimdata[[j]])[1] < (n.sim + n.start) ) stop("\ncgarchsim-->error: external regressor list submatrix has less points than requested simulation length (n.sim + n.start)...", call. = FALSE)
}
} else{
mexsimdata = vector(mode = "list", length = m.sim)
for(j in 1:m.sim) mexsimdata[[i]] = matrix(0, ncol = mxn, nrow = (n.sim + n.start))
}
} else{
mexsimdata = NULL
}
vrmodel = fit@mfit$vrmodel
p = as.numeric(vrmodel$lag)
if(startMethod[1] == "sample"){
if( !is.null(prereturns) ){
if(!is.matrix(prereturns)) stop("\nprereturns must be a matrix\n")
if(dim(prereturns)[1] != p) stop("\nwrong number of rows for prereturns matrix")
if(dim(prereturns)[2] != m) stop("\nwrong number of cols for prereturns matrix")
prereturns = matrix(prereturns, nrow = p, ncol = m)
} else{
prereturns = as.matrix( tail(Data, p) )
}
} else{
prereturns = matrix(0, nrow = p, ncol = m, byrow = TRUE)
}
xlist = vector(mode = "list", length = m.sim)
if( n.sim == 1 && n.start == 0 ){
# SPECIAL ROUTINE TO SPEED UP 1-AHEAD ESTIMATION
tmp = varxsimXX(X = Data, vrmodel$Bcoef, p, m.sim = m.sim, prereturns, resids = t(sapply(simX, FUN = function(x) x)), if(!is.null(mexsimdata)) t(sapply(mexsimdata, FUN = function(x) x)) else NULL)
for(i in 1:m.sim) xlist[[i]] = tmp[i, , drop = FALSE]
simRes = simX
simX = xlist
} else{
if( parallel ){
if( parallel.control$pkg == "multicore" ){
xlist = mclapply(as.list(1:m.sim), FUN = function(i) varxsim(X = Data, vrmodel$Bcoef, p, n.sim,
n.start, prereturns, resids = simX[[i]], mexsimdata[[i]]),
mc.cores = parallel.control$cores)
} else{
# when the GARCH model is without mean equation, simX == simulated residuals
sfExport("Data", "vrmodel", "p", "n.sim", "n.start", "prereturns", "simX", "mexsimdata", local = TRUE)
xlist = sfLapply(as.list(1:m.sim), fun = function(i) rgarch:::varxsim(X = Data, vrmodel$Bcoef, p,
n.sim, n.start, prereturns, resids = simX[[i]], mexsimdata[[i]]))
sfStop()
}
} else{
xlist = lapply(as.list(1:m.sim), FUN = function(i) varxsim(X = Data, vrmodel$Bcoef, p,
n.sim, n.start, prereturns, resids = simX[[i]], mexsimdata[[i]]))
}
simRes = simX
simX = xlist
}
} else{
# reshape
for(j in 1:m.sim) simX[[j]] = tail(simX[[j]], n.sim)
}
msim = list()
msim$simH = simH
msim$simX = simX
msim$simRes = simRes
#msim$Z = z
msim$rseed = rseed
msim$model$Data = Data
#msim$model$garchpars = garchpars
msim$model$n.sim = n.sim
msim$model$m.sim = m.sim
msim$model$n.start = n.start
msim$model$startMethod = startMethod[1]
msim$model$distribution = fit@model$copula.distribution
msim$model$time.varying = FALSE
msim$model$transformation = transformation
ans = new("cGARCHsim",
msim = msim)
return( ans )
}
# time varying copula simulation
.cgarchsim2 = function(fit, n.sim = 1000, n.start = 0, m.sim = 1, startMethod = c("unconditional", "sample"),
presigma = NULL, preresiduals = NULL, prereturns = NULL, preR = NULL, rseed = NULL, mexsimdata = NULL,
vexsimdata = NULL, parallel = FALSE, parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2),
VAR.fit = NULL, ...)
{
# first generate the copula random uniform numbers (static copula)
# ures --> transform --> zres
# pass zres to GARCH simulation
# pass if needed to VAR model
if( is.null(rseed) ){
rseed = as.integer(runif(m.sim, 1, Sys.time()))
} else {
if(length(rseed) == 1) c(rseed, as.integer(runif(m.sim-1, 1, Sys.time()))) else rseed = as.integer( rseed[1:m.sim] )
}
Data = head(fit@mfit$origdata, dim(fit@mfit$origdata)[1] - fit@model$out.sample)
m = dim(Data)[2]
n = dim(Data)[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])
}
}
nsim = n.sim + n.start
spec = fit@uspec
cf = coef(fit, "garch")
cfl = vector(mode = "list", length = m)
for(i in 1:m){
if(spec@type == "equal"){
spec@spec[[i]]@optimization.model$fixed.pars = as.list(cf[,i])
cfl[[i]] = cf[,i]
} else{
spec@spec[[i]]@optimization.model$fixed.pars = as.list(cf[[i]])
cfl[[i]] = cf[[i]]
}
}
if(is.null(preR)){
Rbar = fit@copula$Rt[[length(fit@copula$Rt)]]
} else{
if(dim(preR)[1] != m) stop("\ncgarchsim-->error: wrong dimension for preR\n")
if(dim(preR)[2] != m) stop("\ncgarchsim-->error: wrong dimension for preR\n")
if(any(diag(preR) != 1)) stop("\ncgarchsim-->error: preR diagonals must be 1.\n")
Rbar = preR
}
dist = fit@model$copula.distribution
dx = paste(dist, 1, sep = "")
smplx = .sample.copula(dist = dx, object = fit, Rbar, n.sim, n.start, m.sim, rseed, parallel, parallel.control)
# ures = [0,1] copula random numbers
# now transform back into margins for use in garch
Usim = smplx$Usim
Rt = smplx$simR
transformation = fit@model$transformation
include.skew = sapply(fit@uspec@spec, FUN = function(x) x@distribution.model$include.skew)
include.shape = sapply(fit@uspec@spec, FUN = function(x) x@distribution.model$include.shape)
distx = sapply(fit@uspec@spec, FUN = function(x) x@distribution.model$distribution)
zres = array(NA, dim = c(nsim, m, m.sim))
# parallel:
if( parallel ){
if( parallel.control$pkg == "multicore" ){
mtmp = mclapply(as.list(1:m.sim), FUN = function(i) switch(transformation,
parametric = .qparametric(matrix(Usim[,,i], ncol = m), ucoef = cfl, include.skew, include.shape, dist = distx),
empirical = .qempirical(matrix(Usim[,,i], ncol = m), fit@mfit$zres),
spd = .qspd(matrix(Usim[,,i], ncol = m), sfit = fit@model$sfit)), mc.cores = parallel.control$cores)
for(i in 1:m.sim) zres[,,i] = matrix(mtmp[[i]], ncol = m)
} else{
ssfit = fit@model$sfit
sfInit(parallel = TRUE, cpus = parallel.control$cores)
sfExport("Usim", "cfl", "include.skew", "include.shape", "distx", "ssfit", "transformation", local = TRUE)
mtmp = sfLapply(as.list(1:m.sim), fun = function(i)
switch(transformation,
parametric = rgarch:::.qparametric(matrix(Usim[,,i], ncol = m), ucoef = cfl, include.skew, include.shape, dist = distx),
empirical = rgarch:::.qempirical(matrix(Usim[,,i], ncol = m), fit@mfit$zres),
spd = rgarch:::.qspd(matrix(Usim[,,i], ncol = m), sfit = ssfit)))
for(i in 1:m.sim) zres[,,i] = matrix(mtmp[[i]], ncol = m)
#sfStop()
}
} else{
mtmp = lapply(as.list(1:m.sim), FUN = function(i) switch(transformation,
parametric = .qparametric(matrix(Usim[,,i], ncol = m), ucoef = cfl, include.skew, include.shape, dist = distx),
empirical = .qempirical(matrix(Usim[,,i], ncol = m), fit@mfit$zres),
spd = .qspd(matrix(Usim[,,i], ncol = m), sfit = fit@model$sfit)))
for(i in 1:m.sim) zres[,,i] = matrix(mtmp[[i]], ncol = m)
}
# Now we have the innovations which we feed back into the garch routine to simulate
simRes = simX = simH = simSeries = vector(mode = "list", length = m.sim)
if( is.null(fit@mfit$vrmodel) ){
if(is.null(prereturns)) prereturns = fit@mfit$tailuret
}
if( parallel ){
if( parallel.control$pkg == "multicore" ){
simlist = mclapply(as.list(1:m), FUN = function(i){
maxx = spec@spec[[i]]@optimization.model$maxOrder2;
htmp = ugarchpath(spec@spec[[i]], n.sim = n.sim + n.start, n.start = 0, m.sim = m.sim,
custom.dist = list(name = "sample", distfit = matrix(zres[,i,1:m.sim], ncol = m.sim)),
presigma = if( is.null(presigma) ) tail(fit@mfit$tailusigma[,i], maxx) else tail(presigma[,i], maxx),
preresiduals = if( is.null(preresiduals) ) tail(fit@mfit$tailuresids[,i], maxx) else tail(preresiduals[,i], maxx),
prereturns = if(is.null(fit@mfit$vrmodel)) tail(prereturns[,i], maxx) else NA,
mexsimdata = if( is.null(fit@mfit$vrmodel) ) mexsimdata[[i]] else NULL, vexsimdata = vexsimdata[[i]] )
h = matrix(tail(htmp@path$sigmaSim^2, n.sim), nrow = n.sim)
x = matrix(htmp@path$seriesSim, nrow = n.sim + n.start)
return(list(h = h, x = x))}, mc.cores = parallel.control$cores)
H = array(NA, dim = c(n.sim, m, m.sim))
tmpH = array(NA, dim = c(m, m, n.sim))
for(i in 1:n.sim) H[i,,] = t(sapply(simlist, FUN = function(x) as.numeric(x$h[i,])))
for(i in 1:m.sim){
for(j in 1:n.sim){
tmpH[ , , j] = diag(sqrt( H[j, , i]) ) %*% Rt[[i]][,,j] %*% diag(sqrt( H[j, , i] ) )
}
simH[[i]] = tmpH
}
simxX = array(NA, dim = c(n.sim, m, m.sim))
for(i in 1:n.sim) simxX[i,,] = t(sapply(simlist, FUN = function(x) as.numeric(x$x[i,])))
simX = vector(mode = "list", length = m.sim)
for(i in 1:m.sim) simX[[i]] = matrix(simxX[,,i], nrow = n.sim)
} else{
if(is.null(fit@mfit$vrmodel)) include.vrmodel = FALSE else include.vrmodel = TRUE
sfExport("spec", "n.sim", "n.start", "m.sim", "startMethod", "zres", "presigma",
"preresiduals", "prereturns", "fit", "mexsimdata", "vexsimdata", "include.vrmodel", local = TRUE)
simlist = sfLapply(as.list(1:m), fun = function(i){
maxx = spec@spec[[i]]@optimization.model$maxOrder2;
htmp = rgarch::ugarchpath(spec@spec[[i]], n.sim = n.sim + n.start, n.start = 0, m.sim = m.sim,
custom.dist = list(name = "sample", distfit = matrix(zres[,i,1:m.sim], ncol = m.sim)),
presigma = if( is.null(presigma) ) tail(fit@mfit$tailusigma[,i], maxx) else tail(presigma[,i], maxx),
preresiduals = if( is.null(preresiduals) ) tail(fit@mfit$tailuresids[,i], maxx) else tail(preresiduals[,i], maxx),
prereturns = if(is.null(fit@mfit$vrmodel)) tail(prereturns[,i], maxx) else NA,
mexsimdata = if( is.null(fit@mfit$vrmodel) ) mexsimdata[[i]] else NULL, vexsimdata = vexsimdata[[i]] )
h = matrix(tail(htmp@path$sigmaSim^2, n.sim), nrow = n.sim)
x = matrix(htmp@path$seriesSim, nrow = n.sim + n.start)
return(list(h = h, x = x))})
if( !include.vrmodel ) sfStop()
H = array(NA, dim = c(n.sim, m, m.sim))
tmpH = array(NA, dim = c(m, m, n.sim))
for(i in 1:n.sim) H[i,,] = t(sapply(simlist, FUN = function(x) as.numeric(x$h[i,])))
for(i in 1:m.sim){
for(j in 1:n.sim){
tmpH[ , , j] = diag(sqrt( H[j, , i]) ) %*% Rt[[i]][,,j] %*% diag(sqrt( H[j, , i] ) )
}
simH[[i]] = tmpH
}
simxX = array(NA, dim = c(n.sim, m, m.sim))
for(i in 1:n.sim) simxX[i,,] = t(sapply(simlist, FUN = function(x) as.numeric(x$x[i,])))
simX = vector(mode = "list", length = m.sim)
for(i in 1:m.sim) simX[[i]] = matrix(simxX[,,i], nrow = n.sim)
}
} else{
simlist = lapply(as.list(1:m), FUN = function(i){
maxx = spec@spec[[i]]@optimization.model$maxOrder2;
htmp = ugarchpath(spec@spec[[i]], n.sim = n.sim + n.start, n.start = 0, m.sim = m.sim,
custom.dist = list(name = "sample", distfit = matrix(zres[,i,1:m.sim], ncol = m.sim)),
presigma = if( is.null(presigma) ) tail(fit@mfit$tailusigma[,i], maxx) else tail(presigma[,i], maxx),
preresiduals = if( is.null(preresiduals) ) tail(fit@mfit$tailuresids[,i], maxx) else tail(preresiduals[,i], maxx),
prereturns = if(is.null(fit@mfit$vrmodel)) tail(prereturns[,i], maxx) else NA,
mexsimdata = if( is.null(fit@mfit$vrmodel) ) mexsimdata[[i]] else NULL, vexsimdata = vexsimdata[[i]] )
h = matrix(tail(htmp@path$sigmaSim^2, n.sim), nrow = n.sim)
x = matrix(htmp@path$seriesSim, nrow = n.sim + n.start)
return(list(h = h, x = x))})
H = array(NA, dim = c(n.sim, m, m.sim))
tmpH = array(NA, dim = c(m, m, n.sim))
for(i in 1:n.sim) H[i,,] = t(sapply(simlist, FUN = function(x) as.numeric(x$h[i,])))
for(i in 1:m.sim){
for(j in 1:n.sim){
tmpH[ , , j] = diag(sqrt( H[j, , i]) ) %*% Rt[[i]][,,j] %*% diag(sqrt( H[j, , i] ) )
}
simH[[i]] = tmpH
}
simxX = array(NA, dim = c(n.sim, m, m.sim))
for(i in 1:n.sim) simxX[i,,] = t(sapply(simlist, FUN = function(x) as.numeric(x$x[i,])))
simX = vector(mode = "list", length = m.sim)
for(i in 1:m.sim) simX[[i]] = matrix(simxX[,,i], nrow = n.sim)
}
if( !is.null(fit@mfit$vrmodel) ){
# mexsimdata check
mxn = fit@mfit$vrmodel$mxn
if(mxn > 0 ){
if( !is.null(mexsimdata) ){
if( !is.list(mexsimdata) | length(mexsimdata) != m.sim ) stop("\ncgarchsim-->error: mexsimdata must be a list of length equal to m.sim...", call. = FALSE)
for(j in 1:m.sim){
if( !is.matrix(mexsimdata[[j]]) ) stop("\ncgarchsim-->error: mexsimdata list submatrix must be a matrix...", call. = FALSE)
if( dim(mexsimdata[[j]])[2] != mxn ) stop("\ncgarchsim-->error: wrong number of external regressors in list submatrix!...", call. = FALSE)
if( dim(mexsimdata[[j]])[1] < (n.sim + n.start) ) stop("\ncgarchsim-->error: external regressor list submatrix has less points than requested simulation length (n.sim + n.start)...", call. = FALSE)
}
} else{
mexsimdata = vector(mode = "list", length = m.sim)
for(j in 1:m.sim) mexsimdata[[i]] = matrix(0, ncol = mxn, nrow = (n.sim + n.start))
}
} else{
mexsimdata = NULL
}
vrmodel = fit@mfit$vrmodel
p = as.numeric(vrmodel$lag)
if(startMethod[1] == "sample"){
if( !is.null(prereturns) ){
if(!is.matrix(prereturns)) stop("\nprereturns must be a matrix\n")
if(dim(prereturns)[1] != p) stop("\nwrong number of rows for prereturns matrix")
if(dim(prereturns)[2] != m) stop("\nwrong number of cols for prereturns matrix")
prereturns = matrix(prereturns, nrow = p, ncol = m)
} else{
prereturns = as.matrix( tail(Data, p) )
}
} else{
prereturns = matrix(0, nrow = p, ncol = m, byrow = TRUE)
}
xlist = vector(mode = "list", length = m.sim)
if( n.sim == 1 && n.start == 0 ){
# SPECIAL ROUTINE TO SPEED UP 1-AHEAD ESTIMATION
tmp = varxsimXX(X = Data, vrmodel$Bcoef, p, m.sim = m.sim, prereturns, resids = t(sapply(simX, FUN = function(x) x)), if(!is.null(mexsimdata)) t(sapply(mexsimdata, FUN = function(x) x)) else NULL)
for(i in 1:m.sim) xlist[[i]] = tmp[i, , drop = FALSE]
simRes = simX
simX = xlist
} else{
if( parallel ){
if( parallel.control$pkg == "multicore" ){
xlist = mclapply(as.list(1:m.sim), FUN = function(i) varxsim(X = Data, vrmodel$Bcoef, p, n.sim,
n.start, prereturns, resids = simX[[i]], mexsimdata[[i]]),
mc.cores = parallel.control$cores)
} else{
# when the GARCH model is without mean equation, simX == simulated residuals
sfExport("Data", "vrmodel", "p", "n.sim", "n.start", "prereturns", "simX", "mexsimdata", local = TRUE)
xlist = sfLapply(as.list(1:m.sim), fun = function(i) rgarch:::varxsim(X = Data, vrmodel$Bcoef, p,
n.sim, n.start, prereturns, resids = simX[[i]], mexsimdata[[i]]))
sfStop()
}
} else{
xlist = lapply(as.list(1:m.sim), FUN = function(i) varxsim(X = Data, vrmodel$Bcoef, p,
n.sim, n.start, prereturns, resids = simX[[i]], mexsimdata[[i]]))
}
simRes = simX
simX = xlist
}
} else{
# reshape
for(j in 1:m.sim) simX[[j]] = tail(simX[[j]], n.sim)
}
msim = list()
msim$simH = simH
msim$simR = Rt
msim$simX = simX
msim$simRes = simRes
#msim$Z = z
msim$rseed = rseed
msim$model$Data = Data
#msim$model$garchpars = garchpars
msim$model$n.sim = n.sim
msim$model$m.sim = m.sim
msim$model$n.start = n.start
msim$model$startMethod = startMethod[1]
msim$model$distribution = fit@model$copula.distribution
msim$model$time.varying = TRUE
msim$model$transformation = transformation
ans = new("cGARCHsim",
msim = msim)
return( ans )
}
Try the rgarch package in your browser
Any scripts or data that you put into this service are public.
rgarch documentation built on May 2, 2019, 5:22 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
#################################################################################
##
## R package rgarch by Alexios Ghalanos Copyright (C) 2008, 2009, 2010
## 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.
##
#################################################################################
#################################################################################
# Fit, Filter, Forecast and Simulation
#-------------------------------------
.cgarchspec = function(uspec, VAR = FALSE, VAR.opt = list(robust = FALSE, lag = 1, lag.max = NULL,
lag.criterion = c("AIC", "HQ", "SC", "FPE"), external.regressors = NULL,
robust.control = list("gamma" = 0.25, "delta" = 0.01, "nc" = 10, "ns" = 500)),
dccOrder = c(1,1), distribution.model = list(copula = c("mvnorm", "mvt"), method = c("Kendall", "ML"),
time.varying = FALSE, transformation = c("parametric", "empirical", "spd")),
start.pars = list(), fixed.pars = list())
{
m = length(uspec@spec)
if(is.null(distribution.model$copula)) distribution = "mvnorm" else distribution = tolower(distribution.model$copula)
distribution = distribution[1]
valid.distributions = c("mvnorm", "mvt")
if(!any(distribution == valid.distributions)) stop("\nInvalid Copula Distribution Choice\n", call. = FALSE)
if(is.null(distribution.model$method)) method = "kendall" else method = tolower(distribution.model$method)
method = method[1]
valid.methods = c("kendall", "ml")
if(!any(method == valid.methods)){
if(distribution == "mvt") warning("\nInvalid Rho Method Estimation Choice\n", call. = FALSE)
method = "kendall"
}
if(is.null(distribution.model$time.varying)) timecopula = FALSE else timecopula = as.logical(distribution.model$time.varying)
timecopula = timecopula[1]
if(is.null(distribution.model$transformation)) transformation = "parametric" else transformation = tolower(distribution.model$transformation)
transformation = transformation[1]
valid.transformations = c("parametric", "empirical", "spd")
if(!any(transformation == valid.transformations)) stop("\nInvalid Copula Transformation Choice\n", call. = FALSE)
if(is.null(dccOrder) && timecopula) dccOrder = c(1,1)
if(!timecopula) dccOrder = c(0,0)
if( length(dccOrder) != 2) stop("\nInvalid dccOrder\n", call. = FALSE)
dccOrder = as.integer(dccOrder)
maxdccOrder = max(dccOrder)
maxgarchOrder = max( sapply(uspec@spec, FUN = function(x) max(max(x@mean.model$armaOrder), max(x@variance.model$garchOrder) ) ) )
maxOrder = max( maxdccOrder, maxgarchOrder )
shape.distributions = c("mvt")
skew.distributions = NULL
include.skew = include.shape = FALSE
# we will have to change this eventually to accomodate vector
# skewness and shape parameters for more interesting distributions
if(any(distribution == shape.distributions)) include.shape = TRUE
if(any(distribution == skew.distributions)) include.skew = TRUE
mspec = list()
mspec$timecopula = timecopula
mspec$method = method
mspec$transformation = transformation
mspec$dccOrder = dccOrder
mspec$start.pars = start.pars
mspec$fixed.pars = fixed.pars
mspec$distribution = distribution
mspec$include.shape = include.shape
mspec$include.skew = include.skew
modelnames = c(
if(dccOrder[1] > 0 && timecopula) paste("dcca", 1:dccOrder[1], sep = ""),
if(dccOrder[2] > 0 && timecopula) paste("dccb", 1:dccOrder[2], sep = ""),
if(include.skew) paste("skew"),
if(include.shape) paste("shape"),
if(!timecopula && method == "ml") paste("Cov", 1:((m*m - m)/2), sep = ""))
pos = 1
pos.matrix = matrix(0, ncol = 4, nrow = 5)
colnames(pos.matrix) = c("start", "stop", "include", "npars")
rownames(pos.matrix) = c("dcca","dccb","skew","shape", "Cov")
if(dccOrder[1] > 0 && timecopula){
pos.matrix[1,1:3] = c(pos, pos+dccOrder[1]-1, 1)
pos = max(pos.matrix[1:1, 2]) + 1
}
if(dccOrder[2] > 0 && timecopula){
pos.matrix[2,1:3] = c(pos, pos+dccOrder[2]-1, 1)
pos = max(pos.matrix[1:2, 2]) + 1
}
if(include.skew)
{
pos.matrix[3,1:3] = c(pos, pos, 1)
pos = max(pos.matrix[1:3, 2]) + 1
}
if(include.shape)
{
pos.matrix[4,1:3] = c(pos, pos, 1)
pos = max(pos.matrix[1:4, 2]) + 1
}
if(!timecopula && method == "ml"){
pos.matrix[5,1:3] = c(pos, pos + ((m*m - m)/2) - 1, 1)
}
nn = length(modelnames)
modelmatrix = matrix(0, ncol = 3, nrow = nn)
rownames(modelmatrix) = modelnames
colnames(modelmatrix) = c("opt", "fixed", "start")
fixed.names = names(fixed.pars)
fp = charmatch(fixed.names, modelnames)
if(!is.null(fixed.names) && any(!is.na(fp))){
fixed = fp[!is.na(fp)]
modelmatrix[fixed,2] = 1
fz = charmatch(modelnames, fixed.names)
fz = fz[!is.na(fz)]
fixed.pars = fixed.pars[fz]
names(fixed.pars) = fixed.names[fz]
} else{
fixed.pars = NULL
}
modelmatrix[,1] = 1 - modelmatrix[,2]
start.names = names(start.pars)
sp = charmatch(start.names, modelnames)
if(!is.null(start.names) && any(!is.na(sp))){
start = sp[!is.na(sp)]
modelmatrix[start,3] = 1
sz = charmatch(modelnames, start.names)
sz = sz[!is.na(sz)]
start.pars = start.pars[sz]
} else{
start.pars = NULL
}
opt.pars = vector(mode = "numeric", length = length(modelnames))
names(opt.pars) = modelnames
# optimization model
omodel = list(modelnames = modelnames, opt.pars = opt.pars, start.pars = start.pars,
fixed.pars = fixed.pars, modelmatrix = modelmatrix, maxOrder = maxOrder,
maxgarchOrder = maxgarchOrder, maxdccOrder = maxdccOrder,
pos.matrix = pos.matrix)
mspec$optimization.model = omodel
if( !is(uspec, "uGARCHmultispec") ) stop("\ncgarchspec-->error: uspec must be a uGARCHmultispec object")
mspec$varmodel = list()
if( VAR ){
mspec$varmodel$VAR = TRUE
mspec$varmodel$robust = VAR.opt$robust
mspec$varmodel$lag = VAR.opt$lag
mspec$varmodel$lag.max = VAR.opt$lag.max
mspec$varmodel$lag.criterion = VAR.opt$lag.criterion
mspec$varmodel$external.regressors = VAR.opt$external.regressors
mspec$varmodel$robust.control = VAR.opt$robust.control
if( !is.null(mspec$varmodel$external.regressors) ){
mspec$varmodel$mxn = dim( as.matrix(VAR.opt$external.regressors) )[2]
} else{
mspec$varmodel$mxn = 0
}
# loop and reconstruct the spec to exclude a mean specification
nl = length(uspec@spec)
tmpspec = vector(mode ="list", length = nl)
for(i in 1:nl){
tmpspec[[i]] = ugarchspec(mean.model = list(include.mean = FALSE, armaOrder = c(0,0),
external.regressors = NULL, arfima = FALSE, garchInMean = FALSE),
variance.model = list(model = uspec@spec[[i]]@variance.model$model,
submodel = uspec@spec[[i]]@variance.model$submodel,
garchOrder = uspec@spec[[i]]@variance.model$garchOrder,
external.regressors = uspec@spec[[i]]@variance.model$vexdata,
variance.targeting = !uspec@spec[[i]]@variance.model$include.omega),
distribution.model = uspec@spec[[i]]@distribution.model$distribution,
start.pars = uspec@spec[[i]]@optimization.model$start.pars,
fixed.pars = uspec@spec[[i]]@optimization.model$fixed.pars)
}
uspec = multispec( tmpspec )
} else{
mspec$varmodel$VAR = FALSE
mspec$varmodel$lag = 1
mspec$varmodel$lag.max = 1
mspec$varmodel$lag.criterion = "HQ"
mspec$varmodel$external.regressors = NULL
mspec$varmodel$mxn = 0
}
ans = new("cGARCHspec",
mspec = mspec,
uspec = uspec)
return(ans)
}
setMethod("cgarchspec", signature(uspec = "uGARCHmultispec"), .cgarchspec)
.cgarchfit = function(spec, data, spd.control = list(lower = 0.1, upper = 0.9, type = "pwm", kernel = "epanech"),
fit.control = list(eval.se = TRUE, trace = TRUE, stationarity = TRUE), solver = "solnp", solver.control = list(), out.sample = 0, parallel = FALSE,
parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2), fit = NULL, VAR.fit = NULL, ...)
{
dist = spec@mspec$distribution
tv = ifelse(spec@mspec$timecopula, 1, 0)
dx = paste(dist, tv, sep = "")
ans = switch(dx,
mvt0 = .cgarchfit.student(spec = spec, data = data, spd.control = spd.control, fit.control = fit.control, solver = solver,
solver.control = solver.control, out.sample = out.sample, parallel = parallel,
parallel.control = parallel.control, fit = fit, VAR.fit = VAR.fit, ...),
mvt1 = .cgarchfit.tvstudent(spec = spec, data = data, spd.control = spd.control, fit.control = fit.control, solver = solver,
solver.control = solver.control, out.sample = out.sample, parallel = parallel,
parallel.control = parallel.control, fit = fit, VAR.fit = VAR.fit, ...),
mvnorm0 = .cgarchfit.norm(spec = spec, data = data, spd.control = spd.control, fit.control = fit.control, solver = solver,
solver.control = solver.control, out.sample = out.sample, parallel = parallel,
parallel.control = parallel.control, fit = fit, VAR.fit = VAR.fit, ...),
mvnorm1 = .cgarchfit.tvnorm(spec = spec, data = data, spd.control = spd.control, fit.control = fit.control, solver = solver,
solver.control = solver.control, out.sample = out.sample, parallel = parallel,
parallel.control = parallel.control, fit = fit, VAR.fit = VAR.fit, ...))
return( ans )
}
setMethod("cgarchfit", signature(spec = "cGARCHspec"), .cgarchfit)
.cgarchfit.tvnorm = function(spec, data, spd.control = list(lower = 0.1, upper = 0.9, type = "pwm", kernel = "epanech"),
fit.control = list(eval.se = TRUE, trace = TRUE, stationarity = TRUE), solver = "solnp", solver.control = list(), out.sample = 0, parallel = FALSE,
parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2), fit = NULL, VAR.fit = NULL, ...)
{
tic = Sys.time()
ufit.control = list()
if(is.null(fit.control$stationarity)){
ufit.control$stationarity = TRUE
fit.control$stationarity = TRUE
} else {
ufit.control$stationarity = fit.control$stationarity
fit.control$stationarity = fit.control$stationarity
}
if(is.null(fit.control$scale)){
ufit.control$scale = TRUE
} else{
ufit.control$scale = fit.control$scale
fit.control$scale = NULL
}
if(is.null(fit.control$eval.se)){
eval.se = TRUE
} else{
eval.se = as.logical(fit.control$eval.se)
}
# Parallel Execution Prelim Checks (we need them because of the hessian estimation)
if( parallel ){
os = .Platform$OS.type
if(is.null(parallel.control$pkg)){
if( os == "windows" ) parallel.control$pkg = "snowfall" else parallel.control$pkg = "multicore"
if( is.null(parallel.control$cores) ) parallel.control$cores = 2
} else{
mtype = match(tolower(parallel.control$pkg[1]), c("multicore", "snowfall"))
if(is.na(mtype)) stop("\nParallel Package type not recognized in parallel.control\n")
parallel.control$pkg = tolower(parallel.control$pkg[1])
if( os == "windows" && parallel.control$pkg == "multicore" ) stop("\nmulticore not supported on windows O/S\n")
if( is.null(parallel.control$cores) ) parallel.control$cores = 2 else parallel.control$cores = as.integer(parallel.control$cores[1])
}
}
if( is.null( colnames(data) ) ) cnames = paste("Asset_", 1:k, sep = "") else cnames = colnames(data)
xdata = .extractmdata(data)
if(!is.numeric(out.sample)) stop("\ncgarchfit-->error: out.sample must be numeric\n")
if(as.numeric(out.sample) < 0) stop("\ncgarchfit-->error: out.sample must be positive\n")
n.start = round(out.sample, 0)
n = dim(xdata$data)[1]
k = dim(xdata$data)[2]
if( (n-n.start) < 100) stop("\ncgarchfit-->error: function requires at least 100 data\n points to run\n")
data = xdata$data[1:(n-n.start), , drop = FALSE]
dates = xdata$pos[1:(n-n.start), drop = FALSE]
origdata = xdata$data
origdates = xdata$pos
dformat = xdata$dformat
vrmodel = NULL
if( spec@mspec$varmodel$VAR ){
if(spec@mspec$varmodel$mxn>0){
ex = spec@mspec$varmodel$external.regressors
if( dim(ex)[1] < dim(data)[1] ) stop("\ncgarchfit-->error: external regressor matrix has less points than data!...", call. = FALSE)
orex = ex
ex = ex[1:(n - n.start), , drop = FALSE]
mxn = dim(orex)[2]
} else{
orex = NULL
ex = NULL
mxn = 0
}
if( !is.null(VAR.fit) ){
zdata = VAR.fit$xresiduals
# should be N - p
p = mp = VAR.fit$lag
N = dim(zdata)[1]
K = dim(zdata)[2]
mu = VAR.fit$xfitted
vrmodel = VAR.fit
#spec@mspec$optimization.model$maxOrder = max( mp, spec@mspec$optimization.model$maxdccOrder )
#spec@mspec$optimization.model$maxgarchOrder = mp
} else{
cat("\nestimating VAR model...")
ic = spec@mspec$varmodel$lag.criterion[1]
if( !is.null(spec@mspec$varmodel$lag.max) ){
mp = .varxselect(y = data, lag.max = spec@mspec$varmodel$lag.max, exogen = ex)$selection
mp = as.numeric( mp[which(substr(names(mp), 1, 2) == substr(ic, 1, 2))] )
} else {
mp = spec@mspec$varmodel$lag
}
vrmodel = varxfilter(X = data, p = mp, exogen = ex, Bcoef = NULL, robust = spec@mspec$varmodel$robust,
gamma = spec@mspec$varmodel$robust.control$gamma, delta = spec@mspec$varmodel$robust.control$delta,
nc = spec@mspec$varmodel$robust.control$nc, ns = spec@mspec$varmodel$robust.control$ns)
cat("done!\n")
zdata = vrmodel$xresiduals
# should be N - p
N = dim(zdata)[1]
K = dim(zdata)[2]
mu = vrmodel$xfitted
vrmodel$lag = mp
vrmodel$mxn = mxn
#spec@mspec$optimization.model$maxOrder = max( mp, spec@mspec$optimization.model$maxdccOrder )
#spec@mspec$optimization.model$maxgarchOrder = mp
}
# if we have out of sample data we need to filter
if(out.sample>0){
Bcoef = vrmodel$Bcoef
vrmodel2 = varxfilter(X = origdata, p = mp, exogen = orex, Bcoef = Bcoef)
zorigdata = vrmodel2$xresiduals
} else{
zorigdata = zdata
}
} else{
zdata = data
zorigdata = origdata
orex = NULL
ex = NULL
}
if( is(spec@uspec, "uGARCHmultispec") ){
speclist = spec@uspec
} else{
stop("\ncgarchfit-->error: cGARCHspec must contain multispec object with fixed parameters")
}
maxgarchOrder = spec@mspec$optimization.model$maxgarchOrder
tmpidx = .dccparidx( speclist )
paridx = tmpidx$idx
paridxi = tmpidx$idxi
paridxa = tmpidx$idxa
if( !is.null(fit) && is(fit, "uGARCHmultifit") ){
fitlist = fit
} else{
fitlist = multifit(multispec = speclist, data = zorigdata, out.sample = n.start, solver = solver,
solver.control = solver.control, fit.control = ufit.control, parallel = parallel,
parallel.control = parallel.control)
assign(".fitlist", fitlist, envir = .GlobalEnv)
}
garchpars = as.numeric( unlist(coef(fitlist) ) )
converge = sapply(fitlist@fit, FUN = function(x) x@fit$convergence)
if( any( converge == 1 ) ){
pr = which(converge != 1)
cat("\nNon-Converged:\n")
print(pr)
cat("\ncgarchfit-->error: convergence problem in univariate fit...")
cat("\n...returning uGARCHmultifit object instead...check and resubmit...")
return( fitlist )
}
# create a temporary environment to store values (deleted at end of function)
tempenvir = new.env(hash = TRUE, parent = .GlobalEnv)
if(is.null(fit.control$trace)) trace = FALSE else trace = as.logical(fit.control$trace)
assign("trace", trace, envir = tempenvir)
assign("eval.se", eval.se, envir = tempenvir)
assign("cnames", cnames, envir = tempenvir)
assign("spec", spec, envir = tempenvir)
assign("parallel", parallel, envir = tempenvir)
assign("parallel.control", parallel.control, envir = tempenvir)
assign("paridx", paridx, envir = tempenvir)
assign("paridxi", paridxi, envir = tempenvir)
assign("paridxa", paridxa, envir = tempenvir)
assign("speclist", speclist, envir = tempenvir)
assign("fitlist", fitlist, envir = tempenvir)
omodel = spec@mspec$optimization.model
distribution = spec@mspec$distribution
transformation = spec@mspec$transformation
sig = sigma(fitlist)
res = residuals(fitlist)
#stdresid = res/sig
zres = sapply(fitlist@fit, FUN = function(x) x@fit$z)
ures = matrix(NA, ncol = 5, nrow = dim(zres)[1])
if( maxgarchOrder > 0 ) zres = zres[-(1:maxgarchOrder), , drop = FALSE]
# Transformation of z to [0,1] domain
ans = switch(transformation,
parametric = .pparametric(fitlist, zres),
empirical = .pempirical(zres),
spd = .pspd(zres, spd.control))
if(transformation == "spd"){
ures = ans$ures
sfit = ans$sfit
} else{
ures = ans
sfit = NULL
}
assign("spd.control", spd.control, envir = tempenvir)
# make tail adjustments in order to avoid problem with the quantile functions in
# optimization
if(any(ures > 0.99999)){
xn = which(ures > 0.99999)
ures[xn] = 0.99999
}
if(any(ures < eps)){
xn = which(ures < (1.5*eps))
ures[xn] = eps
}
assign("omodel", omodel, envir = tempenvir)
assign("mspec", spec@mspec, envir = tempenvir)
assign("ures", ures, envir = tempenvir)
assign("solver", solver, envir = tempenvir)
assign("fit.control", fit.control, envir = tempenvir)
Ifn = .dcccon
ILB = 0.0001
IUB = 0.9998
if( solver == "solnp" ) fit.control$stationarity = FALSE else fit.control$stationarity = TRUE
assign("fit.control", fit.control, envir = tempenvir)
ipars = .tvcopulastart(data = ures, garchenv = tempenvir)
pars = ipars$pars
LB = ipars$modelLB
UB = ipars$modelUB
assign("npar", length(pars), envir = tempenvir)
if( any( LB == UB ) )
{
fixid = which(LB == UB)
if( length(fixid) == length(pars) )
{
if( !eval.se )
{
# if all parameters are fixed an no standard erros are to
# be calculated then we return a ugarchfilter object
cat("\ncgarchfit-->warning: all parameters fixed...returning cgarchfilter object instead\n")
return(cgarchfilter(spec = spec, data = ures, out.sample = out.sample,
parallel = parallel, parallel.control = parallel.control, VAR.fit = VAR.fit,
fit = fitlist))
} else{
# if all parameters are fixed but we require standard errors, we
# skip the solver
use.solver = FALSE
fixed = pars[fixid]
pars = pars[-fixid]
}
} else{
# with some parameters fixed we extract them (to be rejoined at end)
# so that they do not enter the solver
use.solver = TRUE
fixed = pars[fixid]
pars = pars[-fixid]
LB = LB[-fixid]
UB = UB[-fixid]
}
} else{
# no fixed parameters, all normal
fixed = NULL
fixid = NULL
use.solver = TRUE
}
assign("fixed", fixed, envir = tempenvir)
assign("fixid", fixid, envir = tempenvir)
assign(".llh", 1, envir = tempenvir)
# get
if( use.solver )
{
solution = .copulasolver(solver, pars, fun = copula.tvnormalLLH, Ifn, ILB,
IUB, gr = NULL, hessian = NULL, control = solver.control,
LB = LB, UB = UB, data = ures, returnType = "llh", garchenv = tempenvir)
sol = solution$sol
hess = solution$hess
opt.pars = sol$par
names(opt.pars) = names(pars)
convergence = sol$convergence
} else{
hess = NULL
timer = Sys.time()-tic
opt.pars = NULL
convergence = 0
sol = list()
sol$message = "all parameters fixed"
}
mfit = list()
mfit$dates = dates
mfit$date.format = dformat
mfit$origdata = origdata
mfit$origdates = origdates
mfit$vrmodel = vrmodel
# convert back from uniform to standard variates
model = list()
model$copula.distribution = distribution
model$out.sample = out.sample
model$transformation = transformation
model$spd.control = spd.control
model$sfit = sfit
model$time.varying = TRUE
# add some tail data of the univariate fits to be used with simulation (we assume that noone is
# going to add more than 25 lags in a univariate fit!
mfit$tailusigma = tail(sigma(fitlist), 25)
mfit$tailuresids = tail(residuals(fitlist), 25)
mfit$tailuret = tail(data, 25)
cfit = copula.tvnormalLLH(opt.pars, data = ures, returnType = "ALL", garchenv = tempenvir)
if(eval.se){
nderiv = .cgarchcvar(f = copula.tvnormalLLH, pars = opt.pars, garchfit = fitlist, model = model,
env = tempenvir, fname = "copula.tvnormalLLH")
mfit$pars = nderiv$pars
mfit$garchpars = coef(fitlist)
# stpars = second stage parameters
mfit$stpars = opt.pars
mfit$scores = nderiv$scores
mfit$A = nderiv$A
mfit$cvar = nderiv$cvar
mfit$matcoef = nderiv$matcoef
} else{
garchNames = vector(mode = "character")
assetNames = cnames
dccNames = spec@mspec$optimization.model$modelnames
m = length(spec@uspec@spec)
for(i in 1:m){
cf = coef(fitlist@fit[[i]])
cfnames = names(cf)
garchNames = c(garchNames, paste(assetNames[i], cfnames, sep = ".") )
}
allnames = c(garchNames, dccNames)
mfit$pars = c(as.numeric(unlist(coef(fitlist))), opt.pars)
mfit$garchpars = coef(fitlist)
mfit$stpars = opt.pars
names(mfit$pars) = allnames
mfit$scores = NULL
mfit$A = NULL
mfit$cvar = NULL
matcoef = matrix(NA, ncol = 4, nrow = length(mfit$pars))
matcoef[,1] = mfit$pars
dimnames(matcoef) = list(allnames, c(" Estimate",
" Std. Error", " t value", "Pr(>|t|)"))
mfit$matcoef = matcoef
}
mfit$paridx = tmpidx
garchllh = sapply(fitlist@fit, FUN = function(x) -x@fit$log.likelihoods)
copllh = cfit$lik
if( maxgarchOrder > 0 ) copllh = c(rep(0, maxgarchOrder), copllh)
mfit$llh = sum( apply(cbind(copllh, garchllh), 1, "sum") )
mfit$timer = Sys.time() - tic
mfit$zres = zres
ans = new("cGARCHfit",
mfit = mfit,
mspec = spec@mspec,
uspec = spec@uspec,
copula = cfit,
model = model)
return(ans)
}
.cgarchfit.norm = function(spec, data, spd.control = list(lower = 0.1, upper = 0.9, type = "pwm", kernel = "epanech"),
fit.control = list(eval.se = TRUE, trace = TRUE, stationarity = TRUE), solver = "solnp", solver.control = list(), out.sample = 0, parallel = FALSE,
parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2), fit = NULL, VAR.fit = NULL, ...)
{
tic = Sys.time()
ufit.control = list()
if(is.null(fit.control$stationarity)){
ufit.control$stationarity = TRUE
} else {
ufit.control$stationarity = fit.control$stationarity
fit.control$stationarity = NULL
}
if(is.null(fit.control$scale)){
ufit.control$scale = TRUE
} else{
ufit.control$scale = fit.control$scale
fit.control$scale = NULL
}
if(is.null(fit.control$eval.se)){
eval.se = TRUE
} else{
eval.se = as.logical(fit.control$eval.se)
}
# Parallel Execution Prelim Checks (we need them because of the hessian estimation)
if( parallel ){
os = .Platform$OS.type
if(is.null(parallel.control$pkg)){
if( os == "windows" ) parallel.control$pkg = "snowfall" else parallel.control$pkg = "multicore"
if( is.null(parallel.control$cores) ) parallel.control$cores = 2
} else{
mtype = match(tolower(parallel.control$pkg[1]), c("multicore", "snowfall"))
if(is.na(mtype)) stop("\nParallel Package type not recognized in parallel.control\n")
parallel.control$pkg = tolower(parallel.control$pkg[1])
if( os == "windows" && parallel.control$pkg == "multicore" ) stop("\nmulticore not supported on windows O/S\n")
if( is.null(parallel.control$cores) ) parallel.control$cores = 2 else parallel.control$cores = as.integer(parallel.control$cores[1])
}
}
if( is.null( colnames(data) ) ) cnames = paste("Asset_", 1:k, sep = "") else cnames = colnames(data)
xdata = .extractmdata(data)
if(!is.numeric(out.sample)) stop("\ncgarchfit-->error: out.sample must be numeric\n")
if(as.numeric(out.sample) < 0) stop("\ncgarchfit-->error: out.sample must be positive\n")
n.start = round(out.sample, 0)
n = dim(xdata$data)[1]
k = dim(xdata$data)[2]
if( (n-n.start) < 100) stop("\ncgarchfit-->error: function requires at least 100 data\n points to run\n")
data = xdata$data[1:(n-n.start), , drop = FALSE]
dates = xdata$pos[1:(n-n.start), drop = FALSE]
origdata = xdata$data
origdates = xdata$pos
dformat = xdata$dformat
vrmodel = NULL
if( spec@mspec$varmodel$VAR ){
if(spec@mspec$varmodel$mxn>0){
ex = spec@mspec$varmodel$external.regressors
if( dim(ex)[1] < dim(data)[1] ) stop("\ncgarchfit-->error: external regressor matrix has less points than data!...", call. = FALSE)
orex = ex
ex = ex[1:(n - n.start), , drop = FALSE]
mxn = dim(orex)[2]
} else{
orex = NULL
ex = NULL
mxn = 0
}
if( !is.null(VAR.fit) ){
zdata = VAR.fit$xresiduals
# should be N - p
p = mp = VAR.fit$lag
N = dim(zdata)[1]
K = dim(zdata)[2]
mu = VAR.fit$xfitted
vrmodel = VAR.fit
#spec@mspec$optimization.model$maxOrder = max( mp, spec@mspec$optimization.model$maxdccOrder )
#spec@mspec$optimization.model$maxgarchOrder = mp
} else{
cat("\nestimating VAR model...")
ic = spec@mspec$varmodel$lag.criterion[1]
if( !is.null(spec@mspec$varmodel$lag.max) ){
mp = .varxselect(y = data, lag.max = spec@mspec$varmodel$lag.max, exogen = ex)$selection
mp = as.numeric( mp[which(substr(names(mp), 1, 2) == substr(ic, 1, 2))] )
} else {
mp = spec@mspec$varmodel$lag
}
vrmodel = varxfilter(X = data, p = mp, exogen = ex, Bcoef = NULL, robust = spec@mspec$varmodel$robust,
gamma = spec@mspec$varmodel$robust.control$gamma, delta = spec@mspec$varmodel$robust.control$delta,
nc = spec@mspec$varmodel$robust.control$nc, ns = spec@mspec$varmodel$robust.control$ns)
cat("done!\n")
zdata = vrmodel$xresiduals
# should be N - p
N = dim(zdata)[1]
K = dim(zdata)[2]
mu = vrmodel$xfitted
vrmodel$lag = mp
vrmodel$mxn = mxn
#spec@mspec$optimization.model$maxOrder = max( mp, spec@mspec$optimization.model$maxdccOrder )
#spec@mspec$optimization.model$maxgarchOrder = mp
}
# if we have out of sample data we need to filter
if(out.sample>0){
Bcoef = vrmodel$Bcoef
vrmodel2 = varxfilter(X = origdata, p = mp, exogen = orex, Bcoef = Bcoef)
zorigdata = vrmodel2$xresiduals
} else{
zorigdata = zdata
}
} else{
zdata = data
zorigdata = origdata
orex = NULL
ex = NULL
}
if( is(spec@uspec, "uGARCHmultispec") ){
speclist = spec@uspec
} else{
stop("\ncgarchfit-->error: cGARCHspec must contain multispec object with fixed parameters")
}
maxgarchOrder = spec@mspec$optimization.model$maxgarchOrder
tmpidx = .dccparidx( speclist )
paridx = tmpidx$idx
paridxi = tmpidx$idxi
paridxa = tmpidx$idxa
if( !is.null(fit) && is(fit, "uGARCHmultifit") ){
fitlist = fit
} else{
fitlist = multifit(multispec = speclist, data = zorigdata, out.sample = n.start, solver = solver,
solver.control = solver.control, fit.control = ufit.control, parallel = parallel,
parallel.control = parallel.control)
assign(".fitlist", fitlist, envir = .GlobalEnv)
}
garchpars = as.numeric( unlist(coef(fitlist) ) )
converge = sapply(fitlist@fit, FUN = function(x) x@fit$convergence)
if( any( converge == 1 ) ){
pr = which(converge != 1)
cat("\nNon-Converged:\n")
print(pr)
cat("\ncgarchfit-->error: convergence problem in univariate fit...")
cat("\n...returning uGARCHmultifit object instead...check and resubmit...")
return( fitlist )
}
# create a temporary environment to store values (deleted at end of function)
tempenvir = new.env(hash = TRUE, parent = .GlobalEnv)
if(is.null(fit.control$trace)) trace = FALSE else trace = as.logical(fit.control$trace)
assign("trace", trace, envir = tempenvir)
assign("method", spec@mspec$method, envir = tempenvir)
assign("eval.se", eval.se, envir = tempenvir)
assign("cnames", cnames, envir = tempenvir)
assign("spec", spec, envir = tempenvir)
assign("parallel", parallel, envir = tempenvir)
assign("parallel.control", parallel.control, envir = tempenvir)
assign("paridx", paridx, envir = tempenvir)
assign("paridxi", paridxi, envir = tempenvir)
assign("paridxa", paridxa, envir = tempenvir)
assign("speclist", speclist, envir = tempenvir)
assign("fitlist", fitlist, envir = tempenvir)
omodel = spec@mspec$optimization.model
distribution = spec@mspec$distribution
transformation = spec@mspec$transformation
sig = sigma(fitlist)
res = residuals(fitlist)
#stdresid = res/sig
assign("transformation", transformation, envir = tempenvir)
zres = sapply(fitlist@fit, FUN = function(x) x@fit$z)
ures = matrix(NA, ncol = 5, nrow = dim(zres)[1])
if( maxgarchOrder > 0 ) zres = zres[-(1:maxgarchOrder), , drop = FALSE]
# Transformation of z to [0,1] domain
ans = switch(transformation,
parametric = .pparametric(fitlist, zres),
empirical = .pempirical(zres),
spd = .pspd(zres, spd.control))
if(transformation == "spd"){
ures = ans$ures
sfit = ans$sfit
} else{
ures = ans
sfit = NULL
}
assign("spd.control", spd.control, envir = tempenvir)
# make tail adjustments in order to avoid problem with the quantile functions in
# optimization
if(any(ures > 0.99999)){
xn = which(ures > 0.99999)
ures[xn] = 0.99999
}
if(any(ures < eps)){
xn = which(ures < (1.5*eps))
ures[xn] = eps
}
assign("omodel", omodel, envir = tempenvir)
assign("mspec", spec@mspec, envir = tempenvir)
assign("ures", ures, envir = tempenvir)
assign("solver", solver, envir = tempenvir)
assign("fit.control", fit.control, envir = tempenvir)
Ifn = NULL
ILB = NULL
IUB = NULL
if( solver == "solnp" ) fit.control$stationarity = FALSE else fit.control$stationarity = TRUE
assign("fit.control", fit.control, envir = tempenvir)
if( spec@mspec$method == "ml" ){
ipars = .copulastart(data = ures, garchenv = tempenvir)
pars = ipars$pars
LB = ipars$modelLB
UB = ipars$modelUB
assign("npar", length(pars), envir = tempenvir)
if( any( LB == UB ) )
{
fixid = which(LB == UB)
if( length(fixid) == length(pars) )
{
if( !eval.se )
{
cat("\ncgarchfit-->warning: all parameters fixed...returning cgarchfilter object instead\n")
spec2 = spec
cf = coef(fitlist)
for(i in 1:k){
if(spec@uspec@type == "equal"){
spec2@uspec@spec[[i]]@optimization.model$fixed.pars = as.list(cf[,i])
} else{
spec2@uspec@spec[[i]]@optimization.model$fixed.pars = as.list(cf[[i]])
}
}
return(cgarchfilter(spec = spec2, filter.control = list(n.old = NULL), spd.control = spd.control,
out.sample = out.sample, parallel = parallel, parallel.control = parallel.control,
VAR.fit = vrmodel))
} else{
# if all parameters are fixed but we require standard errors, we
# skip the solver
use.solver = FALSE
fixed = pars[fixid]
pars = pars[-fixid]
}
} else{
# with some parameters fixed we extract them (to be rejoined at end)
# so that they do not enter the solver
use.solver = TRUE
fixed = pars[fixid]
pars = pars[-fixid]
LB = LB[-fixid]
UB = UB[-fixid]
}
} else{
# no fixed parameters, all normal
fixed = NULL
fixid = NULL
use.solver = TRUE
}
assign("fixed", fixed, envir = tempenvir)
assign("fixid", fixid, envir = tempenvir)
assign(".llh", 1, envir = tempenvir)
# get
if( use.solver )
{
# overide solver control for static R
solution = .copulasolver(solver, pars, fun = copula.normalLLH, Ifn, ILB,
IUB, gr = NULL, hessian = NULL, control = solver.control,
LB = LB, UB = UB, data = ures, returnType = "llh", garchenv = tempenvir)
sol = solution$sol
hess = solution$hess
opt.pars = sol$par
names(opt.pars) = names(pars)
convergence = sol$convergence
} else{
hess = NULL
timer = Sys.time()-tic
opt.pars = NULL
convergence = 0
sol = list()
sol$message = "all parameters fixed"
}
} else{
assign("npar", 0, envir = tempenvir)
fixed = NULL
fixid = NULL
use.solver = FALSE
assign("fixed", fixed, envir = tempenvir)
assign("fixid", fixid, envir = tempenvir)
assign(".llh", 1, envir = tempenvir)
hess = NULL
timer = Sys.time()-tic
opt.pars = NULL
convergence = 0
sol = list()
sol$message = NULL
}
mfit = list()
mfit$dates = dates
mfit$date.format = dformat
mfit$origdata = origdata
mfit$origdates = origdates
mfit$vrmodel = vrmodel
# convert back from uniform to standard variates
model = list()
model$copula.distribution = distribution
model$out.sample = out.sample
model$transformation = transformation
model$spd.control = spd.control
model$sfit = sfit
model$time.varying = FALSE
# add some tail data of the univariate fits to be used with simulation
mfit$tailusigma = tail(sigma(fitlist), 25)
mfit$tailuresids = tail(residuals(fitlist), 25)
mfit$tailuret = tail(data, 25)
cfit = copula.normalLLH(opt.pars, data = ures, returnType = "ALL", garchenv = tempenvir)
if(eval.se && spec@mspec$method == "ml"){
nderiv = .cgarchcvar(f = copula.normalLLH, pars = opt.pars, garchfit = fitlist, model = model,
env = tempenvir, fname = "copula.normalLLH")
mfit$pars = nderiv$pars
mfit$garchpars = coef(fitlist)
# stpars = second stage parameters
mfit$stpars = opt.pars
mfit$scores = nderiv$scores
mfit$A = nderiv$A
mfit$cvar = nderiv$cvar
mfit$matcoef = nderiv$matcoef
} else{
if( spec@mspec$method == "ml" ){
garchNames = vector(mode = "character")
assetNames = cnames
dccNames = spec@mspec$optimization.model$modelnames
m = length(spec@uspec@spec)
for(i in 1:m){
cf = coef(fitlist@fit[[i]])
cfnames = names(cf)
garchNames = c(garchNames, paste(assetNames[i], cfnames, sep = ".") )
}
allnames = c(garchNames, dccNames)
mfit$pars = c(as.numeric(unlist(coef(fitlist))), opt.pars)
mfit$garchpars = coef(fitlist)
# stpars = second stage parameters
mfit$stpars = opt.pars
names(mfit$pars) = allnames
mfit$scores = NULL
mfit$A = NULL
mfit$cvar = NULL
matcoef = matrix(NA, ncol = 4, nrow = length(mfit$pars))
matcoef[,1] = mfit$pars
dimnames(matcoef) = list(allnames, c(" Estimate",
" Std. Error", " t value", "Pr(>|t|)"))
mfit$matcoef = matcoef
} else{
garchNames = vector(mode = "character")
assetNames = cnames
dccNames = NULL
m = length(spec@uspec@spec)
for(i in 1:m){
cf = coef(fitlist@fit[[i]])
cfnames = names(cf)
garchNames = c(garchNames, paste(assetNames[i], cfnames, sep = ".") )
}
allnames = c(garchNames, dccNames)
mfit$pars = c(as.numeric(unlist(coef(fitlist))), opt.pars)
mfit$garchpars = coef(fitlist)
# stpars = second stage parameters
mfit$stpars = NULL
names(mfit$pars) = allnames
matcoef = NULL
mfit$scores = NULL
mfit$A = NULL
mfit$cvar = NULL
for(i in 1:m){
matcoef = rbind(matcoef, fitlist@fit[[i]]@fit$robust.matcoef)
}
dimnames(matcoef) = list(allnames, c(" Estimate",
" Std. Error", " t value", "Pr(>|t|)"))
mfit$matcoef = matcoef
}
}
mfit$paridx = tmpidx
garchllh = sapply(fitlist@fit, FUN = function(x) -x@fit$log.likelihoods)
copllh = cfit$lik
if( maxgarchOrder > 0 ) copllh = c(rep(0, maxgarchOrder), copllh)
mfit$llh = sum( apply(cbind(copllh, garchllh), 1, "sum") )
mfit$zres = zres
mfit$timer = Sys.time() - tic
ans = new("cGARCHfit",
mfit = mfit,
mspec = spec@mspec,
uspec = spec@uspec,
copula = cfit,
model = model)
return(ans)
}
.cgarchfit.tvstudent = function(spec, data, spd.control = list(lower = 0.1, upper = 0.9, type = "pwm", kernel = "epanech"),
fit.control = list(eval.se = TRUE, trace = TRUE, stationarity = TRUE), solver = "solnp", solver.control = list(), out.sample = 0, parallel = FALSE,
parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2), fit = NULL, VAR.fit = NULL, ...)
{
tic = Sys.time()
ufit.control = list()
if(is.null(fit.control$stationarity)){
ufit.control$stationarity = TRUE
fit.control$stationarity = TRUE
} else {
ufit.control$stationarity = fit.control$stationarity
fit.control$stationarity = fit.control$stationarity
}
if(is.null(fit.control$scale)){
ufit.control$scale = TRUE
} else{
ufit.control$scale = fit.control$scale
fit.control$scale = NULL
}
if(is.null(fit.control$eval.se)){
eval.se = TRUE
} else{
eval.se = as.logical(fit.control$eval.se)
}
# Parallel Execution Prelim Checks (we need them because of the hessian estimation)
if( parallel ){
os = .Platform$OS.type
if(is.null(parallel.control$pkg)){
if( os == "windows" ) parallel.control$pkg = "snowfall" else parallel.control$pkg = "multicore"
if( is.null(parallel.control$cores) ) parallel.control$cores = 2
} else{
mtype = match(tolower(parallel.control$pkg[1]), c("multicore", "snowfall"))
if(is.na(mtype)) stop("\nParallel Package type not recognized in parallel.control\n")
parallel.control$pkg = tolower(parallel.control$pkg[1])
if( os == "windows" && parallel.control$pkg == "multicore" ) stop("\nmulticore not supported on windows O/S\n")
if( is.null(parallel.control$cores) ) parallel.control$cores = 2 else parallel.control$cores = as.integer(parallel.control$cores[1])
}
}
if( is.null( colnames(data) ) ) cnames = paste("Asset_", 1:k, sep = "") else cnames = colnames(data)
xdata = .extractmdata(data)
if(!is.numeric(out.sample)) stop("\ncgarchfit-->error: out.sample must be numeric\n")
if(as.numeric(out.sample) < 0) stop("\ncgarchfit-->error: out.sample must be positive\n")
n.start = round(out.sample, 0)
n = dim(xdata$data)[1]
k = dim(xdata$data)[2]
if( (n-n.start) < 100) stop("\ncgarchfit-->error: function requires at least 100 data\n points to run\n")
data = xdata$data[1:(n-n.start), , drop = FALSE]
dates = xdata$pos[1:(n-n.start), drop = FALSE]
origdata = xdata$data
origdates = xdata$pos
dformat = xdata$dformat
vrmodel = NULL
if( spec@mspec$varmodel$VAR ){
if(spec@mspec$varmodel$mxn>0){
ex = spec@mspec$varmodel$external.regressors
if( dim(ex)[1] < dim(data)[1] ) stop("\ncgarchfit-->error: external regressor matrix has less points than data!...", call. = FALSE)
orex = ex
ex = ex[1:(n - n.start), , drop = FALSE]
mxn = dim(orex)[2]
} else{
orex = NULL
ex = NULL
mxn = 0
}
if( !is.null(VAR.fit) ){
zdata = VAR.fit$xresiduals
# should be N - p
p = mp = VAR.fit$lag
N = dim(zdata)[1]
K = dim(zdata)[2]
mu = VAR.fit$xfitted
vrmodel = VAR.fit
#spec@mspec$optimization.model$maxOrder = max( mp, spec@mspec$optimization.model$maxdccOrder )
#spec@mspec$optimization.model$maxgarchOrder = mp
} else{
cat("\nestimating VAR model...")
ic = spec@mspec$varmodel$lag.criterion[1]
if( !is.null(spec@mspec$varmodel$lag.max) ){
mp = .varxselect(y = data, lag.max = spec@mspec$varmodel$lag.max, exogen = ex)$selection
mp = as.numeric( mp[which(substr(names(mp), 1, 2) == substr(ic, 1, 2))] )
} else {
mp = spec@mspec$varmodel$lag
}
vrmodel = varxfilter(X = data, p = mp, exogen = ex, Bcoef = NULL, robust = spec@mspec$varmodel$robust,
gamma = spec@mspec$varmodel$robust.control$gamma, delta = spec@mspec$varmodel$robust.control$delta,
nc = spec@mspec$varmodel$robust.control$nc, ns = spec@mspec$varmodel$robust.control$ns)
cat("done!\n")
zdata = vrmodel$xresiduals
# should be N - p
N = dim(zdata)[1]
K = dim(zdata)[2]
mu = vrmodel$xfitted
vrmodel$lag = mp
vrmodel$mxn = mxn
#spec@mspec$optimization.model$maxOrder = max( mp, spec@mspec$optimization.model$maxdccOrder )
#spec@mspec$optimization.model$maxgarchOrder = mp
}
# if we have out of sample data we need to filter
if(out.sample>0){
Bcoef = vrmodel$Bcoef
vrmodel2 = varxfilter(X = origdata, p = mp, exogen = orex, Bcoef = Bcoef)
zorigdata = vrmodel2$xresiduals
} else{
zorigdata = zdata
}
} else{
zdata = data
zorigdata = origdata
orex = NULL
ex = NULL
}
if( is(spec@uspec, "uGARCHmultispec") ){
speclist = spec@uspec
} else{
stop("\ncgarchfit-->error: cGARCHspec must contain multispec object with fixed parameters")
}
maxgarchOrder = spec@mspec$optimization.model$maxgarchOrder
tmpidx = .dccparidx( speclist )
paridx = tmpidx$idx
paridxi = tmpidx$idxi
paridxa = tmpidx$idxa
if( !is.null(fit) && is(fit, "uGARCHmultifit") ){
fitlist = fit
} else{
fitlist = multifit(multispec = speclist, data = zorigdata, out.sample = n.start, solver = solver,
solver.control = solver.control, fit.control = ufit.control, parallel = parallel,
parallel.control = parallel.control)
assign(".fitlist", fitlist, envir = .GlobalEnv)
}
garchpars = as.numeric( unlist(coef(fitlist) ) )
converge = sapply(fitlist@fit, FUN = function(x) x@fit$convergence)
if( any( converge == 1 ) ){
pr = which(converge != 1)
cat("\nNon-Converged:\n")
print(pr)
cat("\ncgarchfit-->error: convergence problem in univariate fit...")
cat("\n...returning uGARCHmultifit object instead...check and resubmit...")
return( fitlist )
}
# create a temporary environment to store values (deleted at end of function)
tempenvir = new.env(hash = TRUE, parent = .GlobalEnv)
if(is.null(fit.control$trace)) trace = FALSE else trace = as.logical(fit.control$trace)
assign("trace", trace, envir = tempenvir)
assign("eval.se", eval.se, envir = tempenvir)
assign("cnames", cnames, envir = tempenvir)
assign("spec", spec, envir = tempenvir)
assign("parallel", parallel, envir = tempenvir)
assign("parallel.control", parallel.control, envir = tempenvir)
assign("paridx", paridx, envir = tempenvir)
assign("paridxi", paridxi, envir = tempenvir)
assign("paridxa", paridxa, envir = tempenvir)
assign("speclist", speclist, envir = tempenvir)
assign("fitlist", fitlist, envir = tempenvir)
omodel = spec@mspec$optimization.model
distribution = spec@mspec$distribution
transformation = spec@mspec$transformation
sig = sigma(fitlist)
res = residuals(fitlist)
#stdresid = res/sig
assign("transformation", transformation, envir = tempenvir)
zres = sapply(fitlist@fit, FUN = function(x) x@fit$z)
ures = matrix(NA, ncol = 5, nrow = dim(zres)[1])
if( maxgarchOrder > 0 ) zres = zres[-(1:maxgarchOrder), , drop = FALSE]
# Transformation of z to [0,1] domain
ans = switch(transformation,
parametric = .pparametric(fitlist, zres),
empirical = .pempirical(zres),
spd = .pspd(zres, spd.control))
if(transformation == "spd"){
ures = ans$ures
sfit = ans$sfit
} else{
ures = ans
sfit = NULL
}
assign("spd.control", spd.control, envir = tempenvir)
# make tail adjustments in order to avoid problem with the quantile functions in
# optimization
if(any(ures > 0.99999)){
xn = which(ures > 0.99999)
ures[xn] = 0.99999
}
if(any(ures < eps)){
xn = which(ures < (1.5*eps))
ures[xn] = eps
}
assign("omodel", omodel, envir = tempenvir)
assign("mspec", spec@mspec, envir = tempenvir)
assign("ures", ures, envir = tempenvir)
assign("solver", solver, envir = tempenvir)
assign("fit.control", fit.control, envir = tempenvir)
Ifn = .dcccon
ILB = 0.0001
IUB = 0.9998
if( solver == "solnp" ) fit.control$stationarity = FALSE else fit.control$stationarity = TRUE
assign("fit.control", fit.control, envir = tempenvir)
ipars = .tvcopulastart(data = ures, garchenv = tempenvir)
pars = ipars$pars
LB = ipars$modelLB
UB = ipars$modelUB
assign("npar", length(pars), envir = tempenvir)
if( any( LB == UB ) )
{
fixid = which(LB == UB)
if( length(fixid) == length(pars) )
{
if( !eval.se )
{
cat("\ncgarchfit-->warning: all parameters fixed...returning cgarchfilter object instead\n")
spec2 = spec
cf = coef(fitlist)
for(i in 1:k){
if(spec@uspec@type == "equal"){
spec2@uspec@spec[[i]]@optimization.model$fixed.pars = as.list(cf[,i])
} else{
spec2@uspec@spec[[i]]@optimization.model$fixed.pars = as.list(cf[[i]])
}
}
return(cgarchfilter(spec = spec2, filter.control = list(n.old = NULL), spd.control = spd.control,
out.sample = out.sample, parallel = parallel, parallel.control = parallel.control,
VAR.fit = vrmodel))
} else{
# if all parameters are fixed but we require standard errors, we
# skip the solver
use.solver = FALSE
fixed = pars[fixid]
pars = pars[-fixid]
}
} else{
# with some parameters fixed we extract them (to be rejoined at end)
# so that they do not enter the solver
use.solver = TRUE
fixed = pars[fixid]
pars = pars[-fixid]
LB = LB[-fixid]
UB = UB[-fixid]
}
} else{
# no fixed parameters, all normal
fixed = NULL
fixid = NULL
use.solver = TRUE
}
assign("fixed", fixed, envir = tempenvir)
assign("fixid", fixid, envir = tempenvir)
assign(".llh", 1, envir = tempenvir)
# get
if( use.solver )
{
solution = .copulasolver(solver, pars, fun = copula.tvstudentLLH, Ifn, ILB,
IUB, gr = NULL, hessian = NULL, control = solver.control,
LB = LB, UB = UB, data = ures, returnType = "llh", garchenv = tempenvir)
sol = solution$sol
hess = solution$hess
opt.pars = sol$par
names(opt.pars) = names(pars)
convergence = sol$convergence
} else{
hess = NULL
timer = Sys.time()-tic
opt.pars = NULL
convergence = 0
sol = list()
sol$message = "all parameters fixed"
}
mfit = list()
mfit$dates = dates
mfit$date.format = dformat
mfit$origdata = origdata
mfit$origdates = origdates
mfit$vrmodel = vrmodel
# convert back from uniform to standard variates
model = list()
model$copula.distribution = distribution
model$out.sample = out.sample
model$transformation = transformation
model$spd.control = spd.control
model$sfit = sfit
model$time.varying = TRUE
# add some tail data of the univariate fits to be used with simulation
mfit$tailusigma = tail(sigma(fitlist), 25)
mfit$tailuresids = tail(residuals(fitlist), 25)
mfit$tailuret = tail(data, 25)
cfit = copula.tvstudentLLH(opt.pars, data = ures, returnType = "ALL", garchenv = tempenvir)
if(eval.se){
nderiv = .cgarchcvar(f = copula.tvstudentLLH, pars = opt.pars, garchfit = fitlist, model = model,
env = tempenvir, fname = "copula.tvstudentLLH")
mfit$pars = nderiv$pars
mfit$garchpars = coef(fitlist)
# stpars = second stage parameters
mfit$stpars = opt.pars
mfit$scores = nderiv$scores
mfit$A = nderiv$A
mfit$cvar = nderiv$cvar
mfit$matcoef = nderiv$matcoef
} else{
garchNames = vector(mode = "character")
assetNames = cnames
dccNames = spec@mspec$optimization.model$modelnames
m = length(spec@uspec@spec)
for(i in 1:m){
cf = coef(fitlist@fit[[i]])
cfnames = names(cf)
garchNames = c(garchNames, paste(assetNames[i], cfnames, sep = ".") )
}
allnames = c(garchNames, dccNames)
mfit$pars = c(as.numeric(unlist(coef(fitlist))), opt.pars)
mfit$garchpars = coef(fitlist)
# stpars = second stage parameters
mfit$stpars = opt.pars
names(mfit$pars) = allnames
mfit$scores = NULL
mfit$A = NULL
mfit$cvar = NULL
matcoef = matrix(NA, ncol = 4, nrow = length(mfit$pars))
matcoef[,1] = mfit$pars
dimnames(matcoef) = list(allnames, c(" Estimate",
" Std. Error", " t value", "Pr(>|t|)"))
mfit$matcoef = matcoef
}
mfit$paridx = tmpidx
garchllh = sapply(fitlist@fit, FUN = function(x) -x@fit$log.likelihoods)
copllh = cfit$lik
if( maxgarchOrder > 0 ) copllh = c(rep(0, maxgarchOrder), copllh)
mfit$llh = sum( apply(cbind(copllh, garchllh), 1, "sum") )
mfit$zres = zres
mfit$timer = Sys.time() - tic
ans = new("cGARCHfit",
mfit = mfit,
mspec = spec@mspec,
uspec = spec@uspec,
copula = cfit,
model = model)
return(ans)
}
.cgarchfit.student = function(spec, data, spd.control = list(lower = 0.1, upper = 0.9, type = "pwm", kernel = "epanech"),
fit.control = list(eval.se = TRUE, trace = TRUE, stationarity = TRUE), solver = "solnp", solver.control = list(), out.sample = 0, parallel = FALSE,
parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2), fit = NULL, VAR.fit = NULL, ...)
{
tic = Sys.time()
ufit.control = list()
if(is.null(fit.control$stationarity)){
ufit.control$stationarity = TRUE
} else {
ufit.control$stationarity = fit.control$stationarity
fit.control$stationarity = NULL
}
if(is.null(fit.control$scale)){
ufit.control$scale = TRUE
} else{
ufit.control$scale = fit.control$scale
fit.control$scale = NULL
}
if(is.null(fit.control$eval.se)){
eval.se = TRUE
} else{
eval.se = as.logical(fit.control$eval.se)
}
# Parallel Execution Prelim Checks (we need them because of the hessian estimation)
if( parallel ){
os = .Platform$OS.type
if(is.null(parallel.control$pkg)){
if( os == "windows" ) parallel.control$pkg = "snowfall" else parallel.control$pkg = "multicore"
if( is.null(parallel.control$cores) ) parallel.control$cores = 2
} else{
mtype = match(tolower(parallel.control$pkg[1]), c("multicore", "snowfall"))
if(is.na(mtype)) stop("\nParallel Package type not recognized in parallel.control\n")
parallel.control$pkg = tolower(parallel.control$pkg[1])
if( os == "windows" && parallel.control$pkg == "multicore" ) stop("\nmulticore not supported on windows O/S\n")
if( is.null(parallel.control$cores) ) parallel.control$cores = 2 else parallel.control$cores = as.integer(parallel.control$cores[1])
}
}
if( is.null( colnames(data) ) ) cnames = paste("Asset_", 1:k, sep = "") else cnames = colnames(data)
xdata = .extractmdata(data)
if(!is.numeric(out.sample)) stop("\ncgarchfit-->error: out.sample must be numeric\n")
if(as.numeric(out.sample) < 0) stop("\ncgarchfit-->error: out.sample must be positive\n")
n.start = round(out.sample, 0)
n = dim(xdata$data)[1]
k = dim(xdata$data)[2]
if( (n-n.start) < 100) stop("\ncgarchfit-->error: function requires at least 100 data\n points to run\n")
data = xdata$data[1:(n-n.start), , drop = FALSE]
dates = xdata$pos[1:(n-n.start), drop = FALSE]
origdata = xdata$data
origdates = xdata$pos
dformat = xdata$dformat
vrmodel = NULL
if( spec@mspec$varmodel$VAR ){
if(spec@mspec$varmodel$mxn>0){
ex = spec@mspec$varmodel$external.regressors
if( dim(ex)[1] < dim(data)[1] ) stop("\ncgarchfit-->error: external regressor matrix has less points than data!...", call. = FALSE)
orex = ex
ex = ex[1:(n - n.start), , drop = FALSE]
mxn = dim(orex)[2]
} else{
orex = NULL
ex = NULL
mxn = 0
}
if( !is.null(VAR.fit) ){
zdata = VAR.fit$xresiduals
# should be N - p
p = mp = VAR.fit$lag
N = dim(zdata)[1]
K = dim(zdata)[2]
mu = VAR.fit$xfitted
vrmodel = VAR.fit
#spec@mspec$optimization.model$maxOrder = max( mp, spec@mspec$optimization.model$maxdccOrder )
#spec@mspec$optimization.model$maxgarchOrder = mp
} else{
cat("\nestimating VAR model...")
ic = spec@mspec$varmodel$lag.criterion[1]
if( !is.null(spec@mspec$varmodel$lag.max) ){
mp = .varxselect(y = data, lag.max = spec@mspec$varmodel$lag.max, exogen = ex)$selection
mp = as.numeric( mp[which(substr(names(mp), 1, 2) == substr(ic, 1, 2))] )
} else {
mp = spec@mspec$varmodel$lag
}
vrmodel = varxfilter(X = data, p = mp, exogen = ex, Bcoef = NULL, robust = spec@mspec$varmodel$robust,
gamma = spec@mspec$varmodel$robust.control$gamma, delta = spec@mspec$varmodel$robust.control$delta,
nc = spec@mspec$varmodel$robust.control$nc, ns = spec@mspec$varmodel$robust.control$ns)
cat("done!\n")
zdata = vrmodel$xresiduals
# should be N - p
N = dim(zdata)[1]
K = dim(zdata)[2]
mu = vrmodel$xfitted
vrmodel$lag = mp
vrmodel$mxn = mxn
#spec@mspec$optimization.model$maxOrder = max( mp, spec@mspec$optimization.model$maxdccOrder )
#spec@mspec$optimization.model$maxgarchOrder = mp
}
# if we have out of sample data we need to filter
if(out.sample>0){
Bcoef = vrmodel$Bcoef
vrmodel2 = varxfilter(X = origdata, p = mp, exogen = orex, Bcoef = Bcoef)
zorigdata = vrmodel2$xresiduals
} else{
zorigdata = zdata
}
} else{
zdata = data
zorigdata = origdata
orex = NULL
ex = NULL
}
if( is(spec@uspec, "uGARCHmultispec") ){
speclist = spec@uspec
} else{
stop("\ncgarchfit-->error: cGARCHspec must contain multispec object with fixed parameters")
}
maxgarchOrder = spec@mspec$optimization.model$maxgarchOrder
tmpidx = .dccparidx( speclist )
paridx = tmpidx$idx
paridxi = tmpidx$idxi
paridxa = tmpidx$idxa
if( !is.null(fit) && is(fit, "uGARCHmultifit") ){
fitlist = fit
} else{
fitlist = multifit(multispec = speclist, data = zorigdata, out.sample = n.start, solver = solver,
solver.control = solver.control, fit.control = ufit.control, parallel = parallel,
parallel.control = parallel.control)
assign(".fitlist", fitlist, envir = .GlobalEnv)
}
garchpars = as.numeric( unlist(coef(fitlist) ) )
converge = sapply(fitlist@fit, FUN = function(x) x@fit$convergence)
if( any( converge == 1 ) ){
pr = which(converge != 1)
cat("\nNon-Converged:\n")
print(pr)
cat("\ncgarchfit-->error: convergence problem in univariate fit...")
cat("\n...returning uGARCHmultifit object instead...check and resubmit...")
return( fitlist )
}
# create a temporary environment to store values (deleted at end of function)
tempenvir = new.env(hash = TRUE, parent = .GlobalEnv)
if(is.null(fit.control$trace)) trace = FALSE else trace = as.logical(fit.control$trace)
assign("trace", trace, envir = tempenvir)
assign("method", spec@mspec$method, envir = tempenvir)
assign("eval.se", eval.se, envir = tempenvir)
assign("cnames", cnames, envir = tempenvir)
assign("spec", spec, envir = tempenvir)
assign("parallel", parallel, envir = tempenvir)
assign("parallel.control", parallel.control, envir = tempenvir)
assign("paridx", paridx, envir = tempenvir)
assign("paridxi", paridxi, envir = tempenvir)
assign("paridxa", paridxa, envir = tempenvir)
assign("speclist", speclist, envir = tempenvir)
assign("fitlist", fitlist, envir = tempenvir)
omodel = spec@mspec$optimization.model
distribution = spec@mspec$distribution
transformation = spec@mspec$transformation
sig = sigma(fitlist)
res = residuals(fitlist)
#stdresid = res/sig
assign("transformation", transformation, envir = tempenvir)
zres = sapply(fitlist@fit, FUN = function(x) x@fit$z)
ures = matrix(NA, ncol = 5, nrow = dim(zres)[1])
if( maxgarchOrder > 0 ) zres = zres[-(1:maxgarchOrder), , drop = FALSE]
# Transformation of z to [0,1] domain
ans = switch(transformation,
parametric = .pparametric(fitlist, zres),
empirical = .pempirical(zres),
spd = .pspd(zres, spd.control))
if(transformation == "spd"){
ures = ans$ures
sfit = ans$sfit
} else{
ures = ans
sfit = NULL
}
assign("spd.control", spd.control, envir = tempenvir)
# make tail adjustments in order to avoid problem with the quantile functions in
# optimization
if(any(ures > 0.99999)){
xn = which(ures > 0.99999)
ures[xn] = 0.99999
}
if(any(ures < eps)){
xn = which(ures < (1.5*eps))
ures[xn] = eps
}
assign("omodel", omodel, envir = tempenvir)
assign("mspec", spec@mspec, envir = tempenvir)
assign("ures", ures, envir = tempenvir)
assign("solver", solver, envir = tempenvir)
assign("fit.control", fit.control, envir = tempenvir)
Ifn = NULL
ILB = NULL
IUB = NULL
if( solver == "solnp" ) fit.control$stationarity = FALSE else fit.control$stationarity = TRUE
assign("fit.control", fit.control, envir = tempenvir)
ipars = .copulastart(data = ures, garchenv = tempenvir)
pars = ipars$pars
LB = ipars$modelLB
UB = ipars$modelUB
assign("npar", length(pars), envir = tempenvir)
if( any( LB == UB ) )
{
fixid = which(LB == UB)
if( length(fixid) == length(pars) )
{
if( !eval.se )
{
# if all parameters are fixed an no standard erros are to
# be calculated then we return a ugarchfilter object
cat("\ncgarchfit-->warning: all parameters fixed...returning cgarchfilter object instead\n")
spec2 = spec
cf = coef(fitlist)
for(i in 1:k){
if(spec@uspec@type == "equal"){
spec2@uspec@spec[[i]]@optimization.model$fixed.pars = as.list(cf[,i])
} else{
spec2@uspec@spec[[i]]@optimization.model$fixed.pars = as.list(cf[[i]])
}
}
return(cgarchfilter(spec = spec2, filter.control = list(n.old = NULL), spd.control = spd.control,
out.sample = out.sample, parallel = parallel, parallel.control = parallel.control,
VAR.fit = vrmodel))
} else{
# if all parameters are fixed but we require standard errors, we
# skip the solver
use.solver = FALSE
fixed = pars[fixid]
pars = pars[-fixid]
}
} else{
# with some parameters fixed we extract them (to be rejoined at end)
# so that they do not enter the solver
use.solver = TRUE
fixed = pars[fixid]
pars = pars[-fixid]
LB = LB[-fixid]
UB = UB[-fixid]
}
} else{
# no fixed parameters, all normal
fixed = NULL
fixid = NULL
use.solver = TRUE
}
assign("fixed", fixed, envir = tempenvir)
assign("fixid", fixid, envir = tempenvir)
assign(".llh", 1, envir = tempenvir)
# get
if( use.solver )
{
# overide solver control for static R and use plain optim (L-BFGS-B)
if(spec@mspec$method == "kendall"){
solution = .copulasolver(solver = "lbfgs", pars, fun = copula.studentLLH, Ifn, ILB,
IUB, gr = NULL, hessian = NULL, control = list(),
LB = LB, UB = UB, data = ures, returnType = "llh", garchenv = tempenvir)
sol = solution$sol
hess = solution$hess
opt.pars = sol$par
names(opt.pars) = names(pars)
convergence = sol$convergence
} else{
solution = .copulasolver(solver, pars, fun = copula.studentLLH, Ifn, ILB,
IUB, gr = NULL, hessian = NULL, control = solver.control,
LB = LB, UB = UB, data = ures, returnType = "llh", garchenv = tempenvir)
sol = solution$sol
hess = solution$hess
opt.pars = sol$par
names(opt.pars) = names(pars)
convergence = sol$convergence
}
} else{
hess = NULL
timer = Sys.time()-tic
opt.pars = NULL
convergence = 0
sol = list()
sol$message = "all parameters fixed"
}
mfit = list()
mfit$dates = dates
mfit$date.format = dformat
mfit$origdata = origdata
mfit$origdates = origdates
mfit$vrmodel = vrmodel
# convert back from uniform to standard variates
model = list()
model$copula.distribution = distribution
model$out.sample = out.sample
model$transformation = transformation
model$spd.control = spd.control
model$sfit = sfit
model$time.varying = FALSE
# add some tail data of the univariate fits to be used with simulation
mfit$tailusigma = tail(sigma(fitlist), 25)
mfit$tailuresids = tail(residuals(fitlist), 25)
mfit$tailuret = tail(data, 25)
cfit = copula.studentLLH(opt.pars, data = ures, returnType = "ALL", garchenv = tempenvir)
if(eval.se){
nderiv = .cgarchcvar(f = copula.studentLLH, pars = opt.pars, garchfit = fitlist, model = model,
env = tempenvir, fname = "copula.studentLLH")
mfit$pars = nderiv$pars
mfit$garchpars = coef(fitlist)
mfit$stpars = opt.pars
mfit$scores = nderiv$scores
mfit$A = nderiv$A
mfit$cvar = nderiv$cvar
mfit$matcoef = nderiv$matcoef
} else{
garchNames = vector(mode = "character")
assetNames = cnames
dccNames = spec@mspec$optimization.model$modelnames
m = length(spec@uspec@spec)
for(i in 1:m){
cf = coef(fitlist@fit[[i]])
cfnames = names(cf)
garchNames = c(garchNames, paste(assetNames[i], cfnames, sep = ".") )
}
allnames = c(garchNames, dccNames)
mfit$pars = c(as.numeric(unlist(coef(fitlist))), opt.pars)
names(mfit$pars) = allnames
mfit$garchpars = coef(fitlist)
mfit$stpars = opt.pars
mfit$scores = NULL
mfit$A = NULL
mfit$cvar = NULL
matcoef = matrix(NA, ncol = 4, nrow = length(mfit$pars))
matcoef[,1] = mfit$pars
dimnames(matcoef) = list(allnames, c(" Estimate",
" Std. Error", " t value", "Pr(>|t|)"))
mfit$matcoef = matcoef
}
mfit$paridx = tmpidx
garchllh = sapply(fitlist@fit, FUN = function(x) -x@fit$log.likelihoods)
copllh = cfit$lik
if( maxgarchOrder > 0 ) copllh = c(rep(0, maxgarchOrder), copllh)
mfit$llh = sum( apply(cbind(copllh, garchllh), 1, "sum") )
mfit$zres = zres
mfit$timer = Sys.time() - tic
ans = new("cGARCHfit",
mfit = mfit,
mspec = spec@mspec,
uspec = spec@uspec,
copula = cfit,
model = model)
return(ans)
}
.cgarchfilter = function(spec, data, out.sample = 0, filter.control = list(n.old = NULL), spd.control = list(lower = 0.1, upper = 0.9, type = "pwm", kernel = "epanech"),
parallel = FALSE, parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2), VAR.fit = NULL, ...)
{
dist = spec@mspec$distribution
tv = ifelse(spec@mspec$timecopula, 1, 0)
dx = paste(dist, tv, sep = "")
ans = switch(dx,
mvt0 = .cgarchfilter.student(spec = spec, data = data, spd.control = spd.control, out.sample = out.sample, parallel = parallel,
parallel.control = parallel.control, VAR.fit = VAR.fit, ...),
mvt1 = .cgarchfilter.tvstudent(spec = spec, data = data, spd.control = spd.control, out.sample = out.sample, parallel = parallel,
parallel.control = parallel.control, VAR.fit = VAR.fit, ...),
mvnorm0 = .cgarchfilter.norm(spec = spec, data = data, spd.control = spd.control, out.sample = out.sample, parallel = parallel,
parallel.control = parallel.control, VAR.fit = VAR.fit, ...),
mvnorm1 = .cgarchfilter.tvnorm(spec = spec, data = data, spd.control = spd.control, out.sample = out.sample, parallel = parallel,
parallel.control = parallel.control, VAR.fit = VAR.fit, ...))
return( ans )
}
setMethod("cgarchfilter", signature(spec = "cGARCHspec"), .cgarchfilter)
.cgarchfilter.tvstudent = function(spec, data, filter.control = list(n.old = NULL), spd.control = list(lower = 0.1, upper = 0.9, type = "pwm", kernel = "epanech"),
out.sample = 0, parallel = FALSE, parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2),
VAR.fit = NULL, ...)
{
tic = Sys.time()
if( parallel ){
os = .Platform$OS.type
if(is.null(parallel.control$pkg)){
if( os == "windows" ) parallel.control$pkg = "snowfall" else parallel.control$pkg = "multicore"
if( is.null(parallel.control$cores) ) parallel.control$cores = 2
} else{
mtype = match(tolower(parallel.control$pkg[1]), c("multicore", "snowfall"))
if(is.na(mtype)) stop("\nParallel Package type not recognized in parallel.control\n")
parallel.control$pkg = tolower(parallel.control$pkg[1])
if( os == "windows" && parallel.control$pkg == "multicore" ) stop("\nmulticore not supported on windows O/S\n")
if( is.null(parallel.control$cores) ) parallel.control$cores = 2 else parallel.control$cores = as.integer(parallel.control$cores[1])
}
}
n.old = filter.control$n.old
if( is.null( colnames(data) ) ) cnames = paste("Asset_", 1:k, sep = "") else cnames = colnames(data)
xdata = .extractmdata(data)
if(!is.numeric(out.sample)) stop("\ncgarchfilter-->error: out.sample must be numeric\n")
if(as.numeric(out.sample) < 0) stop("\ncgarchfilter-->error: out.sample must be positive\n")
n.start = round(out.sample, 0)
n = dim(xdata$data)[1]
k = dim(xdata$data)[2]
if( (n - n.start) < 100) stop("\ncgarchfilter-->error: function requires at least 100 data\n points to run\n")
data = xdata$data[1:(n - n.start), , drop = FALSE]
dates = xdata$pos[1:(n - n.start)]
origdata = xdata$data
origdates = xdata$pos
dformat = xdata$dformat
vrmodel = NULL
if( spec@mspec$varmodel$VAR ){
if(spec@mspec$varmodel$mxn>0){
ex = spec@mspec$varmodel$external.regressors
if( dim(ex)[1] < dim(data)[1] ) stop("\ndccfilter-->error: external regressor matrix has less points than data!...", call. = FALSE)
orex = ex
ex = ex[1:(n - n.start), , drop = FALSE]
mxn = dim(orex)[2]
} else{
orex = NULL
ex = NULL
mxn = 0
}
if( !is.null(VAR.fit) ){
zdata = VAR.fit$xresiduals
# should be N - p
p = mp = VAR.fit$lag
N = dim(zdata)[1]
K = dim(zdata)[2]
mu = VAR.fit$xfitted
vrmodel = VAR.fit
if(out.sample>0){
Bcoef = VAR.fit$Bcoef
vrmodel2 = varxfilter(X = origdata, p = mp, exogen = orex, Bcoef = Bcoef)
zorigdata = vrmodel2$xresiduals
} else{
zorigdata = zdata
}
} else{
cat("\nestimating VAR model...")
ic = spec@mspec$varmodel$lag.criterion[1]
if( !is.null(spec@mspec$varmodel$lag.max) ){
mp = .varxselect(y = data, lag.max = spec@mspec$varmodel$lag.max, exogen = ex)$selection
mp = as.numeric( mp[which(substr(names(mp), 1, 2) == substr(ic, 1, 2))] )
} else {
mp = spec@mspec$varmodel$lag
}
vrmodel = varxfilter(X = data, p = mp, exogen = ex, Bcoef = NULL, robust = spec@mspec$varmodel$robust,
gamma = spec@mspec$varmodel$robust.control$gamma, delta = spec@mspec$varmodel$robust.control$delta,
nc = spec@mspec$varmodel$robust.control$nc, ns = spec@mspec$varmodel$robust.control$ns)
cat("done!\n")
zdata = vrmodel$xresiduals
# should be N - p
N = dim(zdata)[1]
K = dim(zdata)[2]
mu = vrmodel$xfitted
vrmodel$lag = mp
vrmodel$mxn = mxn
# if we have out of sample data we need to filter
if(out.sample>0){
Bcoef = vrmodel$Bcoef
vrmodel2 = varxfilter(X = origdata, p = mp, exogen = orex, Bcoef = Bcoef)
zorigdata = vrmodel2$xresiduals
} else{
zorigdata = zdata
}
}
} else{
zdata = data
zorigdata = origdata
orex = NULL
ex = NULL
}
if( is(spec@uspec, "uGARCHmultispec") ){
speclist = spec@uspec
} else{
stop("\ncgarchfilter-->error: cGARCHspec must contain multispec object with fixed parameters")
}
maxgarchOrder = spec@mspec$optimization.model$maxgarchOrder
filterlist = multifilter(multifitORspec = spec@uspec, data = as.data.frame(zorigdata), out.sample = out.sample,
parallel = parallel, parallel.control = parallel.control, n.old = n.old)
tmpidx = .dccparidx( speclist )
paridx = tmpidx$idx
paridxi = tmpidx$idxi
paridxa = tmpidx$idxa
garchpars = as.numeric( unlist(coef(filterlist) ) )
# create a temporary environment to store values (deleted at end of function)
tempenvir = new.env(hash = TRUE, parent = .GlobalEnv)
assign("trace", FALSE, envir = tempenvir)
assign("cnames", cnames, envir = tempenvir)
assign("spec", spec, envir = tempenvir)
assign("parallel", parallel, envir = tempenvir)
assign("parallel.control", parallel.control, envir = tempenvir)
assign("paridx", paridx, envir = tempenvir)
assign("paridxi", paridxi, envir = tempenvir)
assign("paridxa", paridxa, envir = tempenvir)
assign("speclist", speclist, envir = tempenvir)
assign("filterlist", filterlist, envir = tempenvir)
omodel = spec@mspec$optimization.model
distribution = spec@mspec$distribution
transformation = spec@mspec$transformation
sig = sigma(filterlist)
res = residuals(filterlist)
#stdresid = res/sig
zres = sapply(filterlist@filter, FUN = function(x) x@filter$z)
ures = matrix(NA, ncol = 5, nrow = dim(zres)[1])
if( maxgarchOrder > 0 ) zres = zres[-(1:maxgarchOrder), , drop = FALSE]
assign("spd.control", spd.control, envir = tempenvir)
# Transformation of z to [0,1] domain
ans = switch(transformation,
parametric = .pparametric.filter(filterlist@filter, zres),
empirical = .pempirical(zres),
spd = .pspd(zres, spd.control))
if(transformation == "spd"){
ures = ans$ures
sfilter = ans$sfit
} else{
ures = ans
sfilter = NULL
}
# make tail adjustments in order to avoid problem with the quantile functions in
# optimization
if(any(ures > 0.99999)){
xn = which(ures > 0.99999)
ures[xn] = 0.99999
}
if(any(ures < eps)){
xn = which(ures < (1.5*eps))
ures[xn] = eps
}
assign("omodel", omodel, envir = tempenvir)
assign("mspec", spec@mspec, envir = tempenvir)
assign("ures", ures, envir = tempenvir)
assign("solver", "solnp", envir = tempenvir)
ipars = .tvcopulastart(data = ures, garchenv = tempenvir)
pars = ipars$pars
assign("npar", length(pars), envir = tempenvir)
fixed = NULL
fixid = NULL
assign("fixed", fixed, envir = tempenvir)
assign("fixid", fixid, envir = tempenvir)
assign(".llh", 1, envir = tempenvir)
mfilter = list()
mfilter$dates = dates
mfilter$date.format = dformat
mfilter$origdata = origdata
mfilter$origdates = origdates
mfilter$vrmodel = vrmodel
# convert back from uniform to standard variates
model = list()
model$copula.distribution = distribution
model$out.sample = out.sample
model$transformation = transformation
model$spd.control = spd.control
model$sfilter = sfilter
model$time.varying = TRUE
# add some tail data of the univariate fits to be used with simulation
mfilter$tailusigma = tail(sigma(filterlist), 25)
mfilter$tailuresids = tail(residuals(filterlist), 25)
mfilter$tailuret = tail(data, 25)
cfilter = copula.tvstudentLLH(pars, data = ures, returnType = "ALL", garchenv = tempenvir)
garchNames = vector(mode = "character")
assetNames = cnames
dccNames = spec@mspec$optimization.model$modelnames
m = length(spec@uspec@spec)
for(i in 1:m){
cf = coef(filterlist@filter[[i]])
cfnames = names(cf)
garchNames = c(garchNames, paste(assetNames[i], cfnames, sep = ".") )
}
allnames = c(garchNames, dccNames)
mfilter$pars = c(as.numeric(coef(filterlist)), pars)
mfilter$garchpars = coef(filterlist)
mfilter$stpars = pars
names(mfilter$pars) = allnames
mfilter$scores = NULL
mfilter$A = NULL
mfilter$cvar = NULL
matcoef = matrix(NA, ncol = 4, nrow = length(mfilter$pars))
matcoef[,1] = mfilter$pars
dimnames(matcoef) = list(allnames, c(" Estimate",
" Std. Error", " t value", "Pr(>|t|)"))
mfilter$matcoef = matcoef
mfilter$paridx = tmpidx
garchllh = sapply(filterlist@filter, FUN = function(x) x@filter$log.likelihoods)
copllh = cfilter$lik
if( maxgarchOrder > 0 ) copllh = c(rep(0, maxgarchOrder), copllh)
mfilter$llh = sum( apply(cbind(copllh, garchllh), 1, "sum") )
mfilter$zres = zres
mfilter$timer = Sys.time() - tic
ans = new("cGARCHfilter",
mfilter = mfilter,
mspec = spec@mspec,
uspec = spec@uspec,
copula = cfilter,
model = model)
return(ans)
}
.cgarchfilter.student = function(spec, data, filter.control = list(n.old = NULL), spd.control = list(lower = 0.1, upper = 0.9, type = "pwm", kernel = "epanech"),
out.sample = 0, parallel = FALSE, parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2),
VAR.fit = NULL, ...)
{
tic = Sys.time()
if( parallel ){
os = .Platform$OS.type
if(is.null(parallel.control$pkg)){
if( os == "windows" ) parallel.control$pkg = "snowfall" else parallel.control$pkg = "multicore"
if( is.null(parallel.control$cores) ) parallel.control$cores = 2
} else{
mtype = match(tolower(parallel.control$pkg[1]), c("multicore", "snowfall"))
if(is.na(mtype)) stop("\nParallel Package type not recognized in parallel.control\n")
parallel.control$pkg = tolower(parallel.control$pkg[1])
if( os == "windows" && parallel.control$pkg == "multicore" ) stop("\nmulticore not supported on windows O/S\n")
if( is.null(parallel.control$cores) ) parallel.control$cores = 2 else parallel.control$cores = as.integer(parallel.control$cores[1])
}
}
n.old = filter.control$n.old
if( is.null( colnames(data) ) ) cnames = paste("Asset_", 1:k, sep = "") else cnames = colnames(data)
xdata = .extractmdata(data)
if(!is.numeric(out.sample)) stop("\ncgarchfilter-->error: out.sample must be numeric\n")
if(as.numeric(out.sample) < 0) stop("\ncgarchfilter-->error: out.sample must be positive\n")
n.start = round(out.sample, 0)
n = dim(xdata$data)[1]
k = dim(xdata$data)[2]
if( (n - n.start) < 100) stop("\ncgarchfilter-->error: function requires at least 100 data\n points to run\n")
data = xdata$data[1:(n - n.start), , drop = FALSE]
dates = xdata$pos[1:(n - n.start)]
origdata = xdata$data
origdates = xdata$pos
dformat = xdata$dformat
vrmodel = NULL
if( spec@mspec$varmodel$VAR ){
if(spec@mspec$varmodel$mxn>0){
ex = spec@mspec$varmodel$external.regressors
if( dim(ex)[1] < dim(data)[1] ) stop("\ndccfilter-->error: external regressor matrix has less points than data!...", call. = FALSE)
orex = ex
ex = ex[1:(n - n.start), , drop = FALSE]
mxn = dim(orex)[2]
} else{
orex = NULL
ex = NULL
mxn = 0
}
if( !is.null(VAR.fit) ){
zdata = VAR.fit$xresiduals
# should be N - p
p = mp = VAR.fit$lag
N = dim(zdata)[1]
K = dim(zdata)[2]
mu = VAR.fit$xfitted
vrmodel = VAR.fit
if(out.sample>0){
Bcoef = VAR.fit$Bcoef
vrmodel2 = varxfilter(X = origdata, p = mp, exogen = orex, Bcoef = Bcoef)
zorigdata = vrmodel2$xresiduals
} else{
zorigdata = zdata
}
} else{
cat("\nestimating VAR model...")
ic = spec@mspec$varmodel$lag.criterion[1]
if( !is.null(spec@mspec$varmodel$lag.max) ){
mp = .varxselect(y = data, lag.max = spec@mspec$varmodel$lag.max, exogen = ex)$selection
mp = as.numeric( mp[which(substr(names(mp), 1, 2) == substr(ic, 1, 2))] )
} else {
mp = spec@mspec$varmodel$lag
}
vrmodel = varxfilter(X = data, p = mp, exogen = ex, Bcoef = NULL, robust = spec@mspec$varmodel$robust,
gamma = spec@mspec$varmodel$robust.control$gamma, delta = spec@mspec$varmodel$robust.control$delta,
nc = spec@mspec$varmodel$robust.control$nc, ns = spec@mspec$varmodel$robust.control$ns)
cat("done!\n")
zdata = vrmodel$xresiduals
# should be N - p
N = dim(zdata)[1]
K = dim(zdata)[2]
mu = vrmodel$xfitted
vrmodel$lag = mp
vrmodel$mxn = mxn
# if we have out of sample data we need to filter
if(out.sample>0){
Bcoef = vrmodel$Bcoef
vrmodel2 = varxfilter(X = origdata, p = mp, exogen = orex, Bcoef = Bcoef)
zorigdata = vrmodel2$xresiduals
} else{
zorigdata = zdata
}
}
} else{
zdata = data
zorigdata = origdata
orex = NULL
ex = NULL
}
if( is(spec@uspec, "uGARCHmultispec") ){
speclist = spec@uspec
} else{
stop("\ncgarchfilter-->error: cGARCHspec must contain multispec object with fixed parameters")
}
maxgarchOrder = spec@mspec$optimization.model$maxgarchOrder
filterlist = multifilter(multifitORspec = spec@uspec, data = as.data.frame(zorigdata), out.sample = out.sample,
parallel = parallel, parallel.control = parallel.control, n.old = n.old)
assign(".filterlist", filterlist, envir = .GlobalEnv)
garchpars = as.numeric( unlist(coef(filterlist) ) )
tempenvir = new.env(hash = TRUE, parent = .GlobalEnv)
assign("trace", FALSE, envir = tempenvir)
assign("method", spec@mspec$method, envir = tempenvir)
assign("cnames", cnames, envir = tempenvir)
assign("spec", spec, envir = tempenvir)
assign("parallel", parallel, envir = tempenvir)
assign("parallel.control", parallel.control, envir = tempenvir)
#assign("paridx", paridx, envir = tempenvir)
#assign("paridxi", paridxi, envir = tempenvir)
#assign("paridxa", paridxa, envir = tempenvir)
assign("speclist", speclist, envir = tempenvir)
assign("filterlist", filterlist, envir = tempenvir)
omodel = spec@mspec$optimization.model
distribution = spec@mspec$distribution
transformation = spec@mspec$transformation
sig = sigma(filterlist)
res = residuals(filterlist)
#stdresid = res/sig
assign("transformation", transformation, envir = tempenvir)
zres = sapply(filterlist@filter, FUN = function(x) x@filter$z)
ures = matrix(NA, ncol = 5, nrow = dim(zres)[1])
if( maxgarchOrder > 0 ) zres = zres[-(1:maxgarchOrder), , drop = FALSE]
assign("spd.control", spd.control, envir = tempenvir)
# Transformation of z to [0,1] domain
ans = switch(transformation,
parametric = .pparametric.filter(filterlist@filter, zres),
empirical = .pempirical(zres),
spd = .pspd(zres, spd.control))
if(transformation == "spd"){
ures = ans$ures
sfilter = ans$sfit
} else{
ures = ans
sfilter = NULL
}
# make tail adjustments in order to avoid problem with the quantile functions in
# optimization
if(any(ures > 0.99999)){
xn = which(ures > 0.99999)
ures[xn] = 0.99999
}
if(any(ures < eps)){
xn = which(ures < (1.5*eps))
ures[xn] = eps
}
assign("omodel", omodel, envir = tempenvir)
assign("mspec", spec@mspec, envir = tempenvir)
assign("ures", ures, envir = tempenvir)
ipars = .copulastart(data = ures, garchenv = tempenvir)
pars = ipars$pars
assign("npar", length(pars), envir = tempenvir)
# no fixed parameters, all normal
fixed = NULL
fixid = NULL
use.solver = TRUE
assign("fixed", fixed, envir = tempenvir)
assign("fixid", fixid, envir = tempenvir)
assign(".llh", 1, envir = tempenvir)
mfilter = list()
mfilter$dates = dates
mfilter$date.format = dformat
mfilter$origdata = origdata
mfilter$origdates = origdates
mfilter$vrmodel = vrmodel
# convert back from uniform to standard variates
model = list()
model$copula.distribution = distribution
model$out.sample = out.sample
model$transformation = transformation
model$spd.control = spd.control
model$sfit = sfilter
model$time.varying = FALSE
# add some tail data of the univariate fits to be used with simulation
mfilter$tailusigma = tail(sigma(filterlist), 25)
mfilter$tailuresids = tail(residuals(filterlist), 25)
mfilter$tailuret = tail(data, 25)
cfilter = copula.studentLLH(pars, data = ures, returnType = "ALL", garchenv = tempenvir)
garchNames = vector(mode = "character")
assetNames = cnames
dccNames = spec@mspec$optimization.model$modelnames
m = length(spec@uspec@spec)
for(i in 1:m){
cf = coef(filterlist@filter[[i]])
cfnames = names(cf)
garchNames = c(garchNames, paste(assetNames[i], cfnames, sep = ".") )
}
allnames = c(garchNames, dccNames)
mfilter$pars = c(as.numeric(coef(filterlist)), pars)
names(mfilter$pars) = allnames
mfilter$garchpars = coef(filterlist)
mfilter$stpars = pars
mfilter$scores = NULL
mfilter$A = NULL
mfilter$cvar = NULL
matcoef = matrix(NA, ncol = 4, nrow = length(mfilter$pars))
matcoef[,1] = mfilter$pars
dimnames(matcoef) = list(allnames, c(" Estimate",
" Std. Error", " t value", "Pr(>|t|)"))
mfilter$matcoef = matcoef
# remember that the filter returns the correct sign for the log.likelihoods (unlike fit which returns
# the negative of.
garchllh = sapply(filterlist@filter, FUN = function(x) x@filter$log.likelihoods)
copllh = cfilter$lik
if( maxgarchOrder > 0 ) copllh = c(rep(0, maxgarchOrder), copllh)
mfilter$llh = sum( apply(cbind(copllh, garchllh), 1, "sum") )
mfilter$zres = zres
mfilter$timer = Sys.time() - tic
ans = new("cGARCHfilter",
mfilter = mfilter,
mspec = spec@mspec,
uspec = spec@uspec,
copula = cfilter,
model = model)
return(ans)
}
.cgarchfilter.tvnorm = function(spec, data, filter.control = list(n.old = NULL), spd.control = list(lower = 0.1, upper = 0.9, type = "pwm", kernel = "epanech"),
out.sample = 0, parallel = FALSE, parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2),
VAR.fit = NULL, ...)
{
tic = Sys.time()
if( parallel ){
os = .Platform$OS.type
if(is.null(parallel.control$pkg)){
if( os == "windows" ) parallel.control$pkg = "snowfall" else parallel.control$pkg = "multicore"
if( is.null(parallel.control$cores) ) parallel.control$cores = 2
} else{
mtype = match(tolower(parallel.control$pkg[1]), c("multicore", "snowfall"))
if(is.na(mtype)) stop("\nParallel Package type not recognized in parallel.control\n")
parallel.control$pkg = tolower(parallel.control$pkg[1])
if( os == "windows" && parallel.control$pkg == "multicore" ) stop("\nmulticore not supported on windows O/S\n")
if( is.null(parallel.control$cores) ) parallel.control$cores = 2 else parallel.control$cores = as.integer(parallel.control$cores[1])
}
}
n.old = filter.control$n.old
if( is.null( colnames(data) ) ) cnames = paste("Asset_", 1:k, sep = "") else cnames = colnames(data)
xdata = .extractmdata(data)
if(!is.numeric(out.sample)) stop("\ncgarchfilter-->error: out.sample must be numeric\n")
if(as.numeric(out.sample) < 0) stop("\ncgarchfilter-->error: out.sample must be positive\n")
n.start = round(out.sample, 0)
n = dim(xdata$data)[1]
k = dim(xdata$data)[2]
if( (n - n.start) < 100) stop("\ncgarchfilter-->error: function requires at least 100 data\n points to run\n")
data = xdata$data[1:(n - n.start), , drop = FALSE]
dates = xdata$pos[1:(n - n.start)]
origdata = xdata$data
origdates = xdata$pos
dformat = xdata$dformat
vrmodel = NULL
if( spec@mspec$varmodel$VAR ){
if(spec@mspec$varmodel$mxn>0){
ex = spec@mspec$varmodel$external.regressors
if( dim(ex)[1] < dim(data)[1] ) stop("\ndccfilter-->error: external regressor matrix has less points than data!...", call. = FALSE)
orex = ex
ex = ex[1:(n - n.start), , drop = FALSE]
mxn = dim(orex)[2]
} else{
orex = NULL
ex = NULL
mxn = 0
}
if( !is.null(VAR.fit) ){
zdata = VAR.fit$xresiduals
# should be N - p
p = mp = VAR.fit$lag
N = dim(zdata)[1]
K = dim(zdata)[2]
mu = VAR.fit$xfitted
vrmodel = VAR.fit
if(out.sample>0){
Bcoef = VAR.fit$Bcoef
vrmodel2 = varxfilter(X = origdata, p = mp, exogen = orex, Bcoef = Bcoef)
zorigdata = vrmodel2$xresiduals
} else{
zorigdata = zdata
}
} else{
cat("\nestimating VAR model...")
ic = spec@mspec$varmodel$lag.criterion[1]
if( !is.null(spec@mspec$varmodel$lag.max) ){
mp = .varxselect(y = data, lag.max = spec@mspec$varmodel$lag.max, exogen = ex)$selection
mp = as.numeric( mp[which(substr(names(mp), 1, 2) == substr(ic, 1, 2))] )
} else {
mp = spec@mspec$varmodel$lag
}
vrmodel = varxfilter(X = data, p = mp, exogen = ex, Bcoef = NULL, robust = spec@mspec$varmodel$robust,
gamma = spec@mspec$varmodel$robust.control$gamma, delta = spec@mspec$varmodel$robust.control$delta,
nc = spec@mspec$varmodel$robust.control$nc, ns = spec@mspec$varmodel$robust.control$ns)
cat("done!\n")
zdata = vrmodel$xresiduals
# should be N - p
N = dim(zdata)[1]
K = dim(zdata)[2]
mu = vrmodel$xfitted
vrmodel$lag = mp
vrmodel$mxn = mxn
# if we have out of sample data we need to filter
if(out.sample>0){
Bcoef = vrmodel$Bcoef
vrmodel2 = varxfilter(X = origdata, p = mp, exogen = orex, Bcoef = Bcoef)
zorigdata = vrmodel2$xresiduals
} else{
zorigdata = zdata
}
}
} else{
zdata = data
zorigdata = origdata
orex = NULL
ex = NULL
}
if( is(spec@uspec, "uGARCHmultispec") ){
speclist = spec@uspec
} else{
stop("\ncgarchfilter-->error: cGARCHspec must contain multispec object with fixed parameters")
}
maxgarchOrder = spec@mspec$optimization.model$maxgarchOrder
filterlist = multifilter(multifitORspec = spec@uspec, data = as.data.frame(zorigdata), out.sample = out.sample,
parallel = parallel, parallel.control = parallel.control, n.old = n.old)
tmpidx = .dccparidx( speclist )
paridx = tmpidx$idx
paridxi = tmpidx$idxi
paridxa = tmpidx$idxa
garchpars = as.numeric( unlist(coef(filterlist) ) )
# create a temporary environment to store values (deleted at end of function)
tempenvir = new.env(hash = TRUE, parent = .GlobalEnv)
assign("trace", FALSE, envir = tempenvir)
assign("cnames", cnames, envir = tempenvir)
assign("spec", spec, envir = tempenvir)
assign("parallel", parallel, envir = tempenvir)
assign("parallel.control", parallel.control, envir = tempenvir)
assign("paridx", paridx, envir = tempenvir)
assign("paridxi", paridxi, envir = tempenvir)
assign("paridxa", paridxa, envir = tempenvir)
assign("speclist", speclist, envir = tempenvir)
assign("filterlist", filterlist, envir = tempenvir)
omodel = spec@mspec$optimization.model
distribution = spec@mspec$distribution
transformation = spec@mspec$transformation
sig = sigma(filterlist)
res = residuals(filterlist)
#stdresid = res/sig
zres = sapply(filterlist@filter, FUN = function(x) x@filter$z)
ures = matrix(NA, ncol = 5, nrow = dim(zres)[1])
if( maxgarchOrder > 0 ) zres = zres[-(1:maxgarchOrder), , drop = FALSE]
assign("spd.control", spd.control, envir = tempenvir)
# Transformation of z to [0,1] domain
ans = switch(transformation,
parametric = .pparametric.filter(filterlist@filter, zres),
empirical = .pempirical(zres),
spd = .pspd(zres, spd.control))
if(transformation == "spd"){
ures = ans$ures
sfilter = ans$sfit
} else{
ures = ans
sfilter = NULL
}
# make tail adjustments in order to avoid problem with the quantile functions in
# optimization
if(any(ures > 0.99999)){
xn = which(ures > 0.99999)
ures[xn] = 0.99999
}
if(any(ures < eps)){
xn = which(ures < (1.5*eps))
ures[xn] = eps
}
assign("omodel", omodel, envir = tempenvir)
assign("mspec", spec@mspec, envir = tempenvir)
assign("ures", ures, envir = tempenvir)
assign("solver", "solnp", envir = tempenvir)
ipars = .tvcopulastart(data = ures, garchenv = tempenvir)
pars = ipars$pars
assign("npar", length(pars), envir = tempenvir)
fixed = NULL
fixid = NULL
assign("fixed", fixed, envir = tempenvir)
assign("fixid", fixid, envir = tempenvir)
assign(".llh", 1, envir = tempenvir)
mfilter = list()
mfilter$dates = dates
mfilter$date.format = dformat
mfilter$origdata = origdata
mfilter$origdates = origdates
mfilter$vrmodel = vrmodel
# convert back from uniform to standard variates
model = list()
model$copula.distribution = distribution
model$out.sample = out.sample
model$transformation = transformation
model$spd.control = spd.control
model$sfilter = sfilter
model$time.varying = TRUE
# add some tail data of the univariate fits to be used with simulation
mfilter$tailusigma = tail(sigma(filterlist), 25)
mfilter$tailuresids = tail(residuals(filterlist), 25)
mfilter$tailuret = tail(data, 25)
cfilter = copula.tvnormalLLH(pars, data = ures, returnType = "ALL", garchenv = tempenvir)
garchNames = vector(mode = "character")
assetNames = cnames
dccNames = spec@mspec$optimization.model$modelnames
m = length(spec@uspec@spec)
for(i in 1:m){
cf = coef(filterlist@filter[[i]])
cfnames = names(cf)
garchNames = c(garchNames, paste(assetNames[i], cfnames, sep = ".") )
}
allnames = c(garchNames, dccNames)
mfilter$pars = c(as.numeric(coef(filterlist)), pars)
mfilter$garchpars = coef(filterlist)
mfilter$stpars = pars
names(mfilter$pars) = allnames
mfilter$scores = NULL
mfilter$A = NULL
mfilter$cvar = NULL
matcoef = matrix(NA, ncol = 4, nrow = length(mfilter$pars))
matcoef[,1] = mfilter$pars
dimnames(matcoef) = list(allnames, c(" Estimate",
" Std. Error", " t value", "Pr(>|t|)"))
mfilter$matcoef = matcoef
mfilter$paridx = tmpidx
garchllh = sapply(filterlist@filter, FUN = function(x) x@filter$log.likelihoods)
copllh = cfilter$lik
if( maxgarchOrder > 0 ) copllh = c(rep(0, maxgarchOrder), copllh)
mfilter$llh = sum( apply(cbind(copllh, garchllh), 1, "sum") )
mfilter$zres = zres
mfilter$timer = Sys.time() - tic
ans = new("cGARCHfilter",
mfilter = mfilter,
mspec = spec@mspec,
uspec = spec@uspec,
copula = cfilter,
model = model)
return(ans)
}
.cgarchfilter.norm = function(spec, data, filter.control = list(n.old = NULL), spd.control = list(lower = 0.1, upper = 0.9, type = "pwm", kernel = "epanech"),
out.sample = 0, parallel = FALSE, parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2),
VAR.fit = NULL, ...)
{
tic = Sys.time()
if( parallel ){
os = .Platform$OS.type
if(is.null(parallel.control$pkg)){
if( os == "windows" ) parallel.control$pkg = "snowfall" else parallel.control$pkg = "multicore"
if( is.null(parallel.control$cores) ) parallel.control$cores = 2
} else{
mtype = match(tolower(parallel.control$pkg[1]), c("multicore", "snowfall"))
if(is.na(mtype)) stop("\nParallel Package type not recognized in parallel.control\n")
parallel.control$pkg = tolower(parallel.control$pkg[1])
if( os == "windows" && parallel.control$pkg == "multicore" ) stop("\nmulticore not supported on windows O/S\n")
if( is.null(parallel.control$cores) ) parallel.control$cores = 2 else parallel.control$cores = as.integer(parallel.control$cores[1])
}
}
n.old = filter.control$n.old
if( is.null( colnames(data) ) ) cnames = paste("Asset_", 1:k, sep = "") else cnames = colnames(data)
xdata = .extractmdata(data)
if(!is.numeric(out.sample)) stop("\ncgarchfilter-->error: out.sample must be numeric\n")
if(as.numeric(out.sample) < 0) stop("\ncgarchfilter-->error: out.sample must be positive\n")
n.start = round(out.sample, 0)
n = dim(xdata$data)[1]
k = dim(xdata$data)[2]
if( (n - n.start) < 100) stop("\ncgarchfilter-->error: function requires at least 100 data\n points to run\n")
data = xdata$data[1:(n - n.start), , drop = FALSE]
dates = xdata$pos[1:(n - n.start)]
origdata = xdata$data
origdates = xdata$pos
dformat = xdata$dformat
vrmodel = NULL
if( spec@mspec$varmodel$VAR ){
if(spec@mspec$varmodel$mxn>0){
ex = spec@mspec$varmodel$external.regressors
if( dim(ex)[1] < dim(data)[1] ) stop("\ndccfilter-->error: external regressor matrix has less points than data!...", call. = FALSE)
orex = ex
ex = ex[1:(n - n.start), , drop = FALSE]
mxn = dim(orex)[2]
} else{
orex = NULL
ex = NULL
mxn = 0
}
if( !is.null(VAR.fit) ){
zdata = VAR.fit$xresiduals
# should be N - p
p = mp = VAR.fit$lag
N = dim(zdata)[1]
K = dim(zdata)[2]
mu = VAR.fit$xfitted
vrmodel = VAR.fit
if(out.sample>0){
Bcoef = VAR.fit$Bcoef
vrmodel2 = varxfilter(X = origdata, p = mp, exogen = orex, Bcoef = Bcoef)
zorigdata = vrmodel2$xresiduals
} else{
zorigdata = zdata
}
} else{
cat("\nestimating VAR model...")
ic = spec@mspec$varmodel$lag.criterion[1]
if( !is.null(spec@mspec$varmodel$lag.max) ){
mp = .varxselect(y = data, lag.max = spec@mspec$varmodel$lag.max, exogen = ex)$selection
mp = as.numeric( mp[which(substr(names(mp), 1, 2) == substr(ic, 1, 2))] )
} else {
mp = spec@mspec$varmodel$lag
}
vrmodel = varxfilter(X = data, p = mp, exogen = ex, Bcoef = NULL, robust = spec@mspec$varmodel$robust,
gamma = spec@mspec$varmodel$robust.control$gamma, delta = spec@mspec$varmodel$robust.control$delta,
nc = spec@mspec$varmodel$robust.control$nc, ns = spec@mspec$varmodel$robust.control$ns)
cat("done!\n")
zdata = vrmodel$xresiduals
# should be N - p
N = dim(zdata)[1]
K = dim(zdata)[2]
mu = vrmodel$xfitted
vrmodel$lag = mp
vrmodel$mxn = mxn
# if we have out of sample data we need to filter
if(out.sample>0){
Bcoef = vrmodel$Bcoef
vrmodel2 = varxfilter(X = origdata, p = mp, exogen = orex, Bcoef = Bcoef)
zorigdata = vrmodel2$xresiduals
} else{
zorigdata = zdata
}
}
} else{
zdata = data
zorigdata = origdata
orex = NULL
ex = NULL
}
if( is(spec@uspec, "uGARCHmultispec") ){
speclist = spec@uspec
} else{
stop("\ncgarchfilter-->error: cGARCHspec must contain multispec object with fixed parameters")
}
maxgarchOrder = spec@mspec$optimization.model$maxgarchOrder
filterlist = multifilter(multifitORspec = spec@uspec, data = as.data.frame(zorigdata), out.sample = out.sample,
parallel = parallel, parallel.control = parallel.control, n.old = n.old)
assign(".filterlist", filterlist, envir = .GlobalEnv)
garchpars = as.numeric( unlist(coef(filterlist) ) )
tempenvir = new.env(hash = TRUE, parent = .GlobalEnv)
assign("trace", FALSE, envir = tempenvir)
assign("method", spec@mspec$method, envir = tempenvir)
assign("cnames", cnames, envir = tempenvir)
assign("spec", spec, envir = tempenvir)
assign("parallel", parallel, envir = tempenvir)
assign("parallel.control", parallel.control, envir = tempenvir)
#assign("paridx", paridx, envir = tempenvir)
#assign("paridxi", paridxi, envir = tempenvir)
#assign("paridxa", paridxa, envir = tempenvir)
assign("speclist", speclist, envir = tempenvir)
assign("filterlist", filterlist, envir = tempenvir)
omodel = spec@mspec$optimization.model
distribution = spec@mspec$distribution
transformation = spec@mspec$transformation
sig = sigma(filterlist)
res = residuals(filterlist)
#stdresid = res/sig
assign("transformation", transformation, envir = tempenvir)
zres = sapply(filterlist@filter, FUN = function(x) x@filter$z)
ures = matrix(NA, ncol = 5, nrow = dim(zres)[1])
if( maxgarchOrder > 0 ) zres = zres[-(1:maxgarchOrder), , drop = FALSE]
assign("spd.control", spd.control, envir = tempenvir)
# Transformation of z to [0,1] domain
ans = switch(transformation,
parametric = .pparametric.filter(filterlist@filter, zres),
empirical = .pempirical(zres),
spd = .pspd(zres, spd.control))
if(transformation == "spd"){
ures = ans$ures
sfilter = ans$sfit
} else{
ures = ans
sfilter = NULL
}
# make tail adjustments in order to avoid problem with the quantile functions in
# optimization
if(any(ures > 0.99999)){
xn = which(ures > 0.99999)
ures[xn] = 0.99999
}
if(any(ures < eps)){
xn = which(ures < (1.5*eps))
ures[xn] = eps
}
assign("omodel", omodel, envir = tempenvir)
assign("mspec", spec@mspec, envir = tempenvir)
assign("ures", ures, envir = tempenvir)
ipars = .copulastart(data = ures, garchenv = tempenvir)
pars = ipars$pars
assign("npar", length(pars), envir = tempenvir)
# no fixed parameters, all normal
fixed = NULL
fixid = NULL
use.solver = TRUE
assign("fixed", fixed, envir = tempenvir)
assign("fixid", fixid, envir = tempenvir)
assign(".llh", 1, envir = tempenvir)
mfilter = list()
mfilter$dates = dates
mfilter$date.format = dformat
mfilter$origdata = origdata
mfilter$origdates = origdates
mfilter$vrmodel = vrmodel
# convert back from uniform to standard variates
model = list()
model$copula.distribution = distribution
model$out.sample = out.sample
model$transformation = transformation
model$spd.control = spd.control
model$sfit = sfilter
model$time.varying = FALSE
# add some tail data of the univariate fits to be used with simulation
mfilter$tailusigma = tail(sigma(filterlist), 25)
mfilter$tailuresids = tail(residuals(filterlist), 25)
mfilter$tailuret = tail(data, 25)
cfilter = copula.normalLLH(pars, data = ures, returnType = "ALL", garchenv = tempenvir)
garchNames = vector(mode = "character")
assetNames = cnames
dccNames = spec@mspec$optimization.model$modelnames
m = length(spec@uspec@spec)
for(i in 1:m){
cf = coef(filterlist@filter[[i]])
cfnames = names(cf)
garchNames = c(garchNames, paste(assetNames[i], cfnames, sep = ".") )
}
allnames = c(garchNames, dccNames)
mfilter$pars = c(as.numeric(coef(filterlist)), pars)
names(mfilter$pars) = allnames
mfilter$garchpars = coef(filterlist)
mfilter$stpars = pars
mfilter$scores = NULL
mfilter$A = NULL
mfilter$cvar = NULL
matcoef = matrix(NA, ncol = 4, nrow = length(mfilter$pars))
matcoef[,1] = mfilter$pars
dimnames(matcoef) = list(allnames, c(" Estimate",
" Std. Error", " t value", "Pr(>|t|)"))
mfilter$matcoef = matcoef
# remember that the filter returns the correct sign for the log.likelihoods (unlike fit which returns
# the negative of.
garchllh = sapply(filterlist@filter, FUN = function(x) x@filter$log.likelihoods)
copllh = cfilter$lik
if( maxgarchOrder > 0 ) copllh = c(rep(0, maxgarchOrder), copllh)
mfilter$llh = sum( apply(cbind(copllh, garchllh), 1, "sum") )
mfilter$zres = zres
mfilter$timer = Sys.time() - tic
ans = new("cGARCHfilter",
mfilter = mfilter,
mspec = spec@mspec,
uspec = spec@uspec,
copula = cfilter,
model = model)
return(ans)
}
.cgarchsim = function(fit, n.sim = 1000, n.start = 0, m.sim = 1, startMethod = c("unconditional", "sample"),
presigma = NULL, preresiduals = NULL, prereturns = NULL, preR = NULL, rseed = NULL, mexsimdata = NULL,
vexsimdata = NULL, parallel = FALSE, parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2),
VAR.fit = NULL, ...)
{
dist = fit@model$copula.distribution
tv = ifelse(fit@model$time.varying, 1, 0)
dx = paste(dist, tv, sep = "")
ans = switch(dx,
mvt0 = .cgarchsim1(fit = fit, n.sim = n.sim, n.start = n.start, m.sim = m.sim, startMethod = startMethod,
presigma = presigma, preresiduals = preresiduals, prereturns = prereturns, preR = preR, rseed = rseed,
mexsimdata = mexsimdata, vexsimdata = vexsimdata, parallel = parallel, parallel.control = parallel.control,
VAR.fit = VAR.fit, ...),
mvt1 = .cgarchsim2(fit = fit, n.sim = n.sim, n.start = n.start, m.sim = m.sim, startMethod = startMethod,
presigma = presigma, preresiduals = preresiduals, prereturns = prereturns, preR = preR, rseed = rseed,
mexsimdata = mexsimdata, vexsimdata = vexsimdata, parallel = parallel, parallel.control = parallel.control,
VAR.fit = VAR.fit, ...),
mvnorm0 = .cgarchsim1(fit = fit, n.sim = n.sim, n.start = n.start, m.sim = m.sim, startMethod = startMethod,
presigma = presigma, preresiduals = preresiduals, prereturns = prereturns, preR = preR, rseed = rseed,
mexsimdata = mexsimdata, vexsimdata = vexsimdata, parallel = parallel, parallel.control = parallel.control,
VAR.fit = VAR.fit, ...),
mvnorm1 = .cgarchsim2(fit = fit, n.sim = n.sim, n.start = n.start, m.sim = m.sim, startMethod = startMethod,
presigma = presigma, preresiduals = preresiduals, prereturns = prereturns, preR = preR, rseed = rseed,
mexsimdata = mexsimdata, vexsimdata = vexsimdata, parallel = parallel, parallel.control = parallel.control,
VAR.fit = VAR.fit, ...))
return( ans )
}
setMethod("cgarchsim", signature(fit = "cGARCHfit"), .cgarchsim)
# non time varying copula simulation
.cgarchsim1 = function(fit, n.sim = 1000, n.start = 0, m.sim = 1, startMethod = c("unconditional", "sample"),
presigma = NULL, preresiduals = NULL, prereturns = NULL, preR = NULL, rseed = NULL, mexsimdata = NULL,
vexsimdata = NULL, parallel = FALSE, parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2),
VAR.fit = NULL, ...)
{
# first generate the copula random uniform numbers (static copula)
# ures --> transform --> zres
# pass zres to GARCH simulation
# pass if needed to VAR model
if( is.null(rseed) ){
rseed = as.integer(runif(m.sim, 1, Sys.time()))
} else {
if(length(rseed) == 1) c(rseed, as.integer(runif(m.sim-1, 1, Sys.time()))) else rseed = as.integer( rseed[1:m.sim] )
}
Data = head(fit@mfit$origdata, dim(fit@mfit$origdata)[1] - fit@model$out.sample)
m = dim(Data)[2]
n = dim(Data)[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])
}
}
nsim = n.sim + n.start
spec = fit@uspec
cf = coef(fit, "garch")
cfl = vector(mode = "list", length = m)
for(i in 1:m){
if(spec@type == "equal"){
spec@spec[[i]]@optimization.model$fixed.pars = as.list(cf[,i])
cfl[[i]] = cf[,i]
} else{
spec@spec[[i]]@optimization.model$fixed.pars = as.list(cf[[i]])
cfl[[i]] = cf[[i]]
}
}
if(is.null(preR)){
Rbar = fit@copula$Rt
} else{
if(dim(preR)[1] != m) stop("\ncgarchsim-->error: wrong dimension for preR\n")
if(dim(preR)[2] != m) stop("\ncgarchsim-->error: wrong dimension for preR\n")
if(any(diag(preR) != 1)) stop("\ncgarchsim-->error: preR diagonals must be 1.\n")
Rbar = preR
}
dist = fit@model$copula.distribution
dx = paste(dist, 0, sep = "")
ures = .sample.copula(fit, Rbar, dist = dx, n.sim = n.sim, n.start = n.start, m.sim = m.sim, rseed = rseed)
# ures = [0,1] copula random numbers
# now transform back into margins for use in garch
transformation = fit@model$transformation
include.skew = sapply(fit@uspec@spec, FUN = function(x) x@distribution.model$include.skew)
include.shape = sapply(fit@uspec@spec, FUN = function(x) x@distribution.model$include.shape)
distx = sapply(fit@uspec@spec, FUN = function(x) x@distribution.model$distribution)
zres = array(NA, dim = c(nsim, m, m.sim))
# parallel:
if( parallel ){
if( parallel.control$pkg == "multicore" ){
mtmp = mclapply(as.list(1:m.sim), FUN = function(i) switch(transformation,
parametric = .qparametric(ures[,,i], ucoef = cfl, include.skew, include.shape, dist = distx),
empirical = .qempirical(matrix(ures[,,i], ncol = m), fit@mfit$zres),
spd = .qspd(ures[,,i], sfit = fit@model$sfit)), mc.cores = parallel.control$cores)
for(i in 1:m.sim) zres[,,i] = mtmp[[i]]
} else{
ssfit = fit@model$sfit
sfInit(parallel = TRUE, cpus = parallel.control$cores)
sfExport("ures", "cfl", "include.skew", "include.shape", "distx", "ssfit", "transformation", local = TRUE)
mtmp = sfLapply(as.list(1:m.sim), fun = function(i)
switch(transformation,
parametric = rgarch:::.qparametric(ures[,,i], ucoef = cfl, include.skew, include.shape, dist = distx),
empirical = rgarch:::.qempirical(matrix(ures[,,i], ncol = m), fit@mfit$zres),
spd = rgarch:::.qspd(matrix(ures[,,i], ncol = m), sfit = ssfit)))
for(i in 1:m.sim) zres[,,i] = mtmp[[i]]
#sfStop()
}
} else{
mtmp = lapply(as.list(1:m.sim), FUN = function(i) switch(transformation,
parametric = .qparametric(matrix(ures[,,i], ncol = m), ucoef = cfl, include.skew, include.shape, dist = distx),
empirical = .qempirical(matrix(ures[,,i], ncol = m), fit@mfit$zres),
spd = .qspd(matrix(ures[,,i], ncol = m), sfit = fit@model$sfit)))
for(i in 1:m.sim) zres[,,i] = mtmp[[i]]
}
# Now we have the innovations which we feed back into the garch routine to simulate
simRes = simX = simR = simQ = simH = simSeries = vector(mode = "list", length = m.sim)
if( is.null(fit@mfit$vrmodel) ){
if(is.null(prereturns)) prereturns = fit@mfit$tailuret
}
if( parallel ){
if( parallel.control$pkg == "multicore" ){
simlist = mclapply(as.list(1:m), FUN = function(i){
maxx = spec@spec[[i]]@optimization.model$maxOrder2;
htmp = ugarchpath(spec@spec[[i]], n.sim = n.sim + n.start, n.start = 0, m.sim = m.sim,
custom.dist = list(name = "sample", distfit = matrix(zres[,i,1:m.sim], ncol = m.sim)),
presigma = if( is.null(presigma) ) tail(fit@mfit$tailusigma[,i], maxx) else tail(presigma[,i], maxx),
preresiduals = if( is.null(preresiduals) ) tail(fit@mfit$tailuresids[,i], maxx) else tail(preresiduals[,i], maxx),
prereturns = if(is.null(fit@mfit$vrmodel)) tail(prereturns[,i], maxx) else NA,
mexsimdata = if( is.null(fit@mfit$vrmodel) ) mexsimdata[[i]] else NULL, vexsimdata = vexsimdata[[i]] )
h = matrix(tail(htmp@path$sigmaSim^2, n.sim), nrow = n.sim)
x = matrix(htmp@path$seriesSim, nrow = n.sim + n.start)
return(list(h = h, x = x))}, mc.cores = parallel.control$cores)
H = array(NA, dim = c(n.sim, m, m.sim))
tmpH = array(NA, dim = c(m, m, n.sim))
for(i in 1:n.sim) H[i,,] = t(sapply(simlist, FUN = function(x) as.numeric(x$h[i,])))
for(i in 1:m.sim){
for(j in 1:n.sim){
tmpH[ , , j] = diag(sqrt( H[j, , i]) ) %*% fit@copula$Rt %*% diag(sqrt( H[j, , i] ) )
}
simH[[i]] = tmpH
}
simxX = array(NA, dim = c(n.sim, m, m.sim))
for(i in 1:n.sim) simxX[i,,] = t(sapply(simlist, FUN = function(x) as.numeric(x$x[i,])))
simX = vector(mode = "list", length = m.sim)
for(i in 1:m.sim) simX[[i]] = matrix(simxX[,,i], nrow = n.sim)
} else{
if(is.null(fit@mfit$vrmodel)) include.vrmodel = FALSE else include.vrmodel = TRUE
sfExport("spec", "n.sim", "n.start", "m.sim", "startMethod", "zres", "presigma",
"preresiduals", "prereturns", "fit", "mexsimdata", "vexsimdata", "include.vrmodel", local = TRUE)
simlist = sfLapply(as.list(1:m), fun = function(i){
maxx = spec@spec[[i]]@optimization.model$maxOrder2;
htmp = rgarch::ugarchpath(spec@spec[[i]], n.sim = n.sim + n.start, n.start = 0, m.sim = m.sim,
custom.dist = list(name = "sample", distfit = matrix(zres[,i,1:m.sim], ncol = m.sim)),
presigma = if( is.null(presigma) ) tail(fit@mfit$tailusigma[,i], maxx) else tail(presigma[,i], maxx),
preresiduals = if( is.null(preresiduals) ) tail(fit@mfit$tailuresids[,i], maxx) else tail(preresiduals[,i], maxx),
prereturns = if(is.null(fit@mfit$vrmodel)) tail(prereturns[,i], maxx) else NA,
mexsimdata = if( is.null(fit@mfit$vrmodel) ) mexsimdata[[i]] else NULL, vexsimdata = vexsimdata[[i]] )
h = matrix(tail(htmp@path$sigmaSim^2, n.sim), nrow = n.sim)
x = matrix(htmp@path$seriesSim, nrow = n.sim + n.start)
return(list(h = h, x = x))})
if( !include.vrmodel ) sfStop()
H = array(NA, dim = c(n.sim, m, m.sim))
tmpH = array(NA, dim = c(m, m, n.sim))
for(i in 1:n.sim) H[i,,] = t(sapply(simlist, FUN = function(x) as.numeric(x$h[i,])))
for(i in 1:m.sim){
for(j in 1:n.sim){
tmpH[ , , j] = diag(sqrt( H[j, , i]) ) %*% fit@copula$Rt %*% diag(sqrt( H[j, , i] ) )
}
simH[[i]] = tmpH
}
simxX = array(NA, dim = c(n.sim, m, m.sim))
for(i in 1:n.sim) simxX[i,,] = t(sapply(simlist, FUN = function(x) as.numeric(x$x[i,])))
simX = vector(mode = "list", length = m.sim)
for(i in 1:m.sim) simX[[i]] = matrix(simxX[,,i], nrow = n.sim)
}
} else{
simlist = lapply(as.list(1:m), FUN = function(i){
maxx = spec@spec[[i]]@optimization.model$maxOrder2;
htmp = ugarchpath(spec@spec[[i]], n.sim = n.sim + n.start, n.start = 0, m.sim = m.sim,
custom.dist = list(name = "sample", distfit = matrix(zres[,i,1:m.sim], ncol = m.sim)),
presigma = if( is.null(presigma) ) tail(fit@mfit$tailusigma[,i], maxx) else tail(presigma[,i], maxx),
preresiduals = if( is.null(preresiduals) ) tail(fit@mfit$tailuresids[,i], maxx) else tail(preresiduals[,i], maxx),
prereturns = if(is.null(fit@mfit$vrmodel)) tail(prereturns[,i], maxx) else NA,
mexsimdata = if( is.null(fit@mfit$vrmodel) ) mexsimdata[[i]] else NULL, vexsimdata = vexsimdata[[i]] )
h = matrix(tail(htmp@path$sigmaSim^2, n.sim), nrow = n.sim)
x = matrix(htmp@path$seriesSim, nrow = n.sim + n.start)
return(list(h = h, x = x))})
H = array(NA, dim = c(n.sim, m, m.sim))
tmpH = array(NA, dim = c(m, m, n.sim))
for(i in 1:n.sim) H[i,,] = t(sapply(simlist, FUN = function(x) as.numeric(x$h[i,])))
for(i in 1:m.sim){
for(j in 1:n.sim){
tmpH[ , , j] = diag(sqrt( H[j, , i]) ) %*% fit@copula$Rt %*% diag(sqrt( H[j, , i] ) )
}
simH[[i]] = tmpH
}
simxX = array(NA, dim = c(n.sim, m, m.sim))
for(i in 1:n.sim) simxX[i,,] = t(sapply(simlist, FUN = function(x) as.numeric(x$x[i,])))
simX = vector(mode = "list", length = m.sim)
for(i in 1:m.sim) simX[[i]] = matrix(simxX[,,i], nrow = n.sim)
}
if( !is.null(fit@mfit$vrmodel) ){
# mexsimdata check
mxn = fit@mfit$vrmodel$mxn
if(mxn > 0 ){
if( !is.null(mexsimdata) ){
if( !is.list(mexsimdata) | length(mexsimdata) != m.sim ) stop("\ncgarchsim-->error: mexsimdata must be a list of length equal to m.sim...", call. = FALSE)
for(j in 1:m.sim){
if( !is.matrix(mexsimdata[[j]]) ) stop("\ncgarchsim-->error: mexsimdata list submatrix must be a matrix...", call. = FALSE)
if( dim(mexsimdata[[j]])[2] != mxn ) stop("\ncgarchsim-->error: wrong number of external regressors in list submatrix!...", call. = FALSE)
if( dim(mexsimdata[[j]])[1] < (n.sim + n.start) ) stop("\ncgarchsim-->error: external regressor list submatrix has less points than requested simulation length (n.sim + n.start)...", call. = FALSE)
}
} else{
mexsimdata = vector(mode = "list", length = m.sim)
for(j in 1:m.sim) mexsimdata[[i]] = matrix(0, ncol = mxn, nrow = (n.sim + n.start))
}
} else{
mexsimdata = NULL
}
vrmodel = fit@mfit$vrmodel
p = as.numeric(vrmodel$lag)
if(startMethod[1] == "sample"){
if( !is.null(prereturns) ){
if(!is.matrix(prereturns)) stop("\nprereturns must be a matrix\n")
if(dim(prereturns)[1] != p) stop("\nwrong number of rows for prereturns matrix")
if(dim(prereturns)[2] != m) stop("\nwrong number of cols for prereturns matrix")
prereturns = matrix(prereturns, nrow = p, ncol = m)
} else{
prereturns = as.matrix( tail(Data, p) )
}
} else{
prereturns = matrix(0, nrow = p, ncol = m, byrow = TRUE)
}
xlist = vector(mode = "list", length = m.sim)
if( n.sim == 1 && n.start == 0 ){
# SPECIAL ROUTINE TO SPEED UP 1-AHEAD ESTIMATION
tmp = varxsimXX(X = Data, vrmodel$Bcoef, p, m.sim = m.sim, prereturns, resids = t(sapply(simX, FUN = function(x) x)), if(!is.null(mexsimdata)) t(sapply(mexsimdata, FUN = function(x) x)) else NULL)
for(i in 1:m.sim) xlist[[i]] = tmp[i, , drop = FALSE]
simRes = simX
simX = xlist
} else{
if( parallel ){
if( parallel.control$pkg == "multicore" ){
xlist = mclapply(as.list(1:m.sim), FUN = function(i) varxsim(X = Data, vrmodel$Bcoef, p, n.sim,
n.start, prereturns, resids = simX[[i]], mexsimdata[[i]]),
mc.cores = parallel.control$cores)
} else{
# when the GARCH model is without mean equation, simX == simulated residuals
sfExport("Data", "vrmodel", "p", "n.sim", "n.start", "prereturns", "simX", "mexsimdata", local = TRUE)
xlist = sfLapply(as.list(1:m.sim), fun = function(i) rgarch:::varxsim(X = Data, vrmodel$Bcoef, p,
n.sim, n.start, prereturns, resids = simX[[i]], mexsimdata[[i]]))
sfStop()
}
} else{
xlist = lapply(as.list(1:m.sim), FUN = function(i) varxsim(X = Data, vrmodel$Bcoef, p,
n.sim, n.start, prereturns, resids = simX[[i]], mexsimdata[[i]]))
}
simRes = simX
simX = xlist
}
} else{
# reshape
for(j in 1:m.sim) simX[[j]] = tail(simX[[j]], n.sim)
}
msim = list()
msim$simH = simH
msim$simX = simX
msim$simRes = simRes
#msim$Z = z
msim$rseed = rseed
msim$model$Data = Data
#msim$model$garchpars = garchpars
msim$model$n.sim = n.sim
msim$model$m.sim = m.sim
msim$model$n.start = n.start
msim$model$startMethod = startMethod[1]
msim$model$distribution = fit@model$copula.distribution
msim$model$time.varying = FALSE
msim$model$transformation = transformation
ans = new("cGARCHsim",
msim = msim)
return( ans )
}
# time varying copula simulation
.cgarchsim2 = function(fit, n.sim = 1000, n.start = 0, m.sim = 1, startMethod = c("unconditional", "sample"),
presigma = NULL, preresiduals = NULL, prereturns = NULL, preR = NULL, rseed = NULL, mexsimdata = NULL,
vexsimdata = NULL, parallel = FALSE, parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2),
VAR.fit = NULL, ...)
{
# first generate the copula random uniform numbers (static copula)
# ures --> transform --> zres
# pass zres to GARCH simulation
# pass if needed to VAR model
if( is.null(rseed) ){
rseed = as.integer(runif(m.sim, 1, Sys.time()))
} else {
if(length(rseed) == 1) c(rseed, as.integer(runif(m.sim-1, 1, Sys.time()))) else rseed = as.integer( rseed[1:m.sim] )
}
Data = head(fit@mfit$origdata, dim(fit@mfit$origdata)[1] - fit@model$out.sample)
m = dim(Data)[2]
n = dim(Data)[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])
}
}
nsim = n.sim + n.start
spec = fit@uspec
cf = coef(fit, "garch")
cfl = vector(mode = "list", length = m)
for(i in 1:m){
if(spec@type == "equal"){
spec@spec[[i]]@optimization.model$fixed.pars = as.list(cf[,i])
cfl[[i]] = cf[,i]
} else{
spec@spec[[i]]@optimization.model$fixed.pars = as.list(cf[[i]])
cfl[[i]] = cf[[i]]
}
}
if(is.null(preR)){
Rbar = fit@copula$Rt[[length(fit@copula$Rt)]]
} else{
if(dim(preR)[1] != m) stop("\ncgarchsim-->error: wrong dimension for preR\n")
if(dim(preR)[2] != m) stop("\ncgarchsim-->error: wrong dimension for preR\n")
if(any(diag(preR) != 1)) stop("\ncgarchsim-->error: preR diagonals must be 1.\n")
Rbar = preR
}
dist = fit@model$copula.distribution
dx = paste(dist, 1, sep = "")
smplx = .sample.copula(dist = dx, object = fit, Rbar, n.sim, n.start, m.sim, rseed, parallel, parallel.control)
# ures = [0,1] copula random numbers
# now transform back into margins for use in garch
Usim = smplx$Usim
Rt = smplx$simR
transformation = fit@model$transformation
include.skew = sapply(fit@uspec@spec, FUN = function(x) x@distribution.model$include.skew)
include.shape = sapply(fit@uspec@spec, FUN = function(x) x@distribution.model$include.shape)
distx = sapply(fit@uspec@spec, FUN = function(x) x@distribution.model$distribution)
zres = array(NA, dim = c(nsim, m, m.sim))
# parallel:
if( parallel ){
if( parallel.control$pkg == "multicore" ){
mtmp = mclapply(as.list(1:m.sim), FUN = function(i) switch(transformation,
parametric = .qparametric(matrix(Usim[,,i], ncol = m), ucoef = cfl, include.skew, include.shape, dist = distx),
empirical = .qempirical(matrix(Usim[,,i], ncol = m), fit@mfit$zres),
spd = .qspd(matrix(Usim[,,i], ncol = m), sfit = fit@model$sfit)), mc.cores = parallel.control$cores)
for(i in 1:m.sim) zres[,,i] = matrix(mtmp[[i]], ncol = m)
} else{
ssfit = fit@model$sfit
sfInit(parallel = TRUE, cpus = parallel.control$cores)
sfExport("Usim", "cfl", "include.skew", "include.shape", "distx", "ssfit", "transformation", local = TRUE)
mtmp = sfLapply(as.list(1:m.sim), fun = function(i)
switch(transformation,
parametric = rgarch:::.qparametric(matrix(Usim[,,i], ncol = m), ucoef = cfl, include.skew, include.shape, dist = distx),
empirical = rgarch:::.qempirical(matrix(Usim[,,i], ncol = m), fit@mfit$zres),
spd = rgarch:::.qspd(matrix(Usim[,,i], ncol = m), sfit = ssfit)))
for(i in 1:m.sim) zres[,,i] = matrix(mtmp[[i]], ncol = m)
#sfStop()
}
} else{
mtmp = lapply(as.list(1:m.sim), FUN = function(i) switch(transformation,
parametric = .qparametric(matrix(Usim[,,i], ncol = m), ucoef = cfl, include.skew, include.shape, dist = distx),
empirical = .qempirical(matrix(Usim[,,i], ncol = m), fit@mfit$zres),
spd = .qspd(matrix(Usim[,,i], ncol = m), sfit = fit@model$sfit)))
for(i in 1:m.sim) zres[,,i] = matrix(mtmp[[i]], ncol = m)
}
# Now we have the innovations which we feed back into the garch routine to simulate
simRes = simX = simH = simSeries = vector(mode = "list", length = m.sim)
if( is.null(fit@mfit$vrmodel) ){
if(is.null(prereturns)) prereturns = fit@mfit$tailuret
}
if( parallel ){
if( parallel.control$pkg == "multicore" ){
simlist = mclapply(as.list(1:m), FUN = function(i){
maxx = spec@spec[[i]]@optimization.model$maxOrder2;
htmp = ugarchpath(spec@spec[[i]], n.sim = n.sim + n.start, n.start = 0, m.sim = m.sim,
custom.dist = list(name = "sample", distfit = matrix(zres[,i,1:m.sim], ncol = m.sim)),
presigma = if( is.null(presigma) ) tail(fit@mfit$tailusigma[,i], maxx) else tail(presigma[,i], maxx),
preresiduals = if( is.null(preresiduals) ) tail(fit@mfit$tailuresids[,i], maxx) else tail(preresiduals[,i], maxx),
prereturns = if(is.null(fit@mfit$vrmodel)) tail(prereturns[,i], maxx) else NA,
mexsimdata = if( is.null(fit@mfit$vrmodel) ) mexsimdata[[i]] else NULL, vexsimdata = vexsimdata[[i]] )
h = matrix(tail(htmp@path$sigmaSim^2, n.sim), nrow = n.sim)
x = matrix(htmp@path$seriesSim, nrow = n.sim + n.start)
return(list(h = h, x = x))}, mc.cores = parallel.control$cores)
H = array(NA, dim = c(n.sim, m, m.sim))
tmpH = array(NA, dim = c(m, m, n.sim))
for(i in 1:n.sim) H[i,,] = t(sapply(simlist, FUN = function(x) as.numeric(x$h[i,])))
for(i in 1:m.sim){
for(j in 1:n.sim){
tmpH[ , , j] = diag(sqrt( H[j, , i]) ) %*% Rt[[i]][,,j] %*% diag(sqrt( H[j, , i] ) )
}
simH[[i]] = tmpH
}
simxX = array(NA, dim = c(n.sim, m, m.sim))
for(i in 1:n.sim) simxX[i,,] = t(sapply(simlist, FUN = function(x) as.numeric(x$x[i,])))
simX = vector(mode = "list", length = m.sim)
for(i in 1:m.sim) simX[[i]] = matrix(simxX[,,i], nrow = n.sim)
} else{
if(is.null(fit@mfit$vrmodel)) include.vrmodel = FALSE else include.vrmodel = TRUE
sfExport("spec", "n.sim", "n.start", "m.sim", "startMethod", "zres", "presigma",
"preresiduals", "prereturns", "fit", "mexsimdata", "vexsimdata", "include.vrmodel", local = TRUE)
simlist = sfLapply(as.list(1:m), fun = function(i){
maxx = spec@spec[[i]]@optimization.model$maxOrder2;
htmp = rgarch::ugarchpath(spec@spec[[i]], n.sim = n.sim + n.start, n.start = 0, m.sim = m.sim,
custom.dist = list(name = "sample", distfit = matrix(zres[,i,1:m.sim], ncol = m.sim)),
presigma = if( is.null(presigma) ) tail(fit@mfit$tailusigma[,i], maxx) else tail(presigma[,i], maxx),
preresiduals = if( is.null(preresiduals) ) tail(fit@mfit$tailuresids[,i], maxx) else tail(preresiduals[,i], maxx),
prereturns = if(is.null(fit@mfit$vrmodel)) tail(prereturns[,i], maxx) else NA,
mexsimdata = if( is.null(fit@mfit$vrmodel) ) mexsimdata[[i]] else NULL, vexsimdata = vexsimdata[[i]] )
h = matrix(tail(htmp@path$sigmaSim^2, n.sim), nrow = n.sim)
x = matrix(htmp@path$seriesSim, nrow = n.sim + n.start)
return(list(h = h, x = x))})
if( !include.vrmodel ) sfStop()
H = array(NA, dim = c(n.sim, m, m.sim))
tmpH = array(NA, dim = c(m, m, n.sim))
for(i in 1:n.sim) H[i,,] = t(sapply(simlist, FUN = function(x) as.numeric(x$h[i,])))
for(i in 1:m.sim){
for(j in 1:n.sim){
tmpH[ , , j] = diag(sqrt( H[j, , i]) ) %*% Rt[[i]][,,j] %*% diag(sqrt( H[j, , i] ) )
}
simH[[i]] = tmpH
}
simxX = array(NA, dim = c(n.sim, m, m.sim))
for(i in 1:n.sim) simxX[i,,] = t(sapply(simlist, FUN = function(x) as.numeric(x$x[i,])))
simX = vector(mode = "list", length = m.sim)
for(i in 1:m.sim) simX[[i]] = matrix(simxX[,,i], nrow = n.sim)
}
} else{
simlist = lapply(as.list(1:m), FUN = function(i){
maxx = spec@spec[[i]]@optimization.model$maxOrder2;
htmp = ugarchpath(spec@spec[[i]], n.sim = n.sim + n.start, n.start = 0, m.sim = m.sim,
custom.dist = list(name = "sample", distfit = matrix(zres[,i,1:m.sim], ncol = m.sim)),
presigma = if( is.null(presigma) ) tail(fit@mfit$tailusigma[,i], maxx) else tail(presigma[,i], maxx),
preresiduals = if( is.null(preresiduals) ) tail(fit@mfit$tailuresids[,i], maxx) else tail(preresiduals[,i], maxx),
prereturns = if(is.null(fit@mfit$vrmodel)) tail(prereturns[,i], maxx) else NA,
mexsimdata = if( is.null(fit@mfit$vrmodel) ) mexsimdata[[i]] else NULL, vexsimdata = vexsimdata[[i]] )
h = matrix(tail(htmp@path$sigmaSim^2, n.sim), nrow = n.sim)
x = matrix(htmp@path$seriesSim, nrow = n.sim + n.start)
return(list(h = h, x = x))})
H = array(NA, dim = c(n.sim, m, m.sim))
tmpH = array(NA, dim = c(m, m, n.sim))
for(i in 1:n.sim) H[i,,] = t(sapply(simlist, FUN = function(x) as.numeric(x$h[i,])))
for(i in 1:m.sim){
for(j in 1:n.sim){
tmpH[ , , j] = diag(sqrt( H[j, , i]) ) %*% Rt[[i]][,,j] %*% diag(sqrt( H[j, , i] ) )
}
simH[[i]] = tmpH
}
simxX = array(NA, dim = c(n.sim, m, m.sim))
for(i in 1:n.sim) simxX[i,,] = t(sapply(simlist, FUN = function(x) as.numeric(x$x[i,])))
simX = vector(mode = "list", length = m.sim)
for(i in 1:m.sim) simX[[i]] = matrix(simxX[,,i], nrow = n.sim)
}
if( !is.null(fit@mfit$vrmodel) ){
# mexsimdata check
mxn = fit@mfit$vrmodel$mxn
if(mxn > 0 ){
if( !is.null(mexsimdata) ){
if( !is.list(mexsimdata) | length(mexsimdata) != m.sim ) stop("\ncgarchsim-->error: mexsimdata must be a list of length equal to m.sim...", call. = FALSE)
for(j in 1:m.sim){
if( !is.matrix(mexsimdata[[j]]) ) stop("\ncgarchsim-->error: mexsimdata list submatrix must be a matrix...", call. = FALSE)
if( dim(mexsimdata[[j]])[2] != mxn ) stop("\ncgarchsim-->error: wrong number of external regressors in list submatrix!...", call. = FALSE)
if( dim(mexsimdata[[j]])[1] < (n.sim + n.start) ) stop("\ncgarchsim-->error: external regressor list submatrix has less points than requested simulation length (n.sim + n.start)...", call. = FALSE)
}
} else{
mexsimdata = vector(mode = "list", length = m.sim)
for(j in 1:m.sim) mexsimdata[[i]] = matrix(0, ncol = mxn, nrow = (n.sim + n.start))
}
} else{
mexsimdata = NULL
}
vrmodel = fit@mfit$vrmodel
p = as.numeric(vrmodel$lag)
if(startMethod[1] == "sample"){
if( !is.null(prereturns) ){
if(!is.matrix(prereturns)) stop("\nprereturns must be a matrix\n")
if(dim(prereturns)[1] != p) stop("\nwrong number of rows for prereturns matrix")
if(dim(prereturns)[2] != m) stop("\nwrong number of cols for prereturns matrix")
prereturns = matrix(prereturns, nrow = p, ncol = m)
} else{
prereturns = as.matrix( tail(Data, p) )
}
} else{
prereturns = matrix(0, nrow = p, ncol = m, byrow = TRUE)
}
xlist = vector(mode = "list", length = m.sim)
if( n.sim == 1 && n.start == 0 ){
# SPECIAL ROUTINE TO SPEED UP 1-AHEAD ESTIMATION
tmp = varxsimXX(X = Data, vrmodel$Bcoef, p, m.sim = m.sim, prereturns, resids = t(sapply(simX, FUN = function(x) x)), if(!is.null(mexsimdata)) t(sapply(mexsimdata, FUN = function(x) x)) else NULL)
for(i in 1:m.sim) xlist[[i]] = tmp[i, , drop = FALSE]
simRes = simX
simX = xlist
} else{
if( parallel ){
if( parallel.control$pkg == "multicore" ){
xlist = mclapply(as.list(1:m.sim), FUN = function(i) varxsim(X = Data, vrmodel$Bcoef, p, n.sim,
n.start, prereturns, resids = simX[[i]], mexsimdata[[i]]),
mc.cores = parallel.control$cores)
} else{
# when the GARCH model is without mean equation, simX == simulated residuals
sfExport("Data", "vrmodel", "p", "n.sim", "n.start", "prereturns", "simX", "mexsimdata", local = TRUE)
xlist = sfLapply(as.list(1:m.sim), fun = function(i) rgarch:::varxsim(X = Data, vrmodel$Bcoef, p,
n.sim, n.start, prereturns, resids = simX[[i]], mexsimdata[[i]]))
sfStop()
}
} else{
xlist = lapply(as.list(1:m.sim), FUN = function(i) varxsim(X = Data, vrmodel$Bcoef, p,
n.sim, n.start, prereturns, resids = simX[[i]], mexsimdata[[i]]))
}
simRes = simX
simX = xlist
}
} else{
# reshape
for(j in 1:m.sim) simX[[j]] = tail(simX[[j]], n.sim)
}
msim = list()
msim$simH = simH
msim$simR = Rt
msim$simX = simX
msim$simRes = simRes
#msim$Z = z
msim$rseed = rseed
msim$model$Data = Data
#msim$model$garchpars = garchpars
msim$model$n.sim = n.sim
msim$model$m.sim = m.sim
msim$model$n.start = n.start
msim$model$startMethod = startMethod[1]
msim$model$distribution = fit@model$copula.distribution
msim$model$time.varying = TRUE
msim$model$transformation = transformation
ans = new("cGARCHsim",
msim = msim)
return( ans )
}
Try the rgarch package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.