Nothing
#################################################################################
##
## R package rgarch by Alexios Ghalanos Copyright (C) 2008, 2009, 2010, 2011
## This file is part of the R package rgarch.
##
## The R package rgarch is free software: you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation, either version 3 of the License, or
## (at your option) any later version.
##
## The R package rgarch is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
## GNU General Public License for more details.
##
#################################################################################
.dccspec = function(uspec, VAR = FALSE, VAR.opt = list(robust = FALSE, lag = 1, lag.max = NULL,
lag.criterion = c("AIC", "HQ", "SC", "FPE"), external.regressors = NULL,
robust.control = list("gamma" = 0.25, "delta" = 0.01, "nc" = 10, "ns" = 500)),
dccOrder = c(1,1), distribution = c("mvnorm", "mvt", "mvlaplace"),
start.pars = list(), fixed.pars = list())
{
if(is.null(distribution)) distribution = "mvnorm"
distribution = distribution[1]
valid.distributions = c("mvnorm", "mvt", "mvlaplace")
if(!any(distribution == valid.distributions)) stop("\nInvalid Distribution Choice\n", call. = FALSE)
if(is.null(dccOrder)) dccOrder =c(1,1)
if(length(dccOrder)!=2) stop("\nInvalid dccOrder\n", call. = FALSE)
dccOrder = as.integer(dccOrder)
maxdccOrder = max(dccOrder)
maxgarchOrder = max( sapply(uspec@spec, FUN = function(x) max(max(x@mean.model$armaOrder), max(x@variance.model$garchOrder) ) ) )
maxOrder = max( maxdccOrder, maxgarchOrder)
shape.distributions = c("mvt")
skew.distributions = NULL
include.skew = include.shape = FALSE
# we will have to change this eventually to accomodate vector
# skewness and shape parameters for more interesting distributions
if(any(distribution == shape.distributions)) include.shape = TRUE
if(any(distribution == skew.distributions)) include.skew = TRUE
mspec = list()
mspec$dccOrder = dccOrder
mspec$start.pars = start.pars
mspec$fixed.pars = fixed.pars
mspec$distribution = distribution
mspec$include.shape = include.shape
mspec$include.skew = include.skew
modelnames = c(
if(dccOrder[1] > 0) paste("dcca", 1:dccOrder[1], sep = ""),
if(dccOrder[2] > 0) paste("dccb", 1:dccOrder[2], sep = ""),
if(include.skew) paste("skew"),
if(include.shape) paste("shape"))
pos = 1
pos.matrix = matrix(0, ncol = 4, nrow = 4)
colnames(pos.matrix) = c("start", "stop", "include", "npars")
rownames(pos.matrix) = c("dcca","dccb","skew","shape")
if(dccOrder[1] > 0){
pos.matrix[1,1:3] = c(pos, pos+dccOrder[1]-1, 1)
pos = max(pos.matrix[1:1, 2]) + 1
}
if(dccOrder[2] > 0){
pos.matrix[2,1:3] = c(pos, pos+dccOrder[2]-1, 1)
pos = max(pos.matrix[1:2, 2]) + 1
}
if(include.skew)
{
pos.matrix[3,1:3] = c(pos, pos, 1)
pos = max(pos.matrix[1:3, 2]) + 1
}
if(include.shape)
{
pos.matrix[4,1:3] = c(pos, pos, 1)
}
nn = length(modelnames)
modelmatrix = matrix(0, ncol = 3, nrow = nn)
rownames(modelmatrix) = modelnames
colnames(modelmatrix) = c("opt", "fixed", "start")
fixed.names = names(fixed.pars)
fp = charmatch(fixed.names, modelnames)
if(!is.null(fixed.names) && any(!is.na(fp))){
fixed = fp[!is.na(fp)]
modelmatrix[fixed,2] = 1
fz = charmatch(modelnames, fixed.names)
fz = fz[!is.na(fz)]
fixed.pars = fixed.pars[fz]
names(fixed.pars) = fixed.names[fz]
} else{
fixed.pars = NULL
}
modelmatrix[,1] = 1 - modelmatrix[,2]
start.names = names(start.pars)
sp = charmatch(start.names, modelnames)
if(!is.null(start.names) && any(!is.na(sp))){
start = sp[!is.na(sp)]
modelmatrix[start,3] = 1
sz = charmatch(modelnames, start.names)
sz = sz[!is.na(sz)]
start.pars = start.pars[sz]
} else{
start.pars = NULL
}
opt.pars = vector(mode = "numeric", length = length(modelnames))
names(opt.pars) = modelnames
# optimization model
omodel = list(modelnames = modelnames, opt.pars = opt.pars, start.pars = start.pars,
fixed.pars = fixed.pars, modelmatrix = modelmatrix, maxOrder = maxOrder,
maxgarchOrder = maxgarchOrder, maxdccOrder = maxdccOrder,
pos.matrix = pos.matrix)
mspec$optimization.model = omodel
if( !is(uspec, "uGARCHmultispec") ) stop("\ndccspec-->error: uspec must be a uGARCHmultispec object")
mspec$varmodel = list()
if( VAR ){
mspec$varmodel$VAR = TRUE
mspec$varmodel$robust = VAR.opt$robust
mspec$varmodel$lag = VAR.opt$lag
mspec$varmodel$lag.max = VAR.opt$lag.max
mspec$varmodel$lag.criterion = VAR.opt$lag.criterion
mspec$varmodel$external.regressors = VAR.opt$external.regressors
mspec$varmodel$robust.control = VAR.opt$robust.control
if( !is.null(mspec$varmodel$external.regressors) ){
mspec$varmodel$mxn = dim( as.matrix(VAR.opt$external.regressors) )[2]
} else{
mspec$varmodel$mxn = 0
}
# loop and reconstruct the spec to exclude a mean specification
nl = length(uspec@spec)
tmpspec = vector(mode ="list", length = nl)
for(i in 1:nl){
tmpspec[[i]] = ugarchspec(mean.model = list(include.mean = FALSE, armaOrder = c(0,0),
external.regressors = NULL, arfima = FALSE, garchInMean = FALSE),
variance.model = list(model = uspec@spec[[i]]@variance.model$model,
submodel = uspec@spec[[i]]@variance.model$submodel,
garchOrder = uspec@spec[[i]]@variance.model$garchOrder,
external.regressors = uspec@spec[[i]]@variance.model$vexdata,
variance.targeting = !uspec@spec[[i]]@variance.model$include.omega),
distribution.model = uspec@spec[[i]]@distribution.model$distribution,
start.pars = uspec@spec[[i]]@optimization.model$start.pars,
fixed.pars = uspec@spec[[i]]@optimization.model$fixed.pars)
}
uspec = multispec( tmpspec )
} else{
mspec$varmodel$VAR = FALSE
mspec$varmodel$lag = 1
mspec$varmodel$lag.max = 1
mspec$varmodel$lag.criterion = "HQ"
mspec$varmodel$external.regressors = NULL
mspec$varmodel$mxn = 0
}
ans = new("DCCspec",
mspec = mspec,
uspec = uspec)
return(ans)
}
dccfit.norm = function(spec, data, out.sample = 0, solver = "solnp", solver.control = list(),
fit.control = list(eval.se = TRUE, stationarity = TRUE, scale = TRUE),
parallel = FALSE, parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2),
fit = NULL, VAR.fit = NULL, ...)
{
# we only do 2-stage estimation
# first get the univariate fits
tic = Sys.time()
ufit.control = list()
if(is.null(fit.control$stationarity)){
ufit.control$stationarity = TRUE
} else {
ufit.control$stationarity = fit.control$stationarity
fit.control$stationarity = NULL
}
if(is.null(fit.control$scale)){
ufit.control$scale = TRUE
} else{
ufit.control$scale = fit.control$scale
fit.control$scale = NULL
}
if(is.null(fit.control$eval.se)) fit.control$eval.se = TRUE
##################
# Parallel Execution Prelim Checks (we need them because of the hessian estimation)
if( parallel ){
os = .Platform$OS.type
if(is.null(parallel.control$pkg)){
if( os == "windows" ) parallel.control$pkg = "snowfall" else parallel.control$pkg = "multicore"
if( is.null(parallel.control$cores) ) parallel.control$cores = 2
} else{
mtype = match(tolower(parallel.control$pkg[1]), c("multicore", "snowfall"))
if(is.na(mtype)) stop("\nParallel Package type not recognized in parallel.control\n")
parallel.control$pkg = tolower(parallel.control$pkg[1])
if( os == "windows" && parallel.control$pkg == "multicore" ) stop("\nmulticore not supported on windows O/S\n")
if( is.null(parallel.control$cores) ) parallel.control$cores = 2 else parallel.control$cores = as.integer(parallel.control$cores[1])
}
}
if( is.null( colnames(data) ) ) cnames = paste("Asset_", 1:k, sep = "") else cnames = colnames(data)
xdata = .extractmdata(data)
if(!is.numeric(out.sample)) stop("\ndccfit-->error: out.sample must be numeric\n")
if(as.numeric(out.sample) < 0) stop("\ndccfit-->error: out.sample must be positive\n")
n.start = round(out.sample, 0)
n = dim(xdata$data)[1]
k = dim(xdata$data)[2]
if( (n-n.start) < 100) stop("\ndccfit-->error: function requires at least 100 data\n points to run\n")
data = xdata$data[1:(n-n.start), , drop = FALSE]
dates = xdata$pos[1:(n-n.start), drop = FALSE]
origdata = xdata$data
origdates = xdata$pos
dformat = xdata$dformat
vrmodel = NULL
if( spec@mspec$varmodel$VAR ){
if(spec@mspec$varmodel$mxn>0){
ex = spec@mspec$varmodel$external.regressors
if( dim(ex)[1] < dim(data)[1] ) stop("\ndccfit-->error: external regressor matrix has less points than data!...", call. = FALSE)
orex = ex
ex = ex[1:(n - n.start), , drop = FALSE]
mxn = dim(orex)[2]
} else{
orex = NULL
ex = NULL
mxn = 0
}
if( !is.null(VAR.fit) ){
zdata = VAR.fit$xresiduals
# should be N - p
p = mp = VAR.fit$lag
N = dim(zdata)[1]
K = dim(zdata)[2]
mu = VAR.fit$xfitted
vrmodel = VAR.fit
#spec@mspec$optimization.model$maxOrder = max( mp, spec@mspec$optimization.model$maxdccOrder )
#spec@mspec$optimization.model$maxgarchOrder = mp
} else{
cat("\nestimating VAR model...")
ic = spec@mspec$varmodel$lag.criterion[1]
if( !is.null(spec@mspec$varmodel$lag.max) ){
mp = .varxselect(y = data, lag.max = spec@mspec$varmodel$lag.max, exogen = ex)$selection
mp = as.numeric( mp[which(substr(names(mp), 1, 2) == substr(ic, 1, 2))] )
} else {
mp = spec@mspec$varmodel$lag
}
vrmodel = varxfilter(X = data, p = mp, exogen = ex, Bcoef = NULL, robust = spec@mspec$varmodel$robust,
gamma = spec@mspec$varmodel$robust.control$gamma, delta = spec@mspec$varmodel$robust.control$delta,
nc = spec@mspec$varmodel$robust.control$nc, ns = spec@mspec$varmodel$robust.control$ns)
cat("done!\n")
zdata = vrmodel$xresiduals
# should be N - p
N = dim(zdata)[1]
K = dim(zdata)[2]
mu = vrmodel$xfitted
vrmodel$lag = mp
vrmodel$mxn = mxn
#spec@mspec$optimization.model$maxOrder = max( mp, spec@mspec$optimization.model$maxdccOrder )
#spec@mspec$optimization.model$maxgarchOrder = mp
}
# if we have out of sample data we need to filter
if(out.sample>0){
Bcoef = vrmodel$Bcoef
vrmodel2 = varxfilter(X = origdata, p = mp, exogen = orex, Bcoef = Bcoef)
zorigdata = vrmodel2$xresiduals
} else{
zorigdata = zdata
}
} else{
zdata = data
zorigdata = origdata
orex = NULL
ex = NULL
}
if( is(spec@uspec, "uGARCHmultispec") ){
speclist = spec@uspec
} else{
stop("\ndccfit-->error: DCCspec must contain multispec object with fixed parameters")
}
maxgarchOrder = spec@mspec$optimization.model$maxgarchOrder
tmpidx = .dccparidx( speclist )
paridx = tmpidx$idx
paridxi = tmpidx$idxi
paridxa = tmpidx$idxa
if( !is.null(fit) && is(fit, "uGARCHmultifit") ){
fitlist = fit
} else{
fitlist = multifit(multispec = speclist, data = zorigdata, out.sample = n.start, solver = solver,
solver.control = solver.control, fit.control = ufit.control, parallel = parallel, parallel.control = parallel.control)
assign(".fitlist", fitlist, envir = .GlobalEnv)
}
garchpars = as.numeric( unlist(coef(fitlist) ) )
converge = sapply(fitlist@fit, FUN = function(x) x@fit$convergence)
if( any( converge == 1 ) ){
pr = which(converge != 1)
cat("\nNon-Converged:\n")
print(pr)
cat("\ndccfit-->error: convergence problem in univariate fit...")
cat("\n...returning uGARCHmultifit object instead...check and resubmit...")
return( fitlist )
}
# create a temporary environment to store values (deleted at end of function)
tempenvir = new.env(hash = TRUE, parent = .GlobalEnv)
assign("parallel", parallel, envir = tempenvir)
assign("parallel.control", parallel.control, envir = tempenvir)
assign("eval.se", fit.control$eval.se, envir = tempenvir)
assign("paridx", paridx, envir = tempenvir)
assign("paridxi", paridxi, envir = tempenvir)
assign("paridxa", paridxa, envir = tempenvir)
assign("speclist", speclist, envir = tempenvir)
assign("fitlist", fitlist, envir = tempenvir)
omodel = spec@mspec$optimization.model
sig = sigma(fitlist)
res = residuals(fitlist)
stdresid = res/sig
if( maxgarchOrder > 0 ) stdresid = stdresid[-(1:maxgarchOrder), , drop = FALSE]
Qbar = cov(stdresid)
H = sig^2
assign("omodel", omodel, envir = tempenvir)
assign("mspec", spec@mspec, envir = tempenvir)
assign("stdresid", stdresid, envir = tempenvir)
assign("Qbar", Qbar, envir = tempenvir)
assign("solver", solver, envir = tempenvir)
assign("fit.control", fit.control, envir = tempenvir)
assign("H", H, envir = tempenvir)
Ifn = .dcccon
ILB = 0.0001
IUB = 0.9998
if( solver == "solnp" ) fit.control$stationarity = FALSE else fit.control$stationarity = TRUE
assign("fit.control", fit.control, envir = tempenvir)
# DCC parameters
ipars = .dccstart(data = zdata, garchenv = tempenvir)
pars = ipars$pars
LB = ipars$modelLB
UB = ipars$modelUB
assign("npar", length(pars), envir = tempenvir)
if( any( LB == UB ) )
{
fixid = which(LB == UB)
if( length(fixid) == length(pars) )
{
if( !fit.control$eval.se )
{
# if all parameters are fixed an no standard erros are to
# be calculated then we return a ugarchfilter object
cat("\nrdccfit-->warning: all parameters fixed...returning dccfilter object instead\n")
return(dccfilter(spec = spec, data = zdata, out.sample = out.sample,
parallel = parallel, parallel.control = parallel.control, VAR.fit = VAR.fit))
} else{
# if all parameters are fixed but we require standard errors, we
# skip the solver
use.solver = FALSE
fixed = pars[fixid]
pars = pars[-fixid]
}
} else{
# with some parameters fixed we extract them (to be rejoined at end)
# so that they do not enter the solver
use.solver = TRUE
fixed = pars[fixid]
pars = pars[-fixid]
LB = LB[-fixid]
UB = UB[-fixid]
}
} else{
# no fixed parameters, all normal
fixed = NULL
fixid = NULL
use.solver = TRUE
}
assign("fixed", fixed, envir = tempenvir)
assign("fixid", fixid, envir = tempenvir)
assign(".llh", 1, envir = tempenvir)
# get
if( use.solver )
{
solution = .dccsolver(solver, pars, fun = normal.dccLLH1, Ifn, ILB,
IUB, gr = NULL, hessian = NULL, control = solver.control,
LB = LB, UB = UB, data = zdata, returnType = "llh", garchenv = tempenvir)
sol = solution$sol
hess = solution$hess
opt.pars = sol$par
names(opt.pars) = names(pars)
convergence = sol$convergence
} else{
hess = NULL
timer = Sys.time()-tic
opt.pars = NULL
convergence = 0
sol = list()
sol$message = "all parameters fixed"
}
fit = list()
model = list()
model$vrmodel = vrmodel
model$asset.names = cnames
model$distribution = spec@mspec$distribution
model$include.skew = spec@mspec$include.skew
model$include.shape = spec@mspec$include.shape
model$dccOrder = spec@mspec$dccOrder
model$out.sample = out.sample
model$pos.matrix = spec@mspec$optimization.model$pos.matrix
model$model.matrix = spec@mspec$optimization.model$modelmatrix
fit$spec = spec
# check convergence else write message/return
if( convergence == 0 ){
if( !is.null(fixed) ){
npar = get("npar", tempenvir)
pall = vector(mode = "numeric", length = npar)
pall[fixid] = fixed
pall[-fixid] = opt.pars
opt.pars = pall
Names = omodel$modelnames
names(opt.pars) = Names
assign("fixed", NULL, envir = tempenvir)
assign("fixid", NULL, envir = tempenvir)
}
fit = .dccmakefitmodel(garchmodel = "dccnorm", f = normal.dccLLH2,
pars = opt.pars, garchpars, zdata, garchenv = tempenvir,
timer = 0, message = sol$message, fname = "normal.dccLLH2")
fit$dates = dates
fit$date.format = dformat
fit$origdata = origdata
fit$origdates = origdates
fit$timer = Sys.time() - tic
} else{
fit$message = sol$message
fit$convergence = 1
}
fit$model = model
# make model list to return some usefule information which
ans = new("DCCfit",
mfit = fit,
ufit = fitlist)
rm(tempenvir)
return(ans)
}
dccfit.student = function(spec, data, out.sample = 0, solver = "solnp", solver.control = list(),
fit.control = list(eval.se = TRUE, stationarity = TRUE, scale = FALSE), parallel = FALSE,
parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2), fit = NULL, VAR.fit = NULL, ...)
{
# we only do 2-stage estimation
# first get the univariate fits
tic = Sys.time()
ufit.control = list()
if(is.null(fit.control$stationarity)){
ufit.control$stationarity = TRUE
} else {
ufit.control$stationarity = fit.control$stationarity
fit.control$stationarity = NULL
}
if(is.null(fit.control$scale)){
ufit.control$scale = TRUE
} else{
ufit.control$scale = fit.control$scale
fit.control$scale = NULL
}
if(is.null(fit.control$eval.se)) fit.control$eval.se = TRUE
##################
# Parallel Execution Prelim Checks (we need them because of the hessian estimation)
if( parallel ){
os = .Platform$OS.type
if(is.null(parallel.control$pkg)){
if( os == "windows" ) parallel.control$pkg = "snowfall" else parallel.control$pkg = "multicore"
if( is.null(parallel.control$cores) ) parallel.control$cores = 2
} else{
mtype = match(tolower(parallel.control$pkg[1]), c("multicore", "snowfall"))
if(is.na(mtype)) stop("\nParallel Package type not recognized in parallel.control\n")
parallel.control$pkg = tolower(parallel.control$pkg[1])
if( os == "windows" && parallel.control$pkg == "multicore" ) stop("\nmulticore not supported on windows O/S\n")
if( is.null(parallel.control$cores) ) parallel.control$cores = 2 else parallel.control$cores = as.integer(parallel.control$cores[1])
}
}
if( is.null( colnames(data) ) ) cnames = paste("Asset_", 1:k, sep = "") else cnames = colnames(data)
xdata = .extractmdata(data)
if(!is.numeric(out.sample)) stop("\ndccfit-->error: out.sample must be numeric\n")
if(as.numeric(out.sample) < 0) stop("\ndccfit-->error: out.sample must be positive\n")
n.start = round(out.sample, 0)
n = dim(xdata$data)[1]
k = dim(xdata$data)[2]
if( (n-n.start) < 100) stop("\nrdcc-->error: function requires at least 100 data\n points to run\n")
data = xdata$data[1:(n-n.start), , drop = FALSE]
dates = xdata$pos[1:(n-n.start), drop = FALSE]
origdata = xdata$data
origdates = xdata$pos
dformat = xdata$dformat
vrmodel = NULL
if( spec@mspec$varmodel$VAR ){
if(spec@mspec$varmodel$mxn>0){
ex = spec@mspec$varmodel$external.regressors
if( dim(ex)[1] < dim(data)[1] ) stop("\ndccfit-->error: external regressor matrix has less points than data!...", call. = FALSE)
orex = ex
ex = ex[1:(n - n.start), , drop = FALSE]
mxn = dim(orex)[2]
} else{
orex = NULL
ex = NULL
mxn = 0
}
if( !is.null(VAR.fit) ){
zdata = VAR.fit$xresiduals
# should be N - p
p = mp = VAR.fit$lag
N = dim(zdata)[1]
K = dim(zdata)[2]
mu = VAR.fit$xfitted
vrmodel = VAR.fit
#spec@mspec$optimization.model$maxOrder = max( mp, spec@mspec$optimization.model$maxdccOrder )
#spec@mspec$optimization.model$maxgarchOrder = mp
} else{
cat("\nestimating VAR model...")
ic = spec@mspec$varmodel$lag.criterion[1]
if( !is.null(spec@mspec$varmodel$lag.max) ){
mp = .varxselect(y = data, lag.max = spec@mspec$varmodel$lag.max, exogen = ex)$selection
mp = as.numeric( mp[which(substr(names(mp), 1, 2) == substr(ic, 1, 2))] )
} else {
mp = spec@mspec$varmodel$lag
}
vrmodel = varxfilter(X = data, p = mp, exogen = ex, Bcoef = NULL, robust = spec@mspec$varmodel$robust,
gamma = spec@mspec$varmodel$robust.control$gamma, delta = spec@mspec$varmodel$robust.control$delta,
nc = spec@mspec$varmodel$robust.control$nc, ns = spec@mspec$varmodel$robust.control$ns)
cat("done!\n")
zdata = vrmodel$xresiduals
# should be N - p
N = dim(zdata)[1]
K = dim(zdata)[2]
mu = vrmodel$xfitted
vrmodel$lag = mp
vrmodel$mxn = mxn
#spec@mspec$optimization.model$maxOrder = max( mp, spec@mspec$optimization.model$maxdccOrder )
#spec@mspec$optimization.model$maxgarchOrder = mp
}
# if we have out of sample data we need to filter
if(out.sample>0){
Bcoef = vrmodel$Bcoef
vrmodel2 = varxfilter(X = origdata, p = mp, exogen = orex, Bcoef = Bcoef)
zorigdata = vrmodel2$xresiduals
} else{
zorigdata = zdata
}
} else{
zdata = data
zorigdata = origdata
orex = NULL
ex = NULL
}
if( is(spec@uspec, "uGARCHmultispec") ){
speclist = spec@uspec
} else{
stop("\ndccfit-->error: DCCspec must contain multispec object with fixed parameters")
}
maxgarchOrder = spec@mspec$optimization.model$maxgarchOrder
if( !is.null(fit) && is(fit, "uGARCHmultifit") ){
fitlist = fit
} else{
fitlist = multifit(multispec = speclist, data = zorigdata, out.sample = n.start, solver = solver,
solver.control = solver.control, fit.control = ufit.control, parallel = parallel, parallel.control = parallel.control)
assign(".fitlist", fitlist, envir = .GlobalEnv)
}
garchpars = as.numeric( unlist(coef(fitlist) ) )
tmpidx = .dccparidx( speclist )
paridx = tmpidx$idx
paridxi = tmpidx$idxi
paridxa = tmpidx$idxa
converge = sapply(fitlist@fit, FUN = function(x) x@fit$convergence)
if( any( converge == 1 ) ){
pr = which(converge != 1)
cat("\nNon-Converged:\n")
print(pr)
cat("\ndccfit-->error: convergence problem in univariate fit...")
cat("\n...returning uGARCHmultifit object instead...check and resubmit...")
return( fitlist )
}
# create a temporary environment to store values (deleted at end of function)
tempenvir = new.env(hash = TRUE, parent = .GlobalEnv)
assign("parallel", parallel, envir = tempenvir)
assign("parallel.control", parallel.control, envir = tempenvir)
assign("eval.se", fit.control$eval.se, envir = tempenvir)
assign("paridx", paridx, envir = tempenvir)
assign("paridxi", paridxi, envir = tempenvir)
assign("paridxa", paridxa, envir = tempenvir)
assign("speclist", speclist, envir = tempenvir)
assign("fitlist", fitlist, envir = tempenvir)
omodel = spec@mspec$optimization.model
sig = sigma(fitlist)
res = residuals(fitlist)
stdresid = res/sig
if(maxgarchOrder > 0 ) stdresid = stdresid[-(1:maxgarchOrder), ]
Qbar = cov(stdresid)
H = sig^2
assign("omodel", omodel, envir = tempenvir)
assign("mspec", spec@mspec, envir = tempenvir)
assign("stdresid", stdresid, envir = tempenvir)
assign("Qbar", Qbar, envir = tempenvir)
assign("solver", solver, envir = tempenvir)
assign("fit.control", fit.control, envir = tempenvir)
assign("H", H, envir = tempenvir)
Ifn = .dcccon
ILB = 0.0001
IUB = 0.9998
if(solver == "solnp") fit.control$stationarity = FALSE else fit.control$stationarity = TRUE
assign("fit.control", fit.control, envir = tempenvir)
# DCC parameters
ipars = .dccstart(data = zdata, garchenv = tempenvir)
pars = ipars$pars
LB = ipars$modelLB
UB = ipars$modelUB
assign("npar", length(pars), envir = tempenvir)
if( any( LB == UB ) )
{
fixid = which(LB == UB)
if( length(fixid) == length(pars) )
{
if( !fit.control$eval.se )
{
# if all parameters are fixed an no standard erros are to
# be calculated then we return a ugarchfilter object
cat("\nrdccfit-->warning: all parameters fixed...returning dccfilter object instead\n")
return(dccfilter(spec = spec, data = zdata, out.sample = out.sample,
parallel = parallel, parallel.control = parallel.control, VAR.fit = VAR.fit))
} else{
# if all parameters are fixed but we require standard errors, we
# skip the solver
use.solver = FALSE
fixed = pars[fixid]
pars = pars[-fixid]
}
} else{
# with some parameters fixed we extract them (to be rejoined at end)
# so that they do not enter the solver
use.solver = TRUE
fixed = pars[fixid]
pars = pars[-fixid]
LB = LB[-fixid]
UB = UB[-fixid]
}
} else{
# no fixed parameters, all normal
fixed = NULL
fixid = NULL
use.solver = TRUE
}
assign("fixed", fixed, envir = tempenvir)
assign("fixid", fixid, envir = tempenvir)
assign(".llh", 1, envir = tempenvir)
# get
if( use.solver )
{
solution = .dccsolver(solver, pars, fun = student.dccLLH1, Ifn, ILB,
IUB, gr = NULL, hessian = NULL, control = solver.control,
LB = LB, UB = UB, data = zdata, returnType = "llh", garchenv = tempenvir)
sol = solution$sol
hess = solution$hess
opt.pars = sol$par
names(opt.pars) = names(pars)
convergence = sol$convergence
} else{
hess = NULL
timer = Sys.time()-tic
opt.pars = NULL
convergence = 0
sol = list()
sol$message = "all parameters fixed"
}
fit = list()
model = list()
model$vrmodel = vrmodel
model$asset.names = cnames
model$distribution = spec@mspec$distribution
model$include.skew = spec@mspec$include.skew
model$include.shape = spec@mspec$include.shape
model$dccOrder = spec@mspec$dccOrder
model$out.sample = out.sample
model$pos.matrix = spec@mspec$optimization.model$pos.matrix
model$model.matrix = spec@mspec$optimization.model$modelmatrix
fit$spec = spec
# check convergence else write message/return
if( convergence == 0 ){
if( !is.null(fixed) ){
npar = get("npar", tempenvir)
pall = vector(mode = "numeric", length = npar)
pall[fixid] = fixed
pall[-fixid] = opt.pars
opt.pars = pall
Names = omodel$modelnames
names(opt.pars) = Names
assign("fixed", NULL, envir = tempenvir)
assign("fixid", NULL, envir = tempenvir)
}
fit = .dccmakefitmodel(garchmodel = "dccstudent", f = student.dccLLH2,
pars = opt.pars, garchpars, zdata, garchenv = tempenvir,
timer = 0, message = sol$message, fname = "student.dccLLH2")
fit$dates = dates
fit$date.format = dformat
fit$origdata = origdata
fit$origdates = origdates
fit$timer = Sys.time() - tic
} else{
fit$message = sol$message
fit$convergence = 1
}
fit$model = model
# make model list to return some usefule information which
ans = new("DCCfit",
mfit = fit,
ufit = fitlist)
rm(tempenvir)
return(ans)
}
dccfit.laplace = function(spec, data, out.sample = 0, solver = "solnp", solver.control = list(),
fit.control = list(eval.se = TRUE, stationarity = TRUE, scale = FALSE),
parallel = FALSE, parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2), fit = NULL, VAR.fit = NULL, ...)
{
# we only do 2-stage estimation
# first get the univariate fits
tic = Sys.time()
ufit.control = list()
if(is.null(fit.control$stationarity)){
ufit.control$stationarity = TRUE
} else {
ufit.control$stationarity = fit.control$stationarity
fit.control$stationarity = NULL
}
if(is.null(fit.control$scale)){
ufit.control$scale = TRUE
} else{
ufit.control$scale = fit.control$scale
fit.control$scale = NULL
}
if(is.null(fit.control$eval.se)) fit.control$eval.se = TRUE
##################
# Parallel Execution Prelim Checks (we need them because of the hessian estimation)
if( parallel ){
os = .Platform$OS.type
if(is.null(parallel.control$pkg)){
if( os == "windows" ) parallel.control$pkg = "snowfall" else parallel.control$pkg = "multicore"
if( is.null(parallel.control$cores) ) parallel.control$cores = 2
} else{
mtype = match(tolower(parallel.control$pkg[1]), c("multicore", "snowfall"))
if(is.na(mtype)) stop("\nParallel Package type not recognized in parallel.control\n")
parallel.control$pkg = tolower(parallel.control$pkg[1])
if( os == "windows" && parallel.control$pkg == "multicore" ) stop("\nmulticore not supported on windows O/S\n")
if( is.null(parallel.control$cores) ) parallel.control$cores = 2 else parallel.control$cores = as.integer(parallel.control$cores[1])
}
}
if( is.null( colnames(data) ) ) cnames = paste("Asset_", 1:k, sep = "") else cnames = colnames(data)
xdata = .extractmdata(data)
if(!is.numeric(out.sample)) stop("\ndccfit-->error: out.sample must be numeric\n")
if(as.numeric(out.sample) < 0) stop("\ndccfit-->error: out.sample must be positive\n")
n.start = round(out.sample, 0)
n = dim(xdata$data)[1]
k = dim(xdata$data)[2]
if( (n-n.start) < 100) stop("\nrdcc-->error: function requires at least 100 data\n points to run\n")
data = xdata$data[1:(n-n.start), , drop = FALSE]
dates = xdata$pos[1:(n-n.start), drop = FALSE]
origdata = xdata$data
origdates = xdata$pos
dformat = xdata$dformat
vrmodel = NULL
if( spec@mspec$varmodel$VAR ){
if(spec@mspec$varmodel$mxn>0){
ex = spec@mspec$varmodel$external.regressors
if( dim(ex)[1] < dim(data)[1] ) stop("\ndccfit-->error: external regressor matrix has less points than data!...", call. = FALSE)
orex = ex
ex = ex[1:(n - n.start), , drop = FALSE]
mxn = dim(orex)[2]
} else{
orex = NULL
ex = NULL
mxn = 0
}
if( !is.null(VAR.fit) ){
zdata = VAR.fit$xresiduals
# should be N - p
p = mp = VAR.fit$lag
N = dim(zdata)[1]
K = dim(zdata)[2]
mu = VAR.fit$xfitted
vrmodel = VAR.fit
#spec@mspec$optimization.model$maxOrder = max( mp, spec@mspec$optimization.model$maxdccOrder )
#spec@mspec$optimization.model$maxgarchOrder = mp
} else{
cat("\nestimating VAR model...")
ic = spec@mspec$varmodel$lag.criterion[1]
if( !is.null(spec@mspec$varmodel$lag.max) ){
mp = .varxselect(y = data, lag.max = spec@mspec$varmodel$lag.max, exogen = ex)$selection
mp = as.numeric( mp[which(substr(names(mp), 1, 2) == substr(ic, 1, 2))] )
} else {
mp = spec@mspec$varmodel$lag
}
vrmodel = varxfilter(X = data, p = mp, exogen = ex, Bcoef = NULL, robust = spec@mspec$varmodel$robust,
gamma = spec@mspec$varmodel$robust.control$gamma, delta = spec@mspec$varmodel$robust.control$delta,
nc = spec@mspec$varmodel$robust.control$nc, ns = spec@mspec$varmodel$robust.control$ns)
cat("done!\n")
zdata = vrmodel$xresiduals
# should be N - p
N = dim(zdata)[1]
K = dim(zdata)[2]
mu = vrmodel$xfitted
vrmodel$lag = mp
vrmodel$mxn = mxn
#spec@mspec$optimization.model$maxOrder = max( mp, spec@mspec$optimization.model$maxdccOrder )
#spec@mspec$optimization.model$maxgarchOrder = mp
}
# if we have out of sample data we need to filter
if(out.sample>0){
Bcoef = vrmodel$Bcoef
vrmodel2 = varxfilter(X = origdata, p = mp, exogen = orex, Bcoef = Bcoef)
zorigdata = vrmodel2$xresiduals
} else{
zorigdata = zdata
}
} else{
zdata = data
zorigdata = origdata
orex = NULL
ex = NULL
}
if( is(spec@uspec, "uGARCHmultispec") ){
speclist = spec@uspec
} else{
stop("\ndccfit-->error: DCCspec must contain multispec object with fixed parameters")
}
maxgarchOrder = spec@mspec$optimization.model$maxgarchOrder
tmpidx = .dccparidx( speclist )
paridx = tmpidx$idx
paridxi = tmpidx$idxi
paridxa = tmpidx$idxa
if( !is.null(fit) && is(fit, "uGARCHmultifit") ){
fitlist = fit
} else{
fitlist = multifit(multispec = speclist, data = zorigdata, out.sample = n.start, solver = solver,
solver.control = solver.control, fit.control = ufit.control, parallel = parallel, parallel.control = parallel.control)
assign(".fitlist", fitlist, envir = .GlobalEnv)
}
garchpars = as.numeric( unlist(coef(fitlist) ) )
converge = sapply(fitlist@fit, FUN = function(x) x@fit$convergence)
if( any( converge == 1 ) ){
pr = which(converge != 1)
cat("\nNon-Converged:\n")
print(pr)
cat("\ndccfit-->error: convergence problem in univariate fit...")
cat("\n...returning uGARCHmultifit object instead...check and resubmit...")
return( fitlist )
}
# create a temporary environment to store values (deleted at end of function)
tempenvir = new.env(hash = TRUE, parent = .GlobalEnv)
assign("parallel", parallel, envir = tempenvir)
assign("parallel.control", parallel.control, envir = tempenvir)
assign("eval.se", fit.control$eval.se, envir = tempenvir)
assign("paridx", paridx, envir = tempenvir)
assign("paridxi", paridxi, envir = tempenvir)
assign("paridxa", paridxa, envir = tempenvir)
assign("speclist", speclist, envir = tempenvir)
assign("fitlist", fitlist, envir = tempenvir)
omodel = spec@mspec$optimization.model
sig = sigma(fitlist)
res = residuals(fitlist)
stdresid = res/sig
if( maxgarchOrder > 0 ) stdresid = stdresid[-(1:maxgarchOrder), , drop = FALSE]
Qbar = cov(stdresid)
H = sig^2
assign("omodel", omodel, envir = tempenvir)
assign("mspec", spec@mspec, envir = tempenvir)
assign("stdresid", stdresid, envir = tempenvir)
assign("Qbar", Qbar, envir = tempenvir)
assign("solver", solver, envir = tempenvir)
assign("fit.control", fit.control, envir = tempenvir)
assign("H", H, envir = tempenvir)
Ifn = .dcccon
ILB = 0.0001
IUB = 0.9998
if(solver == "solnp") fit.control$stationarity = FALSE else fit.control$stationarity = TRUE
assign("fit.control", fit.control, envir = tempenvir)
# DCC parameters
ipars = .dccstart(data = zdata, garchenv = tempenvir)
pars = ipars$pars
LB = ipars$modelLB
UB = ipars$modelUB
assign("npar", length(pars), envir = tempenvir)
if( any( LB == UB ) )
{
fixid = which(LB == UB)
if( length(fixid) == length(pars) )
{
if( !fit.control$eval.se )
{
# if all parameters are fixed an no standard erros are to
# be calculated then we return a ugarchfilter object
cat("\nrdccfit-->warning: all parameters fixed...returning dccfilter object instead\n")
return(dccfilter(spec = spec, data = zdata, out.sample = out.sample,
parallel = parallel, parallel.control = parallel.control, VAR.fit = VAR.fit))
} else{
# if all parameters are fixed but we require standard errors, we
# skip the solver
use.solver = FALSE
fixed = pars[fixid]
pars = pars[-fixid]
}
} else{
# with some parameters fixed we extract them (to be rejoined at end)
# so that they do not enter the solver
use.solver = TRUE
fixed = pars[fixid]
pars = pars[-fixid]
LB = LB[-fixid]
UB = UB[-fixid]
}
} else{
# no fixed parameters, all normal
fixed = NULL
fixid = NULL
use.solver = TRUE
}
assign("fixed", fixed, envir = tempenvir)
assign("fixid", fixid, envir = tempenvir)
assign(".llh", 1, envir = tempenvir)
# get
if( use.solver )
{
solution = .dccsolver(solver, pars, fun = laplace.dccLLH1, Ifn, ILB,
IUB, gr = NULL, hessian = NULL, control = solver.control,
LB = LB, UB = UB, data = zdata, returnType = "llh", garchenv = tempenvir)
sol = solution$sol
hess = solution$hess
opt.pars = sol$par
names(opt.pars) = names(pars)
convergence = sol$convergence
# normal.dccLLH2(opt.pars, garchpars, data, returnType = "all", garchenv = tempenvir)$llh
} else{
hess = NULL
timer = Sys.time()-tic
opt.pars = NULL
convergence = 0
sol = list()
sol$message = "all parameters fixed"
}
fit = list()
model = list()
model$vrmodel = vrmodel
model$asset.names = cnames
model$distribution = spec@mspec$distribution
model$include.skew = spec@mspec$include.skew
model$include.shape = spec@mspec$include.shape
model$dccOrder = spec@mspec$dccOrder
model$out.sample = out.sample
model$pos.matrix = spec@mspec$optimization.model$pos.matrix
model$model.matrix = spec@mspec$optimization.model$modelmatrix
fit$spec = spec
# check convergence else write message/return
if( convergence == 0 ){
if( !is.null(fixed) ){
npar = get("npar", tempenvir)
pall = vector(mode = "numeric", length = npar)
pall[fixid] = fixed
pall[-fixid] = opt.pars
opt.pars = pall
Names = omodel$modelnames
names(opt.pars) = Names
assign("fixed", NULL, envir = tempenvir)
assign("fixid", NULL, envir = tempenvir)
}
fit = .dccmakefitmodel(garchmodel = "dcclaplace", f = laplace.dccLLH2,
pars = opt.pars, garchpars, zdata, garchenv = tempenvir,
timer = 0, message = sol$message, fname = "laplace.dccLLH2")
fit$dates = dates
fit$date.format = dformat
fit$origdata = origdata
fit$origdates = origdates
fit$timer = Sys.time() - tic
} else{
fit$message = sol$message
fit$convergence = 1
}
fit$model = model
# make model list to return some usefule information which
ans = new("DCCfit",
mfit = fit,
ufit = fitlist)
rm(tempenvir)
return(ans)
}
dccfilter.norm = function(spec, data, out.sample = 0, filter.control = list(n.old = NULL),
parallel = FALSE, parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2), VAR.fit = NULL, ...)
{
tic = Sys.time()
if( parallel ){
os = .Platform$OS.type
if(is.null(parallel.control$pkg)){
if( os == "windows" ) parallel.control$pkg = "snowfall" else parallel.control$pkg = "multicore"
if( is.null(parallel.control$cores) ) parallel.control$cores = 2
} else{
mtype = match(tolower(parallel.control$pkg[1]), c("multicore", "snowfall"))
if(is.na(mtype)) stop("\nParallel Package type not recognized in parallel.control\n")
parallel.control$pkg = tolower(parallel.control$pkg[1])
if( os == "windows" && parallel.control$pkg == "multicore" ) stop("\nmulticore not supported on windows O/S\n")
if( is.null(parallel.control$cores) ) parallel.control$cores = 2 else parallel.control$cores = as.integer(parallel.control$cores[1])
}
}
n.old = filter.control$n.old
if( is.null( colnames(data) ) ) cnames = paste("Asset_", 1:k, sep = "") else cnames = colnames(data)
xdata = .extractmdata(data)
if(!is.numeric(out.sample)) stop("\ndccfilter-->error: out.sample must be numeric\n")
if(as.numeric(out.sample) < 0) stop("\ndccfilter-->error: out.sample must be positive\n")
n.start = round(out.sample, 0)
n = dim(xdata$data)[1]
k = dim(xdata$data)[2]
if( (n - n.start) < 100) stop("\ndccfilter-->error: function requires at least 100 data\n points to run\n")
data = xdata$data[1:(n - n.start), , drop = FALSE]
dates = xdata$pos[1:(n - n.start)]
origdata = xdata$data
origdates = xdata$pos
dformat = xdata$dformat
vrmodel = NULL
if( spec@mspec$varmodel$VAR ){
if(spec@mspec$varmodel$mxn>0){
ex = spec@mspec$varmodel$external.regressors
if( dim(ex)[1] < dim(data)[1] ) stop("\ndccfilter-->error: external regressor matrix has less points than data!...", call. = FALSE)
orex = ex
ex = ex[1:(n - n.start), , drop = FALSE]
mxn = dim(orex)[2]
} else{
orex = NULL
ex = NULL
mxn = 0
}
if( !is.null(VAR.fit) ){
zdata = VAR.fit$xresiduals
# should be N - p
p = mp = VAR.fit$lag
N = dim(zdata)[1]
K = dim(zdata)[2]
mu = VAR.fit$xfitted
vrmodel = VAR.fit
if(out.sample>0){
Bcoef = VAR.fit$Bcoef
vrmodel2 = varxfilter(X = origdata, p = mp, exogen = orex, Bcoef = Bcoef)
zorigdata = vrmodel2$xresiduals
} else{
zorigdata = zdata
}
} else{
cat("\nestimating VAR model...")
ic = spec@mspec$varmodel$lag.criterion[1]
if( !is.null(spec@mspec$varmodel$lag.max) ){
mp = .varxselect(y = data, lag.max = spec@mspec$varmodel$lag.max, exogen = ex)$selection
mp = as.numeric( mp[which(substr(names(mp), 1, 2) == substr(ic, 1, 2))] )
} else {
mp = spec@mspec$varmodel$lag
}
vrmodel = varxfilter(X = data, p = mp, exogen = ex, Bcoef = NULL, robust = spec@mspec$varmodel$robust,
gamma = spec@mspec$varmodel$robust.control$gamma, delta = spec@mspec$varmodel$robust.control$delta,
nc = spec@mspec$varmodel$robust.control$nc, ns = spec@mspec$varmodel$robust.control$ns)
cat("done!\n")
zdata = vrmodel$xresiduals
# should be N - p
N = dim(zdata)[1]
K = dim(zdata)[2]
mu = vrmodel$xfitted
vrmodel$lag = mp
vrmodel$mxn = mxn
# if we have out of sample data we need to filter
if(out.sample>0){
Bcoef = vrmodel$Bcoef
vrmodel2 = varxfilter(X = origdata, p = mp, exogen = orex, Bcoef = Bcoef)
zorigdata = vrmodel2$xresiduals
} else{
zorigdata = zdata
}
}
} else{
zdata = data
zorigdata = origdata
orex = NULL
ex = NULL
}
if( is(spec@uspec, "uGARCHmultispec") ){
speclist = spec@uspec
} else{
stop("\ndccfilter-->error: DCCspec must contain multispec object with fixed parameters")
}
maxgarchOrder = spec@mspec$optimization.model$maxgarchOrder
filterlist = multifilter(multifitORspec = spec@uspec, data = as.data.frame(zorigdata), out.sample = out.sample,
parallel = parallel, parallel.control = parallel.control, n.old = n.old)
tmpidx = .dccparidx( spec@uspec )
paridx = tmpidx$idx
paridxi = tmpidx$idxi
paridxa = tmpidx$idxa
garchpars = as.numeric( unlist(coef(filterlist) ) )
tempenvir = new.env(hash = TRUE, parent = .GlobalEnv)
fit.control = list(stationarity = TRUE, eval.se = FALSE)
assign("fixed", NULL, envir = tempenvir)
assign("fixid", NULL, envir = tempenvir)
assign("parallel", parallel, envir = tempenvir)
assign("parallel.control", parallel.control, envir = tempenvir)
assign("eval.se", FALSE, envir = tempenvir)
assign("mspec", spec@mspec, envir = tempenvir)
assign("paridx", paridx, envir = tempenvir)
assign("paridxi", paridxi, envir = tempenvir)
assign("paridxa", paridxa, envir = tempenvir)
assign("speclist", speclist, envir = tempenvir)
assign("filterlist", filterlist, envir = tempenvir)
assign("fit.control", fit.control, envir = tempenvir)
omodel = spec@mspec$optimization.model
sig = sigma(filterlist)
res = residuals(filterlist)
stdresid = res/sig
if(maxgarchOrder > 0 ) stdresid = stdresid[-(1:maxgarchOrder), ]
Qbar = cov(stdresid)
H = sig^2
assign("omodel", omodel, envir = tempenvir)
assign("mspec", spec@mspec, envir = tempenvir)
assign("stdresid", stdresid, envir = tempenvir)
assign("Qbar", Qbar, envir = tempenvir)
assign("H", H, envir = tempenvir)
dccpars = as.numeric( unlist(spec@mspec$optimization.model$fixed.pars) )
assign("npar", length(dccpars), envir = tempenvir)
filt = .dccmakefiltermodel(garchmodel = "dccnorm", f = normal.dccLLH2,
pars = dccpars, garchpars, zdata, garchenv = tempenvir,
timer = 0, message = 0, fname = "normal.dccLLH2")
model = list()
model$vrmodel = vrmodel
model$asset.names = cnames
model$distribution = spec@mspec$distribution
model$include.skew = spec@mspec$include.skew
model$include.shape = spec@mspec$include.shape
model$dccOrder = spec@mspec$dccOrder
model$out.sample = out.sample
filt$dates = dates
filt$date.format = dformat
filt$origdata = origdata
filt$origdates = origdates
filt$model = model
filt$timer = Sys.time() - tic
ans = new("DCCfilter",
mfilter = filt,
ufilter = filterlist)
rm(tempenvir)
return(ans)
}
dccfilter.student = function(spec, data, out.sample = 0, filter.control = list(n.old = NULL),
parallel = FALSE, parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2), VAR.fit = NULL, ...)
{
tic = Sys.time()
if( parallel ){
os = .Platform$OS.type
if(is.null(parallel.control$pkg)){
if( os == "windows" ) parallel.control$pkg = "snowfall" else parallel.control$pkg = "multicore"
if( is.null(parallel.control$cores) ) parallel.control$cores = 2
} else{
mtype = match(tolower(parallel.control$pkg[1]), c("multicore", "snowfall"))
if(is.na(mtype)) stop("\nParallel Package type not recognized in parallel.control\n")
parallel.control$pkg = tolower(parallel.control$pkg[1])
if( os == "windows" && parallel.control$pkg == "multicore" ) stop("\nmulticore not supported on windows O/S\n")
if( is.null(parallel.control$cores) ) parallel.control$cores = 2 else parallel.control$cores = as.integer(parallel.control$cores[1])
}
}
n.old = filter.control$n.old
if( is.null( colnames(data) ) ) cnames = paste("Asset_", 1:k, sep = "") else cnames = colnames(data)
xdata = .extractmdata(data)
if(!is.numeric(out.sample)) stop("\ndccfilter-->error: out.sample must be numeric\n")
if(as.numeric(out.sample) < 0) stop("\ndccfilter-->error: out.sample must be positive\n")
n.start = round(out.sample, 0)
n = dim(xdata$data)[1]
k = dim(xdata$data)[2]
if( (n - n.start) < 100) stop("\ndccfilter-->error: function requires at least 100 data\n points to run\n")
data = xdata$data[1:(n - n.start), , drop = FALSE]
dates = xdata$pos[1:(n - n.start)]
origdata = xdata$data
origdates = xdata$pos
dformat = xdata$dformat
vrmodel = NULL
if( spec@mspec$varmodel$VAR ){
if(spec@mspec$varmodel$mxn>0){
ex = spec@mspec$varmodel$external.regressors
if( dim(ex)[1] < dim(data)[1] ) stop("\ndccfilter-->error: external regressor matrix has less points than data!...", call. = FALSE)
orex = ex
ex = ex[1:(n - n.start), , drop = FALSE]
mxn = dim(orex)[2]
} else{
orex = NULL
ex = NULL
mxn = 0
}
if( !is.null(VAR.fit) ){
zdata = VAR.fit$xresiduals
# should be N - p
p = mp = VAR.fit$lag
N = dim(zdata)[1]
K = dim(zdata)[2]
mu = VAR.fit$xfitted
vrmodel = VAR.fit
if(out.sample>0){
Bcoef = VAR.fit$Bcoef
vrmodel2 = varxfilter(X = origdata, p = mp, exogen = orex, Bcoef = Bcoef)
zorigdata = vrmodel2$xresiduals
} else{
zorigdata = zdata
}
} else{
cat("\nestimating VAR model...")
ic = spec@mspec$varmodel$lag.criterion[1]
if( !is.null(spec@mspec$varmodel$lag.max) ){
mp = .varxselect(y = data, lag.max = spec@mspec$varmodel$lag.max, exogen = ex)$selection
mp = as.numeric( mp[which(substr(names(mp), 1, 2) == substr(ic, 1, 2))] )
} else {
mp = spec@mspec$varmodel$lag
}
vrmodel = varxfilter(X = data, p = mp, exogen = ex, Bcoef = NULL, robust = spec@mspec$varmodel$robust,
gamma = spec@mspec$varmodel$robust.control$gamma, delta = spec@mspec$varmodel$robust.control$delta,
nc = spec@mspec$varmodel$robust.control$nc, ns = spec@mspec$varmodel$robust.control$ns)
cat("done!\n")
zdata = vrmodel$xresiduals
# should be N - p
N = dim(zdata)[1]
K = dim(zdata)[2]
mu = vrmodel$xfitted
vrmodel$lag = mp
vrmodel$mxn = mxn
# if we have out of sample data we need to filter
if(out.sample>0){
Bcoef = vrmodel$Bcoef
vrmodel2 = varxfilter(X = origdata, p = mp, exogen = orex, Bcoef = Bcoef)
zorigdata = vrmodel2$xresiduals
} else{
zorigdata = zdata
}
}
} else{
zdata = data
zorigdata = origdata
orex = NULL
ex = NULL
}
if( is(spec@uspec, "uGARCHmultispec") ){
speclist = spec@uspec
} else{
stop("\ndccfilter-->error: DCCspec must contain multispec object with fixed parameters")
}
maxgarchOrder = spec@mspec$optimization.model$maxgarchOrder
filterlist = multifilter(multifitORspec = spec@uspec, data = as.data.frame(zorigdata), out.sample = out.sample,
parallel = parallel, parallel.control = parallel.control, n.old = n.old)
tmpidx = .dccparidx( spec@uspec )
paridx = tmpidx$idx
paridxi = tmpidx$idxi
paridxa = tmpidx$idxa
garchpars = as.numeric( unlist(coef(filterlist) ) )
tempenvir = new.env(hash = TRUE, parent = .GlobalEnv)
fit.control = list(stationarity = TRUE, eval.se = FALSE)
assign("fixed", NULL, envir = tempenvir)
assign("fixid", NULL, envir = tempenvir)
assign("parallel", parallel, envir = tempenvir)
assign("parallel.control", parallel.control, envir = tempenvir)
assign("eval.se", FALSE, envir = tempenvir)
assign("mspec", spec@mspec, envir = tempenvir)
assign("paridx", paridx, envir = tempenvir)
assign("paridxi", paridxi, envir = tempenvir)
assign("paridxa", paridxa, envir = tempenvir)
assign("speclist", speclist, envir = tempenvir)
assign("filterlist", filterlist, envir = tempenvir)
assign("fit.control", fit.control, envir = tempenvir)
omodel = spec@mspec$optimization.model
sig = sigma(filterlist)
res = residuals(filterlist)
stdresid = res/sig
if(maxgarchOrder > 0 ) stdresid = stdresid[-(1:maxgarchOrder), ]
Qbar = cov(stdresid)
H = sig^2
assign("omodel", omodel, envir = tempenvir)
assign("mspec", spec@mspec, envir = tempenvir)
assign("stdresid", stdresid, envir = tempenvir)
assign("Qbar", Qbar, envir = tempenvir)
assign("H", H, envir = tempenvir)
dccpars = as.numeric( unlist(spec@mspec$optimization.model$fixed.pars) )
assign("npar", length(dccpars), envir = tempenvir)
filt = .dccmakefiltermodel(garchmodel = "dccstudent", f = student.dccLLH2,
pars = dccpars, garchpars, zdata, garchenv = tempenvir,
timer = 0, message = 0, fname = "student.dccLLH2")
model = list()
model$vrmodel = vrmodel
model$asset.names = cnames
model$distribution = spec@mspec$distribution
model$include.skew = spec@mspec$include.skew
model$include.shape = spec@mspec$include.shape
model$dccOrder = spec@mspec$dccOrder
model$out.sample = out.sample
filt$dates = dates
filt$date.format = dformat
filt$origdata = origdata
filt$origdates = origdates
filt$model = model
filt$timer = Sys.time() - tic
ans = new("DCCfilter",
mfilter = filt,
ufilter = filterlist)
rm(tempenvir)
return(ans)
}
dccfilter.laplace = function(spec, data, out.sample = 0, filter.control = list(n.old = NULL),
parallel = FALSE, parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2), VAR.fit = NULL, ...)
{
tic = Sys.time()
if( parallel ){
os = .Platform$OS.type
if(is.null(parallel.control$pkg)){
if( os == "windows" ) parallel.control$pkg = "snowfall" else parallel.control$pkg = "multicore"
if( is.null(parallel.control$cores) ) parallel.control$cores = 2
} else{
mtype = match(tolower(parallel.control$pkg[1]), c("multicore", "snowfall"))
if(is.na(mtype)) stop("\nParallel Package type not recognized in parallel.control\n")
parallel.control$pkg = tolower(parallel.control$pkg[1])
if( os == "windows" && parallel.control$pkg == "multicore" ) stop("\nmulticore not supported on windows O/S\n")
if( is.null(parallel.control$cores) ) parallel.control$cores = 2 else parallel.control$cores = as.integer(parallel.control$cores[1])
}
}
n.old = filter.control$n.old
if( is.null( colnames(data) ) ) cnames = paste("Asset_", 1:k, sep = "") else cnames = colnames(data)
xdata = .extractmdata(data)
if(!is.numeric(out.sample)) stop("\ndccfilter-->error: out.sample must be numeric\n")
if(as.numeric(out.sample) < 0) stop("\ndccfilter-->error: out.sample must be positive\n")
n.start = round(out.sample, 0)
n = dim(xdata$data)[1]
k = dim(xdata$data)[2]
if( (n - n.start) < 100) stop("\ndccfilter-->error: function requires at least 100 data\n points to run\n")
data = xdata$data[1:(n - n.start), , drop = FALSE]
dates = xdata$pos[1:(n - n.start)]
origdata = xdata$data
origdates = xdata$pos
dformat = xdata$dformat
vrmodel = NULL
if( spec@mspec$varmodel$VAR ){
if(spec@mspec$varmodel$mxn>0){
ex = spec@mspec$varmodel$external.regressors
if( dim(ex)[1] < dim(data)[1] ) stop("\ndccfilter-->error: external regressor matrix has less points than data!...", call. = FALSE)
orex = ex
ex = ex[1:(n - n.start), , drop = FALSE]
mxn = dim(orex)[2]
} else{
orex = NULL
ex = NULL
mxn = 0
}
if( !is.null(VAR.fit) ){
zdata = VAR.fit$xresiduals
# should be N - p
p = mp = VAR.fit$lag
N = dim(zdata)[1]
K = dim(zdata)[2]
mu = VAR.fit$xfitted
vrmodel = VAR.fit
if(out.sample>0){
Bcoef = VAR.fit$Bcoef
vrmodel2 = varxfilter(X = origdata, p = mp, exogen = orex, Bcoef = Bcoef)
zorigdata = vrmodel2$xresiduals
} else{
zorigdata = zdata
}
} else{
cat("\nestimating VAR model...")
ic = spec@mspec$varmodel$lag.criterion[1]
if( !is.null(spec@mspec$varmodel$lag.max) ){
mp = .varxselect(y = data, lag.max = spec@mspec$varmodel$lag.max, exogen = ex)$selection
mp = as.numeric( mp[which(substr(names(mp), 1, 2) == substr(ic, 1, 2))] )
} else {
mp = spec@mspec$varmodel$lag
}
vrmodel = varxfilter(X = data, p = mp, exogen = ex, Bcoef = NULL, robust = spec@mspec$varmodel$robust,
gamma = spec@mspec$varmodel$robust.control$gamma, delta = spec@mspec$varmodel$robust.control$delta,
nc = spec@mspec$varmodel$robust.control$nc, ns = spec@mspec$varmodel$robust.control$ns)
cat("done!\n")
zdata = vrmodel$xresiduals
# should be N - p
N = dim(zdata)[1]
K = dim(zdata)[2]
mu = vrmodel$xfitted
vrmodel$lag = mp
vrmodel$mxn = mxn
# if we have out of sample data we need to filter
if(out.sample>0){
Bcoef = vrmodel$Bcoef
vrmodel2 = varxfilter(X = origdata, p = mp, exogen = orex, Bcoef = Bcoef)
zorigdata = vrmodel2$xresiduals
} else{
zorigdata = zdata
}
}
} else{
zdata = data
zorigdata = origdata
orex = NULL
ex = NULL
}
if( is(spec@uspec, "uGARCHmultispec") ){
speclist = spec@uspec
} else{
stop("\ndccfilter-->error: DCCspec must contain multispec object with fixed parameters")
}
maxgarchOrder = spec@mspec$optimization.model$maxgarchOrder
filterlist = multifilter(multifitORspec = spec@uspec, data = as.data.frame(zorigdata), out.sample = out.sample,
parallel = parallel, parallel.control = parallel.control, n.old = n.old)
tmpidx = .dccparidx( spec@uspec )
paridx = tmpidx$idx
paridxi = tmpidx$idxi
paridxa = tmpidx$idxa
garchpars = as.numeric( unlist(coef(filterlist) ) )
tempenvir = new.env(hash = TRUE, parent = .GlobalEnv)
fit.control = list(stationarity = TRUE, eval.se = FALSE)
assign("fixed", NULL, envir = tempenvir)
assign("fixid", NULL, envir = tempenvir)
assign("parallel", parallel, envir = tempenvir)
assign("parallel.control", parallel.control, envir = tempenvir)
assign("eval.se", FALSE, envir = tempenvir)
assign("mspec", spec@mspec, envir = tempenvir)
assign("paridx", paridx, envir = tempenvir)
assign("paridxi", paridxi, envir = tempenvir)
assign("paridxa", paridxa, envir = tempenvir)
assign("speclist", speclist, envir = tempenvir)
assign("filterlist", filterlist, envir = tempenvir)
assign("fit.control", fit.control, envir = tempenvir)
omodel = spec@mspec$optimization.model
sig = sigma(filterlist)
res = residuals(filterlist)
stdresid = res/sig
if(maxgarchOrder > 0 ) stdresid = stdresid[-(1:maxgarchOrder), ]
Qbar = cov(stdresid)
H = sig^2
assign("omodel", omodel, envir = tempenvir)
assign("mspec", spec@mspec, envir = tempenvir)
assign("stdresid", stdresid, envir = tempenvir)
assign("Qbar", Qbar, envir = tempenvir)
assign("H", H, envir = tempenvir)
dccpars = as.numeric( unlist(spec@mspec$optimization.model$fixed.pars) )
assign("npar", length(dccpars), envir = tempenvir)
filt = .dccmakefiltermodel(garchmodel = "dcclaplace", f = laplace.dccLLH2,
pars = dccpars, garchpars, zdata, garchenv = tempenvir,
timer = 0, message = 0, fname = "laplace.dccLLH2")
model = list()
model$vrmodel = vrmodel
model$asset.names = cnames
model$distribution = spec@mspec$distribution
model$include.skew = spec@mspec$include.skew
model$include.shape = spec@mspec$include.shape
model$dccOrder = spec@mspec$dccOrder
model$out.sample = out.sample
filt$dates = dates
filt$date.format = dformat
filt$origdata = origdata
filt$origdates = origdates
filt$model = model
filt$timer = Sys.time() - tic
ans = new("DCCfilter",
mfilter = filt,
ufilter = filterlist)
rm(tempenvir)
return(ans)
}
.dccforecast = function(fit, n.ahead = 1, n.roll = 0, external.forecasts = list(mregfor = NULL, vregfor = NULL),
parallel = FALSE, parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2), ...)
{
# checks for out.sample wrt n.roll
ns = fit@mfit$model$out.sample
if( n.roll > ns ) stop("n.roll must not be greater than out.sample!")
if(n.roll>1 && n.ahead>1) stop("\ngogarchforecast-->error: n.ahead must be equal to 1 when using n.roll\n")
# checks for external forecasts
tf = n.ahead + n.roll
if( !is.null( external.forecasts$mregfor ) ){
mregfor = external.forecasts$mregfor
if( !is.matrix(mregfor) ) stop("\nmregfor must be a matrix.")
if( dim(mregfor)[1] < tf ) stop("\nmregfor must have at least n.ahead + n.roll observations to be used")
mregfor = mregfor[1:tf, , drop = FALSE]
} else{
mregfor = NULL
}
if( !is.null( external.forecasts$vregfor ) ){
if( !is.matrix(vregfor) ) stop("\nvregfor must be a matrix.")
if( dim(vregfor)[1] < tf ) stop("\nvregfor must have at least n.ahead + n.roll observations to be used")
vregfor = vregfor[1:tf, , drop = FALSE]
}
vrmodel = fit@mfit$model$vrmodel
if( !is.null(vrmodel) ){
mxn = vrmodel$mxn
if( mxn > 0 ){
if( is.null(external.forecasts$mregfor ) ){
warning("\nExternal Regressor Forecasts Matrix NULL...setting to zero...\n")
mregfor = matrix(0, ncol = mxn, nrow = (n.roll + n.ahead) )
} else{
mxn = vrmodel$mxn
if( dim(mregfor)[2] != mxn ) stop("\ndccforecast-->error: wrong number of external regressors!...", call. = FALSE)
if( dim(mregfor)[1] < (n.roll + n.ahead) ) stop("\ndccforecast-->error: external regressor matrix has less points than requested forecast length (1+n.roll) x n.ahead!...", call. = FALSE)
}
} else{
mregfor = NULL
}
if(n.roll>1 && n.ahead>1) stop("\ndccforecast-->error: n.ahead must be equal to 1 when using n.roll\n")
if( n.ahead == 1 && (n.roll > ns) ) stop("\ndccforecast-->error: n.roll greater than out.sample!", call. = FALSE)
#VARf = .varpredict(fit, n.ahead, n.roll, mregfor)$Mu
VARf = varxforecast(X = fit@mfit$origdata, Bcoef = vrmodel$Bcoef, p = vrmodel$lag, out.sample = ns,
n.ahead, n.roll, mregfor)
} else{
VARf = NULL
}
exf = external.forecasts
if( !is.null(fit@mfit$model$vrmodel) ){
exf$mregfor = NULL
}
ans = .dccforecastm(fit, n.ahead = n.ahead, n.roll = n.roll, external.forecasts = exf,
parallel = parallel, parallel.control = parallel.control, ...)
model = list()
model$n.roll = n.roll
model$n.ahead = n.ahead
model$asset.names = fit@mfit$model$asset.names
model$distribution = fit@mfit$model$distribution
model$include.skew = fit@mfit$model$include.skew
model$include.shape = fit@mfit$model$include.shape
model$dccOrder = fit@mfit$model$dccOrder
model$out.sample = ns
model$dccpars = coef(fit, type = "dcc")
model$garchpars = coef(fit, type = "garch")
model$origdata = fit@mfit$origdata
model$origdates = fit@mfit$origdates
# keep last 50 points of R and H for plotting
model$R = last( rcor(fit), min(length(fit@mfit$R), 50) )
model$H = last( rcov(fit), min(length(fit@mfit$R), 50) )
model$fitted = tail(fitted(fit), min(length(fit@mfit$R), 50) )
model$sigma = tail(sigma(fit), min(length(fit@mfit$R), 50) )
mforecast = list( H = ans$H, R = ans$R, Q = ans$Q, Rbar = ans$Rbar, VARf = VARf )
mforecast$model = model
uforecast = ans$ufor
mforecast$shape = mforecast$skew = NULL
if( mforecast$model$include.shape ) mforecast$shape = rshape(fit)
if( mforecast$model$include.skew ) mforecast$skew = rskew(fit)
ans = new("DCCforecast",
mforecast = mforecast,
uforecast = uforecast)
return( ans )
}
dccsim1.norm = function(fitORspec, n.sim = 1000, n.start = 0, m.sim = 1, startMethod = c("unconditional", "sample"),
presigma = NULL, preresiduals = NULL, prereturns = NULL, preR = NULL, preH = NULL, rseed = NULL, mexsimdata = NULL,
vexsimdata = NULL, parallel = FALSE, parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2), VAR.fit = NULL, ...)
{
fit = fitORspec
Data = head(fit@mfit$origdata, dim(fit@mfit$origdata)[1] - fit@mfit$model$out.sample)
m = dim(Data)[2]
n = dim(Data)[1]
dccOrder = fit@mfit$model$dccOrder
mo = max( dccOrder )
if( is.null(rseed) ){
rseed = as.integer(runif(1, 1, Sys.time()))
} else {
if(length(rseed) == 1) rseed = as.integer(rseed[1]) else rseed = as.integer( rseed[1:m.sim] )
}
# we use the GED distribution with nu = 1 which corresponds to the Laplace case.
if(length(rseed) == 1){
tmp = rnorm(m * (n.sim + n.start + mo) * m.sim, 0, 1)
z = array(tmp, dim = c(n.sim + n.start + mo, m, m.sim))
} else{
z = array(NA, dim = c(n.sim + n.start + mo, m, m.sim))
for(i in 1:m.sim){
set.seed( rseed[i] )
z[,,i] = matrix(rnorm(m * (n.sim + n.start + mo), 0, 1), nrow = n.sim + n.start + mo, ncol = m)
}
}
# ok now to expand rseed
if(length(rseed) == 1){
rseed = c(rseed, as.integer(runif(m.sim, 1, Sys.time())))
}
if( is.null(preR) ){
if( startMethod[1] == "sample" ){
preR = last(rcor(fit))[,,1]
diag(preR) = 1
} else{
preR = cor(as.matrix(Data))
}
} else{
chk = dcc.symcheck(preR, m, d = rep(1, m))
}
preQ = preR
uncv = sapply(fit@ufit@fit, FUN = function(x) uncvariance(x))
if( is.null(preH) ){
preH = diag(sqrt(uncv)) %*% preR %*% diag(sqrt(uncv))
} else{
chk = dcc.symcheck(preH, m, d = NULL)
}
if( !is.null(presigma) ){
if( !is.matrix(presigma) ) stop("\ndccsim-->error: presigma must be a matrix.")
if( dim(presigma)[2] != m ) stop("\ndccsim-->error: wrong column dimension for presigma.")
}
if( !is.null(preresiduals) ){
if( !is.matrix(preresiduals) ) stop("\ndccsim-->error: preresiduals must be a matrix.")
if( dim(preresiduals)[2] != m ) stop("\ndccsim-->error: wrong column dimension for preresiduals.")
}
if( !is.null(prereturns) ){
if( !is.matrix(prereturns) ) stop("\ndccsim-->error: prereturns must be a matrix.")
if( dim(prereturns)[2] != m ) stop("\ndccsim-->error: wrong column dimension for prereturns.")
}
garchpars = coef(fit, type = "garch")
dccpars = coef(fit, type = "dcc")
dcca = dccpars[1:dccOrder[1]]
dccb = dccpars[(dccOrder[1]+1):sum(dccOrder[1:2])]
simRes = simX = simR = simQ = simH = simSeries = vector(mode = "list", length = m.sim)
if( parallel ){
os = .Platform$OS.type
if(is.null(parallel.control$pkg)){
if( os == "windows" ) parallel.control$pkg = "snowfall" else parallel.control$pkg = "multicore"
if( is.null(parallel.control$cores) ) parallel.control$cores = 2
} else{
mtype = match(tolower(parallel.control$pkg[1]), c("multicore", "snowfall"))
if(is.na(mtype)) stop("\nParallel Package type not recognized in parallel.control\n")
parallel.control$pkg = tolower(parallel.control$pkg[1])
if( os == "windows" && parallel.control$pkg == "multicore" ) stop("\nmulticore not supported on windows O/S\n")
if( is.null(parallel.control$cores) ) parallel.control$cores = 2 else parallel.control$cores = as.integer(parallel.control$cores[1])
}
if( parallel.control$pkg == "multicore" ){
simH = vector(mode = "list", length = m.sim)
simX = vector(mode = "list", length = m.sim)
mtmp = mclapply(as.list(1:m.sim), FUN = function(j)
.dccsimf(Z = z[,,j], Qbar = preQ, Rbar = preR, dcca, dccb, mo, n.sim, n.start, m, rseed[j]),
mc.cores = parallel.control$cores)
# need to pass m.sim and matrix of simulated z's to ugarchsim (speedup)
simR = lapply(mtmp, FUN = function(x) if(is.matrix(x$R)) array(x$R, dim = c(m, m, n.sim)) else last(x$R, n.sim))
simQ = lapply(mtmp, FUN = function(x) if(is.matrix(x$Q)) array(x$Q, dim = c(m, m, n.sim)) else last(x$Q, n.sim))
stdresid = vector(mode = "list", length = m)
for(i in 1:m) stdresid[[i]] = sapply(mtmp, FUN = function(x) x$stdresid[,i])
xtmp = mclapply(as.list(1:m), FUN = function(j){
maxx = fit@ufit@fit[[j]]@model$maxOrder2;
htmp = ugarchsim(fit@ufit@fit[[j]], n.sim = n.sim + n.start, n.start = 0, m.sim = m.sim,
startMethod = startMethod[1], custom.dist = list(name = "sample",
distfit = matrix(stdresid[[j]][-(1:mo), ], ncol = m.sim)),
presigma = if( is.null(presigma) ) NA else tail(presigma[,j], maxx),
preresiduals = if( is.null(preresiduals) ) NA else tail(preresiduals[,j], maxx),
prereturns = if( is.null(prereturns) | !is.null(fit@mfit$model$vrmodel) ) NA else tail(prereturns[,j], maxx),
mexsimdata = if( is.null(fit@mfit$model$vrmodel) ) mexsimdata[[j]] else NULL, vexsimdata = vexsimdata[[j]] )
h = matrix(tail(htmp@simulation$sigmaSim^2, n.sim), nrow = n.sim)
x = matrix(htmp@simulation$seriesSim, nrow = n.sim + n.start)
return(list(h = h, x = x))}, mc.cores = parallel.control$cores)
H = array(NA, dim = c(n.sim, m, m.sim))
tmpH = array(NA, dim = c(m, m, n.sim))
for(i in 1:n.sim) H[i,,] = t(sapply(xtmp, FUN = function(x) as.numeric(x$h[i,])))
for(i in 1:m.sim){
for(j in 1:n.sim){
tmpH[ , , j] = diag(sqrt( H[j, , i]) ) %*% simR[[i]][ , , j] %*% diag(sqrt( H[j, , i] ) )
}
simH[[i]] = tmpH
}
simxX = array(NA, dim = c(n.sim, m, m.sim))
for(i in 1:n.sim) simxX[i,,] = t(sapply(xtmp, FUN = function(x) as.numeric(x$x[i,])))
simX = vector(mode = "list", length = m.sim)
for(i in 1:m.sim) simX[[i]] = matrix(simxX[,,i], nrow = n.sim)
} else{
simH = vector(mode = "list", length = m.sim)
simX = vector(mode = "list", length = m.sim)
sfInit(parallel = TRUE, cpus = parallel.control$cores)
sfExport("z", "preQ", "preR", "dcca", "dccb", "mo", "n.sim", "n.start", "m", "rseed", local = TRUE)
mtmp = sfLapply(as.list(1:m.sim), fun = function(j)
rgarch:::.dccsimf(Z = z[,,j], Qbar = preQ, Rbar = preR, dcca, dccb, mo, n.sim, n.start, m, rseed[j]))
# need to pass m.sim and matrix of simulated z's to ugarchsim (speedup)
simR = lapply(mtmp, FUN = function(x) if(is.matrix(x$R)) array(x$R, dim = c(m, m, n.sim)) else last(x$R, n.sim))
simQ = lapply(mtmp, FUN = function(x) if(is.matrix(x$Q)) array(x$Q, dim = c(m, m, n.sim)) else last(x$Q, n.sim))
stdresid = vector(mode = "list", length = m)
for(i in 1:m) stdresid[[i]] = sapply(mtmp, FUN = function(x) x$stdresid[,i])
sfExport("fit", "n.sim", "n.start", "m.sim", "startMethod", "stdresid", "presigma",
"preresiduals", "prereturns", "mexsimdata", "vexsimdata", local = TRUE)
xtmp = sfLapply(as.list(1:m), fun = function(j){
maxx = fit@ufit@fit[[j]]@model$maxOrder2;
htmp = rgarch::ugarchsim(fit@ufit@fit[[j]], n.sim = n.sim + n.start, n.start = 0, m.sim = m.sim,
startMethod = startMethod[1], custom.dist = list(name = "sample",
distfit = matrix(stdresid[[j]][-(1:mo), ], ncol = m.sim)),
presigma = if( is.null(presigma) ) NA else tail(presigma[,j], maxx),
preresiduals = if( is.null(preresiduals) ) NA else tail(preresiduals[,j], maxx),
prereturns = if( is.null(prereturns) | !is.null(fit@mfit$model$vrmodel) ) NA else tail(prereturns[,j], maxx),
mexsimdata = if( is.null(fit@mfit$model$vrmodel) ) mexsimdata[[j]] else NULL, vexsimdata = vexsimdata[[j]] )
h = matrix(tail(htmp@simulation$sigmaSim^2, n.sim), nrow = n.sim)
x = matrix(htmp@simulation$seriesSim, nrow = n.sim + n.start)
return(list(h = h, x = x))})
if( is.null(fit@mfit$model$vrmodel) ) sfStop()
H = array(NA, dim = c(n.sim, m, m.sim))
tmpH = array(NA, dim = c(m, m, n.sim))
for(i in 1:n.sim) H[i,,] = t(sapply(xtmp, FUN = function(x) as.numeric(x$h[i,])))
for(i in 1:m.sim){
for(j in 1:n.sim){
tmpH[ , , j] = diag(sqrt( H[j, , i]) ) %*% simR[[i]][ , , j] %*% diag(sqrt( H[j, , i] ) )
}
simH[[i]] = tmpH
}
simxX = array(NA, dim = c(n.sim, m, m.sim))
for(i in 1:n.sim) simxX[i,,] = t(sapply(xtmp, FUN = function(x) as.numeric(x$x[i,])))
simX = vector(mode = "list", length = m.sim)
for(i in 1:m.sim) simX[[i]] = matrix(simxX[,,i], nrow = n.sim)
}
} else{
simH = vector(mode = "list", length = m.sim)
simX = vector(mode = "list", length = m.sim)
mtmp = lapply(as.list(1:m.sim), FUN = function(j) .dccsimf(Z = z[,,j], Qbar = preQ, Rbar = preR, dcca, dccb, mo, n.sim, n.start, m, rseed[j]))
# need to pass m.sim and matrix of simulated z's to ugarchsim (speedup)
simR = lapply(mtmp, FUN = function(x) if(is.matrix(x$R)) array(x$R, dim = c(m, m, n.sim)) else last(x$R, n.sim))
simQ = lapply(mtmp, FUN = function(x) if(is.matrix(x$Q)) array(x$Q, dim = c(m, m, n.sim)) else last(x$Q, n.sim))
stdresid = vector(mode = "list", length = m)
for(i in 1:m) stdresid[[i]] = sapply(mtmp, FUN = function(x) x$stdresid[,i])
xtmp = lapply(as.list(1:m), FUN = function(j){
maxx = fit@ufit@fit[[j]]@model$maxOrder2;
htmp = ugarchsim(fit@ufit@fit[[j]], n.sim = n.sim + n.start, n.start = 0, m.sim = m.sim,
startMethod = startMethod[1], custom.dist = list(name = "sample",
distfit = matrix(stdresid[[j]][-(1:mo), ], ncol = m.sim)),
presigma = if( is.null(presigma) ) NA else tail(presigma[,j], maxx),
preresiduals = if( is.null(preresiduals) ) NA else tail(preresiduals[,j], maxx),
prereturns = if( is.null(prereturns) | !is.null(fit@mfit$model$vrmodel) ) NA else tail(prereturns[,j], maxx),
mexsimdata = if( is.null(fit@mfit$model$vrmodel) ) mexsimdata[[j]] else NULL, vexsimdata = vexsimdata[[j]] )
h = matrix(tail(htmp@simulation$sigmaSim^2, n.sim), nrow = n.sim)
x = matrix(htmp@simulation$seriesSim, nrow = n.sim + n.start)
return(list(h = h, x = x))})
H = array(NA, dim = c(n.sim, m, m.sim))
tmpH = array(NA, dim = c(m, m, n.sim))
for(i in 1:n.sim) H[i,,] = t(sapply(xtmp, FUN = function(x) as.numeric(x$h[i,])))
for(i in 1:m.sim){
for(j in 1:n.sim){
tmpH[ , , j] = diag(sqrt( H[j, , i]) ) %*% simR[[i]][ , , j] %*% diag(sqrt( H[j, , i] ) )
}
simH[[i]] = tmpH
}
simxX = array(NA, dim = c(n.sim, m, m.sim))
for(i in 1:n.sim) simxX[i,,] = t(sapply(xtmp, FUN = function(x) as.numeric(x$x[i,])))
simX = vector(mode = "list", length = m.sim)
for(i in 1:m.sim) simX[[i]] = matrix(simxX[,,i], nrow = n.sim)
}
if( !is.null(fit@mfit$model$vrmodel) ){
# mexsimdata check
mxn = fit@mfit$model$vrmodel$mxn
if(mxn > 0 ){
if( !is.null(mexsimdata) ){
if( !is.list(mexsimdata) | length(mexsimdata) != m.sim ) stop("\ndccsim-->error: mexsimdata must be a list of length equal to m.sim...", call. = FALSE)
for(j in 1:m.sim){
if( !is.matrix(mexsimdata[[j]]) ) stop("\ndccsim-->error: mexsimdata list submatrix must be a matrix...", call. = FALSE)
if( dim(mexsimdata[[j]])[2] != mxn ) stop("\ndccsim-->error: wrong number of external regressors in list submatrix!...", call. = FALSE)
if( dim(mexsimdata[[j]])[1] < (n.sim + n.start) ) stop("\ndccsim-->error: external regressor list submatrix has less points than requested simulation length (n.sim + n.start)...", call. = FALSE)
}
} else{
mexsimdata = vector(mode = "list", length = m.sim)
for(j in 1:m.sim) mexsimdata[[i]] = matrix(0, ncol = mxn, nrow = (n.sim + n.start))
}
} else{
mexsimdata = NULL
}
vrmodel = fit@mfit$model$vrmodel
p = as.numeric(vrmodel$lag)
if(startMethod == "sample"){
if( !is.null(prereturns) ){
if(!is.matrix(prereturns)) stop("\nprereturns must be a matrix\n")
if(dim(prereturns)[1] != p) stop("\nwrong number of rows for prereturns matrix")
if(dim(prereturns)[2] != m) stop("\nwrong number of cols for prereturns matrix")
prereturns = matrix(prereturns, nrow = p, ncol = m)
} else{
prereturns = as.matrix( tail(Data, p) )
}
} else{
prereturns = matrix(0, nrow = p, ncol = m, byrow = TRUE)
}
xlist = vector(mode = "list", length = m.sim)
if( n.sim == 1 && n.start == 0 ){
# SPECIAL ROUTINE TO SPEED UP 1-AHEAD ESTIMATION
tmp = varxsimXX(X = Data, vrmodel$Bcoef, p, m.sim = m.sim, prereturns, resids = t(sapply(simX, FUN = function(x) x)), if(!is.null(mexsimdata)) t(sapply(mexsimdata, FUN = function(x) x)) else NULL)
for(i in 1:m.sim) xlist[[i]] = tmp[i, , drop = FALSE]
simRes = simX
simX = xlist
} else{
if( parallel ){
if( parallel.control$pkg == "multicore" ){
xlist = mclapply(as.list(1:m.sim), FUN = function(i) varxsim(X = Data, vrmodel$Bcoef, p, n.sim,
n.start, prereturns, resids = simX[[i]], mexsimdata[[i]]),
mc.cores = parallel.control$cores)
} else{
# when the GARCH model is without mean equation, simX == simulated residuals
sfExport("Data", "vrmodel", "p", "n.sim", "n.start", "prereturns", "simX", "mexsimdata", local = TRUE)
xlist = sfLapply(as.list(1:m.sim), fun = function(i) rgarch:::varxsim(X = Data, vrmodel$Bcoef, p,
n.sim, n.start, prereturns, resids = simX[[i]], mexsimdata[[i]]))
sfStop()
}
} else{
xlist = lapply(as.list(1:m.sim), FUN = function(i) varxsim(X = Data, vrmodel$Bcoef, p,
n.sim, n.start, prereturns, resids = simX[[i]], mexsimdata[[i]]))
}
simRes = simX
simX = xlist
}
} else{
# reshape
for(j in 1:m.sim) simX[[j]] = tail(simX[[j]], n.sim)
}
msim = list()
msim$varmodel = fit@mfit$vrmodel
msim$simH = simH
msim$simR = simR
msim$simQ = simQ
msim$simX = simX
msim$simRes = simRes
msim$Z = z
msim$rseed = rseed
msim$model$Data = Data
msim$model$garchpars = garchpars
msim$model$dccpars = dccpars
msim$model$n.sim = n.sim
msim$model$m.sim = m.sim
msim$model$n.start = n.start
msim$model$startMethod = startMethod[1]
msim$model$dccOrder = dccOrder
msim$model$distribution = "mvnorm"
ans = new("DCCsim",
msim = msim)
return( ans )
}
dccsim2.norm = function(fitORspec, n.sim = 1000, n.start = 0, m.sim = 1, startMethod = c("unconditional",
"sample"), presigma = NULL, preresiduals = NULL,
prereturns = NULL, preR = NULL, preH = NULL, rseed = NULL, mexsimdata = NULL, vexsimdata = NULL,
parallel = FALSE, parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2), VAR.fit = NULL, ...)
{
spec = fitORspec
# check that VAR.fit is included if used
if( spec@mspec$varmodel$VAR ){
if( is.null(VAR.fit) ) stop("\ndccsim-->error: VAR.fit must not be NULL for VAR method when calling dccsim using spec!", call. = FALSE)
}
dccOrder = spec@mspec$dccOrder
mo = max( dccOrder )
m = length(spec@uspec@spec)
if( is.null(rseed) ){
rseed = as.integer(runif(1, 1, Sys.time()))
} else {
if(length(rseed) == 1) rseed = as.integer(rseed[1]) else rseed = as.integer( rseed[1:m.sim] )
}
# we use the GED distribution with nu = 1 which corresponds to the Laplace case.
if(length(rseed) == 1){
tmp = rnorm(m * (n.sim + n.start + mo) * m.sim, 0, 1)
z = array(tmp, dim = c(n.sim + n.start + mo, m, m.sim))
} else{
z = array(NA, dim = c(n.sim + n.start + mo, m, m.sim))
for(i in 1:m.sim){
set.seed( rseed[i] )
z[,,i] = matrix(rnorm(m * (n.sim + n.start + mo), 0, 1), nrow = n.sim + n.start + mo, ncol = m)
}
}
if(length(rseed) == 1){
rseed = c(rseed, as.integer(runif(m.sim, 1, Sys.time())))
}
if( is.null(preR) ){
stop("\ndccsim-->error: preR cannot be NULL when method uses only DCC specification!")
} else{
chk = dcc.symcheck(preR, m, d = rep(1, m))
}
preQ = preR
# uncvariance has uGARCHfit AND uGARCHspec dispatch methods
uncv = sapply(spec@uspec@spec, FUN = function(x) uncvariance(x))
if( is.null(preH) ){
preH = diag(sqrt(uncv)) %*% preR %*% diag(sqrt(uncv))
} else{
chk = dcc.symcheck(preH, m, d = NULL)
}
if( !is.null(presigma) ){
if( !is.matrix(presigma) ) stop("\ndccsim-->error: presigma must be a matrix.")
if( dim(presigma)[2] != m ) stop("\ndccsim-->error: wrong column dimension for presigma.")
}
if( !is.null(preresiduals) ){
if( !is.matrix(preresiduals) ) stop("\ndccsim-->error: preresiduals must be a matrix.")
if( dim(preresiduals)[2] != m ) stop("\ndccsim-->error: wrong column dimension for preresiduals.")
}
if( !is.null(prereturns) ){
if( !is.matrix(prereturns) ) stop("\ndccsim-->error: prereturns must be a matrix.")
if( dim(prereturns)[2] != m ) stop("\ndccsim-->error: wrong column dimension for prereturns.")
}
garchpars = sapply(spec@uspec@spec, FUN = function(x) unlist(x@optimization.model$fixed.pars))
dccpars = unlist(spec@mspec$optimization.model$fixed.pars)
dcca = dccpars[1:dccOrder[1]]
dccb = dccpars[(dccOrder[1]+1):sum(dccOrder[1:2])]
simRes = simX = simR = simQ = simH = simSeries = vector(mode = "list", length = m.sim)
if( parallel ){
os = .Platform$OS.type
if(is.null(parallel.control$pkg)){
if( os == "windows" ) parallel.control$pkg = "snowfall" else parallel.control$pkg = "multicore"
if( is.null(parallel.control$cores) ) parallel.control$cores = 2
} else{
mtype = match(tolower(parallel.control$pkg[1]), c("multicore", "snowfall"))
if(is.na(mtype)) stop("\nParallel Package type not recognized in parallel.control\n")
parallel.control$pkg = tolower(parallel.control$pkg[1])
if( os == "windows" && parallel.control$pkg == "multicore" ) stop("\nmulticore not supported on windows O/S\n")
if( is.null(parallel.control$cores) ) parallel.control$cores = 2 else parallel.control$cores = as.integer(parallel.control$cores[1])
}
if( parallel.control$pkg == "multicore" ){
simH = vector(mode = "list", length = m.sim)
simX = vector(mode = "list", length = m.sim)
mtmp = mclapply(as.list(1:m.sim), FUN = function(j)
.dccsimf(Z = z[,,j], Qbar = preQ, Rbar = preR, dcca, dccb, mo, n.sim, n.start, m, rseed[j]),
mc.cores = parallel.control$cores)
simR = lapply(mtmp, FUN = function(x) if(is.matrix(x$R)) array(x$R, dim = c(m, m, n.sim)) else last(x$R, n.sim))
simQ = lapply(mtmp, FUN = function(x) if(is.matrix(x$Q)) array(x$Q, dim = c(m, m, n.sim)) else last(x$Q, n.sim))
stdresid = vector(mode = "list", length = m)
for(i in 1:m) stdresid[[i]] = sapply(mtmp, FUN = function(x) x$stdresid[,i])
xtmp = mclapply(as.list(1:m), FUN = function(j){
maxx = spec@uspec@spec[[j]]@optimization.model$maxOrder2;
htmp = ugarchpath(spec@uspec@spec[[j]], n.sim = n.sim + n.start, n.start = 0, m.sim = m.sim,
custom.dist = list(name = "sample", distfit = matrix(stdresid[[j]][-(1:mo), ], ncol = m.sim)),
presigma = if( is.null(presigma) ) NA else tail(presigma[,j], maxx),
preresiduals = if( is.null(preresiduals) ) NA else tail(preresiduals[,j], maxx),
prereturns = if( is.null(prereturns) | spec@mspec$varmodel$VAR ) NA else tail(prereturns[,j], maxx),
mexsimdata = if( !spec@mspec$varmodel$VAR ) mexsimdata[[j]] else NULL, vexsimdata = vexsimdata[[j]] )
h = matrix(tail(htmp@path$sigmaSim^2, n.sim), nrow = n.sim)
x = matrix(htmp@path$seriesSim, nrow = n.sim + n.start)
return(list(h = h, x = x))}, mc.cores = parallel.control$cores)
H = array(NA, dim = c(n.sim, m, m.sim))
tmpH = array(NA, dim = c(m, m, n.sim))
for(i in 1:n.sim) H[i,,] = t(sapply(xtmp, FUN = function(x) as.numeric(x$h[i,])))
for(i in 1:m.sim){
for(j in 1:n.sim){
tmpH[ , , j] = diag(sqrt( H[j, , i]) ) %*% simR[[i]][ , , j] %*% diag(sqrt( H[j, , i] ) )
}
simH[[i]] = tmpH
}
simxX = array(NA, dim = c(n.sim, m, m.sim))
for(i in 1:n.sim) simxX[i,,] = t(sapply(xtmp, FUN = function(x) as.numeric(x$x[i,])))
simX = vector(mode = "list", length = m.sim)
for(i in 1:m.sim) simX[[i]] = matrix(simxX[,,i], nrow = n.sim)
} else{
sfInit(parallel = TRUE, cpus = parallel.control$cores)
sfExport("spec", "z", "preQ", "preR", "dcca", "dccb", "mo", "n.sim", "n.start",
"m", "prereturns", "preresiduals", "presigma", "mexsimdata", "vexsimdata",
"rseed", local = TRUE)
simH = vector(mode = "list", length = m.sim)
simX = vector(mode = "list", length = m.sim)
mtmp = sfLapply(as.list(1:m.sim), fun = function(j) .dccsimf(Z = z[,,j], Qbar = preQ, Rbar = preR, dcca, dccb, mo, n.sim, n.start, m, rseed[j]))
simR = lapply(mtmp, FUN = function(x) if(is.matrix(x$R)) array(x$R, dim = c(m, m, n.sim)) else last(x$R, n.sim))
simQ = lapply(mtmp, FUN = function(x) if(is.matrix(x$Q)) array(x$Q, dim = c(m, m, n.sim)) else last(x$Q, n.sim))
stdresid = vector(mode = "list", length = m)
for(i in 1:m) stdresid[[i]] = sapply(mtmp, FUN = function(x) x$stdresid[,i])
xtmp = sfLapply(as.list(1:m), fun = function(j){
maxx = spec@uspec@spec[[j]]@optimization.model$maxOrder2;
htmp = rgarch::ugarchpath(spec@uspec@spec[[j]], n.sim = n.sim + n.start, n.start = 0, m.sim = m.sim,
custom.dist = list(name = "sample", distfit = matrix(stdresid[[j]][-(1:mo), ], ncol = m.sim)),
presigma = if( is.null(presigma) ) NA else tail(presigma[,j], maxx),
preresiduals = if( is.null(preresiduals) ) NA else tail(preresiduals[,j], maxx),
prereturns = if( is.null(prereturns) | spec@mspec$varmodel$VAR ) NA else tail(prereturns[,j], maxx),
mexsimdata = if( !spec@mspec$varmodel$VAR ) mexsimdata[[j]] else NULL, vexsimdata = vexsimdata[[j]] )
h = matrix(tail(htmp@path$sigmaSim^2, n.sim), nrow = n.sim)
x = matrix(htmp@path$seriesSim, nrow = n.sim + n.start)
return(list(h = h, x = x))})
if( !spec@mspec$varmodel$VAR ) sfStop()
H = array(NA, dim = c(n.sim, m, m.sim))
tmpH = array(NA, dim = c(m, m, n.sim))
for(i in 1:n.sim) H[i,,] = t(sapply(xtmp, FUN = function(x) as.numeric(x$h[i,])))
for(i in 1:m.sim){
for(j in 1:n.sim){
tmpH[ , , j] = diag(sqrt( H[j, , i]) ) %*% simR[[i]][ , , j] %*% diag(sqrt( H[j, , i] ) )
}
simH[[i]] = tmpH
}
simxX = array(NA, dim = c(n.sim, m, m.sim))
for(i in 1:n.sim) simxX[i,,] = t(sapply(xtmp, FUN = function(x) as.numeric(x$x[i,])))
simX = vector(mode = "list", length = m.sim)
for(i in 1:m.sim) simX[[i]] = matrix(simxX[,,i], nrow = n.sim)
}
} else{
simH = vector(mode = "list", length = m.sim)
simX = vector(mode = "list", length = m.sim)
mtmp = lapply(as.list(1:m.sim), FUN = function(j) .dccsimf(Z = z[,,j], Qbar = preQ, Rbar = preR, dcca, dccb, mo, n.sim, n.start, m, rseed[j]))
simR = lapply(mtmp, FUN = function(x) if(is.matrix(x$R)) array(x$R, dim = c(m, m, n.sim)) else last(x$R, n.sim))
simQ = lapply(mtmp, FUN = function(x) if(is.matrix(x$Q)) array(x$Q, dim = c(m, m, n.sim)) else last(x$Q, n.sim))
stdresid = vector(mode = "list", length = m)
for(i in 1:m) stdresid[[i]] = sapply(mtmp, FUN = function(x) x$stdresid[,i])
xtmp = lapply(as.list(1:m), FUN = function(j){
maxx = spec@uspec@spec[[j]]@optimization.model$maxOrder2;
htmp = ugarchpath(spec@uspec@spec[[j]], n.sim = n.sim + n.start, n.start = 0, m.sim = m.sim,
custom.dist = list(name = "sample", distfit = matrix(stdresid[[j]][-(1:mo), ], ncol = m.sim)),
presigma = if( is.null(presigma) ) NA else tail(presigma[,j], maxx),
preresiduals = if( is.null(preresiduals) ) NA else tail(preresiduals[,j], maxx),
prereturns = if( is.null(prereturns) | spec@mspec$varmodel$VAR ) NA else tail(prereturns[,j], maxx),
mexsimdata = if( !spec@mspec$varmodel$VAR ) mexsimdata[[j]] else NULL, vexsimdata = vexsimdata[[j]] )
h = matrix(tail(htmp@path$sigmaSim^2, n.sim), nrow = n.sim)
x = matrix(htmp@path$seriesSim, nrow = n.sim + n.start)
return(list(h = h, x = x))})
H = array(NA, dim = c(n.sim, m, m.sim))
tmpH = array(NA, dim = c(m, m, n.sim))
for(i in 1:n.sim) H[i,,] = t(sapply(xtmp, FUN = function(x) as.numeric(x$h[i,])))
for(i in 1:m.sim){
for(j in 1:n.sim){
tmpH[ , , j] = diag(sqrt( H[j, , i]) ) %*% simR[[i]][ , , j] %*% diag(sqrt( H[j, , i] ) )
}
simH[[i]] = tmpH
}
simxX = array(NA, dim = c(n.sim, m, m.sim))
for(i in 1:n.sim) simxX[i,,] = t(sapply(xtmp, FUN = function(x) as.numeric(x$x[i,])))
simX = vector(mode = "list", length = m.sim)
for(i in 1:m.sim) simX[[i]] = matrix(simxX[,,i], nrow = n.sim)
}
if( spec@mspec$varmodel$VAR ){
# mexsimdata check
mxn = VAR.fit$mxn
Data = VAR.fit$xfitted
if(mxn > 0 ){
if( !is.null(mexsimdata) ){
if( !is.list(mexsimdata) | length(mexsimdata) != m.sim ) stop("\ndccsim-->error: mexsimdata must be a list of length equal to m.sim...", call. = FALSE)
for(j in 1:m.sim){
if( !is.matrix(mexsimdata[[j]]) ) stop("\ndccsim-->error: mexsimdata list submatrix must be a matrix...", call. = FALSE)
if( dim(mexsimdata[[j]])[2] != mxn ) stop("\ndccsim-->error: wrong number of external regressors in list submatrix!...", call. = FALSE)
if( dim(mexsimdata[[j]])[1] < (n.sim + n.start) ) stop("\ndccsim-->error: external regressor list submatrix has less points than requested simulation length (n.sim + n.start)...", call. = FALSE)
}
} else{
mexsimdata = vector(mode = "list", length = m.sim)
for(j in 1:m.sim) mexsimdata[[i]] = matrix(0, ncol = mxn, nrow = (n.sim + n.start))
}
} else{
mexsimdata = NULL
}
p = VAR.fit$lag
if( !is.null(prereturns) ){
if(!is.matrix(prereturns)) stop("\nprereturns must be a matrix\n")
if(dim(prereturns)[1] != p) stop("\nwrong number of rows for prereturns matrix")
if(dim(prereturns)[2] != m) stop("\nwrong number of cols for prereturns matrix")
prereturns = matrix(prereturns, nrow = p, ncol = m)
} else{
prereturns = matrix(0, nrow = p, ncol = m, byrow = TRUE)
}
xlist = vector(mode = "list", length = m.sim)
if( n.sim == 1 && n.start == 0 ){
# SPECIAL ROUTINE TO SPEED UP 1-AHEAD ESTIMATION
tmp = varxsimXX(X = Data, VAR.fit$Bcoef, p, m.sim = m.sim, prereturns, resids = t(sapply(simX, FUN = function(x) x)), if(!is.null(mexsimdata)) t(sapply(mexsimdata, FUN = function(x) x)) else NULL)
for(i in 1:m.sim) xlist[[i]] = tmp[i, , drop = FALSE]
simRes = simX
simX = xlist
} else{
if( parallel ){
if( parallel.control$pkg == "multicore" ){
xlist = mclapply(as.list(1:m.sim), FUN = function(i) varxsim(X = Data, VAR.fit$Bcoef, p, n.sim,
n.start, prereturns, resids = simX[[i]], mexsimdata[[i]]),
mc.cores = parallel.control$cores)
} else{
# when the GARCH model is without mean equation, simX == simulated residuals
sfExport("Data", "VAR.fit", "p", "n.sim", "n.start", "prereturns", "simX", "mexsimdata", local = TRUE)
xlist = sfLapply(as.list(1:m.sim), fun = function(i) rgarch:::varxsim(X = Data, VAR.fit$Bcoef, p,
n.sim, n.start, prereturns, resids = simX[[i]], mexsimdata[[i]]))
sfStop()
}
} else{
xlist = lapply(as.list(1:m.sim), FUN = function(i) varxsim(X = Data, VAR.fit$Bcoef, p,
n.sim, n.start, prereturns, resids = simX[[i]], mexsimdata[[i]]))
}
simRes = simX
simX = xlist
}
} else{
# reshape
for(j in 1:m.sim) simX[[j]] = tail(simX[[j]], n.sim)
}
msim = list()
msim$simH = simH
msim$simR = simR
msim$simQ = simQ
msim$simX = simX
msim$simRes = simRes
msim$Z = z
msim$rseed = rseed
msim$model$Data = NULL
msim$model$garchpars = garchpars
msim$model$dccpars = dccpars
msim$model$n.sim = n.sim
msim$model$m.sim = m.sim
msim$model$n.start = n.start
msim$model$startMethod = "unconditional"
msim$model$dccOrder = dccOrder
msim$model$distribution = "mvnorm"
ans = new("DCCsim",
msim = msim)
return( ans )
}
dccsim1.student = function(fitORspec, n.sim = 1000, n.start = 0, m.sim = 1, startMethod = c("unconditional", "sample"),
presigma = NULL, preresiduals = NULL, prereturns = NULL, preR = NULL, preH = NULL, rseed = NULL, mexsimdata = NULL,
vexsimdata = NULL, parallel = FALSE, parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2), VAR.fit = NULL, ...)
{
fit = fitORspec
Data = head(fit@mfit$origdata, dim(fit@mfit$origdata)[1] - fit@mfit$model$out.sample)
m = dim(Data)[2]
n = dim(Data)[1]
dccOrder = fit@mfit$model$dccOrder
shape = rshape(fit)
mo = max( dccOrder )
if( is.null(rseed) ){
rseed = as.integer(runif(1, 1, Sys.time()))
} else {
if(length(rseed) == 1) rseed = as.integer(rseed[1]) else rseed = as.integer( rseed[1:m.sim] )
}
# we use the GED distribution with nu = 1 which corresponds to the Laplace case.
if(length(rseed) == 1){
tmp = rstd(m * (n.sim + n.start + mo) * m.sim, 0, 1, nu = shape)
z = array(tmp, dim = c(n.sim + n.start + mo, m, m.sim))
} else{
z = array(NA, dim = c(n.sim + n.start + mo, m, m.sim))
for(i in 1:m.sim){
set.seed( rseed[i] )
z[,,i] = matrix(rstd(m * (n.sim + n.start + mo), 0, 1, nu = shape), nrow = n.sim + n.start + mo, ncol = m)
}
}
if(length(rseed) == 1){
rseed = c(rseed, as.integer(runif(m.sim, 1, Sys.time())))
}
if( is.null(preR) ){
if( startMethod[1] == "sample" ){
preR = last(rcor(fit))[,,1]
diag(preR) = 1
} else{
preR = cor(as.matrix(Data))
}
} else{
chk = dcc.symcheck(preR, m, d = rep(1, m))
}
preQ = preR
uncv = sapply(fit@ufit@fit, FUN = function(x) uncvariance(x))
if( is.null(preH) ){
preH = diag(sqrt(uncv)) %*% preR %*% diag(sqrt(uncv))
} else{
chk = dcc.symcheck(preH, m, d = NULL)
}
if( !is.null(presigma) ){
if( !is.matrix(presigma) ) stop("\ndccsim-->error: presigma must be a matrix.")
if( dim(presigma)[2] != m ) stop("\ndccsim-->error: wrong column dimension for presigma.")
}
if( !is.null(preresiduals) ){
if( !is.matrix(preresiduals) ) stop("\ndccsim-->error: preresiduals must be a matrix.")
if( dim(preresiduals)[2] != m ) stop("\ndccsim-->error: wrong column dimension for preresiduals.")
}
if( !is.null(prereturns) ){
if( !is.matrix(prereturns) ) stop("\ndccsim-->error: prereturns must be a matrix.")
if( dim(prereturns)[2] != m ) stop("\ndccsim-->error: wrong column dimension for prereturns.")
}
garchpars = coef(fit, type = "garch")
dccpars = coef(fit, type = "dcc")
dcca = dccpars[1:dccOrder[1]]
dccb = dccpars[(dccOrder[1]+1):sum(dccOrder[1:2])]
simRes = simX = simR = simQ = simH = simSeries = vector(mode = "list", length = m.sim)
if( parallel ){
os = .Platform$OS.type
if(is.null(parallel.control$pkg)){
if( os == "windows" ) parallel.control$pkg = "snowfall" else parallel.control$pkg = "multicore"
if( is.null(parallel.control$cores) ) parallel.control$cores = 2
} else{
mtype = match(tolower(parallel.control$pkg[1]), c("multicore", "snowfall"))
if(is.na(mtype)) stop("\nParallel Package type not recognized in parallel.control\n")
parallel.control$pkg = tolower(parallel.control$pkg[1])
if( os == "windows" && parallel.control$pkg == "multicore" ) stop("\nmulticore not supported on windows O/S\n")
if( is.null(parallel.control$cores) ) parallel.control$cores = 2 else parallel.control$cores = as.integer(parallel.control$cores[1])
}
if( parallel.control$pkg == "multicore" ){
simH = vector(mode = "list", length = m.sim)
simX = vector(mode = "list", length = m.sim)
mtmp = mclapply(as.list(1:m.sim), FUN = function(j)
.dccsimf(Z = z[,,j], Qbar = preQ, Rbar = preR, dcca, dccb, mo, n.sim, n.start, m, rseed[j]),
mc.cores = parallel.control$cores)
# need to pass m.sim and matrix of simulated z's to ugarchsim (speedup)
simR = lapply(mtmp, FUN = function(x) if(is.matrix(x$R)) array(x$R, dim = c(m, m, n.sim)) else last(x$R, n.sim))
simQ = lapply(mtmp, FUN = function(x) if(is.matrix(x$Q)) array(x$Q, dim = c(m, m, n.sim)) else last(x$Q, n.sim))
stdresid = vector(mode = "list", length = m)
for(i in 1:m) stdresid[[i]] = sapply(mtmp, FUN = function(x) x$stdresid[,i])
xtmp = mclapply(as.list(1:m), FUN = function(j){
maxx = fit@ufit@fit[[j]]@model$maxOrder2;
htmp = ugarchsim(fit@ufit@fit[[j]], n.sim = n.sim + n.start, n.start = 0, m.sim = m.sim,
startMethod = startMethod[1], custom.dist = list(name = "sample",
distfit = matrix(stdresid[[j]][-(1:mo), ], ncol = m.sim)),
presigma = if( is.null(presigma) ) NA else tail(presigma[,j], maxx),
preresiduals = if( is.null(preresiduals) ) NA else tail(preresiduals[,j], maxx),
prereturns = if( is.null(prereturns) | !is.null(fit@mfit$model$vrmodel) ) NA else tail(prereturns[,j], maxx),
mexsimdata = if( is.null(fit@mfit$model$vrmodel) ) mexsimdata[[j]] else NULL, vexsimdata = vexsimdata[[j]] )
h = matrix(tail(htmp@simulation$sigmaSim^2, n.sim), nrow = n.sim)
x = matrix(htmp@simulation$seriesSim, nrow = n.sim + n.start)
return(list(h = h, x = x))}, mc.cores = parallel.control$cores)
H = array(NA, dim = c(n.sim, m, m.sim))
tmpH = array(NA, dim = c(m, m, n.sim))
for(i in 1:n.sim) H[i,,] = t(sapply(xtmp, FUN = function(x) as.numeric(x$h[i,])))
for(i in 1:m.sim){
for(j in 1:n.sim){
tmpH[ , , j] = diag(sqrt( H[j, , i]) ) %*% simR[[i]][ , , j] %*% diag(sqrt( H[j, , i] ) )
}
simH[[i]] = tmpH
}
simxX = array(NA, dim = c(n.sim, m, m.sim))
for(i in 1:n.sim) simxX[i,,] = t(sapply(xtmp, FUN = function(x) as.numeric(x$x[i,])))
simX = vector(mode = "list", length = m.sim)
for(i in 1:m.sim) simX[[i]] = matrix(simxX[,,i], nrow = n.sim)
} else{
simH = vector(mode = "list", length = m.sim)
simX = vector(mode = "list", length = m.sim)
sfInit(parallel = TRUE, cpus = parallel.control$cores)
sfExport("z", "preQ", "preR", "dcca", "dccb", "mo", "n.sim", "n.start", "m", "rseed", local = TRUE)
mtmp = sfLapply(as.list(1:m.sim), fun = function(j)
rgarch:::.dccsimf(Z = z[,,j], Qbar = preQ, Rbar = preR, dcca, dccb, mo, n.sim, n.start, m, rseed[j]))
# need to pass m.sim and matrix of simulated z's to ugarchsim (speedup)
simR = lapply(mtmp, FUN = function(x) if(is.matrix(x$R)) array(x$R, dim = c(m, m, n.sim)) else last(x$R, n.sim))
simQ = lapply(mtmp, FUN = function(x) if(is.matrix(x$Q)) array(x$Q, dim = c(m, m, n.sim)) else last(x$Q, n.sim))
stdresid = vector(mode = "list", length = m)
for(i in 1:m) stdresid[[i]] = sapply(mtmp, FUN = function(x) x$stdresid[,i])
sfExport("fit", "n.sim", "n.start", "m.sim", "startMethod", "stdresid", "presigma",
"preresiduals", "prereturns", "mexsimdata", "vexsimdata", local = TRUE)
xtmp = sfLapply(as.list(1:m), fun = function(j){
maxx = fit@ufit@fit[[j]]@model$maxOrder2;
htmp = rgarch::ugarchsim(fit@ufit@fit[[j]], n.sim = n.sim + n.start, n.start = 0, m.sim = m.sim,
startMethod = startMethod[1], custom.dist = list(name = "sample",
distfit = matrix(stdresid[[j]][-(1:mo), ], ncol = m.sim)),
presigma = if( is.null(presigma) ) NA else tail(presigma[,j], maxx),
preresiduals = if( is.null(preresiduals) ) NA else tail(preresiduals[,j], maxx),
prereturns = if( is.null(prereturns) | !is.null(fit@mfit$model$vrmodel) ) NA else tail(prereturns[,j], maxx),
mexsimdata = if( is.null(fit@mfit$model$vrmodel) ) mexsimdata[[j]] else NULL, vexsimdata = vexsimdata[[j]] )
h = matrix(tail(htmp@simulation$sigmaSim^2, n.sim), nrow = n.sim)
x = matrix(htmp@simulation$seriesSim, nrow = n.sim + n.start)
return(list(h = h, x = x))})
if( is.null(fit@mfit$model$vrmodel) ) sfStop()
H = array(NA, dim = c(n.sim, m, m.sim))
tmpH = array(NA, dim = c(m, m, n.sim))
for(i in 1:n.sim) H[i,,] = t(sapply(xtmp, FUN = function(x) as.numeric(x$h[i,])))
for(i in 1:m.sim){
for(j in 1:n.sim){
tmpH[ , , j] = diag(sqrt( H[j, , i]) ) %*% simR[[i]][ , , j] %*% diag(sqrt( H[j, , i] ) )
}
simH[[i]] = tmpH
}
simxX = array(NA, dim = c(n.sim, m, m.sim))
for(i in 1:n.sim) simxX[i,,] = t(sapply(xtmp, FUN = function(x) as.numeric(x$x[i,])))
simX = vector(mode = "list", length = m.sim)
for(i in 1:m.sim) simX[[i]] = matrix(simxX[,,i], nrow = n.sim)
}
} else{
simH = vector(mode = "list", length = m.sim)
simX = vector(mode = "list", length = m.sim)
mtmp = lapply(as.list(1:m.sim), FUN = function(j) .dccsimf(Z = z[,,j], Qbar = preQ, Rbar = preR, dcca, dccb, mo, n.sim, n.start, m, rseed[j]))
# need to pass m.sim and matrix of simulated z's to ugarchsim (speedup)
simR = lapply(mtmp, FUN = function(x) if(is.matrix(x$R)) array(x$R, dim = c(m, m, n.sim)) else last(x$R, n.sim))
simQ = lapply(mtmp, FUN = function(x) if(is.matrix(x$Q)) array(x$Q, dim = c(m, m, n.sim)) else last(x$Q, n.sim))
stdresid = vector(mode = "list", length = m)
for(i in 1:m) stdresid[[i]] = sapply(mtmp, FUN = function(x) x$stdresid[,i])
xtmp = lapply(as.list(1:m), FUN = function(j){
maxx = fit@ufit@fit[[j]]@model$maxOrder2;
htmp = ugarchsim(fit@ufit@fit[[j]], n.sim = n.sim + n.start, n.start = 0, m.sim = m.sim,
startMethod = startMethod[1], custom.dist = list(name = "sample",
distfit = matrix(stdresid[[j]][-(1:mo), ], ncol = m.sim)),
presigma = if( is.null(presigma) ) NA else tail(presigma[,j], maxx),
preresiduals = if( is.null(preresiduals) ) NA else tail(preresiduals[,j], maxx),
prereturns = if( is.null(prereturns) | !is.null(fit@mfit$model$vrmodel) ) NA else tail(prereturns[,j], maxx),
mexsimdata = if( is.null(fit@mfit$model$vrmodel) ) mexsimdata[[j]] else NULL, vexsimdata = vexsimdata[[j]] )
h = matrix(tail(htmp@simulation$sigmaSim^2, n.sim), nrow = n.sim)
x = matrix(htmp@simulation$seriesSim, nrow = n.sim + n.start)
return(list(h = h, x = x))})
H = array(NA, dim = c(n.sim, m, m.sim))
tmpH = array(NA, dim = c(m, m, n.sim))
for(i in 1:n.sim) H[i,,] = t(sapply(xtmp, FUN = function(x) as.numeric(x$h[i,])))
for(i in 1:m.sim){
for(j in 1:n.sim){
tmpH[ , , j] = diag(sqrt( H[j, , i]) ) %*% simR[[i]][ , , j] %*% diag(sqrt( H[j, , i] ) )
}
simH[[i]] = tmpH
}
simxX = array(NA, dim = c(n.sim, m, m.sim))
for(i in 1:n.sim) simxX[i,,] = t(sapply(xtmp, FUN = function(x) as.numeric(x$x[i,])))
simX = vector(mode = "list", length = m.sim)
for(i in 1:m.sim) simX[[i]] = matrix(simxX[,,i], nrow = n.sim)
}
if( !is.null(fit@mfit$model$vrmodel) ){
# mexsimdata check
mxn = fit@mfit$model$vrmodel$mxn
if(mxn > 0 ){
if( !is.null(mexsimdata) ){
if( !is.list(mexsimdata) | length(mexsimdata) != m.sim ) stop("\ndccsim-->error: mexsimdata must be a list of length equal to m.sim...", call. = FALSE)
for(j in 1:m.sim){
if( !is.matrix(mexsimdata[[j]]) ) stop("\ndccsim-->error: mexsimdata list submatrix must be a matrix...", call. = FALSE)
if( dim(mexsimdata[[j]])[2] != mxn ) stop("\ndccsim-->error: wrong number of external regressors in list submatrix!...", call. = FALSE)
if( dim(mexsimdata[[j]])[1] < (n.sim + n.start) ) stop("\ndccsim-->error: external regressor list submatrix has less points than requested simulation length (n.sim + n.start)...", call. = FALSE)
}
} else{
mexsimdata = vector(mode = "list", length = m.sim)
for(j in 1:m.sim) mexsimdata[[i]] = matrix(0, ncol = mxn, nrow = (n.sim + n.start))
}
} else{
mexsimdata = NULL
}
vrmodel = fit@mfit$model$vrmodel
p = as.numeric(vrmodel$lag)
if(startMethod == "sample"){
if( !is.null(prereturns) ){
if(!is.matrix(prereturns)) stop("\nprereturns must be a matrix\n")
if(dim(prereturns)[1] != p) stop("\nwrong number of rows for prereturns matrix")
if(dim(prereturns)[2] != m) stop("\nwrong number of cols for prereturns matrix")
prereturns = matrix(prereturns, nrow = p, ncol = m)
} else{
prereturns = as.matrix( tail(Data, p) )
}
} else{
prereturns = matrix(0, nrow = p, ncol = m, byrow = TRUE)
}
xlist = vector(mode = "list", length = m.sim)
if( n.sim == 1 && n.start == 0 ){
# SPECIAL ROUTINE TO SPEED UP 1-AHEAD ESTIMATION
tmp = varxsimXX(X = Data, vrmodel$Bcoef, p, m.sim = m.sim, prereturns, resids = t(sapply(simX, FUN = function(x) x)), if(!is.null(mexsimdata)) t(sapply(mexsimdata, FUN = function(x) x)) else NULL)
for(i in 1:m.sim) xlist[[i]] = tmp[i, , drop = FALSE]
simRes = simX
simX = xlist
} else{
if( parallel ){
if( parallel.control$pkg == "multicore" ){
xlist = mclapply(as.list(1:m.sim), FUN = function(i) varxsim(X = Data, vrmodel$Bcoef, p, n.sim,
n.start, prereturns, resids = simX[[i]], mexsimdata[[i]]),
mc.cores = parallel.control$cores)
} else{
# when the GARCH model is without mean equation, simX == simulated residuals
sfExport("Data", "vrmodel", "p", "n.sim", "n.start", "prereturns", "simX", "mexsimdata", local = TRUE)
xlist = sfLapply(as.list(1:m.sim), fun = function(i) rgarch:::varxsim(X = Data, vrmodel$Bcoef, p,
n.sim, n.start, prereturns, resids = simX[[i]], mexsimdata[[i]]))
sfStop()
}
} else{
xlist = lapply(as.list(1:m.sim), FUN = function(i) varxsim(X = Data, vrmodel$Bcoef, p,
n.sim, n.start, prereturns, resids = simX[[i]], mexsimdata[[i]]))
}
simRes = simX
simX = xlist
}
} else{
# reshape
for(j in 1:m.sim) simX[[j]] = tail(simX[[j]], n.sim)
}
msim = list()
msim$simH = simH
msim$simR = simR
msim$simQ = simQ
msim$simX = simX
msim$simRes = simRes
msim$Z = z
msim$rseed = rseed
msim$model$Data = Data
msim$model$garchpars = garchpars
msim$model$dccpars = dccpars
msim$model$n.sim = n.sim
msim$model$m.sim = m.sim
msim$model$n.start = n.start
msim$model$startMethod = startMethod[1]
msim$model$dccOrder = dccOrder
msim$model$distribution = "mvt"
ans = new("DCCsim",
msim = msim)
return( ans )
}
dccsim2.student = function(fitORspec, n.sim = 1000, n.start = 0, m.sim = 1, startMethod = c("unconditional",
"sample"), presigma = NULL, preresiduals = NULL,
prereturns = NULL, preR = NULL, preH = NULL, rseed = NULL, mexsimdata = NULL,
vexsimdata = NULL, parallel = FALSE, parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2), VAR.fit = NULL, ...)
{
spec = fitORspec
# check that VAR.fit is included if used
if( spec@mspec$varmodel$VAR ){
if( is.null(VAR.fit) ) stop("\ndccsim-->error: VAR.fit must not be NULL for VAR method when calling dccsim using spec!", call. = FALSE)
}
dccOrder = spec@mspec$dccOrder
mo = max( dccOrder )
m = length(spec@uspec@spec)
shape = spec@mspec$optimization.model$fixed.pars["shape"]
if( is.null(rseed) ){
rseed = as.integer(runif(1, 1, Sys.time()))
} else {
if(length(rseed) == 1) rseed = as.integer(rseed[1]) else rseed = as.integer( rseed[1:m.sim] )
}
# we use the GED distribution with nu = 1 which corresponds to the Laplace case.
if(length(rseed) == 1){
tmp = rstd(m * (n.sim + n.start + mo) * m.sim, 0, 1, nu = shape)
z = array(tmp, dim = c(n.sim + n.start + mo, m, m.sim))
} else{
z = array(NA, dim = c(n.sim + n.start + mo, m, m.sim))
for(i in 1:m.sim){
set.seed( rseed[i] )
z[,,i] = matrix(rstd(m * (n.sim + n.start + mo), 0, 1, nu = shape), nrow = n.sim + n.start + mo, ncol = m)
}
}
if(length(rseed) == 1){
rseed = c(rseed, as.integer(runif(m.sim, 1, Sys.time())))
}
if( is.null(preR) ){
stop("\ndccsim-->error: preR cannot be NULL when method uses only DCC specification!\n")
} else{
chk = dcc.symcheck(preR, m, d = rep(1, m))
}
preQ = preR
# uncvariance has uGARCHfit AND uGARCHspec dispatch methods
uncv = sapply(spec@uspec@spec, FUN = function(x) uncvariance(x))
if( is.null(preH) ){
preH = diag(sqrt(uncv)) %*% preR %*% diag(sqrt(uncv))
} else{
chk = dcc.symcheck(preH, m, d = NULL)
}
if( !is.null(presigma) ){
if( !is.matrix(presigma) ) stop("\ndccsim-->error: presigma must be a matrix.")
if( dim(presigma)[2] != m ) stop("\ndccsim-->error: wrong column dimension for presigma.")
}
if( !is.null(preresiduals) ){
if( !is.matrix(preresiduals) ) stop("\ndccsim-->error: preresiduals must be a matrix.")
if( dim(preresiduals)[2] != m ) stop("\ndccsim-->error: wrong column dimension for preresiduals.")
}
if( !is.null(prereturns) ){
if( !is.matrix(prereturns) ) stop("\ndccsim-->error: prereturns must be a matrix.")
if( dim(prereturns)[2] != m ) stop("\ndccsim-->error: wrong column dimension for prereturns.")
}
garchpars = sapply(spec@uspec@spec, FUN = function(x) unlist(x@optimization.model$fixed.pars))
dccpars = unlist(spec@mspec$optimization.model$fixed.pars)
dcca = dccpars[1:dccOrder[1]]
dccb = dccpars[(dccOrder[1]+1):sum(dccOrder[1:2])]
simRes = simX = simR = simQ = simH = simSeries = vector(mode = "list", length = m.sim)
if( parallel ){
os = .Platform$OS.type
if(is.null(parallel.control$pkg)){
if( os == "windows" ) parallel.control$pkg = "snowfall" else parallel.control$pkg = "multicore"
if( is.null(parallel.control$cores) ) parallel.control$cores = 2
} else{
mtype = match(tolower(parallel.control$pkg[1]), c("multicore", "snowfall"))
if(is.na(mtype)) stop("\nParallel Package type not recognized in parallel.control\n")
parallel.control$pkg = tolower(parallel.control$pkg[1])
if( os == "windows" && parallel.control$pkg == "multicore" ) stop("\nmulticore not supported on windows O/S\n")
if( is.null(parallel.control$cores) ) parallel.control$cores = 2 else parallel.control$cores = as.integer(parallel.control$cores[1])
}
if( parallel.control$pkg == "multicore" ){
simH = vector(mode = "list", length = m.sim)
simX = vector(mode = "list", length = m.sim)
mtmp = mclapply(as.list(1:m.sim), FUN = function(j)
.dccsimf(Z = z[,,j], Qbar = preQ, Rbar = preR, dcca, dccb, mo, n.sim, n.start, m, rseed[j]),
mc.cores = parallel.control$cores)
simR = lapply(mtmp, FUN = function(x) if(is.matrix(x$R)) array(x$R, dim = c(m, m, n.sim)) else last(x$R, n.sim))
simQ = lapply(mtmp, FUN = function(x) if(is.matrix(x$Q)) array(x$Q, dim = c(m, m, n.sim)) else last(x$Q, n.sim))
stdresid = vector(mode = "list", length = m)
for(i in 1:m) stdresid[[i]] = sapply(mtmp, FUN = function(x) x$stdresid[,i])
xtmp = mclapply(as.list(1:m), FUN = function(j){
maxx = spec@uspec@spec[[j]]@optimization.model$maxOrder2;
htmp = ugarchpath(spec@uspec@spec[[j]], n.sim = n.sim + n.start, n.start = 0, m.sim = m.sim,
custom.dist = list(name = "sample", distfit = matrix(stdresid[[j]][-(1:mo), ], ncol = m.sim)),
presigma = if( is.null(presigma) ) NA else tail(presigma[,j], maxx),
preresiduals = if( is.null(preresiduals) ) NA else tail(preresiduals[,j], maxx),
prereturns = if( is.null(prereturns) | spec@mspec$varmodel$VAR ) NA else tail(prereturns[,j], maxx),
mexsimdata = if( !spec@mspec$varmodel$VAR ) mexsimdata[[j]] else NULL, vexsimdata = vexsimdata[[j]] )
h = matrix(tail(htmp@path$sigmaSim^2, n.sim), nrow = n.sim)
x = matrix(htmp@path$seriesSim, nrow = n.sim + n.start)
return(list(h = h, x = x))}, mc.cores = parallel.control$cores)
H = array(NA, dim = c(n.sim, m, m.sim))
tmpH = array(NA, dim = c(m, m, n.sim))
for(i in 1:n.sim) H[i,,] = t(sapply(xtmp, FUN = function(x) as.numeric(x$h[i,])))
for(i in 1:m.sim){
for(j in 1:n.sim){
tmpH[ , , j] = diag(sqrt( H[j, , i]) ) %*% simR[[i]][ , , j] %*% diag(sqrt( H[j, , i] ) )
}
simH[[i]] = tmpH
}
simxX = array(NA, dim = c(n.sim, m, m.sim))
for(i in 1:n.sim) simxX[i,,] = t(sapply(xtmp, FUN = function(x) as.numeric(x$x[i,])))
simX = vector(mode = "list", length = m.sim)
for(i in 1:m.sim) simX[[i]] = matrix(simxX[,,i], nrow = n.sim)
} else{
sfInit(parallel = TRUE, cpus = parallel.control$cores)
sfExport("spec", "z", "preQ", "preR", "dcca", "dccb", "mo", "n.sim", "n.start",
"m", "prereturns", "preresiduals", "presigma", "mexsimdata", "vexsimdata",
"rseed", local = TRUE)
simH = vector(mode = "list", length = m.sim)
simX = vector(mode = "list", length = m.sim)
mtmp = sfLapply(as.list(1:m.sim), fun = function(j) .dccsimf(Z = z[,,j], Qbar = preQ, Rbar = preR, dcca, dccb, mo, n.sim, n.start, m, rseed[j]))
simR = lapply(mtmp, FUN = function(x) if(is.matrix(x$R)) array(x$R, dim = c(m, m, n.sim)) else last(x$R, n.sim))
simQ = lapply(mtmp, FUN = function(x) if(is.matrix(x$Q)) array(x$Q, dim = c(m, m, n.sim)) else last(x$Q, n.sim))
stdresid = vector(mode = "list", length = m)
for(i in 1:m) stdresid[[i]] = sapply(mtmp, FUN = function(x) x$stdresid[,i])
xtmp = sfLapply(as.list(1:m), fun = function(j){
maxx = spec@uspec@spec[[j]]@optimization.model$maxOrder2;
htmp = rgarch::ugarchpath(spec@uspec@spec[[j]], n.sim = n.sim + n.start, n.start = 0, m.sim = m.sim,
custom.dist = list(name = "sample", distfit = matrix(stdresid[[j]][-(1:mo), ], ncol = m.sim)),
presigma = if( is.null(presigma) ) NA else tail(presigma[,j], maxx),
preresiduals = if( is.null(preresiduals) ) NA else tail(preresiduals[,j], maxx),
prereturns = if( is.null(prereturns) | spec@mspec$varmodel$VAR ) NA else tail(prereturns[,j], maxx),
mexsimdata = if( !spec@mspec$varmodel$VAR ) mexsimdata[[j]] else NULL, vexsimdata = vexsimdata[[j]] )
h = matrix(tail(htmp@path$sigmaSim^2, n.sim), nrow = n.sim)
x = matrix(htmp@path$seriesSim, nrow = n.sim + n.start)
return(list(h = h, x = x))})
if( !spec@mspec$varmodel$VAR ) sfStop()
H = array(NA, dim = c(n.sim, m, m.sim))
tmpH = array(NA, dim = c(m, m, n.sim))
for(i in 1:n.sim) H[i,,] = t(sapply(xtmp, FUN = function(x) as.numeric(x$h[i,])))
for(i in 1:m.sim){
for(j in 1:n.sim){
tmpH[ , , j] = diag(sqrt( H[j, , i]) ) %*% simR[[i]][ , , j] %*% diag(sqrt( H[j, , i] ) )
}
simH[[i]] = tmpH
}
simxX = array(NA, dim = c(n.sim, m, m.sim))
for(i in 1:n.sim) simxX[i,,] = t(sapply(xtmp, FUN = function(x) as.numeric(x$x[i,])))
simX = vector(mode = "list", length = m.sim)
for(i in 1:m.sim) simX[[i]] = matrix(simxX[,,i], nrow = n.sim)
}
} else{
simH = vector(mode = "list", length = m.sim)
simX = vector(mode = "list", length = m.sim)
mtmp = lapply(as.list(1:m.sim), FUN = function(j) .dccsimf(Z = z[,,j], Qbar = preQ, Rbar = preR, dcca, dccb, mo, n.sim, n.start, m, rseed[j]))
simR = lapply(mtmp, FUN = function(x) if(is.matrix(x$R)) array(x$R, dim = c(m, m, n.sim)) else last(x$R, n.sim))
simQ = lapply(mtmp, FUN = function(x) if(is.matrix(x$Q)) array(x$Q, dim = c(m, m, n.sim)) else last(x$Q, n.sim))
stdresid = vector(mode = "list", length = m)
for(i in 1:m) stdresid[[i]] = sapply(mtmp, FUN = function(x) x$stdresid[,i])
xtmp = lapply(as.list(1:m), FUN = function(j){
maxx = spec@uspec@spec[[j]]@optimization.model$maxOrder2;
htmp = ugarchpath(spec@uspec@spec[[j]], n.sim = n.sim + n.start, n.start = 0, m.sim = m.sim,
custom.dist = list(name = "sample", distfit = matrix(stdresid[[j]][-(1:mo), ], ncol = m.sim)),
presigma = if( is.null(presigma) ) NA else tail(presigma[,j], maxx),
preresiduals = if( is.null(preresiduals) ) NA else tail(preresiduals[,j], maxx),
prereturns = if( is.null(prereturns) | spec@mspec$varmodel$VAR ) NA else tail(prereturns[,j], maxx),
mexsimdata = if( !spec@mspec$varmodel$VAR ) mexsimdata[[j]] else NULL, vexsimdata = vexsimdata[[j]] )
h = matrix(tail(htmp@path$sigmaSim^2, n.sim), nrow = n.sim)
x = matrix(htmp@path$seriesSim, nrow = n.sim + n.start)
return(list(h = h, x = x))})
H = array(NA, dim = c(n.sim, m, m.sim))
tmpH = array(NA, dim = c(m, m, n.sim))
for(i in 1:n.sim) H[i,,] = t(sapply(xtmp, FUN = function(x) as.numeric(x$h[i,])))
for(i in 1:m.sim){
for(j in 1:n.sim){
tmpH[ , , j] = diag(sqrt( H[j, , i]) ) %*% simR[[i]][ , , j] %*% diag(sqrt( H[j, , i] ) )
}
simH[[i]] = tmpH
}
simxX = array(NA, dim = c(n.sim, m, m.sim))
for(i in 1:n.sim) simxX[i,,] = t(sapply(xtmp, FUN = function(x) as.numeric(x$x[i,])))
simX = vector(mode = "list", length = m.sim)
for(i in 1:m.sim) simX[[i]] = matrix(simxX[,,i], nrow = n.sim)
}
if( spec@mspec$varmodel$VAR ){
# mexsimdata check
mxn = VAR.fit$mxn
Data = VAR.fit$xfitted
if(mxn > 0 ){
if( !is.null(mexsimdata) ){
if( !is.list(mexsimdata) | length(mexsimdata) != m.sim ) stop("\ndccsim-->error: mexsimdata must be a list of length equal to m.sim...", call. = FALSE)
for(j in 1:m.sim){
if( !is.matrix(mexsimdata[[j]]) ) stop("\ndccsim-->error: mexsimdata list submatrix must be a matrix...", call. = FALSE)
if( dim(mexsimdata[[j]])[2] != mxn ) stop("\ndccsim-->error: wrong number of external regressors in list submatrix!...", call. = FALSE)
if( dim(mexsimdata[[j]])[1] < (n.sim + n.start) ) stop("\ndccsim-->error: external regressor list submatrix has less points than requested simulation length (n.sim + n.start)...", call. = FALSE)
}
} else{
mexsimdata = vector(mode = "list", length = m.sim)
for(j in 1:m.sim) mexsimdata[[i]] = matrix(0, ncol = mxn, nrow = (n.sim + n.start))
}
} else{
mexsimdata = NULL
}
p = VAR.fit$lag
if( !is.null(prereturns) ){
if(!is.matrix(prereturns)) stop("\nprereturns must be a matrix\n")
if(dim(prereturns)[1] != p) stop("\nwrong number of rows for prereturns matrix")
if(dim(prereturns)[2] != m) stop("\nwrong number of cols for prereturns matrix")
prereturns = matrix(prereturns, nrow = p, ncol = m)
} else{
prereturns = matrix(0, nrow = p, ncol = m, byrow = TRUE)
}
xlist = vector(mode = "list", length = m.sim)
if( n.sim == 1 && n.start == 0 ){
# SPECIAL ROUTINE TO SPEED UP 1-AHEAD ESTIMATION
tmp = varxsimXX(X = Data, VAR.fit$Bcoef, p, m.sim = m.sim, prereturns, resids = t(sapply(simX, FUN = function(x) x)), if(!is.null(mexsimdata)) t(sapply(mexsimdata, FUN = function(x) x)) else NULL)
for(i in 1:m.sim) xlist[[i]] = tmp[i, , drop = FALSE]
simRes = simX
simX = xlist
} else{
if( parallel ){
if( parallel.control$pkg == "multicore" ){
xlist = mclapply(as.list(1:m.sim), FUN = function(i) varxsim(X = Data, VAR.fit$Bcoef, p, n.sim,
n.start, prereturns, resids = simX[[i]], mexsimdata[[i]]),
mc.cores = parallel.control$cores)
} else{
# when the GARCH model is without mean equation, simX == simulated residuals
sfExport("Data", "VAR.fit", "p", "n.sim", "n.start", "prereturns", "simX", "mexsimdata", local = TRUE)
xlist = sfLapply(as.list(1:m.sim), fun = function(i) rgarch:::varxsim(X = Data, VAR.fit$Bcoef, p,
n.sim, n.start, prereturns, resids = simX[[i]], mexsimdata[[i]]))
sfStop()
}
} else{
xlist = lapply(as.list(1:m.sim), FUN = function(i) varxsim(X = Data, VAR.fit$Bcoef, p,
n.sim, n.start, prereturns, resids = simX[[i]], mexsimdata[[i]]))
}
simRes = simX
simX = xlist
}
} else{
# reshape
for(j in 1:m.sim) simX[[j]] = tail(simX[[j]], n.sim)
}
msim = list()
msim$simH = simH
msim$simR = simR
msim$simQ = simQ
msim$simX = simX
msim$simRes = simRes
msim$Z = z
msim$rseed = rseed
msim$model$Data = NULL
msim$model$garchpars = garchpars
msim$model$dccpars = dccpars
msim$model$n.sim = n.sim
msim$model$m.sim = m.sim
msim$model$n.start = n.start
msim$model$startMethod = "unconditional"
msim$model$dccOrder = dccOrder
msim$model$distribution = "mvt"
ans = new("DCCsim",
msim = msim)
return( ans )
}
dccsim1.laplace = function(fitORspec, n.sim = 1000, n.start = 0, m.sim = 1, startMethod = c("unconditional", "sample"),
presigma = NULL, preresiduals = NULL, prereturns = NULL, preR = NULL, preH = NULL, rseed = NULL, mexsimdata = NULL,
vexsimdata = NULL, parallel = FALSE, parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2), VAR.fit = NULL, ...)
{
fit = fitORspec
Data = head(fit@mfit$origdata, dim(fit@mfit$origdata)[1] - fit@mfit$model$out.sample)
m = dim(Data)[2]
n = dim(Data)[1]
dccOrder = fit@mfit$model$dccOrder
mo = max( dccOrder )
if( is.null(rseed) ){
rseed = as.integer(runif(1, 1, Sys.time()))
} else {
if(length(rseed) == 1) rseed = as.integer(rseed[1]) else rseed = as.integer( rseed[1:m.sim] )
}
# we use the GED distribution with nu = 1 which corresponds to the Laplace case.
if(length(rseed) == 1){
tmp = rged(m * (n.sim + n.start + mo) * m.sim, 0, 1, nu = 1)
z = array(tmp, dim = c(n.sim + n.start + mo, m, m.sim))
} else{
z = array(NA, dim = c(n.sim + n.start + mo, m, m.sim))
for(i in 1:m.sim){
set.seed( rseed[i] )
z[,,i] = matrix(rged(m * (n.sim + n.start + mo), 0, 1, nu = 1), nrow = n.sim + n.start + mo, ncol = m)
}
}
if(length(rseed) == 1){
rseed = c(rseed, as.integer(runif(m.sim, 1, Sys.time())))
}
if( is.null(preR) ){
if( startMethod[1] == "sample" ){
preR = last(rcor(fit))[,,1]
diag(preR) = 1
} else{
preR = cor(as.matrix(Data))
}
} else{
chk = dcc.symcheck(preR, m, d = rep(1, m))
}
preQ = preR
uncv = sapply(fit@ufit@fit, FUN = function(x) uncvariance(x))
if( is.null(preH) ){
preH = diag(sqrt(uncv)) %*% preR %*% diag(sqrt(uncv))
} else{
chk = dcc.symcheck(preH, m, d = NULL)
}
if( !is.null(presigma) ){
if( !is.matrix(presigma) ) stop("\ndccsim-->error: presigma must be a matrix.")
if( dim(presigma)[2] != m ) stop("\ndccsim-->error: wrong column dimension for presigma.")
}
if( !is.null(preresiduals) ){
if( !is.matrix(preresiduals) ) stop("\ndccsim-->error: preresiduals must be a matrix.")
if( dim(preresiduals)[2] != m ) stop("\ndccsim-->error: wrong column dimension for preresiduals.")
}
if( !is.null(prereturns) ){
if( !is.matrix(prereturns) ) stop("\ndccsim-->error: prereturns must be a matrix.")
if( dim(prereturns)[2] != m ) stop("\ndccsim-->error: wrong column dimension for prereturns.")
}
garchpars = coef(fit, type = "garch")
dccpars = coef(fit, type = "dcc")
dcca = dccpars[1:dccOrder[1]]
dccb = dccpars[(dccOrder[1]+1):sum(dccOrder[1:2])]
simRes = simX = simR = simQ = simH = simSeries = vector(mode = "list", length = m.sim)
if( parallel ){
os = .Platform$OS.type
if(is.null(parallel.control$pkg)){
if( os == "windows" ) parallel.control$pkg = "snowfall" else parallel.control$pkg = "multicore"
if( is.null(parallel.control$cores) ) parallel.control$cores = 2
} else{
mtype = match(tolower(parallel.control$pkg[1]), c("multicore", "snowfall"))
if(is.na(mtype)) stop("\nParallel Package type not recognized in parallel.control\n")
parallel.control$pkg = tolower(parallel.control$pkg[1])
if( os == "windows" && parallel.control$pkg == "multicore" ) stop("\nmulticore not supported on windows O/S\n")
if( is.null(parallel.control$cores) ) parallel.control$cores = 2 else parallel.control$cores = as.integer(parallel.control$cores[1])
}
if( parallel.control$pkg == "multicore" ){
simH = vector(mode = "list", length = m.sim)
simX = vector(mode = "list", length = m.sim)
mtmp = mclapply(as.list(1:m.sim), FUN = function(j)
.dccsimf(Z = z[,,j], Qbar = preQ, Rbar = preR, dcca, dccb, mo, n.sim, n.start, m, rseed[j]),
mc.cores = parallel.control$cores)
# need to pass m.sim and matrix of simulated z's to ugarchsim (speedup)
simR = lapply(mtmp, FUN = function(x) if(is.matrix(x$R)) array(x$R, dim = c(m, m, n.sim)) else last(x$R, n.sim))
simQ = lapply(mtmp, FUN = function(x) if(is.matrix(x$Q)) array(x$Q, dim = c(m, m, n.sim)) else last(x$Q, n.sim))
stdresid = vector(mode = "list", length = m)
for(i in 1:m) stdresid[[i]] = sapply(mtmp, FUN = function(x) x$stdresid[,i])
xtmp = mclapply(as.list(1:m), FUN = function(j){
maxx = fit@ufit@fit[[j]]@model$maxOrder2;
htmp = ugarchsim(fit@ufit@fit[[j]], n.sim = n.sim + n.start, n.start = 0, m.sim = m.sim,
startMethod = startMethod[1], custom.dist = list(name = "sample",
distfit = matrix(stdresid[[j]][-(1:mo), ], ncol = m.sim)),
presigma = if( is.null(presigma) ) NA else tail(presigma[,j], maxx),
preresiduals = if( is.null(preresiduals) ) NA else tail(preresiduals[,j], maxx),
prereturns = if( is.null(prereturns) | !is.null(fit@mfit$model$vrmodel) ) NA else tail(prereturns[,j], maxx),
mexsimdata = if( is.null(fit@mfit$model$vrmodel) ) mexsimdata[[j]] else NULL, vexsimdata = vexsimdata[[j]] )
h = matrix(tail(htmp@simulation$sigmaSim^2, n.sim), nrow = n.sim)
x = matrix(htmp@simulation$seriesSim, nrow = n.sim + n.start)
return(list(h = h, x = x))}, mc.cores = parallel.control$cores)
H = array(NA, dim = c(n.sim, m, m.sim))
tmpH = array(NA, dim = c(m, m, n.sim))
for(i in 1:n.sim) H[i,,] = t(sapply(xtmp, FUN = function(x) as.numeric(x$h[i,])))
for(i in 1:m.sim){
for(j in 1:n.sim){
tmpH[ , , j] = diag(sqrt( H[j, , i]) ) %*% simR[[i]][ , , j] %*% diag(sqrt( H[j, , i] ) )
}
simH[[i]] = tmpH
}
simxX = array(NA, dim = c(n.sim, m, m.sim))
for(i in 1:n.sim) simxX[i,,] = t(sapply(xtmp, FUN = function(x) as.numeric(x$x[i,])))
simX = vector(mode = "list", length = m.sim)
for(i in 1:m.sim) simX[[i]] = matrix(simxX[,,i], nrow = n.sim)
} else{
simH = vector(mode = "list", length = m.sim)
simX = vector(mode = "list", length = m.sim)
sfInit(parallel = TRUE, cpus = parallel.control$cores)
sfExport("z", "preQ", "preR", "dcca", "dccb", "mo", "n.sim", "n.start", "m", "rseed", local = TRUE)
mtmp = sfLapply(as.list(1:m.sim), fun = function(j)
rgarch:::.dccsimf(Z = z[,,j], Qbar = preQ, Rbar = preR, dcca, dccb, mo, n.sim, n.start, m, rseed[j]))
# need to pass m.sim and matrix of simulated z's to ugarchsim (speedup)
simR = lapply(mtmp, FUN = function(x) if(is.matrix(x$R)) array(x$R, dim = c(m, m, n.sim)) else last(x$R, n.sim))
simQ = lapply(mtmp, FUN = function(x) if(is.matrix(x$Q)) array(x$Q, dim = c(m, m, n.sim)) else last(x$Q, n.sim))
stdresid = vector(mode = "list", length = m)
for(i in 1:m) stdresid[[i]] = sapply(mtmp, FUN = function(x) x$stdresid[,i])
sfExport("fit", "n.sim", "n.start", "m.sim", "startMethod", "stdresid", "presigma",
"preresiduals", "prereturns", "mexsimdata", "vexsimdata", local = TRUE)
xtmp = sfLapply(as.list(1:m), fun = function(j){
maxx = fit@ufit@fit[[j]]@model$maxOrder2;
htmp = rgarch::ugarchsim(fit@ufit@fit[[j]], n.sim = n.sim + n.start, n.start = 0, m.sim = m.sim,
startMethod = startMethod[1], custom.dist = list(name = "sample",
distfit = matrix(stdresid[[j]][-(1:mo), ], ncol = m.sim)),
presigma = if( is.null(presigma) ) NA else tail(presigma[,j], maxx),
preresiduals = if( is.null(preresiduals) ) NA else tail(preresiduals[,j], maxx),
prereturns = if( is.null(prereturns) | !is.null(fit@mfit$model$vrmodel) ) NA else tail(prereturns[,j], maxx),
mexsimdata = if( is.null(fit@mfit$model$vrmodel) ) mexsimdata[[j]] else NULL, vexsimdata = vexsimdata[[j]] )
h = matrix(tail(htmp@simulation$sigmaSim^2, n.sim), nrow = n.sim)
x = matrix(htmp@simulation$seriesSim, nrow = n.sim + n.start)
return(list(h = h, x = x))})
if( is.null(fit@mfit$model$vrmodel) ) sfStop()
H = array(NA, dim = c(n.sim, m, m.sim))
tmpH = array(NA, dim = c(m, m, n.sim))
for(i in 1:n.sim) H[i,,] = t(sapply(xtmp, FUN = function(x) as.numeric(x$h[i,])))
for(i in 1:m.sim){
for(j in 1:n.sim){
tmpH[ , , j] = diag(sqrt( H[j, , i]) ) %*% simR[[i]][ , , j] %*% diag(sqrt( H[j, , i] ) )
}
simH[[i]] = tmpH
}
simxX = array(NA, dim = c(n.sim, m, m.sim))
for(i in 1:n.sim) simxX[i,,] = t(sapply(xtmp, FUN = function(x) as.numeric(x$x[i,])))
simX = vector(mode = "list", length = m.sim)
for(i in 1:m.sim) simX[[i]] = matrix(simxX[,,i], nrow = n.sim)
}
} else{
simH = vector(mode = "list", length = m.sim)
simX = vector(mode = "list", length = m.sim)
mtmp = lapply(as.list(1:m.sim), FUN = function(j) .dccsimf(Z = z[,,j], Qbar = preQ, Rbar = preR, dcca, dccb, mo, n.sim, n.start, m, rseed[j]))
# need to pass m.sim and matrix of simulated z's to ugarchsim (speedup)
simR = lapply(mtmp, FUN = function(x) if(is.matrix(x$R)) array(x$R, dim = c(m, m, n.sim)) else last(x$R, n.sim))
simQ = lapply(mtmp, FUN = function(x) if(is.matrix(x$Q)) array(x$Q, dim = c(m, m, n.sim)) else last(x$Q, n.sim))
stdresid = vector(mode = "list", length = m)
for(i in 1:m) stdresid[[i]] = sapply(mtmp, FUN = function(x) x$stdresid[,i])
xtmp = lapply(as.list(1:m), FUN = function(j){
maxx = fit@ufit@fit[[j]]@model$maxOrder2;
htmp = ugarchsim(fit@ufit@fit[[j]], n.sim = n.sim + n.start, n.start = 0, m.sim = m.sim,
startMethod = startMethod[1], custom.dist = list(name = "sample",
distfit = matrix(stdresid[[j]][-(1:mo), ], ncol = m.sim)),
presigma = if( is.null(presigma) ) NA else tail(presigma[,j], maxx),
preresiduals = if( is.null(preresiduals) ) NA else tail(preresiduals[,j], maxx),
prereturns = if( is.null(prereturns) | !is.null(fit@mfit$model$vrmodel) ) NA else tail(prereturns[,j], maxx),
mexsimdata = if( is.null(fit@mfit$model$vrmodel) ) mexsimdata[[j]] else NULL, vexsimdata = vexsimdata[[j]] )
h = matrix(tail(htmp@simulation$sigmaSim^2, n.sim), nrow = n.sim)
x = matrix(htmp@simulation$seriesSim, nrow = n.sim + n.start)
return(list(h = h, x = x))})
H = array(NA, dim = c(n.sim, m, m.sim))
tmpH = array(NA, dim = c(m, m, n.sim))
for(i in 1:n.sim) H[i,,] = t(sapply(xtmp, FUN = function(x) as.numeric(x$h[i,])))
for(i in 1:m.sim){
for(j in 1:n.sim){
tmpH[ , , j] = diag(sqrt( H[j, , i]) ) %*% simR[[i]][ , , j] %*% diag(sqrt( H[j, , i] ) )
}
simH[[i]] = tmpH
}
simxX = array(NA, dim = c(n.sim, m, m.sim))
for(i in 1:n.sim) simxX[i,,] = t(sapply(xtmp, FUN = function(x) as.numeric(x$x[i,])))
simX = vector(mode = "list", length = m.sim)
for(i in 1:m.sim) simX[[i]] = matrix(simxX[,,i], nrow = n.sim)
}
if( !is.null(fit@mfit$model$vrmodel) ){
# mexsimdata check
mxn = fit@mfit$model$vrmodel$mxn
if(mxn > 0 ){
if( !is.null(mexsimdata) ){
if( !is.list(mexsimdata) | length(mexsimdata) != m.sim ) stop("\ndccsim-->error: mexsimdata must be a list of length equal to m.sim...", call. = FALSE)
for(j in 1:m.sim){
if( !is.matrix(mexsimdata[[j]]) ) stop("\ndccsim-->error: mexsimdata list submatrix must be a matrix...", call. = FALSE)
if( dim(mexsimdata[[j]])[2] != mxn ) stop("\ndccsim-->error: wrong number of external regressors in list submatrix!...", call. = FALSE)
if( dim(mexsimdata[[j]])[1] < (n.sim + n.start) ) stop("\ndccsim-->error: external regressor list submatrix has less points than requested simulation length (n.sim + n.start)...", call. = FALSE)
}
} else{
mexsimdata = vector(mode = "list", length = m.sim)
for(j in 1:m.sim) mexsimdata[[i]] = matrix(0, ncol = mxn, nrow = (n.sim + n.start))
}
} else{
mexsimdata = NULL
}
vrmodel = fit@mfit$model$vrmodel
p = as.numeric(vrmodel$lag)
if(startMethod == "sample"){
if( !is.null(prereturns) ){
if(!is.matrix(prereturns)) stop("\nprereturns must be a matrix\n")
if(dim(prereturns)[1] != p) stop("\nwrong number of rows for prereturns matrix")
if(dim(prereturns)[2] != m) stop("\nwrong number of cols for prereturns matrix")
prereturns = matrix(prereturns, nrow = p, ncol = m)
} else{
prereturns = as.matrix( tail(Data, p) )
}
} else{
prereturns = matrix(0, nrow = p, ncol = m, byrow = TRUE)
}
xlist = vector(mode = "list", length = m.sim)
if( n.sim == 1 && n.start == 0 ){
# SPECIAL ROUTINE TO SPEED UP 1-AHEAD ESTIMATION
tmp = varxsimXX(X = Data, vrmodel$Bcoef, p, m.sim = m.sim, prereturns, resids = t(sapply(simX, FUN = function(x) x)), if(!is.null(mexsimdata)) t(sapply(mexsimdata, FUN = function(x) x)) else NULL)
for(i in 1:m.sim) xlist[[i]] = tmp[i, , drop = FALSE]
simRes = simX
simX = xlist
} else{
if( parallel ){
if( parallel.control$pkg == "multicore" ){
xlist = mclapply(as.list(1:m.sim), FUN = function(i) varxsim(X = Data, vrmodel$Bcoef, p, n.sim,
n.start, prereturns, resids = simX[[i]], mexsimdata[[i]]),
mc.cores = parallel.control$cores)
} else{
# when the GARCH model is without mean equation, simX == simulated residuals
sfExport("Data", "vrmodel", "p", "n.sim", "n.start", "prereturns", "simX", "mexsimdata", local = TRUE)
xlist = sfLapply(as.list(1:m.sim), fun = function(i) rgarch:::varxsim(X = Data, vrmodel$Bcoef, p,
n.sim, n.start, prereturns, resids = simX[[i]], mexsimdata[[i]]))
sfStop()
}
} else{
xlist = lapply(as.list(1:m.sim), FUN = function(i) varxsim(X = Data, vrmodel$Bcoef, p,
n.sim, n.start, prereturns, resids = simX[[i]], mexsimdata[[i]]))
}
simRes = simX
simX = xlist
}
} else{
# reshape
for(j in 1:m.sim) simX[[j]] = tail(simX[[j]], n.sim)
}
msim = list()
msim$simH = simH
msim$simR = simR
msim$simQ = simQ
msim$simX = simX
msim$simRes = simRes
msim$Z = z
msim$rseed = rseed
msim$model$Data = Data
msim$model$garchpars = garchpars
msim$model$dccpars = dccpars
msim$model$n.sim = n.sim
msim$model$m.sim = m.sim
msim$model$n.start = n.start
msim$model$startMethod = startMethod[1]
msim$model$dccOrder = dccOrder
msim$model$distribution = "mvlaplace"
ans = new("DCCsim",
msim = msim)
return( ans )
}
dccsim2.laplace = function(fitORspec, n.sim = 1000, n.start = 0, m.sim = 1, startMethod = c("unconditional",
"sample"), presigma = NULL, preresiduals = NULL,
prereturns = NULL, preR = NULL, preH = NULL, rseed = NULL, mexsimdata = NULL, vexsimdata = NULL,
parallel = FALSE, parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2), VAR.fit = NULL, ...)
{
spec = fitORspec
# check that VAR.fit is included if used
if( spec@mspec$varmodel$VAR ){
if( is.null(VAR.fit) ) stop("\ndccsim-->error: VAR.fit must not be NULL for VAR method when calling dccsim using spec!", call. = FALSE)
}
dccOrder = spec@mspec$dccOrder
mo = max( dccOrder )
m = length(spec@uspec@spec)
if( is.null(rseed) ){
rseed = as.integer(runif(1, 1, Sys.time()))
} else {
if(length(rseed) == 1) rseed = as.integer(rseed[1]) else rseed = as.integer( rseed[1:m.sim] )
}
# we use the GED distribution with nu = 1 which corresponds to the Laplace case.
if(length(rseed) == 1){
tmp = rged(m * (n.sim + n.start + mo) * m.sim, 0, 1, nu = 1)
z = array(tmp, dim = c(n.sim + n.start + mo, m, m.sim))
} else{
z = array(NA, dim = c(n.sim + n.start + mo, m, m.sim))
for(i in 1:m.sim){
set.seed( rseed[i] )
z[,,i] = matrix(rged(m * (n.sim + n.start + mo), 0, 1, nu = 1), nrow = n.sim + n.start + mo, ncol = m)
}
}
if(length(rseed) == 1){
rseed = c(rseed, as.integer(runif(m.sim, 1, Sys.time())))
}
if( is.null(preR) ){
stop("\ndccsim-->error: preR cannot be NULL when method uses only DCC specification!")
} else{
chk = dcc.symcheck(preR, m, d = rep(1, m))
}
preQ = preR
# uncvariance has uGARCHfit AND uGARCHspec dispatch methods
uncv = sapply(spec@uspec@spec, FUN = function(x) uncvariance(x))
if( is.null(preH) ){
preH = diag(sqrt(uncv)) %*% preR %*% diag(sqrt(uncv))
} else{
chk = dcc.symcheck(preH, m, d = NULL)
}
if( !is.null(presigma) ){
if( !is.matrix(presigma) ) stop("\ndccsim-->error: presigma must be a matrix.")
if( dim(presigma)[2] != m ) stop("\ndccsim-->error: wrong column dimension for presigma.")
}
if( !is.null(preresiduals) ){
if( !is.matrix(preresiduals) ) stop("\ndccsim-->error: preresiduals must be a matrix.")
if( dim(preresiduals)[2] != m ) stop("\ndccsim-->error: wrong column dimension for preresiduals.")
}
if( !is.null(prereturns) ){
if( !is.matrix(prereturns) ) stop("\ndccsim-->error: prereturns must be a matrix.")
if( dim(prereturns)[2] != m ) stop("\ndccsim-->error: wrong column dimension for prereturns.")
}
garchpars = sapply(spec@uspec@spec, FUN = function(x) unlist(x@optimization.model$fixed.pars))
dccpars = unlist(spec@mspec$optimization.model$fixed.pars)
dcca = dccpars[1:dccOrder[1]]
dccb = dccpars[(dccOrder[1]+1):sum(dccOrder[1:2])]
simRes = simX = simR = simQ = simH = simSeries = vector(mode = "list", length = m.sim)
if( parallel ){
os = .Platform$OS.type
if(is.null(parallel.control$pkg)){
if( os == "windows" ) parallel.control$pkg = "snowfall" else parallel.control$pkg = "multicore"
if( is.null(parallel.control$cores) ) parallel.control$cores = 2
} else{
mtype = match(tolower(parallel.control$pkg[1]), c("multicore", "snowfall"))
if(is.na(mtype)) stop("\nParallel Package type not recognized in parallel.control\n")
parallel.control$pkg = tolower(parallel.control$pkg[1])
if( os == "windows" && parallel.control$pkg == "multicore" ) stop("\nmulticore not supported on windows O/S\n")
if( is.null(parallel.control$cores) ) parallel.control$cores = 2 else parallel.control$cores = as.integer(parallel.control$cores[1])
}
if( parallel.control$pkg == "multicore" ){
simH = vector(mode = "list", length = m.sim)
simX = vector(mode = "list", length = m.sim)
mtmp = mclapply(as.list(1:m.sim), FUN = function(j)
.dccsimf(Z = z[,,j], Qbar = preQ, Rbar = preR, dcca, dccb, mo, n.sim, n.start, m, rseed[j]),
mc.cores = parallel.control$cores)
simR = lapply(mtmp, FUN = function(x) if(is.matrix(x$R)) array(x$R, dim = c(m, m, n.sim)) else last(x$R, n.sim))
simQ = lapply(mtmp, FUN = function(x) if(is.matrix(x$Q)) array(x$Q, dim = c(m, m, n.sim)) else last(x$Q, n.sim))
stdresid = vector(mode = "list", length = m)
for(i in 1:m) stdresid[[i]] = sapply(mtmp, FUN = function(x) x$stdresid[,i])
xtmp = mclapply(as.list(1:m), FUN = function(j){
maxx = spec@uspec@spec[[j]]@optimization.model$maxOrder2;
htmp = ugarchpath(spec@uspec@spec[[j]], n.sim = n.sim + n.start, n.start = 0, m.sim = m.sim,
custom.dist = list(name = "sample", distfit = matrix(stdresid[[j]][-(1:mo), ], ncol = m.sim)),
presigma = if( is.null(presigma) ) NA else tail(presigma[,j], maxx),
preresiduals = if( is.null(preresiduals) ) NA else tail(preresiduals[,j], maxx),
prereturns = if( is.null(prereturns) | spec@mspec$varmodel$VAR ) NA else tail(prereturns[,j], maxx),
mexsimdata = if( !spec@mspec$varmodel$VAR ) mexsimdata[[j]] else NULL, vexsimdata = vexsimdata[[j]] )
h = matrix(tail(htmp@path$sigmaSim^2, n.sim), nrow = n.sim)
x = matrix(htmp@path$seriesSim, nrow = n.sim + n.start)
return(list(h = h, x = x))}, mc.cores = parallel.control$cores)
H = array(NA, dim = c(n.sim, m, m.sim))
tmpH = array(NA, dim = c(m, m, n.sim))
for(i in 1:n.sim) H[i,,] = t(sapply(xtmp, FUN = function(x) as.numeric(x$h[i,])))
for(i in 1:m.sim){
for(j in 1:n.sim){
tmpH[ , , j] = diag(sqrt( H[j, , i]) ) %*% simR[[i]][ , , j] %*% diag(sqrt( H[j, , i] ) )
}
simH[[i]] = tmpH
}
simxX = array(NA, dim = c(n.sim, m, m.sim))
for(i in 1:n.sim) simxX[i,,] = t(sapply(xtmp, FUN = function(x) as.numeric(x$x[i,])))
simX = vector(mode = "list", length = m.sim)
for(i in 1:m.sim) simX[[i]] = matrix(simxX[,,i], nrow = n.sim)
} else{
sfInit(parallel = TRUE, cpus = parallel.control$cores)
sfExport("spec", "z", "preQ", "preR", "dcca", "dccb", "mo", "n.sim", "n.start",
"m", "prereturns", "preresiduals", "presigma", "mexsimdata", "vexsimdata",
"rseed", local = TRUE)
simH = vector(mode = "list", length = m.sim)
simX = vector(mode = "list", length = m.sim)
mtmp = sfLapply(as.list(1:m.sim), fun = function(j) .dccsimf(Z = z[,,j], Qbar = preQ, Rbar = preR, dcca, dccb, mo, n.sim, n.start, m, rseed[j]))
simR = lapply(mtmp, FUN = function(x) if(is.matrix(x$R)) array(x$R, dim = c(m, m, n.sim)) else last(x$R, n.sim))
simQ = lapply(mtmp, FUN = function(x) if(is.matrix(x$Q)) array(x$Q, dim = c(m, m, n.sim)) else last(x$Q, n.sim))
stdresid = vector(mode = "list", length = m)
for(i in 1:m) stdresid[[i]] = sapply(mtmp, FUN = function(x) x$stdresid[,i])
xtmp = sfLapply(as.list(1:m), fun = function(j){
maxx = spec@uspec@spec[[j]]@optimization.model$maxOrder2;
htmp = rgarch::ugarchpath(spec@uspec@spec[[j]], n.sim = n.sim + n.start, n.start = 0, m.sim = m.sim,
custom.dist = list(name = "sample", distfit = matrix(stdresid[[j]][-(1:mo), ], ncol = m.sim)),
presigma = if( is.null(presigma) ) NA else tail(presigma[,j], maxx),
preresiduals = if( is.null(preresiduals) ) NA else tail(preresiduals[,j], maxx),
prereturns = if( is.null(prereturns) | spec@mspec$varmodel$VAR ) NA else tail(prereturns[,j], maxx),
mexsimdata = if( !spec@mspec$varmodel$VAR ) mexsimdata[[j]] else NULL, vexsimdata = vexsimdata[[j]] )
h = matrix(tail(htmp@path$sigmaSim^2, n.sim), nrow = n.sim)
x = matrix(htmp@path$seriesSim, nrow = n.sim + n.start)
return(list(h = h, x = x))})
if( !spec@mspec$varmodel$VAR ) sfStop()
H = array(NA, dim = c(n.sim, m, m.sim))
tmpH = array(NA, dim = c(m, m, n.sim))
for(i in 1:n.sim) H[i,,] = t(sapply(xtmp, FUN = function(x) as.numeric(x$h[i,])))
for(i in 1:m.sim){
for(j in 1:n.sim){
tmpH[ , , j] = diag(sqrt( H[j, , i]) ) %*% simR[[i]][ , , j] %*% diag(sqrt( H[j, , i] ) )
}
simH[[i]] = tmpH
}
simxX = array(NA, dim = c(n.sim, m, m.sim))
for(i in 1:n.sim) simxX[i,,] = t(sapply(xtmp, FUN = function(x) as.numeric(x$x[i,])))
simX = vector(mode = "list", length = m.sim)
for(i in 1:m.sim) simX[[i]] = matrix(simxX[,,i], nrow = n.sim)
}
} else{
simH = vector(mode = "list", length = m.sim)
simX = vector(mode = "list", length = m.sim)
mtmp = lapply(as.list(1:m.sim), FUN = function(j) .dccsimf(Z = z[,,j], Qbar = preQ, Rbar = preR, dcca, dccb, mo, n.sim, n.start, m, rseed[j]))
simR = lapply(mtmp, FUN = function(x) if(is.matrix(x$R)) array(x$R, dim = c(m, m, n.sim)) else last(x$R, n.sim))
simQ = lapply(mtmp, FUN = function(x) if(is.matrix(x$Q)) array(x$Q, dim = c(m, m, n.sim)) else last(x$Q, n.sim))
stdresid = vector(mode = "list", length = m)
for(i in 1:m) stdresid[[i]] = sapply(mtmp, FUN = function(x) x$stdresid[,i])
xtmp = lapply(as.list(1:m), FUN = function(j){
maxx = spec@uspec@spec[[j]]@optimization.model$maxOrder2;
htmp = ugarchpath(spec@uspec@spec[[j]], n.sim = n.sim + n.start, n.start = 0, m.sim = m.sim,
custom.dist = list(name = "sample", distfit = matrix(stdresid[[j]][-(1:mo), ], ncol = m.sim)),
presigma = if( is.null(presigma) ) NA else tail(presigma[,j], maxx),
preresiduals = if( is.null(preresiduals) ) NA else tail(preresiduals[,j], maxx),
prereturns = if( is.null(prereturns) | spec@mspec$varmodel$VAR ) NA else tail(prereturns[,j], maxx),
mexsimdata = if( !spec@mspec$varmodel$VAR ) mexsimdata[[j]] else NULL, vexsimdata = vexsimdata[[j]] )
h = matrix(tail(htmp@path$sigmaSim^2, n.sim), nrow = n.sim)
x = matrix(htmp@path$seriesSim, nrow = n.sim + n.start)
return(list(h = h, x = x))})
H = array(NA, dim = c(n.sim, m, m.sim))
tmpH = array(NA, dim = c(m, m, n.sim))
for(i in 1:n.sim) H[i,,] = t(sapply(xtmp, FUN = function(x) as.numeric(x$h[i,])))
for(i in 1:m.sim){
for(j in 1:n.sim){
tmpH[ , , j] = diag(sqrt( H[j, , i]) ) %*% simR[[i]][ , , j] %*% diag(sqrt( H[j, , i] ) )
}
simH[[i]] = tmpH
}
simxX = array(NA, dim = c(n.sim, m, m.sim))
for(i in 1:n.sim) simxX[i,,] = t(sapply(xtmp, FUN = function(x) as.numeric(x$x[i,])))
simX = vector(mode = "list", length = m.sim)
for(i in 1:m.sim) simX[[i]] = matrix(simxX[,,i], nrow = n.sim)
}
if( spec@mspec$varmodel$VAR ){
# mexsimdata check
mxn = VAR.fit$mxn
Data = VAR.fit$xfitted
if(mxn > 0 ){
if( !is.null(mexsimdata) ){
if( !is.list(mexsimdata) | length(mexsimdata) != m.sim ) stop("\ndccsim-->error: mexsimdata must be a list of length equal to m.sim...", call. = FALSE)
for(j in 1:m.sim){
if( !is.matrix(mexsimdata[[j]]) ) stop("\ndccsim-->error: mexsimdata list submatrix must be a matrix...", call. = FALSE)
if( dim(mexsimdata[[j]])[2] != mxn ) stop("\ndccsim-->error: wrong number of external regressors in list submatrix!...", call. = FALSE)
if( dim(mexsimdata[[j]])[1] < (n.sim + n.start) ) stop("\ndccsim-->error: external regressor list submatrix has less points than requested simulation length (n.sim + n.start)...", call. = FALSE)
}
} else{
mexsimdata = vector(mode = "list", length = m.sim)
for(j in 1:m.sim) mexsimdata[[i]] = matrix(0, ncol = mxn, nrow = (n.sim + n.start))
}
} else{
mexsimdata = NULL
}
p = VAR.fit$lag
if( !is.null(prereturns) ){
if(!is.matrix(prereturns)) stop("\nprereturns must be a matrix\n")
if(dim(prereturns)[1] != p) stop("\nwrong number of rows for prereturns matrix")
if(dim(prereturns)[2] != m) stop("\nwrong number of cols for prereturns matrix")
prereturns = matrix(prereturns, nrow = p, ncol = m)
} else{
prereturns = matrix(0, nrow = p, ncol = m, byrow = TRUE)
}
xlist = vector(mode = "list", length = m.sim)
if( n.sim == 1 && n.start == 0 ){
# SPECIAL ROUTINE TO SPEED UP 1-AHEAD ESTIMATION
tmp = varxsimXX(X = Data, VAR.fit$Bcoef, p, m.sim = m.sim, prereturns, resids = t(sapply(simX, FUN = function(x) x)), if(!is.null(mexsimdata)) t(sapply(mexsimdata, FUN = function(x) x)) else NULL)
for(i in 1:m.sim) xlist[[i]] = tmp[i, , drop = FALSE]
simRes = simX
simX = xlist
} else{
if( parallel ){
if( parallel.control$pkg == "multicore" ){
xlist = mclapply(as.list(1:m.sim), FUN = function(i) varxsim(X = Data, VAR.fit$Bcoef, p, n.sim,
n.start, prereturns, resids = simX[[i]], mexsimdata[[i]]),
mc.cores = parallel.control$cores)
} else{
# when the GARCH model is without mean equation, simX == simulated residuals
sfExport("Data", "VAR.fit", "p", "n.sim", "n.start", "prereturns", "simX", "mexsimdata", local = TRUE)
xlist = sfLapply(as.list(1:m.sim), fun = function(i) rgarch:::varxsim(X = Data, VAR.fit$Bcoef, p,
n.sim, n.start, prereturns, resids = simX[[i]], mexsimdata[[i]]))
sfStop()
}
} else{
xlist = lapply(as.list(1:m.sim), FUN = function(i) varxsim(X = Data, VAR.fit$Bcoef, p,
n.sim, n.start, prereturns, resids = simX[[i]], mexsimdata[[i]]))
}
simRes = simX
simX = xlist
}
} else{
# reshape
for(j in 1:m.sim) simX[[j]] = tail(simX[[j]], n.sim)
}
msim = list()
msim$simH = simH
msim$simR = simR
msim$simQ = simQ
msim$simX = simX
msim$simRes = simRes
msim$Z = z
msim$rseed = rseed
msim$model$Data = NULL
msim$model$garchpars = garchpars
msim$model$dccpars = dccpars
msim$model$n.sim = n.sim
msim$model$m.sim = m.sim
msim$model$n.start = n.start
msim$model$startMethod = "unconditional"
msim$model$dccOrder = dccOrder
msim$model$distribution = "mvlaplace"
ans = new("DCCsim",
msim = msim)
return( ans )
}
###################################################################
# Functions for parallel evaluation based on large m.sim
.mcsim = function(fit, z, preQ, preR, dcca, dccb, mo, n.sim, n.start, m, startMethod, prereturns, preresiduals,
presigma, mexsimdata, vexsimdata, rseed)
{
tmp = .dccsimf(Z = z, Qbar = preQ, Rbar = preR, dcca, dccb, mo, n.sim, n.start, m, rseed)
simR = if(is.matrix(tmp$R)) array(tmp$R, dim = c(m, m, n.sim)) else last(tmp$R, n.sim)
simQ = if(is.matrix(tmp$Q)) array(tmp$Q, dim = c(m, m, n.sim)) else last(tmp$Q, n.sim)
stdresid = tmp$stdresid
X = H = matrix(NA, nrow = n.sim + n.start, ncol = m)
simH = array(NA, dim = c(m, m, n.sim))
for(i in 1:m){
maxx = fit@ufit@fit[[i]]@model$maxOrder2
htmp = ugarchsim(fit@ufit@fit[[i]], n.sim = n.sim + n.start, n.start = 0, m.sim = 1,
startMethod = startMethod[1], custom.dist = list(name = "sample",
distfit = matrix(stdresid[-(1:mo), i], ncol = 1)),
presigma = if( is.null(presigma) ) NA else tail(presigma[,i], maxx),
preresiduals = if( is.null(preresiduals) ) NA else tail(preresiduals[,i], maxx),
prereturns = if( is.null(prereturns) | !is.null(fit@mfit$model$vrmodel) ) NA else tail(prereturns[,i], maxx),
mexsimdata = if( is.null(fit@mfit$model$vrmodel) ) mexsimdata[[i]] else NULL, vexsimdata = vexsimdata[[i]] )
H[, i] = htmp@simulation$sigmaSim^2
X[, i] = htmp@simulation$seriesSim
}
H = tail(H, n.sim)
for(i in 1:n.sim){
simH[ , , i] = diag(sqrt( H[i, ]) ) %*% simR[ , , i] %*% diag(sqrt( H[i, ]) )
}
simX = X
return( list(simH = simH, simX = simX, simR = simR, simQ = simQ) )
}
.mcpath = function(spec, z, preQ, preR, dcca, dccb, mo, n.sim, n.start, m, prereturns, preresiduals,
presigma, mexsimdata, vexsimdata, rseed)
{
tmp = .dccsimf(Z = z, Qbar = preQ, Rbar = preR, dcca, dccb, mo, n.sim, n.start, m, rseed)
simR = if(is.matrix(tmp$R)) array(tmp$R, dim = c(m, m, n.sim)) else last(tmp$R, n.sim)
simQ = if(is.matrix(tmp$Q)) array(tmp$Q, dim = c(m, m, n.sim)) else last(tmp$Q, n.sim)
stdresid = tmp$stdresid
X = H = matrix(NA, nrow = n.sim + n.start, ncol = m)
simH = array(NA, dim = c(m, m, n.sim))
for(i in 1:m){
maxx = spec@uspec@spec[[i]]@optimization.model$maxOrder2
htmp = ugarchpath(spec@uspec@spec[[i]], n.sim = n.sim + n.start, n.start = 0, m.sim = 1,
custom.dist = list(name = "sample", distfit = matrix(stdresid[-(1:mo), i], ncol = 1)),
presigma = if( is.null(presigma) ) NA else tail(presigma[,i], maxx),
preresiduals = if( is.null(preresiduals) ) NA else tail(preresiduals[,i], maxx),
prereturns = if( is.null(prereturns) | spec@mspec$varmodel$VAR ) NA else tail(prereturns[,i], maxx),
mexsimdata = if( !spec@mspec$varmodel$VAR ) mexsimdata[[i]] else NULL, vexsimdata = vexsimdata[[i]] )
H[, i] = htmp@path$sigmaSim^2
X[, i] = htmp@path$seriesSim
}
H = tail(H, n.sim)
for(i in 1:n.sim){
simH[ , , i] = diag(sqrt( H[i, ]) ) %*% simR[ , , i] %*% diag(sqrt( H[i, ]) )
}
simX = X
return( list(simH = simH, simX = simX, simR = simR, simQ = simQ) )
}
.rolldcc = function(spec, data, n.ahead = 1, forecast.length = 50, refit.every = 25,
refit.window = c("recursive", "moving"), solver = "solnp",
fit.control = list(eval.se = TRUE),
parallel = FALSE, parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2),
solver.control = list(), save.fit = FALSE, save.wdir = NULL, trace = FALSE, ...)
{
forecast.length = floor(forecast.length/refit.every) * refit.every
WD <- getwd()
if(is.null(save.wdir)){
if (!is.null(WD)) setwd(WD)
} else{
ND = save.wdir
if (!is.null(ND)) setwd(ND)
}
# Setup Indexed Data for Fitting
cnames = colnames(data)
if( is.null(cnames) ) cnames = paste("Asset_", 1:dim(data)[2], sep = "")
xdata = .extractmdata(data)
data = as.matrix(data)
dates = xdata$pos
if(is.null(fit.control$eval.se)) fit.control$eval.se = FALSE
T = dim(data)[1]
startT = T - forecast.length
# forecast window index
fwindex = t(.embed((T-forecast.length - n.ahead + 2):T, refit.every, by = refit.every, TRUE ) )
fitindx.start = c(fwindex[1,] - 1)
fwindex.end = fwindex[refit.every,]
nf = length(fwindex.end)
fitdata = vector(mode = "list", length = nf)
mcfit = vector(mode = "list", length = nf)
#gofit = vector(mode = "list", length = nf)
outdata = vector(mode = "list", length = nf)
distribution = spec@mspec$distribution
forecastlist = cf = lik = vector(mode = "list", length = nf)
n.data = rep(0, nf)
mi = 0
for(i in 1:nf){
if(refit.window[1]=="recursive"){
fitdata = data[1:(fitindx.start[i]+refit.every),]
} else{
fitdata = data[(i*refit.every):(fitindx.start[i]+refit.every),]
}
mcfit = dccfit(spec, data = fitdata, out.sample = refit.every, solver = solver,
fit.control = fit.control, solver.control = solver.control, parallel = parallel,
parallel.control = parallel.control)
cf[[i]] = coef(mcfit)
lik[[i]] = likelihood(mcfit)
n.data[i] = dim(fitdata)[1]
# solver.control = list(tol=1e-6, delta=1e-5)
# mcfit = dccfit(spec, data = fitdata, out.sample = refit.every, solver = solver,
# fit.control = fit.control, solver.control = solver.control)
forecastlist[[i]] = dccforecast(mcfit, n.ahead = n.ahead, n.roll = refit.every - 1,
parallel = parallel, parallel.control = parallel.control)
if(save.fit){
eval(parse(text = paste("mcfitroll",i,"=mcfit",sep = "")))
eval(parse(text = paste("save(mcfitroll",i,",file='mcfitroll",i,".rda')",sep = "")))
}
if(trace) print(i)
mi = mi + 1
}
mspec = list()
mspec$refit.every = refit.every
mspec$forecast.length = forecast.length
mspec$asset.names = cnames
mspec$dccspec = spec
mspec$refit.window = refit.window
mspec$coef = cf
mspec$likelihood = lik
mspec$n.data = n.data
ans = new("DCCroll",
forecast = forecastlist,
spec = mspec)
# return to the current directory
setwd(WD)
return(ans)
}
.dccparidx = function( speclist )
{
m = length(speclist@spec)
type = speclist@type
idx.s = idx.e = vector(mode = "numeric", length = m)
idxi.s = idxi.e = vector(mode = "numeric", length = m)
idxa.s = idxa.e = vector(mode = "numeric", length = m)
pars = lapply( speclist@spec, FUN = function(x) x@optimization.model$modelnames )
for(i in 1:m){
n = length( pars[[i]] )
nx = 0
if( any(substr(pars[[i]], 1, 2) == "mu") ){
nx = c( nx, which( substr(pars[[i]], 1, 2) == "mu" ) )
}
if( any(substr(pars[[i]], 1, 2) == "ar") ){
nx = c( nx, which( substr(pars[[i]], 1, 2) == "ar" ) )
}
if( any(substr(pars[[i]], 1, 2) == "ma") ){
nx = c( nx, which( substr(pars[[i]], 1, 2) == "ma" ) )
}
if( any(substr(pars[[i]], 1, 6) == "inmean") ){
nx = c( nx, which( substr(pars[[i]], 1, 6) == "inmean" ) )
}
if( any(substr(pars[[i]], 1, 7) == "darfima") ){
nx = c( nx, which( substr(pars[[i]], 1, 7) == "darfima" ) )
}
if( any(substr(pars[[i]], 1, 5) == "mxreg") ){
nx = c( nx, which( substr(pars[[i]], 1, 5) == "mxreg" ) )
}
if( i == 1){
idx.s[i] = max(nx) + 1
idx.e[i] = length(pars[[i]])
idxa.s[i] = 1
idxa.e[i] = length(pars[[i]])
} else{
idx.s[i] = idx.e[i - 1] + max(nx) + 1
idx.e[i] = idx.e[i - 1] + n
idxa.s[i] = idxa.e[i - 1] + 1
idxa.e[i] = idxa.e[i - 1] + n
}
idxi.s[i] = max(nx) + 1
idxi.e[i] = n
}
idx = cbind(idx.s, idx.e)
idxi = cbind(idxi.s, idxi.e)
idxa = cbind(idxa.s, idxa.e)
return( list( idx = idx, idxi = idxi, idxa = idxa ) )
}
.dccfilterspec1 = function(spec, dccpars, fitlist)
{
uspec = spec@uspec
m = length(spec@uspec@spec)
newspec = vector(mode = "list", length = m)
for(i in 1:m){
newspec[[i]] = uspec@spec[[i]]
newspec[[i]]@optimization.model$fixed.pars = as.list(coef(fitlist@fit[[i]]))
}
uspec = multispec( newspec )
dspec = dccspec(uspec, dccOrder = spec@mspec$dccOrder,
distribution = spec@mspec$distribution, fixed.pars = as.list(dccpars))
return( dspec )
}
.dccfilterspec2 = function(spec, dccpars, garchpars)
{
uspec = spec@uspec
idx = .dccparidx( uspec )$idxa
m = length(spec)
newspec = vector(mode = "list", length = m)
for(i in 1:m){
newspec[[i]] = uspec[[i]]
gnames = uspec[[i]]@optimization.model$modelnames
gpars = garchpars[idx[i,1]:idx[i,2]]
names(gpars) = gnames
newspec[[i]]@optimization.model$fixed.pars = as.list(gpars)
}
uspec = multispec( newspec )
dspec = dccspec(uspec, dccOrder = spec@mspec$dccOrder,
distribution = spec@mspec$distribution, fixed.pars = as.list(dccpars))
return( dspec )
}
# n-ahead forecast based on the approximation method described in:
# ENGLE, R. and SHEPPARD (2001), K. Theoretical and empirical properties of
# dynamic conditional correlation multivariate garch. NBER Working Papers, No. 8554
.dccforecastm = function(fit, n.ahead = 1, n.roll = 10, external.forecasts = list(mregfor = NULL, vregfor = NULL),
parallel = FALSE, parallel.control = list(pkg = c("multicore", "snowfall"), cores = 2), ...)
{
filterlist = multifilter(multifitORspec = fit@ufit, data = fit@mfit$origdata, out.sample = 0, n.old = NULL,
parallel = parallel, parallel.control = parallel.control)
n.roll = n.roll + 1
m = length(fit@ufit@fit)
out.sample = fit@mfit$out.sample
pos = fit@mfit$model$pos.matrix
# we only allow n.ahead = 1 when using n.roll (for simlicity and avoidance of arrays
fitlist = fit@ufit
mo = max(fit@mfit$model$dccOrder)
# we augmented n.roll before, now subtract
# forclist = multiforecast(multifitORspec = fitlist, n.ahead = n.ahead, n.roll = n.roll - 1)
forclist = multiforecast(multifitORspec = fitlist, n.ahead = n.ahead, n.roll = n.roll - 1,
external.forecasts = external.forecasts, parallel = parallel,
parallel.control = parallel.control, ...)
sig = sigma(filterlist)
resid = residuals(filterlist)
stdresid = resid/sig
T = dim(fit@mfit$H)[3]
pars = coef(fit, type = "dcc")
dcca = dccb = skew = shape = 0
if(pos[1,3] == 1) dcca = pars[pos[1,1]:pos[1,2]]
if(pos[2,3] == 1) dccb = pars[pos[2,1]:pos[2,2]]
## call .nsdccforecast with rolling window
Rbar = Rtfor = Htfor = Qtfor = vector(mode = "list", length = n.roll)
Qstart = last( rcor(fit, type = "Q"), mo )
Rstart = last( rcor(fit, type = "R"), mo )
Hstart = last( rcov(fit), mo)
for(i in 1:n.roll){
xQbar = cov(stdresid[1:(T + i - 1), ])
xstdresids = stdresid[(T - mo + i ):(T + i - 1), , drop = FALSE]
xfsig = matrix( sapply( forclist@forecast, FUN = function(x) unlist(x@forecast$forecasts[[i]]['sigma']) ), ncol = m)
ans = .dccforecastn(Qstart, Rstart, Hstart, xQbar, xstdresids, xfsig, dcca, dccb, n.ahead, mo)
Rtfor[[i]] = ans$Rtfor
Qtfor[[i]] = ans$Qtfor
Htfor[[i]] = ans$Htfor
Rbar[[i]] = ans$Rbar
Qstart = last( .abind(Qstart, ans$Qtfor[, , 1]), mo )
Rstart = last( .abind(Rstart, ans$Rtfor[, , 1]), mo )
Hstart = last( .abind(Hstart, ans$Htfor[, , 1]), mo )
}
forc = list( H = Htfor, R = Rtfor, Q = Qtfor, Rbar = Rbar, ufor = forclist )
return(forc)
}
.dccforecastn = function(Qstart, Rstart, Hstart, Qbar, stdresids, fsig, dcca, dccb, n.ahead, mo)
{
m = dim(Qbar)[1]
Qtfor = Rtfor = Htfor = array(NA, dim = c(m, m, n.ahead + mo))
Qtfor[ , , 1:mo] = Qstart[, , 1:mo]
Rtfor[ , , 1:mo] = Rstart[, , 1:mo]
Htfor[ , , 1:mo] = Hstart[, , 1:mo]
dccsum = sum(dcca) + sum(dccb)
Qt_1 = (1 - dccsum) * Qbar
for(i in 1:n.ahead){
Qtfor[, , mo + i] = Qt_1
if( i == 1 ){
for(j in 1:length(dcca)){
Qtfor[ , , mo + 1] = Qtfor[ , , mo + 1] + dcca[j] * (stdresids[(mo + 1 - j), ] %*% t(stdresids[(mo + 1 - j), ]))
}
for(j in 1:length(dccb)){
Qtfor[ , , mo + 1] = Qtfor[ , , mo + 1] + dccb[j] * Qtfor[ , , mo + 1 - j]
}
Qtmp = diag( 1/sqrt( diag(Qtfor[ , , mo + 1]) ) , m, m)
Rtfor[ , , mo + 1] = Qtmp %*% Qtfor[ , , mo + 1] %*% t(Qtmp)
Dtmp = diag(fsig[1, ], m, m)
Htfor[ , , mo + 1] = Dtmp %*% Rtfor[ , , mo + 1] %*% Dtmp
# now the unconditional calculations
Qt_1star = diag( 1/sqrt( diag(Qtfor[, , mo + 1]) ) , m, m)
ER_1 = Qt_1star %*% Qtfor[, , mo + 1] %*% t(Qt_1star)
Qbarstar = diag( 1/sqrt( diag(Qbar) ) , m, m)
Rbar = Qbarstar %*% Qbar %*% t(Qbarstar)
} else{
Rtfor[, , mo + i] = (1 - dccsum^(i - 1) ) * Rbar + dccsum^(i - 1) * ER_1
Dtmp = diag(fsig[i, ], m, m)
Htfor[, , mo + i] = Dtmp %*% Rtfor[, , mo + i] %*% Dtmp
Qtfor[, , mo + i] = Qtfor[, , mo + 1]
}
}
return( list( Rtfor = Rtfor[, , -(1:mo), drop = FALSE], Htfor = Htfor[, , -(1:mo), drop = FALSE],
Qtfor = Qtfor[, , -(1:mo), drop = FALSE], Rbar = Rbar, ER_1 = ER_1) )
}
.dccsimf = function(Z, Qbar, Rbar, dcca, dccb, mo, n.sim, n.start, m, rseed)
{
n = n.sim + n.start + mo
set.seed(rseed[1]+1)
stdresid = matrix(rnorm(m * (n.sim+n.start+mo)), nrow = n.sim + n.start + mo, ncol = m)
sumdcca = sum(dcca)
sumdccb = sum(dccb)
dccsumx = sumdcca + sumdccb
ldcca = length(dcca)
ldccb = length(dccb)
# ToDo: Assumes that no order is zero!
dccorderx = c(ldcca, ldccb, mo, n)
res = .Call(name = "dccsim", Qbar = as.matrix(Qbar), Rbar = as.matrix(Rbar),
Z = as.matrix(Z), Res = as.matrix(stdresid), dcca = as.numeric(dcca),
dccb = as.numeric(dccb), dccorder = as.numeric(dccorderx), dccsum = as.numeric(dccsumx))
Q = array(NA, dim = c(m, m, n.sim + n.start + mo))
R = array(NA, dim = c(m, m, n.sim + n.start + mo))
for(i in 1:(n.sim + n.start + mo)){
R[,,i] = res[[2]][[i]]
Q[,,i] = res[[3]][[i]]
}
ans = list( Q = Q, R = R, stdresid = res[[3]])
return( ans )
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.