Nothing
R/fdcc-main.R
In rmgarch: Multivariate GARCH Models
Defines functions
.fdcccon
.fullinc_fdcc
.fdccsimf
.fdccsim.spec
.fdccsim.fit
.fdccforecastn2
.fdccforecastn
.fdccforecastm
.fdccforecast
.fdccfilter
getfdccpars
.fdccfit
.fdccspec
#################################################################################
##
## R package rmgarch by Alexios Galanos Copyright (C) 2008-2022.
## This file is part of the R package rmgarch.
##
## The R package rmgarch 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 rmgarch 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.
##
#################################################################################
.fdccspec = function(uspec, VAR = FALSE, 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),
fdccOrder=c(1,1), fdccindex = c(1,1,1,2,2,2,3,3,3),
distribution = c("mvnorm", "mvt", "mvlaplace"),
start.pars = list(), fixed.pars = list())
{
# VAR checks
VAR.opt = list()
if(is.null(VAR)) VAR.opt$VAR = FALSE else VAR.opt$VAR = as.logical(VAR)
if(is.null(robust)) VAR.opt$robust = FALSE else VAR.opt$robust = as.logical(robust)
if(is.null(lag)) VAR.opt$lag = 1 else VAR.opt$lag = as.integer(lag)
if(is.null(lag.max)) VAR.opt$lag.max = NULL else VAR.opt$lag.max = as.integer(min(1, lag.max))
if(is.null(lag.criterion)) VAR.opt$lag.criterion = "AIC" else VAR.opt$lag.criterion = lag.criterion[1]
if(is.null(external.regressors)) VAR.opt$external.regressors = NULL else VAR.opt$external.regressors = external.regressors
rc = list("gamma" = 0.25, "delta" = 0.01, "nc" = 10, "ns" = 500)
rcmatch = match(names(robust.control), c("gamma", "delta", "nc", "ns"))
if(length(rcmatch[!is.na(rcmatch)]) > 0){
rx = which(!is.na(rcmatch))
rc[rcmatch[!is.na(rcmatch)]] = robust.control[rx]
}
VAR.opt$robust.control = rc
.eps = .Machine$double.eps
modeldata = list()
modeldesc = list()
m = length(uspec@spec)
## distribution checks
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)
## fdcc checks
fdccOrder = as.integer(abs(fdccOrder))
if(max(fdccOrder)>1) stop("\nFDCC model only allows (1,1) DCC dynamics and not higher.")
fdccindex = as.integer(abs(fdccindex))
if(any(fdccindex)==0) stop("\nfdccindex cannot contain zeros...index starts at 1!")
lenF = length(unique(fdccindex))
if(length(fdccindex)!=length(uspec@spec)) stop("\nfdccindex length not equal to length of uspec.")
# allow for case of just 1 index
if(!all(fdccindex==1) && any(diff(unique(fdccindex)))!=1) stop("\nfdccindex must use contiguous indexing, starting from one and incrementing by 1.")
modelinc = rep(0, 10)
names(modelinc) = c("var", "mvmxreg", "fdccC", "fdccA", "fdccB", "mshape", "mskew", "aux", "aux", "aux")
if(distribution == "mvt") modelinc[6] = 1
modelinc[3] = lenF
modelinc[4] = fdccOrder[1]
modelinc[5] = fdccOrder[2]
# indicator matrix of groups
sidx = matrix(0, ncol = lenF, nrow = length(fdccindex))
for(i in 1:lenF) sidx[fdccindex==i,i] = 1
# sidx %*% a
# ai_j...i=sector, j = order
# bi_j...
# ci
if( VAR ){
if(is.null(VAR.opt$lag)) modelinc[1] = 1 else modelinc[1] = as.integer( VAR.opt$lag )
if(!is.null(VAR.opt$external.regressors)){
if(!is.matrix(VAR.opt$external.regressors)) stop("\nexternal.regressors must be a matrix.")
modelinc[2] = dim(VAR.opt$external.regressors)[2]
modeldata$mexdata = VAR.opt$external.regressors
} else{
modeldata$mexdata = NULL
}
}
modelinc[10] = which(c("mvnorm", "mvt", "mvlaplace") == distribution)
modeldesc$distribution = distribution
if( !is(uspec, "uGARCHmultispec") ) stop("\ndccspec-->error: uspec must be a uGARCHmultispec object")
varmodel = list()
umodel = vector(mode ="list")
if( modelinc[1]>0 ){
varmodel$robust = VAR.opt$robust
varmodel$lag.max = VAR.opt$lag.max
varmodel$lag.criterion = VAR.opt$lag.criterion
varmodel$robust.control = VAR.opt$robust.control
umodel$modelinc = matrix(0, ncol = m, nrow = 21)
rownames(umodel$modelinc) = names(uspec@spec[[1]]@model$modelinc[1:21])
umodel$modeldesc = list()
umodel$vt = sapply(uspec@spec, FUN = function(x) x@model$modelinc[22])
umodel$modeldesc$vmodel = vector(mode = "character", length = m)
umodel$modeldesc$vsubmodel = vector(mode = "character", length = m)
umodel$start.pars = umodel$fixed.pars = vector(mode = "list", length = m)
umodel$modeldesc$distribution = vector(mode = "character", length = m)
umodel$modeldata = list()
umodel$modeldata$vexdata = vector(mode = "list", length = m)
for(i in 1:m){
# zero the mean equation since we are using VAR
umodel$modeldesc$vmodel[i] = uspec@spec[[i]]@model$modeldesc$vmodel
umodel$modeldesc$vsubmodel[i] = ifelse(is.null(uspec@spec[[i]]@model$modeldesc$vsubmodel),"GARCH",uspec@spec[[i]]@model$modeldesc$vsubmodel)
umodel$modeldesc$distribution[i] = uspec@spec[[i]]@model$modeldesc$distribution
umodel$modelinc[,i] = uspec@spec[[i]]@model$modelinc[1:21]
umodel$modelinc[1:6,i] = 0
umodel$modeldata$vexdata[[i]] = if(is.null(uspec@spec[[i]]@model$modeldata$vexdata)) NA else uspec@spec[[i]]@model$modeldata$vexdata
umodel$start.pars[[i]] = if(is.null(uspec@spec[[i]]@model$start.pars)) NA else uspec@spec[[i]]@model$start.pars
umodel$fixed.pars[[i]] = if(is.null(uspec@spec[[i]]@model$fixed.pars)) NA else uspec@spec[[i]]@model$fixed.pars
umodel$modeldata$mexdata[[i]] = NA
}
} else{
varmodel$lag.max = 1
varmodel$lag.criterion = "HQ"
umodel$modelinc = matrix(0, ncol = m, nrow = 21)
rownames(umodel$modelinc) = names(uspec@spec[[1]]@model$modelinc[1:21])
umodel$modeldesc = list()
umodel$vt = sapply(uspec@spec, FUN = function(x) x@model$modelinc[22])
umodel$modeldesc$vmodel = vector(mode = "character", length = m)
umodel$modeldesc$vsubmodel = vector(mode = "character", length = m)
umodel$start.pars = umodel$fixed.pars = vector(mode = "list", length = m)
umodel$modeldesc$distribution = vector(mode = "character", length = m)
umodel$modeldata = list()
umodel$modeldata$mexdata = vector(mode = "list", length = m)
umodel$modeldata$vexdata = vector(mode = "list", length = m)
for(i in 1:m){
umodel$modeldesc$vmodel[i] = uspec@spec[[i]]@model$modeldesc$vmodel
umodel$modeldesc$vsubmodel[i] = ifelse(is.null(uspec@spec[[i]]@model$modeldesc$vsubmodel),"GARCH",uspec@spec[[i]]@model$modeldesc$vsubmodel)
umodel$modeldesc$distribution[i] = uspec@spec[[i]]@model$modeldesc$distribution
umodel$modelinc[,i] = uspec@spec[[i]]@model$modelinc[1:21]
umodel$modeldata$mexdata[[i]] = if(is.null(uspec@spec[[i]]@model$modeldata$mexdata)) NA else uspec@spec[[i]]@model$modeldata$mexdata
umodel$modeldata$vexdata[[i]] = if(is.null(uspec@spec[[i]]@model$modeldata$vexdata)) NA else uspec@spec[[i]]@model$modeldata$vexdata
umodel$start.pars[[i]] = if(is.null(uspec@spec[[i]]@model$start.pars)) NA else uspec@spec[[i]]@model$start.pars
umodel$fixed.pars[[i]] = if(is.null(uspec@spec[[i]]@model$fixed.pars)) NA else uspec@spec[[i]]@model$fixed.pars
}
}
maxgarchOrder = max( sapply(uspec@spec, FUN = function(x) x@model$maxOrder) )
# Note: We run varxfit with the postpad option using "constant" so that
# it resembles the ARMA-GARCH type functionality
if(modelinc[1]>0){
maxgarchOrder = max(c(maxgarchOrder, modelinc[1]))
}
################################################################
# pars list
pos = 1
pos.matrix = matrix(0, ncol = 3, nrow = 4)
colnames(pos.matrix) = c("start", "stop", "include")
rownames(pos.matrix) = c("fdcca", "fdccb", "mshape", "mskew")
# c("var", "mvmxreg", "fdccC", "fdccA", "fdccB", "mshape", "mskew", "aux", "aux", "aux")
if(modelinc[4]>0) pos.matrix[1,1:3] = c(pos, pos+(modelinc[4]*modelinc[3])-1, 1)
pos = pos + (modelinc[4]*modelinc[3])
if(modelinc[5]>0) pos.matrix[2,1:3] = c(pos, pos+(modelinc[5]*modelinc[3])-1, 1)
pos = pos + (modelinc[5]*modelinc[3])
if(modelinc[6]>0) pos.matrix[3,1:3] = c(pos, pos, 1)
pos = pos + modelinc[6]
if(modelinc[7]>0) pos.matrix[4,1:3] = c(pos, pos, 1)
# minimum parameter matrix = modelinc[3]+dcc1+dccb+mshape+mskew = modelinc[3]+4
mm = max(1, (modelinc[4]*modelinc[3]))+max(1, (modelinc[5]*modelinc[3])) + max(1, modelinc[6]) + max(1, modelinc[6])
pars = matrix(0, ncol = 6, nrow = mm)
colnames(pars) = c("Level", "Fixed", "Include", "Estimate", "LB", "UB")
pidx = matrix(NA, nrow = 4, ncol = 2)
colnames(pidx) = c("begin", "end")
rownames(pidx) = c("fdcca", "fdccb", "mshape", "mskew")
fixed.names = names(fixed.pars)
start.names = names(start.pars)
fixed.pars = unlist(fixed.pars)
start.pars = unlist(start.pars)
pn = 1
pnames = NULL
nx = 0
pidx[1,1] = 1
if(pos.matrix[1,3] == 1){
pn = length( seq(pos.matrix[1,1], pos.matrix[1,2], by = 1) )
xid = 0
for(i in 1:modelinc[4]){
for(j in 1:modelinc[3]){
nnx = paste("fdcca[g=", j,",P=", i, "]", sep="")
pars[(nx+xid+j), 1] = 0.01
sp = na.omit(match(nnx, start.names))
if(length(sp)==1) pars[(nx+xid+j), 1] = start.pars[sp]
pars[(nx+xid+j), 3] = 1
pars[(nx+xid+j), 5] = 0
pars[(nx+xid+j), 6] = 0.999
sp = na.omit(match(nnx, fixed.names))
if(length(sp)==1){
pars[(nx+xid+j), 1] = fixed.pars[sp]
pars[(nx+xid+j), 2] = 1
} else{
pars[(nx+xid+j), 4] = 1
}
pnames = c(pnames, nnx)
}
xid = xid + modelinc[3]
}
} else{
pnames = c(pnames, "fdcca")
}
pidx[1,2] = nx + pn
nx = nx + pn
pn = 1
pidx[2,1] = nx+1
if(pos.matrix[2,3] == 1){
pn = length( seq(pos.matrix[2,1], pos.matrix[2,2], by = 1) )
xid = 0
for(i in 1:modelinc[5]){
for(j in 1:modelinc[3]){
nnx = paste("fdccb[g=", j,",Q=", i, "]", sep="")
pars[(nx+xid+j), 1] = 0.9
sp = na.omit(match(nnx, start.names))
if(length(sp)==1) pars[(nx+xid+j), 1] = start.pars[sp]
pars[(nx+xid+j), 3] = 1
pars[(nx+xid+j), 5] = 0
pars[(nx+xid+j), 6] = 1
sp = na.omit(match(nnx, fixed.names))
if(length(sp)==1){
pars[(nx+xid+j), 1] = fixed.pars[sp]
pars[(nx+xid+j), 2] = 1
} else{
pars[(nx+xid+j), 4] = 1
}
pnames = c(pnames, nnx)
}
xid = xid+modelinc[3]
}
} else{
pnames = c(pnames, "fdccb")
}
pidx[2,2] = nx+pn
nx = nx + pn
pn = 1
pidx[3,1] = nx+1
if(modelinc[6]<=1){
if(pos.matrix[3,3]==1){
pars[nx+pn, 3] = 1
pars[nx+pn, 1] = 5
pars[(nx+pn), 5] = 4
pars[(nx+pn), 6] = 50
if(any(substr(start.names, 1, 6) == "mshape")) pars[nx+pn, 1] = start.pars["mshape"]
if(any(substr(fixed.names, 1, 6) == "mshape")){
pars[nx+pn,2] = 1
pars[nx+pn, 1] = fixed.pars["mshape"]
} else{
pars[nx+pn,4] = 1
}
}
pnames = c(pnames, "mshape")
} else{
if(pos.matrix[3,3] == 1){
pn = length( seq(pos.matrix[3,1], pos.matrix[3,2], by = 1) )
for(i in 1:pn){
nnx = paste("mshape", i, sep="")
pars[(nx+i), 1] = 5
if(any(substr(start.names, 1, nchar(nnx))==nnx)){
nix = which(substr(start.names, 1, nchar(nnx)) == nnx)
pars[(nx+i), 1] = start.pars[nix]
}
pars[(nx+i), 3] = 1
pars[(nx+i), 5] = 4
pars[(nx+i), 6] = 50
if(any(substr(fixed.names, 1, nchar(nnx))==nnx)){
nix = which(substr(fixed.names, 1, nchar(nnx)) == nnx)
pars[(nx+i), 1] = fixed.pars[nix]
pars[(nx+i), 2] = 1
} else{
pars[(nx+i), 4] = 1
}
pnames = c(pnames, nnx)
}
} else{
pnames = c(pnames, "mshape")
}
}
pidx[3,2] = nx+pn
nx = nx + pn
pn = 1
pidx[4,1] = nx+1
if(modelinc[7]<=1){
if(pos.matrix[4,3]==1){
pars[nx+pn, 3] = 1
pars[nx+pn, 1] = 0.5
pars[(nx+pn), 5] = -1
pars[(nx+pn), 6] = 1
if(any(substr(start.names, 1, 5) == "mskew")) pars[nx+pn, 1] = start.pars["mskew"]
if(any(substr(fixed.names, 1, 5) == "mskew")){
pars[nx+pn,2] = 1
pars[nx+pn, 1] = fixed.pars["mskew"]
} else{
pars[nx+pn,4] = 1
}
}
pnames = c(pnames, "mskew")
} else{
if(pos.matrix[4,3] == 1){
pn = length( seq(pos.matrix[4,1], pos.matrix[4,2], by = 1) )
for(i in 1:pn){
nnx = paste("mskew", i, sep="")
pars[(nx+i), 1] = 3
pars[(nx+i), 5] = -1
pars[(nx+i), 6] = 1
if(any(substr(start.names, 1, nchar(nnx))==nnx)){
nix = which(substr(start.names, 1, nchar(nnx)) == nnx)
pars[(nx+i), 1] = start.pars[nix]
}
if(any(substr(fixed.names, 1, nchar(nnx))==nnx)){
nix = which(substr(fixed.names, 1, nchar(nnx)) == nnx)
pars[(nx+i), 1] = fixed.pars[nix]
pars[(nx+i), 2] = 1
} else{
pars[(nx+i), 4] = 1
}
pnames = c(pnames, nnx)
}
} else{
pnames = c(pnames, "mskew")
}
}
pidx[4,2] = nx+pn
rownames(pars) = pnames
modeldesc$type = "2-step"
model = list(modelinc = modelinc, modeldesc = modeldesc, modeldata = modeldata, varmodel = varmodel,
pars = pars, start.pars = start.pars, fixed.pars = fixed.pars, maxgarchOrder = maxgarchOrder,
sidx = sidx, fdccindex = fdccindex, pos.matrix = pos.matrix, pidx = pidx)
model$DCC = "FDCC"
ans = new("DCCspec",
model = model,
umodel = umodel)
return(ans)
}
.fdccfit = function(spec, data, out.sample = 0, solver = "solnp",
solver.control = list(), fit.control = list(eval.se = TRUE,
stationarity = TRUE, scale = FALSE),
cluster = NULL, fit = NULL, VAR.fit = NULL, verbose = FALSE,
realizedVol = NULL, ...)
{
tic = Sys.time()
.eps = .Machine$double.eps
model = spec@model
umodel = spec@umodel
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
# allow first and second stage solvers to be defined
if(length(solver)==2){
garch.solver = solver[1]
solver = solver[2]
} else{
garch.solver = solver[1]
}
solver = match.arg(tolower(solver)[1], c("solnp", "nlminb", "lbfgs","gosolnp"))
#-----------------------------------------------------------------------------------
# Data Extraction
m = dim(data)[2]
if( is.null( colnames(data) ) ) cnames = paste("Asset_", 1:m, sep = "") else cnames = colnames(data)
colnames(umodel$modelinc) = cnames
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]
if( (n-n.start) < 100)
stop("\ndccfit-->error: function requires at least 100 data\n points to run\n")
data = xdata$data
index = xdata$index
period = xdata$period
# save the data to the model spec
model$modeldata$data = data
model$modeldata$index = index
model$modeldata$period = period
T = model$modeldata$T = n - n.start
model$modeldata$n.start = n.start
model$modeldata$asset.names = cnames
#-----------------------------------------------------------------------------------
# VAR model
if( model$modelinc[1]>0 ){
tmp = mvmean.varfit(model = model, data = data, VAR.fit = VAR.fit, T = T,
out.sample = out.sample, cluster = cluster)
model = tmp$model
zdata = tmp$zdata
mu = tmp$mu
varcoef = tmp$varcoef
p = tmp$p
N = tmp$N
} else{
zdata = data
ex = NULL
}
T = dim(zdata)[1] - out.sample
#-----------------------------------------------------------------------------------
# Univariate GARCH fit
# create a multispec list
mspec = .makemultispec(umodel$modelinc, umodel$modeldesc$vmodel, umodel$modeldesc$vsubmodel,
umodel$modeldata$mexdata, umodel$modeldata$vexdata, umodel$start.pars,
umodel$fixed.pars, umodel$vt)
if( !is.null(fit) && is(fit, "uGARCHmultifit") ){
# check VAR and fit:
if(model$modelinc[1]>0){
for(i in 1:m){
if(sum(fit@fit[[i]]@model$modelinc[1:6])>0)
stop("\nThe user supplied fit object has a non-null mean specification but VAR already chosen for mean filtration!!!")
}
}
fitlist = fit
if(spec@model$modelinc[1]>0) model$mu = mu else model$mu = fitted(fitlist)
model$residuals = res = residuals(fitlist)
model$sigma = sig = sigma(fitlist)
} else{
fitlist = multifit(multispec = mspec, data = xts(zdata, index), out.sample = n.start,
solver = garch.solver, solver.control = solver.control,
fit.control = ufit.control, cluster = cluster, realizedVol = realizedVol)
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 )
}
if(spec@model$modelinc[1]>0) model$mu = mu else model$mu = fitted(fitlist)
model$residuals = res = residuals(fitlist)
model$sigma = sig = sigma(fitlist)
}
stdresid = res/sig
Qbar = (t(stdresid) %*% stdresid)/nrow(stdresid)
# In the VAR model we postpadded the first lag values with a constant.
# In an ARMA model, the first lag values are also constant.
#-----------------------------------------------------------------------------------
# Create the Full Model Pars and Indices
modelinc = model$modelinc
# create full par matrix
tmp = .fullinc_fdcc(modelinc, umodel)
midx = tmp$midx
fdccAidx = tmp$fdccAidx
fdccBidx = tmp$fdccBidx
mpars = midx*0
# This is the estimated parameter index
eidx = .estindfn(midx, mspec, model$pars)
unipars = lapply(fitlist@fit, FUN = function(x) x@fit$ipars[x@fit$ipars[,3]==1,1])
if(is.list(unipars)){
for(i in 1:length(unipars)){
uninames = names(unipars[[i]])
mpars[uninames, i] = unipars[[i]]
}
} else{
uninames = rownames(unipars)
mpars[uninames, 1:NCOL(unipars)] = unipars
}
# add include pars from DCC spec (includes the fixed pars)
mpars[which(midx[,m+1]==1), m+1] = as.numeric( model$pars[model$pars[,3]==1,1] )
# DCC parameters
ipars = model$pars
LB = ipars[,5]
UB = ipars[,6]
estidx = as.logical( ipars[,4] )
npars = sum(estidx)
#-----------------------------------------------------------------------------------
H = sig^2
#-----------------------------------------------------------------------------------
# create a temporary environment to store values (deleted at end of function)
mgarchenv = new.env(hash = TRUE)
arglist = list()
arglist$mgarchenv = mgarchenv
arglist$verbose = verbose
arglist$cluster = cluster
arglist$eval.se = fit.control$eval.se
arglist$solver = solver
arglist$fit.control = fit.control
arglist$cnames = cnames
arglist$m = m
arglist$T = T
arglist$data = zdata
arglist$index = index
arglist$realizedVol = realizedVol
arglist$model = model
arglist$fitlist = fitlist
arglist$umodel = umodel
arglist$midx = midx
arglist$eidx = eidx
arglist$mpars = mpars
arglist$ipars = ipars
arglist$estidx = estidx
arglist$dccN = npars
arglist$stdresid = stdresid
arglist$H = H
arglist$Qbar = Qbar
# FDCC specific
arglist$fdccAidx = fdccAidx
arglist$fdccBidx = fdccBidx
arglist$sidx = model$sidx
#-----------------------------------------------------------------------------------
# Check for fixed parameters in DCC
if(any(ipars[,2]==1)){
if(npars == 0){
if(fit.control$eval.se==0) {
warning("\nfdccfit-->warning: all parameters fixed...returning dccfilter object instead\n")
xspex = spec
for(i in 1:m) xspex@umodel$fixed.pars[[i]] = as.list(fitlist@fit[[i]]@model$pars[fitlist@fit[[i]]@model$pars[,3]==1,1])
return(.fdccfilter(spec = xspex, data = xts(data, index), out.sample = out.sample, cluster = cluster, VAR.fit = VAR.fit,
realizedVol = realizedVol))
} else{
# if all parameters are fixed but we require standard errors, we
# skip the solver
use.solver = 0
ipars[ipars[,2]==1, 4] = 1
ipars[ipars[,2]==1, 2] = 0
arglist$pars = ipars
estidx = as.logical( ipars[,4] )
arglist$estidx = estidx
}
} else{
# with some parameters fixed we extract them (to be rejoined at end)
# so that they do not enter the solver
use.solver = 1
}
} else{
use.solver = 1
}
# start counter
assign("rmgarch_llh", 1, envir = mgarchenv)
# get
# calculate constraint length:
if(modelinc[3]==1){
ILB = 0
IUB = 1
} else{
if(modelinc[4]>0) xx = matrix(ipars[arglist$model$pidx["fdcca", 1]:arglist$model$pidx["fdcca", 2], 1], ncol = modelinc[4]) else xx = matrix(0, nrow = modelinc[3], ncol = modelinc[5])
if(modelinc[5]>0) yy = matrix(ipars[arglist$model$pidx["fdccb", 1]:arglist$model$pidx["fdccb", 2], 1], ncol = modelinc[5]) else yy = matrix(0, nrow = modelinc[3], ncol = modelinc[4])
# each pairwise combination
nc = NROW(expand.grid(xx, xx, yy, yy))
ILB = rep(0, nc)
IUB = rep(1, nc)
}
# Asymmetric Spec has different constraints
Ifn = .fdcccon
if( solver == "solnp" | solver == "gosolnp") fit.control$stationarity = FALSE else fit.control$stationarity = TRUE
arglist$fit.control = fit.control
if( use.solver )
{
arglist$returnType = "llh"
solution = switch(model$modeldesc$distribution,
mvnorm = .dccsolver(solver, pars = ipars[estidx, 1], fun = normal.fdccLLH1, Ifn, ILB,
IUB, gr = NULL, hessian = NULL, control = solver.control,
LB = ipars[estidx, 5], UB = ipars[estidx, 6], arglist = arglist))
#mvlaplace = .dccsolver(solver, pars = ipars[estidx, 1], fun = laplace.fdccLLH1, Ifn, ILB,
# IUB, gr = NULL, hessian = NULL, control = solver.control,
# LB = ipars[estidx, 5], UB = ipars[estidx, 6], arglist = arglist),
#mvt = .dccsolver(solver, pars = ipars[estidx, 1], fun = student.fdccLLH1, Ifn, ILB,
# IUB, gr = NULL, hessian = NULL, control = solver.control,
# LB = ipars[estidx, 5], UB = ipars[estidx, 6], arglist = arglist))
sol = solution$sol
hess = solution$hess
timer = Sys.time()-tic
convergence = sol$convergence
mpars[which(eidx[,(m+1)]==1, arr.ind = TRUE),m+1] = sol$pars
ipars[estidx, 1] = sol$pars
arglist$mpars = mpars
arglist$ipars = ipars
} else{
hess = NULL
timer = Sys.time()-tic
convergence = 0
sol = list()
sol$message = "all parameters fixed"
}
fit = list()
# check convergence else write message/return
if( convergence == 0 ){
fit = switch(model$modeldesc$distribution,
mvnorm = .dccmakefitmodel(garchmodel = "fdccnorm", f = normal.fdccLLH2,
arglist = arglist, timer = 0, message = sol$message, fname = "normal.fdccLLH2"))
#mvlaplace = .dccmakefitmodel(garchmodel = "fdcclaplace", f = laplace.fdccLLH2,
# arglist = arglist, timer = 0, message = sol$message, fname = "laplace.fdccLLH2"),
#mvt =.dccmakefitmodel(garchmodel = "fdccstudent", f = student.fdccLLH2,
# arglist = arglist, timer = 0, message = sol$message, fname = "student.fdccLLH2"))
fit$timer = Sys.time() - tic
} else{
fit$message = sol$message
fit$convergence = 1
}
# make model list to return some usefule information which
model$mpars = mpars
model$ipars = ipars
model$pars[,1] = ipars[,1]
model$midx = midx
model$eidx = eidx
model$umodel = umodel
fit$Qbar = Qbar
fit$realizedVol = realizedVol
ans = new("DCCfit",
mfit = fit,
model = model)
return(ans)
}
getfdccpars = function(pars, model){
modelinc = model$modelinc
pidx = model$pidx
sidx = model$sidx
if(modelinc[4]>0){
A = matrix(pars[pidx["fdcca",1]:pidx["fdcca",2],1], ncol = modelinc[4], byrow=FALSE)
A = apply(A, 2, function(x) sidx %*% x)
} else{
A = matrix(0)
}
if(modelinc[5]>0){
B = matrix(pars[pidx["fdccb",1]:pidx["fdccb",2],1], ncol = modelinc[5], byrow=FALSE)
B = apply(B, 2, function(x) sidx %*% x)
} else{
B = matrix(0)
}
# Because only the DCC (1,1) model is allowed, this is simplified:
BB = (B %*% t(B))
AA = (A %*% t(A))
C = matrix(1, NROW(BB), NCOL(BB)) - BB - AA
return(list(A = A, B = B, C = C))
}
.fdccfilter = function(spec, data, out.sample = 0, filter.control = list(n.old = NULL),
cluster = NULL, varcoef = NULL, realizedVol = NULL, ...)
{
tic = Sys.time()
model = spec@model
umodel = spec@umodel
n.old = filter.control$n.old
#-----------------------------------------------------------------------------------
# Data Extraction
m = dim(data)[2]
if( is.null( colnames(data) ) ) cnames = paste("Asset_", 1:m, sep = "") else cnames = colnames(data)
colnames(umodel$modelinc) = cnames
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]
if( (n-n.start) < 100)
stop("\ndccfilter-->error: function requires at least 100 data\n points to run\n")
data = xdata$data
index = xdata$index
period = xdata$period
# save the data to the model spec
model$modeldata$data = data
model$modeldata$index = index
model$modeldata$period = period
T = model$modeldata$T = n - n.start
model$modeldata$n.start = n.start
model$modeldata$asset.names = cnames
#-----------------------------------------------------------------------------------
# VAR model
if( spec@model$modelinc[1]>0 ){
tmp = mvmean.varfilter(model = model, data = data, varcoef = varcoef,
T = T, out.sample = out.sample)
model = tmp$model
zdata = tmp$zdata
mu = tmp$mu
p = tmp$p
N = tmp$N
} else{
zdata = data
ex = NULL
}
T = dim(zdata)[1] - out.sample
#-----------------------------------------------------------------------------------
if(is.null(filter.control$n.old)) n.old = T
#-----------------------------------------------------------------------------------
# Univariate GARCH filter
if(model$modelinc[1]>0){
for(i in 1:m){
if(sum(umodel$modelinc[1:6,i])>0)
stop("\nThe user supplied univariate spec object has a non-null mean specification but VAR already chosen for mean filtration!!!")
}
}
mspec = .makemultispec(umodel$modelinc, umodel$modeldesc$vmodel, umodel$modeldesc$vsubmodel,
umodel$modeldata$mexdata, umodel$modeldata$vexdata, umodel$start.pars,
umodel$fixed.pars, NULL)
filterlist = multifilter(multifitORspec = mspec, data = xts(zdata, index[1:nrow(zdata)]),
out.sample = out.sample, cluster = cluster, n.old = n.old,
realizedVol = realizedVol)
if(spec@model$modelinc[1]>0) model$mu = mu else model$mu = fitted(filterlist)
model$residuals = res = residuals(filterlist)
model$sigma = sig = sigma(filterlist)
stdresid = res/sig
if(is.null(filter.control$n.old)) dcc.old = dim(stdresid)[1] else dcc.old = n.old
#-----------------------------------------------------------------------------------
# Create the Full Model Pars and Indices
modelinc = model$modelinc
# create full par matrix
tmp = .fullinc_fdcc(modelinc, umodel)
midx = tmp$midx
fdccAidx = tmp$fdccAidx
fdccBidx = tmp$fdccBidx
mpars = midx*0
eidx = midx
# This is the estimated parameter index
unipars = sapply(filterlist@filter, FUN = function(x) x@filter$ipars[x@filter$ipars[,3]==1,1])
if(is.list(unipars)){
for(i in 1:length(unipars)){
uninames = names(unipars[[i]])
mpars[uninames, i] = unipars[[i]]
}
} else{
uninames = rownames(unipars)
mpars[uninames, 1:NCOL(unipars)] = unipars
}
# add include pars from DCC spec (includes the fixed pars)
mpars[which(midx[,m+1]==1, arr.ind = TRUE), m+1] = as.numeric( model$pars[model$pars[,3]==1,1] )
# DCC parameters
ipars = model$pars
estidx = as.logical( ipars[,3] )
npars = sum(estidx)
#-----------------------------------------------------------------------------------
Qbar = cov(stdresid[1:dcc.old, ])
H = sig^2
#-----------------------------------------------------------------------------------
mgarchenv = new.env(hash = TRUE)
arglist = list()
arglist$mgarchenv = mgarchenv
arglist$verbose = FALSE
arglist$cluster = cluster
arglist$filter.control = filter.control
arglist$cnames = cnames
arglist$m = m
arglist$T = T
arglist$n.old = n.old
arglist$dcc.old = dcc.old
arglist$data = zdata
arglist$index = index
arglist$realizedVol = realizedVol
arglist$model = model
arglist$filterlist = filterlist
arglist$umodel = umodel
arglist$midx = midx
arglist$eidx = eidx
arglist$mpars = mpars
arglist$ipars = ipars
arglist$estidx = estidx
arglist$dccN = npars
arglist$stdresid = stdresid
arglist$H = H
arglist$Qbar = Qbar
# FDCC specific
arglist$fdccAidx = fdccAidx
arglist$fdccBidx = fdccBidx
arglist$sidx = model$sidx
assign("rmgarch_llh", 1, envir = mgarchenv)
#-----------------------------------------------------------------------------------
filt = switch(model$modeldesc$distribution,
mvnorm = .dccmakefiltermodel(garchmodel = "fdccnorm", f = normalfilter.fdccLLH2,
arglist = arglist, timer = 0, message = 0, fname = "normalfilter.fdccLLH2"))
model$mpars = mpars
model$ipars = ipars
model$pars[,1] = ipars[,1]
model$midx = midx
model$eidx = eidx
model$umodel = umodel
filt$Qbar = Qbar
filt$realizedVol = realizedVol
filt$timer = Sys.time() - tic
ans = new("DCCfilter",
mfilter = filt,
model = model)
return(ans)
}
.fdccforecast = function(fit, n.ahead = 1, n.roll = 0,
external.forecasts = list(mregfor = NULL, vregfor = NULL), cluster = NULL, ...)
{
# checks for out.sample wrt n.roll
model = fit@model
modelinc = model$modelinc
ns = model$modeldata$n.start
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]
}
if( modelinc[1]>0 ){
if( modelinc[2] > 0 ){
if( is.null(external.forecasts$mregfor ) ){
warning("\nExternal Regressor Forecasts Matrix NULL...setting to zero...\n")
mregfor = matrix(0, ncol = modelinc[2], nrow = (n.roll + n.ahead) )
} else{
if( dim(mregfor)[2] != modelinc[2] ) 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.ahead!=1) stop("\ndccforecast-->error: only n.ahead==1 supported\n")
if( n.roll > ns ) stop("\ndccforecast-->error: n.roll greater than out.sample!", call. = FALSE)
#VARf = .varpredict(fit, n.ahead, n.roll, mregfor)$Mu
mu = varxforecast(X = fit@model$modeldata$data, Bcoef = model$varcoef, p = modelinc[1],
out.sample = ns, n.ahead = n.ahead, n.roll = n.roll, mregfor = mregfor)
} else{
mu = NULL
}
exf = external.forecasts
if( modelinc[1] > 0 ){
exf$mregfor = NULL
}
ans = .fdccforecastm(fit, n.ahead = n.ahead, n.roll = n.roll, external.forecasts = exf, cluster = cluster, ...)
if(modelinc[1]==0) mu = ans$mu
model$n.roll = n.roll
model$n.ahead = n.ahead
# 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) )
mforecast = list( H = ans$H, R = ans$R, Q = ans$Q, mu = mu )
ans = new("DCCforecast",
mforecast = mforecast,
model = model)
return( ans )
}
# 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
.fdccforecastm = function(fit, n.ahead = 1, n.roll = 10,
external.forecasts = list(mregfor = NULL, vregfor = NULL),
cluster = NULL, ...)
{
model = fit@model
modelinc = model$modelinc
umodel = model$umodel
m = dim(umodel$modelinc)[2]
ns = fit@model$modeldata$n.start
# first check whether the data needs to be filtered by VAR:
Data = fit@model$modeldata$data
if(modelinc[1]>0){
zdata = varxfilter(Data, p = model$modelinc[1], Bcoef = model$varcoef,
exogen = fit@model$modeldata$mexdata, postpad = c("constant"))$xresiduals
} else{
zdata = Data
}
fpars = lapply(1:m, FUN = function(i) fit@model$mpars[fit@model$midx[,i]==1,i])
mspec = .makemultispec(umodel$modelinc, umodel$modeldesc$vmodel, umodel$modeldesc$vsubmodel,
umodel$modeldata$mexdata, umodel$modeldata$vexdata, umodel$start.pars,
fpars, NULL)
filterlist = multifilter(multifitORspec = mspec, data = xts(zdata, fit@model$modeldata$index[1:nrow(zdata)]), out.sample = 0,
n.old = fit@model$modeldata$T, cluster = cluster, realizedVol = fit@mfit$realizedVol)
# all.equal(head(sigma(filterlist)), head(fit@model$sigma) )
n.roll = n.roll + 1
m = length(mspec@spec)
out.sample = fit@model$modeldata$n.start
mo = max(modelinc[4:5])
# we augmented n.roll before, now subtract
# forclist = multiforecast(multifitORspec = fitlist, n.ahead = n.ahead, n.roll = n.roll - 1)
forclist = multiforecast(multifitORspec = mspec, data = xts(zdata, fit@model$modeldata$index[1:nrow(zdata)]), n.ahead = n.ahead,
out.sample = ns, n.roll = n.roll - 1, external.forecasts = external.forecasts,
cluster = cluster, realizedVol = fit@mfit$realizedVol, ...)
# ToDo:
if(modelinc[1] == 0){
mu = array(NA, dim=c(n.ahead, m, n.roll))
f = lapply(forclist@forecast, function(x) fitted(x))
for(i in 1:n.roll) mu[,,i] = matrix( sapply( f, function(x) x[,i] ), ncol = m)
} else{
mu = NULL
}
sig = sigma(filterlist)
resid = residuals(filterlist)
stdresid = resid/sig
T = dim(fit@mfit$H)[3]
## 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)
f = lapply(forclist@forecast, function(x) sigma(x))
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( f, function(x) x[,i] ), ncol = m)
ans = .fdccforecastn(model, Qstart, Rstart, Hstart, xQbar, xstdresids, xfsig, n.ahead, mo)
Rtfor[[i]] = ans$Rtfor
Qtfor[[i]] = ans$Qtfor
Htfor[[i]] = ans$Htfor
Qstart = last( .abind_safe(Qstart, ans$Qtfor[, , 1]), mo )
Rstart = last( .abind_safe(Rstart, ans$Rtfor[, , 1]), mo )
Hstart = last( .abind_safe(Hstart, ans$Htfor[, , 1]), mo )
}
forc = list( H = Htfor, R = Rtfor, Q = Qtfor, mu = mu )
return(forc)
}
.fdccforecastn = function(model, Qstart, Rstart, Hstart, Qbar, stdresids, fsig, n.ahead, mo)
{
m = dim(Qbar)[1]
modelinc = model$modelinc
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]
tmp = getfdccpars(model$ipars, model)
A = tmp$A
B = tmp$B
BB = (B %*% t(B))
AA = (A %*% t(A))
C = matrix(1, NROW(BB), NCOL(BB)) - BB - AA
Q = C * Qbar
Qstar = diag( 1/sqrt( diag(Q) ) , m, m)
Rbar = Qstar %*% Q %*% t(Qstar)
for(i in 1:n.ahead){
if(i==1){
Qtfor[, , mo + i] = Q
if(modelinc[4]>0){
for(j in 1:modelinc[4]){
Qtfor[ , , mo + i] = Qtfor[ , , mo + i] + (A[,j]%*%t(A[,j])) * (stdresids[(mo + i - j), ] %*% t(stdresids[(mo + i - j), ]))
}
}
if(modelinc[5]>0){
for(j in 1:modelinc[5]){
Qtfor[ , , mo + i] = Qtfor[ , , mo + i] + (B[,j]%*%t(B[,j])) * Qtfor[ , , mo + i - j]
}
}
Qtmp = diag( 1/sqrt( diag(Qtfor[ , , mo + i]) ) , m, m)
Rtfor[ , , mo + i] = Qtmp %*% Qtfor[ , , mo + i] %*% t(Qtmp)
Dtmp = diag(fsig[i, ], m, m)
Htfor[ , , mo + i] = Dtmp %*% Rtfor[ , , mo + i] %*% Dtmp
} else{
Rtfor[ , , mo + i] = Rbar
if(modelinc[4]>0){
for(j in 1:modelinc[4]){
Rtfor[ , , mo + i] = Rtfor[ , , mo + i] + (A[,j]%*%t(A[,j])) * (Rtfor[ , , mo + i - j]-Rbar)
}
}
if(modelinc[5]>0){
for(j in 1:modelinc[5]){
Rtfor[ , , mo + i] = Rtfor[ , , mo + i] + (B[,j]%*%t(B[,j])) * (Rtfor[ , , mo + i - j]-Rbar)
}
}
Dtmp = diag(fsig[i, ], m, m)
Htfor[ , , mo + i] = Dtmp %*% Rtfor[ , , mo + i] %*% Dtmp
}
}
return( list( Rtfor = Rtfor[, , -(1:mo), drop = FALSE], Htfor = Htfor[, , -(1:mo), drop = FALSE],
Qtfor = Qtfor[, , -(1:mo), drop = FALSE] ) )
}
.fdccforecastn2 = function(model, Qstart, Rstart, Hstart, Qbar, stdresids, fsig, n.ahead, mo)
{
m = dim(Qbar)[1]
modelinc = model$modelinc
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]
tmp = getfdccpars(model$ipars, model)
A = tmp$A
B = tmp$B
BB = (B %*% t(B))
AA = (A %*% t(A))
C = matrix(1, NROW(BB), NCOL(BB)) - BB - AA
Q = C * Qbar
for(i in 1:n.ahead){
if(i==1){
Qtfor[, , mo + i] = Q
if(modelinc[4]>0){
for(j in 1:modelinc[4]){
Qtfor[ , , mo + i] = Qtfor[ , , mo + i] + (A[,j]%*%t(A[,j])) * (stdresids[(mo + i - j), ] %*% t(stdresids[(mo + i - j), ]))
}
}
if(modelinc[5]>0){
for(j in 1:modelinc[5]){
Qtfor[ , , mo + i] = Qtfor[ , , mo + i] + (B[,j]%*%t(B[,j])) * Qtfor[ , , mo + i - j]
}
}
Qtmp = diag( 1/sqrt( diag(Qtfor[ , , mo + i]) ) , m, m)
Rtfor[ , , mo + i] = Qtmp %*% Qtfor[ , , mo + i] %*% t(Qtmp)
Dtmp = diag(fsig[i, ], m, m)
Htfor[ , , mo + i] = Dtmp %*% Rtfor[ , , mo + i] %*% Dtmp
} else{
Qtfor[, , mo + i] = Q
if(modelinc[4]>0){
for(j in 1:modelinc[4]){
Qtfor[ , , mo + i] = Qtfor[ , , mo + i] + (A[,j]%*%t(A[,j])) * (Qtfor[ , , mo + i - j])
}
}
if(modelinc[5]>0){
for(j in 1:modelinc[5]){
Qtfor[ , , mo + i] = Qtfor[ , , mo + i] + (B[,j]%*%t(B[,j])) * (Qtfor[ , , mo + i - j])
}
}
Qtmp = diag( 1/sqrt( diag(Qtfor[ , , mo + i]) ) , m, m)
Rtfor[ , , mo + i] = Qtmp %*% Qtfor[ , , mo + i] %*% t(Qtmp)
Dtmp = diag(fsig[i, ], m, m)
Htfor[ , , mo + i] = Dtmp %*% Rtfor[ , , mo + i] %*% Dtmp
}
}
return( list( Rtfor = Rtfor[, , -(1:mo), drop = FALSE], Htfor = Htfor[, , -(1:mo), drop = FALSE],
Qtfor = Qtfor[, , -(1:mo), drop = FALSE] ) )
}
.fdccsim.fit = function(fitORspec, n.sim = 1000, n.start = 0, m.sim = 1,
startMethod = c("unconditional", "sample"), presigma = NULL,
preresiduals = NULL, prereturns = NULL, preQ = NULL, preZ = NULL,
Qbar = NULL, rseed = NULL, mexsimdata = NULL,
vexsimdata = NULL, cluster = NULL, VAR.fit = NULL, prerealized = NULL, ...)
{
fit = fitORspec
T = fit@model$modeldata$T
Data = fit@model$modeldata$data[1:T,]
m = dim(Data)[2]
mo = max(fit@model$modelinc[4:5])
mg = fit@model$maxgarchOrder
startMethod = startMethod[1]
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] )
}
model = fit@model
umodel = model$umodel
modelinc = model$modelinc
fpars = lapply(1:m, FUN = function(i) fit@model$mpars[fit@model$midx[,i]==1,i])
mspec = .makemultispec(umodel$modelinc, umodel$modeldesc$vmodel, umodel$modeldesc$vsubmodel,
umodel$modeldata$mexdata, umodel$modeldata$vexdata, umodel$start.pars,
fpars, NULL)
# Technically, we should pass an array of size (m, m, mo) for the
# sample but how many use a dccOrder > 1? (and it would mean modifying
# the C++ code to deal with arrays/lists for this input).
if(startMethod == "sample"){
if(is.null(preZ)){
preZ = matrix(tail(residuals(fit)/sigma(fit), mo), ncol = m)
} else{
preZ = matrix(tail(preZ, 1), ncol = m, nrow = mo, byrow = TRUE)
}
if(is.null(preQ)){
preQ = fit@mfit$Q[[length(fit@mfit$Q)]]
} else{
dcc.symcheck(preQ, m, d = NULL)
}
Rbar = preQ/(sqrt(diag(preQ)) %*% t(sqrt(diag(preQ))) )
} else{
if(is.null(preZ)){
preZ = matrix(0, ncol = m, nrow = mo)
} else{
preZ = matrix(tail(preZ, 1), ncol = m, nrow = mo, byrow = TRUE)
}
Rbar = cor(Data)
if(is.null(preQ)){
preQ = Rbar
} else{
dcc.symcheck(preQ, m, d = NULL)
Rbar = preQ/(sqrt(diag(preQ)) %*% t(sqrt(diag(preQ))) )
}
}
if(is.null(Qbar)){
Qbar = fit@mfit$Qbar
} else{
dcc.symcheck(Qbar, m, d = NULL)
}
uncv = sapply(mspec@spec, FUN = function(x) uncvariance(x))
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( dim(presigma)[1] != mg )
stop(paste("\ndccsim-->error: wrong row dimension for presigma (need ", mg, " rows.", sep = ""))
} else{
if(startMethod == "sample"){
mx = max(sapply(mspec@spec, FUN = function(x) x@model$maxOrder))
presigma = matrix(NA, ncol = m, nrow = mx)
tmp = last(fit@mfit$H, mx)
for(i in 1:mx) presigma[i,] = sqrt(diag(tmp[,,i]))
}
}
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( dim(preresiduals)[1] != mg )
stop(paste("\ndccsim-->error: wrong row dimension for preresiduals (need ", mg, " rows.", sep = ""))
} else{
if(startMethod == "sample"){
mx = max(sapply(mspec@spec, FUN = function(x) x@model$maxOrder))
preresiduals = matrix(NA, ncol = m, nrow = mx)
tmp = tail(fit@model$residuals, mx)
for(i in 1:mx) preresiduals[i,] = tmp[i,]
}
}
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.")
if( dim(prereturns)[1] != mg )
stop(paste("\ndccsim-->error: wrong row dimension for prereturns (need ", mg, " rows.", sep = ""))
} else{
if(startMethod == "sample"){
mx = max(sapply(mspec@spec, FUN = function(x) x@model$maxOrder))
prereturns = matrix(NA, ncol = m, nrow = mx)
tmp = tail(Data, mx)
for(i in 1:mx) prereturns[i,] = tmp[i,]
}
}
if(fit@model$umodel$modeldesc$vmodel[1]=="realGARCH"){
if( !is.null(prerealized) ){
if( !is.matrix(prerealized) )
stop("\ndccsim-->error: prerealized must be a matrix.")
if( dim(prerealized)[2] != m )
stop("\ndccsim-->error: wrong column dimension for prerealized.")
if( dim(prerealized)[1] != mg )
stop(paste("\ndccsim-->error: wrong row dimension for prerealized (need ", mg, " rows.", sep = ""))
} else{
if(startMethod == "sample"){
mx = max(sapply(mspec@spec, FUN = function(x) x@model$maxOrder))
prerealized = matrix(NA, ncol = m, nrow = mx)
tmp = tail(fit@mfit$realizedVol[1:T,], mx)
for(i in 1:mx) prerealized[i,] = tmp[i,]
}
}
} else{
mx = max(sapply(mspec@spec, FUN = function(x) x@model$maxOrder))
prerealized = matrix(NA, ncol = m, nrow = mx)
}
# switch distributions
if(fit@model$modeldesc$distribution == "mvnorm"){
if(length(rseed) == 1){
set.seed( rseed )
tmp = matrix(rnorm(m * (n.sim + n.start) * m.sim, 0, 1), ncol = m, nrow = n.sim+n.start)
z = array(NA, dim = c(n.sim + n.start + mo, m, m.sim))
for(i in 1:m.sim) z[,,i] = rbind(preZ, tmp)
} 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] = rbind(preZ, matrix(rnorm(m * (n.sim + n.start), 0, 1), nrow = n.sim + n.start, ncol = m))
}
}
} else if(fit@model$modeldesc$distribution == "mvlaplace"){
if(length(rseed) == 1){
set.seed( rseed )
tmp = matrix(rdist(distribution = "ged", n = m * (n.sim + n.start) * m.sim, 0, 1, shape = 1), ncol = m, nrow = n.sim+n.start)
z = array(NA, dim = c(n.sim + n.start + mo, m, m.sim))
for(i in 1:m.sim) z[,,i] = rbind(preZ, tmp)
} 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] = rbind(preZ, matrix(rdist(distribution = "ged", n = m * (n.sim + n.start), 0, 1, shape = 1), nrow = n.sim + n.start, ncol = m))
}
}
} else{
if(length(rseed) == 1){
set.seed( rseed )
tmp = matrix(rdist(distribution = "std", n = m * (n.sim + n.start) * m.sim, 0, 1, shape = rshape(fit)), ncol = m, nrow = n.sim+n.start)
z = array(NA, dim = c(n.sim + n.start + mo, m, m.sim))
for(i in 1:m.sim) z[,,i] = rbind(preZ, tmp)
} 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] = rbind(preZ, matrix(rdist(distribution = "std", n = m * (n.sim + n.start), 0, 1, shape = rshape(fit)), nrow = n.sim + n.start, ncol = m))
}
}
}
# ok now to expand rseed (but still allow for replication given initial seed)
if(length(rseed) == 1){
rseed = c(rseed, (1:m.sim)*(rseed+1))
}
tmp = getfdccpars(model$ipars, model)
A = tmp$A
B = tmp$B
BB = (B %*% t(B))
AA = (A %*% t(A))
C = matrix(1, NROW(BB), NCOL(BB)) - BB - AA
# model, Z, A, B, C, Qbar, preQ, Rbar, mo, n.sim, n.start, m, rseed
simRes = simX = simR = simQ = simH = simSeries = vector(mode = "list", length = m.sim)
if( !is.null(cluster) ){
simH = vector(mode = "list", length = m.sim)
simX = vector(mode = "list", length = m.sim)
clusterEvalQ(cluster, require(rmgarch))
clusterExport(cluster, c("model", "z", "A", "B", "C",
"preQ", "Rbar", "Qbar", "mo", "n.sim", "n.start", "m",
"rseed",".fdccsimf"), envir = environment())
mtmp = parLapply(cluster, as.list(1:m.sim), fun = function(j){
.fdccsimf(model, Z = z[,,j], A = A, B = B, C = C,
Qbar = Qbar, preQ = preQ, Rbar = Rbar, mo = 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))
#simZ = lapply(mtmp, FUN = function(x) x$Z)
simZ = vector(mode = "list", length = m)
for(i in 1:m) simZ[[i]] = sapply(mtmp, FUN = function(x) x$Z[,i])
clusterExport(cluster, c("fit", "n.sim", "n.start", "m.sim",
"startMethod", "simZ", "presigma", "preresiduals",
"prereturns", "mexsimdata", "vexsimdata","prerealized"),
envir = environment())
xtmp = parLapply(cluster, as.list(1:m), fun = function(j){
maxx = mspec@spec[[j]]@model$maxOrder;
htmp = ugarchpath(mspec@spec[[j]], n.sim = n.sim + n.start, n.start = 0, m.sim = m.sim,
custom.dist = list(name = "sample", distfit = matrix(simZ[[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) | model$modelinc[1]>0 ) NA else tail(prereturns[,j], maxx),
mexsimdata = if( model$modelinc[1]==0 ) mexsimdata[[j]] else NULL,
vexsimdata = vexsimdata[[j]], prerealized = tail(prerealized[,j], maxx));
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))
})
} 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){
.fdccsimf(model, Z = z[,,j], A = A, B = B, C = C,
Qbar = Qbar, preQ = preQ, Rbar = Rbar, mo = 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))
#simZ = lapply(mtmp, FUN = function(x) x$Z)
simZ = vector(mode = "list", length = m)
for(i in 1:m) simZ[[i]] = sapply(mtmp, FUN = function(x) x$Z[,i])
xtmp = lapply(as.list(1:m), FUN = function(j){
maxx = mspec@spec[[j]]@model$maxOrder;
htmp = ugarchpath(mspec@spec[[j]], n.sim = n.sim + n.start, n.start = 0, m.sim = m.sim,
custom.dist = list(name = "sample", distfit = matrix(simZ[[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) | model$modelinc[1]>0 ) NA else tail(prereturns[,j], maxx),
mexsimdata = if( model$modelinc[1]==0 ) mexsimdata[[j]] else NULL,
vexsimdata = vexsimdata[[j]], prerealized = tail(prerealized[,j], maxx));
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
}
if(model$modelinc[1]>0){
simxX = array(NA, dim = c(n.sim+n.start, m, m.sim))
for(i in 1:m.sim) simxX[,,i] = 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+n.start)
} else{
simxX = array(NA, dim = c(n.sim+n.start, m, m.sim))
for(i in 1:m.sim) simxX[,,i] = 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(tail(matrix(simxX[,,i], ncol = m), n.sim), nrow = n.sim)
}
if( model$modelinc[1]>0 ){
simRes = simX
simX = mvmean.varsim(model = model, Data = Data, res = simX,
mexsimdata = mexsimdata, prereturns = prereturns, m.sim = m.sim,
n.sim = n.sim, n.start = n.start, startMethod = startMethod,
cluster = cluster)
} 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$simZ = simZ
msim$simRes = simRes
#msim$Z = stdresid
msim$rseed = rseed
model$n.sim = n.sim
model$m.sim = m.sim
model$n.start = n.start
model$startMethod = startMethod[1]
ans = new("DCCsim",
msim = msim,
model = model)
return( ans )
}
.fdccsim.spec = function(fitORspec, n.sim = 1000, n.start = 0, m.sim = 1,
startMethod = c("unconditional", "sample"), presigma = NULL,
preresiduals = NULL, prereturns = NULL, preQ = NULL, preZ = NULL,
Qbar = NULL, rseed = NULL, mexsimdata = NULL,
vexsimdata = NULL, cluster = NULL, VAR.fit = NULL,
prerealized = NULL, ...)
{
spec = fitORspec
startMethod = startMethod[1]
# check that VAR.fit is included if used
if( spec@model$modelinc[1]>0 ){
if( is.null(VAR.fit) ) stop("\ndccsim-->error: VAR.fit must not be NULL for VAR method when calling dccsim using spec!", call. = FALSE)
}
model = spec@model
model$umodel = spec@umodel
m = dim(spec@umodel$modelinc)[2]
mo = max(spec@model$modelinc[4:5])
mg = spec@model$maxgarchOrder
model$modeldata$asset.names = paste("Asset", 1:m, sep = "")
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] )
}
if(is.null(preZ)){
preZ = matrix(0, ncol = m, nrow = mo)
} else{
preZ = matrix(tail(preZ, 1), ncol = m, nrow = mo, byrow = TRUE)
}
if(is.null(preQ)){
stop("\ndccsim-->error: preQ cannot be NULL when method uses spec!")
} else{
dcc.symcheck(preQ, m, d = NULL)
Rbar = preQ/(sqrt(diag(preQ)) %*% t(sqrt(diag(preQ))) )
}
if(is.null(Qbar)){
stop("\ndccsim-->error: Qbar cannot be NULL when method uses spec!")
} else{
dcc.symcheck(Qbar, m, d = NULL)
}
model = spec@model
umodel = spec@umodel
modelinc = model$modelinc
midx = .fullinc(modelinc, umodel)
mspec = .makemultispec(umodel$modelinc, umodel$modeldesc$vmodel, umodel$modeldesc$vsubmodel,
umodel$modeldata$mexdata, umodel$modeldata$vexdata, umodel$start.pars,
umodel$fixed.pars, NULL)
# uncvariance has uGARCHfit AND uGARCHspec dispatch methods
uncv = sapply(mspec@spec, FUN = function(x) uncvariance(x))
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( dim(presigma)[1] != mg )
stop(paste("\ndccsim-->error: wrong row dimension for presigma (need ", mg, " rows.", sep = ""))
}
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( dim(preresiduals)[1] != mg )
stop(paste("\ndccsim-->error: wrong row dimension for preresiduals (need ", mg, " rows.", sep = ""))
}
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.")
if( dim(prereturns)[1] != mg )
stop(paste("\ndccsim-->error: wrong row dimension for prereturns (need ", mg, " rows.", sep = ""))
}
if(spec@umodel$modeldesc$vmodel[1]=="realGARCH"){
if( !is.null(prerealized) ){
if( !is.matrix(prerealized) )
stop("\ndccsim-->error: prerealized must be a matrix.")
if( dim(prerealized)[2] != m )
stop("\ndccsim-->error: wrong column dimension for prerealized.")
if( dim(prerealized)[1] != mg )
stop(paste("\ndccsim-->error: wrong row dimension for prerealized (need ", mg, " rows.", sep = ""))
}
} else{
prerealized = matrix(NA, ncol = m, nrow = mg)
}
# we use the GED distribution with nu = 1 which corresponds to the Laplace case.
if(model$modeldesc$distribution == "mvnorm"){
if(length(rseed) == 1){
set.seed( rseed )
tmp = matrix(rnorm(m * (n.sim + n.start) * m.sim, 0, 1), ncol = m, nrow = n.sim+n.start)
z = array(NA, dim = c(n.sim + n.start + mo, m, m.sim))
for(i in 1:m.sim) z[,,i] = rbind(preZ, tmp)
} 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] = rbind(preZ, matrix(rnorm(m * (n.sim + n.start), 0, 1), nrow = n.sim + n.start, ncol = m))
}
}
} else if(model$modeldesc$distribution == "mvlaplace"){
if(length(rseed) == 1){
set.seed( rseed )
tmp = rdist(distribution = "ged", n = m * (n.sim + n.start) * m.sim, 0, 1, shape = 1)
z = array(NA, dim = c(n.sim + n.start + mo, m, m.sim))
for(i in 1:m.sim) z[,,i] = rbind(preZ, tmp)
} 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] = rbind(preZ, matrix(rdist(distribution = "ged", n = m * (n.sim + n.start), 0, 1, shape = 1), nrow = n.sim + n.start, ncol = m))
}
}
} else{
if(length(rseed) == 1){
set.seed( rseed )
tmp = rdist(distribution = "std", n = m * (n.sim + n.start) * m.sim, 0, 1, shape = model$mpars["mshape",1])
z = array(NA, dim = c(n.sim + n.start + mo, m, m.sim))
for(i in 1:m.sim) z[,,i] = rbind(preZ, tmp)
} 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] = rbind(preZ, matrix(rdist(distribution = "std", n = m * (n.sim + n.start), 0, 1, shape = model$mpars["mshape",1]), nrow = n.sim + n.start, ncol = m))
}
}
}
# ok now to expand rseed
if(length(rseed) == 1){
rseed = c(rseed, as.integer(runif(m.sim, 1, Sys.time())))
}
tmp = getfdccpars(model$ipars, model)
A = tmp$A
B = tmp$B
BB = (B %*% t(B))
AA = (A %*% t(A))
C = matrix(1, NROW(BB), NCOL(BB)) - BB - AA
simRes = simX = simR = simQ = simH = simSeries = vector(mode = "list", length = m.sim)
if( !is.null(cluster) ){
simH = vector(mode = "list", length = m.sim)
simX = vector(mode = "list", length = m.sim)
clusterEvalQ(cluster, require(rmgarch))
clusterExport(cluster, c("model", "z", "A", "B", "C", "preQ",
"Rbar", "Qbar", "mo", "n.sim", "n.start", "m",
"rseed",".fdccsimf"), envir = environment())
mtmp = parLapply(cluster, as.list(1:m.sim), fun = function(j){
.fdccsimf(model, Z = z[,,j], A = A, B = B, C = C,
Qbar = Qbar, preQ = preQ, Rbar = Rbar, mo = 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))
#simZ = lapply(mtmp, FUN = function(x) x$Z)
simZ = vector(mode = "list", length = m)
for(i in 1:m) simZ[[i]] = sapply(mtmp, FUN = function(x) x$Z[,i])
clusterExport(cluster, c("mspec", "n.sim", "n.start", "m.sim",
"startMethod", "simZ", "presigma", "preresiduals",
"prereturns", "mexsimdata", "vexsimdata","prerealized"),
envir = environment())
xtmp = parLapply(cluster, as.list(1:m), fun = function(j){
maxx = mspec@spec[[j]]@model$maxOrder;
htmp = ugarchpath(mspec@spec[[j]], n.sim = n.sim + n.start, n.start = 0, m.sim = m.sim,
custom.dist = list(name = "sample", distfit = matrix(simZ[[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) | model$modelinc[1]>0 ) NA else tail(prereturns[,j], maxx),
mexsimdata = if( model$modelinc[1]==0 ) mexsimdata[[j]] else NULL,
vexsimdata = vexsimdata[[j]], prerealized = tail(prerealized[,j], maxx))
h = matrix(tail(htmp@path$sigmaSim^2, n.sim), nrow = n.sim)
x = matrix(htmp@path$seriesSim, nrow = n.sim + n.start)
xres = matrix(htmp@path$residSim, nrow = n.sim + n.start)
return(list(h = h, x = x, xres = xres))
})
} 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){
.fdccsimf(model, Z = z[,,j], A = A, B = B, C = C,
Qbar = Qbar, preQ = preQ, Rbar = Rbar, mo = 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))
#simZ = lapply(mtmp, FUN = function(x) x$Z)
simZ = vector(mode = "list", length = m)
for(i in 1:m) simZ[[i]] = sapply(mtmp, FUN = function(x) x$Z[,i])
xtmp = lapply(as.list(1:m), FUN = function(j){
maxx = mspec@spec[[j]]@model$maxOrder;
htmp = ugarchpath(mspec@spec[[j]], n.sim = n.sim + n.start, n.start = 0, m.sim = m.sim,
custom.dist = list(name = "sample", distfit = matrix(simZ[[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) | model$modelinc[1]>0 ) NA else tail(prereturns[,j], maxx),
mexsimdata = if( model$modelinc[1]==0 ) mexsimdata[[j]] else NULL, vexsimdata = vexsimdata[[j]],
prerealized = tail(prerealized[,j], maxx));
h = matrix(tail(htmp@path$sigmaSim^2, n.sim), nrow = n.sim);
x = matrix(htmp@path$seriesSim, nrow = n.sim + n.start);
xres = matrix(htmp@path$residSim, nrow = n.sim + n.start);
return(list(h = h, x = x, xres = xres))
})
}
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
}
if(model$modelinc[1]>0){
simxX = array(NA, dim = c(n.sim+n.start, m, m.sim))
for(i in 1:m.sim) simxX[,,i] = 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+n.start)
} else{
simxX = array(NA, dim = c(n.sim+n.start, m, m.sim))
for(i in 1:m.sim) simxX[,,i] = 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(tail(matrix(simxX[,,i], ncol = m), n.sim), nrow = n.sim)
}
simxRes = array(NA, dim = c(n.sim, m, m.sim))
for(i in 1:n.sim) simxRes[i,,] = t(sapply(xtmp, FUN = function(x) as.numeric(x$xres[i,])))
simRes = vector(mode = "list", length = m.sim)
for(i in 1:m.sim) simRes[[i]] = matrix(simxRes[,,i], nrow = n.sim)
if( model$modelinc[1]>0 ){
model$varcoef = VAR.fit$Bcoef
Data = VAR.fit$xfitted
simRes = simX
simX = mvmean.varsim(model = model, Data = Data, res = simX,
mexsimdata = mexsimdata, prereturns = prereturns, m.sim = m.sim,
n.sim = n.sim, n.start = n.start, startMethod = startMethod,
cluster = cluster)
} 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$simZ = simZ
msim$simRes = simRes
#msim$Z = stdresid
msim$rseed = rseed
msim$model$Data = NULL
model$n.sim = n.sim
model$m.sim = m.sim
model$n.start = n.start
model$startMethod = "unconditional"
ans = new("DCCsim",
msim = msim,
model = model)
return( ans )
}
.fdccsimf = function(model, Z, A, B, C, Qbar, preQ, Rbar, mo, n.sim, n.start, m, rseed)
{
modelinc = model$modelinc
ipars = model$pars
idx = model$pidx
n = n.sim + n.start + mo
set.seed(rseed[1]+1)
# All distributions here make use of this since they are mixtures of normals
NZ = matrix(rnorm(m * (n.sim+n.start+mo)), nrow = n.sim + n.start + mo, ncol = m)
# SEXP model, SEXP A, SEXP B, SEXP C, SEXP Qbar, SEXP preQ, SEXP Rbar, SEXP Z, SEXP NZ
res = switch(model$modeldesc$distribution,
mvnorm = .Call("fdccsimmvn",
model = as.integer(modelinc),
A = A,
B = B,
C = C,
Qbar = as.matrix(Qbar),
preQ = as.matrix(preQ),
Rbar = as.matrix(Rbar),
Z = as.matrix(Z),
NZ = as.matrix(NZ),
PACKAGE = "rmgarch")
)
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[[1]][[i]]
}
ans = list( Q = Q, R = R, Z = res[[3]])
return( ans )
}
.fullinc_fdcc = function(modelinc, umodel){
m = dim(umodel$modelinc)[2]
# the par matrix will have the max parameter vis Asset (if an asset does not have
# the parameter it will be zero...for matrix estimation.
# the index places the solver estimates into their place in the matrix
###################################################################################
# Asset1 .... Asset2 .... AssetN ... ........ Joint
# mu
# ar
# ma
# arfima
# archm
# mxreg
# omega
# alpha
# beta
# gamma
# delta
# lambda
# eta1
# eta2
# vxreg
# skew
# shape
# ghlambda
# xi
# fdccC
# fdccA
# fdccB
# mshape
# mskew
vecmax = rep(0, 19)
names(vecmax) = rownames(umodel$modelinc[1:19,])
vecmax = apply(umodel$modelinc, 1, FUN = function(x) max(x) )
maxOrder = apply(umodel$modelinc, 2, FUN = function(x) max(c(x[2], x[3], x[8], x[9])))
sumv = 19 + sum(pmax(1, vecmax[c(2,3,6,8,9,10,11,12,13,15,16)])) - 11
tmpmat = matrix(0, ncol = m+1, nrow = sumv)
nx = 0
pnames = NULL
# mu
if(vecmax[1]>0){
tmpmat[1, 1:m] = umodel$modelinc[1,]
}
nx = nx + max(1, vecmax[1])
pnames = c(pnames, "mu")
# ar
if(vecmax[2]>0){
for(i in 1:vecmax[2]){
tmpmat[nx+i, 1:m] = as.integer( umodel$modelinc[2,] >= i)
pnames = c(pnames, paste("ar", i, sep = ""))
}
} else{
pnames = c(pnames, "ar")
}
nx = nx + max(1, vecmax[2])
# ma
if(vecmax[3]>0){
for(i in 1:vecmax[3]){
tmpmat[nx+i, 1:m] = as.integer( umodel$modelinc[3,] >= i)
pnames = c(pnames, paste("ma", i, sep = ""))
}
} else{
pnames = c(pnames, "ma")
}
nx = nx + max(1, vecmax[3])
# arfima
if(vecmax[4]>0){
tmpmat[nx+1, 1:m] = umodel$modelinc[4, ]
}
nx = nx + max(1, vecmax[4])
pnames = c(pnames, "arfima")
# archm no allowed
nx = nx + max(1, vecmax[5])
pnames = c(pnames, "archm")
# mxreg
if(vecmax[6]>0){
for(i in 1:vecmax[6]){
tmpmat[nx+i, 1:m] = as.integer(umodel$modelinc[6, ] >= i)
pnames = c(pnames, paste("mxreg", i, sep = ""))
}
} else{
pnames = c(pnames, "mxreg")
}
nx = nx + max(1, vecmax[6])
# omega
if(vecmax[7]>0){
tmpmat[nx+1, 1:m] = umodel$modelinc[7, ]
}
nx = nx + max(1, vecmax[7])
pnames = c(pnames, "omega")
# alpha
if(vecmax[8]>0){
for(i in 1:vecmax[8]){
tmpmat[nx+i, 1:m] = as.integer( umodel$modelinc[8, ] >= i)
pnames = c(pnames, paste("alpha", i, sep = ""))
}
} else{
pnames = c(pnames, "alpha")
}
nx = nx + max(1, vecmax[8])
# beta
if(vecmax[9]>0){
for(i in 1:vecmax[9]){
tmpmat[nx+i, 1:m] = as.integer( umodel$modelinc[9, ] >= i)
pnames = c(pnames, paste("beta", i, sep = ""))
}
} else{
pnames = c(pnames, "beta")
}
nx = nx + max(1, vecmax[9])
# gamma
if(vecmax[10]>0){
for(i in 1:vecmax[10]){
tmpmat[nx+i, 1:m] = as.integer( umodel$modelinc[10, ] >= i)
pnames = c(pnames, paste("gamma", i, sep = ""))
}
} else{
pnames = c(pnames, "gamma")
}
nx = nx + max(1, vecmax[10])
# eta1
if(vecmax[11]>0){
for(i in 1:vecmax[11]){
tmpmat[nx+i, 1:m] = as.integer( umodel$modelinc[11, ] >= i)
pnames = c(pnames, paste("eta1", i, sep = ""))
}
} else{
pnames = c(pnames, "eta1")
}
nx = nx + max(1, vecmax[11])
# eta2
if(vecmax[12]>0){
for(i in 1:vecmax[12]){
tmpmat[nx+i, 1:m] = as.integer( umodel$modelinc[12, ] >= i)
pnames = c(pnames, paste("eta2", i, sep = ""))
}
} else{
pnames = c(pnames, "eta2")
}
nx = nx + max(1, vecmax[12])
# delta
if(vecmax[13]>0){
tmpmat[nx+1, 1:m] = umodel$modelinc[13, ]
}
nx = nx + max(1, vecmax[13])
pnames = c(pnames, "delta")
# lambda
if(vecmax[14]>0){
tmpmat[nx+1, 1:m] = umodel$modelinc[14, ]
}
nx = nx + max(1, vecmax[14])
pnames = c(pnames, "lambda")
# vxreg
if(vecmax[15]>0){
for(i in 1:vecmax[15]){
tmpmat[nx+i, 1:m] = as.integer( umodel$modelinc[15, ] >= i)
pnames = c(pnames, paste("vxreg", i, sep = ""))
}
} else{
pnames = c(pnames, "vxreg")
}
nx = nx + max(1, vecmax[15])
# skew
if(vecmax[16]>0){
tmpmat[nx+1, 1:m] = umodel$modelinc[16, ]
}
nx = nx + max(1, vecmax[16])
pnames = c(pnames, "skew")
# shape
if(vecmax[17]>0){
tmpmat[nx+1, 1:m] = umodel$modelinc[17, ]
}
nx = nx + max(1, vecmax[17])
pnames = c(pnames, "shape")
# skew
if(vecmax[18]>0){
tmpmat[nx+1, 1:m] = umodel$modelinc[18, ]
}
nx = nx + max(1, vecmax[18])
pnames = c(pnames, "ghlambda")
if(vecmax[19]>0){
tmpmat[nx+1, 1:m] = umodel$modelinc[19, ]
}
nx = nx + max(1, vecmax[19])
pnames = c(pnames, "xi")
sumdcc = max(1, (modelinc[4]*modelinc[3]))+max(1, (modelinc[5]*modelinc[3])) +
max(1, modelinc[6]) + max(1, modelinc[6])
tmpmat = rbind(tmpmat, matrix(0, ncol = m+1, nrow = sumdcc))
# intercept
if(modelinc[4]>0){
xid = 0
for(i in 1:modelinc[4]){
for(j in 1:modelinc[3]){
tmpmat[(nx+xid+j), m+1] = 1
pnames = c(pnames, paste("fdcca[g=", j,",P=", i, "]", sep=""))
}
xid = xid + modelinc[3]
}
} else{
pnames = c(pnames, "fdcca")
}
fdccAidx = (nx+1):max( (nx+1), (nx+modelinc[3]*modelinc[4]) )
nx = nx + max(1, modelinc[4]*modelinc[3])
if(modelinc[5]>0){
xid = 0
for(i in 1:modelinc[5]){
for(j in 1:modelinc[3]){
tmpmat[(nx+xid+j), m+1] = 1
pnames = c(pnames, paste("fdccb[g=", j,",Q=", i, "]", sep=""))
}
xid = xid + modelinc[3]
}
} else{
pnames = c(pnames, "fdccb")
}
fdccBidx = (nx+1):max( (nx+1), (nx+modelinc[3]*modelinc[5]) )
nx = nx + max(1, modelinc[5]*modelinc[3])
# we allow for possibility of vector valued shape and skew parameters for future expansion
if(modelinc[6]>0){
for(i in 1:modelinc[6]){
tmpmat[nx+i, m+1] = 1
if(modelinc[6]>1) pnames = c(pnames, paste("mshape", i, sep = "")) else pnames = c(pnames, "mshape")
}
} else{
pnames = c(pnames, "mshape")
}
nx = nx + max(1, modelinc[6])
if(modelinc[7]>0){
for(i in 1:modelinc[7]){
tmpmat[nx+i, m+1] = 1
if(modelinc[7]>1) pnames = c(pnames, paste("mskew", i, sep = "")) else pnames = c(pnames, "mskew")
}
} else{
pnames = c(pnames, "mskew")
}
# NOW we have and indicator matrix
colnames(tmpmat) = c(paste("Asset", 1:m, sep = ""), "Joint")
rownames(tmpmat) = pnames
return(list(midx = tmpmat, fdccAidx = fdccAidx, fdccBidx = fdccBidx))
}
.fdcccon = function(pars, arglist){
ipars = arglist$ipars
estidx = arglist$estidx
idx = arglist$model$pidx
ipars[estidx, 1] = pars
modelinc = arglist$model$modelinc
if(modelinc[4]>0) fdccA = matrix(ipars[idx["fdcca", 1]:idx["fdcca", 2], 1], ncol = modelinc[4]) else fdccA = matrix(0, nrow = modelinc[3], ncol = modelinc[5])
if(modelinc[5]>0) fdccB = matrix(ipars[idx["fdccb", 1]:idx["fdccb", 2], 1], ncol = modelinc[5]) else fdccB = matrix(0, nrow = modelinc[3], ncol = modelinc[4])
if(modelinc[3]==1){
ans = fdccA + fdccB
} else{
# each pairwise combination (only allow DCC(1,1))
ans = apply(expand.grid(fdccA, fdccA, fdccB, fdccB), 1, function(x) x[1]*x[2] + x[3]*x[4])
}
return( as.numeric(ans) )
}
Try the rmgarch package in your browser
Any scripts or data that you put into this service are public.
rmgarch documentation built on Feb. 5, 2022, 1:07 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions .fdcccon .fullinc_fdcc .fdccsimf .fdccsim.spec .fdccsim.fit .fdccforecastn2 .fdccforecastn .fdccforecastm .fdccforecast .fdccfilter getfdccpars .fdccfit .fdccspec
#################################################################################
##
## R package rmgarch by Alexios Galanos Copyright (C) 2008-2022.
## This file is part of the R package rmgarch.
##
## The R package rmgarch 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 rmgarch 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.
##
#################################################################################
.fdccspec = function(uspec, VAR = FALSE, 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),
fdccOrder=c(1,1), fdccindex = c(1,1,1,2,2,2,3,3,3),
distribution = c("mvnorm", "mvt", "mvlaplace"),
start.pars = list(), fixed.pars = list())
{
# VAR checks
VAR.opt = list()
if(is.null(VAR)) VAR.opt$VAR = FALSE else VAR.opt$VAR = as.logical(VAR)
if(is.null(robust)) VAR.opt$robust = FALSE else VAR.opt$robust = as.logical(robust)
if(is.null(lag)) VAR.opt$lag = 1 else VAR.opt$lag = as.integer(lag)
if(is.null(lag.max)) VAR.opt$lag.max = NULL else VAR.opt$lag.max = as.integer(min(1, lag.max))
if(is.null(lag.criterion)) VAR.opt$lag.criterion = "AIC" else VAR.opt$lag.criterion = lag.criterion[1]
if(is.null(external.regressors)) VAR.opt$external.regressors = NULL else VAR.opt$external.regressors = external.regressors
rc = list("gamma" = 0.25, "delta" = 0.01, "nc" = 10, "ns" = 500)
rcmatch = match(names(robust.control), c("gamma", "delta", "nc", "ns"))
if(length(rcmatch[!is.na(rcmatch)]) > 0){
rx = which(!is.na(rcmatch))
rc[rcmatch[!is.na(rcmatch)]] = robust.control[rx]
}
VAR.opt$robust.control = rc
.eps = .Machine$double.eps
modeldata = list()
modeldesc = list()
m = length(uspec@spec)
## distribution checks
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)
## fdcc checks
fdccOrder = as.integer(abs(fdccOrder))
if(max(fdccOrder)>1) stop("\nFDCC model only allows (1,1) DCC dynamics and not higher.")
fdccindex = as.integer(abs(fdccindex))
if(any(fdccindex)==0) stop("\nfdccindex cannot contain zeros...index starts at 1!")
lenF = length(unique(fdccindex))
if(length(fdccindex)!=length(uspec@spec)) stop("\nfdccindex length not equal to length of uspec.")
# allow for case of just 1 index
if(!all(fdccindex==1) && any(diff(unique(fdccindex)))!=1) stop("\nfdccindex must use contiguous indexing, starting from one and incrementing by 1.")
modelinc = rep(0, 10)
names(modelinc) = c("var", "mvmxreg", "fdccC", "fdccA", "fdccB", "mshape", "mskew", "aux", "aux", "aux")
if(distribution == "mvt") modelinc[6] = 1
modelinc[3] = lenF
modelinc[4] = fdccOrder[1]
modelinc[5] = fdccOrder[2]
# indicator matrix of groups
sidx = matrix(0, ncol = lenF, nrow = length(fdccindex))
for(i in 1:lenF) sidx[fdccindex==i,i] = 1
# sidx %*% a
# ai_j...i=sector, j = order
# bi_j...
# ci
if( VAR ){
if(is.null(VAR.opt$lag)) modelinc[1] = 1 else modelinc[1] = as.integer( VAR.opt$lag )
if(!is.null(VAR.opt$external.regressors)){
if(!is.matrix(VAR.opt$external.regressors)) stop("\nexternal.regressors must be a matrix.")
modelinc[2] = dim(VAR.opt$external.regressors)[2]
modeldata$mexdata = VAR.opt$external.regressors
} else{
modeldata$mexdata = NULL
}
}
modelinc[10] = which(c("mvnorm", "mvt", "mvlaplace") == distribution)
modeldesc$distribution = distribution
if( !is(uspec, "uGARCHmultispec") ) stop("\ndccspec-->error: uspec must be a uGARCHmultispec object")
varmodel = list()
umodel = vector(mode ="list")
if( modelinc[1]>0 ){
varmodel$robust = VAR.opt$robust
varmodel$lag.max = VAR.opt$lag.max
varmodel$lag.criterion = VAR.opt$lag.criterion
varmodel$robust.control = VAR.opt$robust.control
umodel$modelinc = matrix(0, ncol = m, nrow = 21)
rownames(umodel$modelinc) = names(uspec@spec[[1]]@model$modelinc[1:21])
umodel$modeldesc = list()
umodel$vt = sapply(uspec@spec, FUN = function(x) x@model$modelinc[22])
umodel$modeldesc$vmodel = vector(mode = "character", length = m)
umodel$modeldesc$vsubmodel = vector(mode = "character", length = m)
umodel$start.pars = umodel$fixed.pars = vector(mode = "list", length = m)
umodel$modeldesc$distribution = vector(mode = "character", length = m)
umodel$modeldata = list()
umodel$modeldata$vexdata = vector(mode = "list", length = m)
for(i in 1:m){
# zero the mean equation since we are using VAR
umodel$modeldesc$vmodel[i] = uspec@spec[[i]]@model$modeldesc$vmodel
umodel$modeldesc$vsubmodel[i] = ifelse(is.null(uspec@spec[[i]]@model$modeldesc$vsubmodel),"GARCH",uspec@spec[[i]]@model$modeldesc$vsubmodel)
umodel$modeldesc$distribution[i] = uspec@spec[[i]]@model$modeldesc$distribution
umodel$modelinc[,i] = uspec@spec[[i]]@model$modelinc[1:21]
umodel$modelinc[1:6,i] = 0
umodel$modeldata$vexdata[[i]] = if(is.null(uspec@spec[[i]]@model$modeldata$vexdata)) NA else uspec@spec[[i]]@model$modeldata$vexdata
umodel$start.pars[[i]] = if(is.null(uspec@spec[[i]]@model$start.pars)) NA else uspec@spec[[i]]@model$start.pars
umodel$fixed.pars[[i]] = if(is.null(uspec@spec[[i]]@model$fixed.pars)) NA else uspec@spec[[i]]@model$fixed.pars
umodel$modeldata$mexdata[[i]] = NA
}
} else{
varmodel$lag.max = 1
varmodel$lag.criterion = "HQ"
umodel$modelinc = matrix(0, ncol = m, nrow = 21)
rownames(umodel$modelinc) = names(uspec@spec[[1]]@model$modelinc[1:21])
umodel$modeldesc = list()
umodel$vt = sapply(uspec@spec, FUN = function(x) x@model$modelinc[22])
umodel$modeldesc$vmodel = vector(mode = "character", length = m)
umodel$modeldesc$vsubmodel = vector(mode = "character", length = m)
umodel$start.pars = umodel$fixed.pars = vector(mode = "list", length = m)
umodel$modeldesc$distribution = vector(mode = "character", length = m)
umodel$modeldata = list()
umodel$modeldata$mexdata = vector(mode = "list", length = m)
umodel$modeldata$vexdata = vector(mode = "list", length = m)
for(i in 1:m){
umodel$modeldesc$vmodel[i] = uspec@spec[[i]]@model$modeldesc$vmodel
umodel$modeldesc$vsubmodel[i] = ifelse(is.null(uspec@spec[[i]]@model$modeldesc$vsubmodel),"GARCH",uspec@spec[[i]]@model$modeldesc$vsubmodel)
umodel$modeldesc$distribution[i] = uspec@spec[[i]]@model$modeldesc$distribution
umodel$modelinc[,i] = uspec@spec[[i]]@model$modelinc[1:21]
umodel$modeldata$mexdata[[i]] = if(is.null(uspec@spec[[i]]@model$modeldata$mexdata)) NA else uspec@spec[[i]]@model$modeldata$mexdata
umodel$modeldata$vexdata[[i]] = if(is.null(uspec@spec[[i]]@model$modeldata$vexdata)) NA else uspec@spec[[i]]@model$modeldata$vexdata
umodel$start.pars[[i]] = if(is.null(uspec@spec[[i]]@model$start.pars)) NA else uspec@spec[[i]]@model$start.pars
umodel$fixed.pars[[i]] = if(is.null(uspec@spec[[i]]@model$fixed.pars)) NA else uspec@spec[[i]]@model$fixed.pars
}
}
maxgarchOrder = max( sapply(uspec@spec, FUN = function(x) x@model$maxOrder) )
# Note: We run varxfit with the postpad option using "constant" so that
# it resembles the ARMA-GARCH type functionality
if(modelinc[1]>0){
maxgarchOrder = max(c(maxgarchOrder, modelinc[1]))
}
################################################################
# pars list
pos = 1
pos.matrix = matrix(0, ncol = 3, nrow = 4)
colnames(pos.matrix) = c("start", "stop", "include")
rownames(pos.matrix) = c("fdcca", "fdccb", "mshape", "mskew")
# c("var", "mvmxreg", "fdccC", "fdccA", "fdccB", "mshape", "mskew", "aux", "aux", "aux")
if(modelinc[4]>0) pos.matrix[1,1:3] = c(pos, pos+(modelinc[4]*modelinc[3])-1, 1)
pos = pos + (modelinc[4]*modelinc[3])
if(modelinc[5]>0) pos.matrix[2,1:3] = c(pos, pos+(modelinc[5]*modelinc[3])-1, 1)
pos = pos + (modelinc[5]*modelinc[3])
if(modelinc[6]>0) pos.matrix[3,1:3] = c(pos, pos, 1)
pos = pos + modelinc[6]
if(modelinc[7]>0) pos.matrix[4,1:3] = c(pos, pos, 1)
# minimum parameter matrix = modelinc[3]+dcc1+dccb+mshape+mskew = modelinc[3]+4
mm = max(1, (modelinc[4]*modelinc[3]))+max(1, (modelinc[5]*modelinc[3])) + max(1, modelinc[6]) + max(1, modelinc[6])
pars = matrix(0, ncol = 6, nrow = mm)
colnames(pars) = c("Level", "Fixed", "Include", "Estimate", "LB", "UB")
pidx = matrix(NA, nrow = 4, ncol = 2)
colnames(pidx) = c("begin", "end")
rownames(pidx) = c("fdcca", "fdccb", "mshape", "mskew")
fixed.names = names(fixed.pars)
start.names = names(start.pars)
fixed.pars = unlist(fixed.pars)
start.pars = unlist(start.pars)
pn = 1
pnames = NULL
nx = 0
pidx[1,1] = 1
if(pos.matrix[1,3] == 1){
pn = length( seq(pos.matrix[1,1], pos.matrix[1,2], by = 1) )
xid = 0
for(i in 1:modelinc[4]){
for(j in 1:modelinc[3]){
nnx = paste("fdcca[g=", j,",P=", i, "]", sep="")
pars[(nx+xid+j), 1] = 0.01
sp = na.omit(match(nnx, start.names))
if(length(sp)==1) pars[(nx+xid+j), 1] = start.pars[sp]
pars[(nx+xid+j), 3] = 1
pars[(nx+xid+j), 5] = 0
pars[(nx+xid+j), 6] = 0.999
sp = na.omit(match(nnx, fixed.names))
if(length(sp)==1){
pars[(nx+xid+j), 1] = fixed.pars[sp]
pars[(nx+xid+j), 2] = 1
} else{
pars[(nx+xid+j), 4] = 1
}
pnames = c(pnames, nnx)
}
xid = xid + modelinc[3]
}
} else{
pnames = c(pnames, "fdcca")
}
pidx[1,2] = nx + pn
nx = nx + pn
pn = 1
pidx[2,1] = nx+1
if(pos.matrix[2,3] == 1){
pn = length( seq(pos.matrix[2,1], pos.matrix[2,2], by = 1) )
xid = 0
for(i in 1:modelinc[5]){
for(j in 1:modelinc[3]){
nnx = paste("fdccb[g=", j,",Q=", i, "]", sep="")
pars[(nx+xid+j), 1] = 0.9
sp = na.omit(match(nnx, start.names))
if(length(sp)==1) pars[(nx+xid+j), 1] = start.pars[sp]
pars[(nx+xid+j), 3] = 1
pars[(nx+xid+j), 5] = 0
pars[(nx+xid+j), 6] = 1
sp = na.omit(match(nnx, fixed.names))
if(length(sp)==1){
pars[(nx+xid+j), 1] = fixed.pars[sp]
pars[(nx+xid+j), 2] = 1
} else{
pars[(nx+xid+j), 4] = 1
}
pnames = c(pnames, nnx)
}
xid = xid+modelinc[3]
}
} else{
pnames = c(pnames, "fdccb")
}
pidx[2,2] = nx+pn
nx = nx + pn
pn = 1
pidx[3,1] = nx+1
if(modelinc[6]<=1){
if(pos.matrix[3,3]==1){
pars[nx+pn, 3] = 1
pars[nx+pn, 1] = 5
pars[(nx+pn), 5] = 4
pars[(nx+pn), 6] = 50
if(any(substr(start.names, 1, 6) == "mshape")) pars[nx+pn, 1] = start.pars["mshape"]
if(any(substr(fixed.names, 1, 6) == "mshape")){
pars[nx+pn,2] = 1
pars[nx+pn, 1] = fixed.pars["mshape"]
} else{
pars[nx+pn,4] = 1
}
}
pnames = c(pnames, "mshape")
} else{
if(pos.matrix[3,3] == 1){
pn = length( seq(pos.matrix[3,1], pos.matrix[3,2], by = 1) )
for(i in 1:pn){
nnx = paste("mshape", i, sep="")
pars[(nx+i), 1] = 5
if(any(substr(start.names, 1, nchar(nnx))==nnx)){
nix = which(substr(start.names, 1, nchar(nnx)) == nnx)
pars[(nx+i), 1] = start.pars[nix]
}
pars[(nx+i), 3] = 1
pars[(nx+i), 5] = 4
pars[(nx+i), 6] = 50
if(any(substr(fixed.names, 1, nchar(nnx))==nnx)){
nix = which(substr(fixed.names, 1, nchar(nnx)) == nnx)
pars[(nx+i), 1] = fixed.pars[nix]
pars[(nx+i), 2] = 1
} else{
pars[(nx+i), 4] = 1
}
pnames = c(pnames, nnx)
}
} else{
pnames = c(pnames, "mshape")
}
}
pidx[3,2] = nx+pn
nx = nx + pn
pn = 1
pidx[4,1] = nx+1
if(modelinc[7]<=1){
if(pos.matrix[4,3]==1){
pars[nx+pn, 3] = 1
pars[nx+pn, 1] = 0.5
pars[(nx+pn), 5] = -1
pars[(nx+pn), 6] = 1
if(any(substr(start.names, 1, 5) == "mskew")) pars[nx+pn, 1] = start.pars["mskew"]
if(any(substr(fixed.names, 1, 5) == "mskew")){
pars[nx+pn,2] = 1
pars[nx+pn, 1] = fixed.pars["mskew"]
} else{
pars[nx+pn,4] = 1
}
}
pnames = c(pnames, "mskew")
} else{
if(pos.matrix[4,3] == 1){
pn = length( seq(pos.matrix[4,1], pos.matrix[4,2], by = 1) )
for(i in 1:pn){
nnx = paste("mskew", i, sep="")
pars[(nx+i), 1] = 3
pars[(nx+i), 5] = -1
pars[(nx+i), 6] = 1
if(any(substr(start.names, 1, nchar(nnx))==nnx)){
nix = which(substr(start.names, 1, nchar(nnx)) == nnx)
pars[(nx+i), 1] = start.pars[nix]
}
if(any(substr(fixed.names, 1, nchar(nnx))==nnx)){
nix = which(substr(fixed.names, 1, nchar(nnx)) == nnx)
pars[(nx+i), 1] = fixed.pars[nix]
pars[(nx+i), 2] = 1
} else{
pars[(nx+i), 4] = 1
}
pnames = c(pnames, nnx)
}
} else{
pnames = c(pnames, "mskew")
}
}
pidx[4,2] = nx+pn
rownames(pars) = pnames
modeldesc$type = "2-step"
model = list(modelinc = modelinc, modeldesc = modeldesc, modeldata = modeldata, varmodel = varmodel,
pars = pars, start.pars = start.pars, fixed.pars = fixed.pars, maxgarchOrder = maxgarchOrder,
sidx = sidx, fdccindex = fdccindex, pos.matrix = pos.matrix, pidx = pidx)
model$DCC = "FDCC"
ans = new("DCCspec",
model = model,
umodel = umodel)
return(ans)
}
.fdccfit = function(spec, data, out.sample = 0, solver = "solnp",
solver.control = list(), fit.control = list(eval.se = TRUE,
stationarity = TRUE, scale = FALSE),
cluster = NULL, fit = NULL, VAR.fit = NULL, verbose = FALSE,
realizedVol = NULL, ...)
{
tic = Sys.time()
.eps = .Machine$double.eps
model = spec@model
umodel = spec@umodel
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
# allow first and second stage solvers to be defined
if(length(solver)==2){
garch.solver = solver[1]
solver = solver[2]
} else{
garch.solver = solver[1]
}
solver = match.arg(tolower(solver)[1], c("solnp", "nlminb", "lbfgs","gosolnp"))
#-----------------------------------------------------------------------------------
# Data Extraction
m = dim(data)[2]
if( is.null( colnames(data) ) ) cnames = paste("Asset_", 1:m, sep = "") else cnames = colnames(data)
colnames(umodel$modelinc) = cnames
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]
if( (n-n.start) < 100)
stop("\ndccfit-->error: function requires at least 100 data\n points to run\n")
data = xdata$data
index = xdata$index
period = xdata$period
# save the data to the model spec
model$modeldata$data = data
model$modeldata$index = index
model$modeldata$period = period
T = model$modeldata$T = n - n.start
model$modeldata$n.start = n.start
model$modeldata$asset.names = cnames
#-----------------------------------------------------------------------------------
# VAR model
if( model$modelinc[1]>0 ){
tmp = mvmean.varfit(model = model, data = data, VAR.fit = VAR.fit, T = T,
out.sample = out.sample, cluster = cluster)
model = tmp$model
zdata = tmp$zdata
mu = tmp$mu
varcoef = tmp$varcoef
p = tmp$p
N = tmp$N
} else{
zdata = data
ex = NULL
}
T = dim(zdata)[1] - out.sample
#-----------------------------------------------------------------------------------
# Univariate GARCH fit
# create a multispec list
mspec = .makemultispec(umodel$modelinc, umodel$modeldesc$vmodel, umodel$modeldesc$vsubmodel,
umodel$modeldata$mexdata, umodel$modeldata$vexdata, umodel$start.pars,
umodel$fixed.pars, umodel$vt)
if( !is.null(fit) && is(fit, "uGARCHmultifit") ){
# check VAR and fit:
if(model$modelinc[1]>0){
for(i in 1:m){
if(sum(fit@fit[[i]]@model$modelinc[1:6])>0)
stop("\nThe user supplied fit object has a non-null mean specification but VAR already chosen for mean filtration!!!")
}
}
fitlist = fit
if(spec@model$modelinc[1]>0) model$mu = mu else model$mu = fitted(fitlist)
model$residuals = res = residuals(fitlist)
model$sigma = sig = sigma(fitlist)
} else{
fitlist = multifit(multispec = mspec, data = xts(zdata, index), out.sample = n.start,
solver = garch.solver, solver.control = solver.control,
fit.control = ufit.control, cluster = cluster, realizedVol = realizedVol)
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 )
}
if(spec@model$modelinc[1]>0) model$mu = mu else model$mu = fitted(fitlist)
model$residuals = res = residuals(fitlist)
model$sigma = sig = sigma(fitlist)
}
stdresid = res/sig
Qbar = (t(stdresid) %*% stdresid)/nrow(stdresid)
# In the VAR model we postpadded the first lag values with a constant.
# In an ARMA model, the first lag values are also constant.
#-----------------------------------------------------------------------------------
# Create the Full Model Pars and Indices
modelinc = model$modelinc
# create full par matrix
tmp = .fullinc_fdcc(modelinc, umodel)
midx = tmp$midx
fdccAidx = tmp$fdccAidx
fdccBidx = tmp$fdccBidx
mpars = midx*0
# This is the estimated parameter index
eidx = .estindfn(midx, mspec, model$pars)
unipars = lapply(fitlist@fit, FUN = function(x) x@fit$ipars[x@fit$ipars[,3]==1,1])
if(is.list(unipars)){
for(i in 1:length(unipars)){
uninames = names(unipars[[i]])
mpars[uninames, i] = unipars[[i]]
}
} else{
uninames = rownames(unipars)
mpars[uninames, 1:NCOL(unipars)] = unipars
}
# add include pars from DCC spec (includes the fixed pars)
mpars[which(midx[,m+1]==1), m+1] = as.numeric( model$pars[model$pars[,3]==1,1] )
# DCC parameters
ipars = model$pars
LB = ipars[,5]
UB = ipars[,6]
estidx = as.logical( ipars[,4] )
npars = sum(estidx)
#-----------------------------------------------------------------------------------
H = sig^2
#-----------------------------------------------------------------------------------
# create a temporary environment to store values (deleted at end of function)
mgarchenv = new.env(hash = TRUE)
arglist = list()
arglist$mgarchenv = mgarchenv
arglist$verbose = verbose
arglist$cluster = cluster
arglist$eval.se = fit.control$eval.se
arglist$solver = solver
arglist$fit.control = fit.control
arglist$cnames = cnames
arglist$m = m
arglist$T = T
arglist$data = zdata
arglist$index = index
arglist$realizedVol = realizedVol
arglist$model = model
arglist$fitlist = fitlist
arglist$umodel = umodel
arglist$midx = midx
arglist$eidx = eidx
arglist$mpars = mpars
arglist$ipars = ipars
arglist$estidx = estidx
arglist$dccN = npars
arglist$stdresid = stdresid
arglist$H = H
arglist$Qbar = Qbar
# FDCC specific
arglist$fdccAidx = fdccAidx
arglist$fdccBidx = fdccBidx
arglist$sidx = model$sidx
#-----------------------------------------------------------------------------------
# Check for fixed parameters in DCC
if(any(ipars[,2]==1)){
if(npars == 0){
if(fit.control$eval.se==0) {
warning("\nfdccfit-->warning: all parameters fixed...returning dccfilter object instead\n")
xspex = spec
for(i in 1:m) xspex@umodel$fixed.pars[[i]] = as.list(fitlist@fit[[i]]@model$pars[fitlist@fit[[i]]@model$pars[,3]==1,1])
return(.fdccfilter(spec = xspex, data = xts(data, index), out.sample = out.sample, cluster = cluster, VAR.fit = VAR.fit,
realizedVol = realizedVol))
} else{
# if all parameters are fixed but we require standard errors, we
# skip the solver
use.solver = 0
ipars[ipars[,2]==1, 4] = 1
ipars[ipars[,2]==1, 2] = 0
arglist$pars = ipars
estidx = as.logical( ipars[,4] )
arglist$estidx = estidx
}
} else{
# with some parameters fixed we extract them (to be rejoined at end)
# so that they do not enter the solver
use.solver = 1
}
} else{
use.solver = 1
}
# start counter
assign("rmgarch_llh", 1, envir = mgarchenv)
# get
# calculate constraint length:
if(modelinc[3]==1){
ILB = 0
IUB = 1
} else{
if(modelinc[4]>0) xx = matrix(ipars[arglist$model$pidx["fdcca", 1]:arglist$model$pidx["fdcca", 2], 1], ncol = modelinc[4]) else xx = matrix(0, nrow = modelinc[3], ncol = modelinc[5])
if(modelinc[5]>0) yy = matrix(ipars[arglist$model$pidx["fdccb", 1]:arglist$model$pidx["fdccb", 2], 1], ncol = modelinc[5]) else yy = matrix(0, nrow = modelinc[3], ncol = modelinc[4])
# each pairwise combination
nc = NROW(expand.grid(xx, xx, yy, yy))
ILB = rep(0, nc)
IUB = rep(1, nc)
}
# Asymmetric Spec has different constraints
Ifn = .fdcccon
if( solver == "solnp" | solver == "gosolnp") fit.control$stationarity = FALSE else fit.control$stationarity = TRUE
arglist$fit.control = fit.control
if( use.solver )
{
arglist$returnType = "llh"
solution = switch(model$modeldesc$distribution,
mvnorm = .dccsolver(solver, pars = ipars[estidx, 1], fun = normal.fdccLLH1, Ifn, ILB,
IUB, gr = NULL, hessian = NULL, control = solver.control,
LB = ipars[estidx, 5], UB = ipars[estidx, 6], arglist = arglist))
#mvlaplace = .dccsolver(solver, pars = ipars[estidx, 1], fun = laplace.fdccLLH1, Ifn, ILB,
# IUB, gr = NULL, hessian = NULL, control = solver.control,
# LB = ipars[estidx, 5], UB = ipars[estidx, 6], arglist = arglist),
#mvt = .dccsolver(solver, pars = ipars[estidx, 1], fun = student.fdccLLH1, Ifn, ILB,
# IUB, gr = NULL, hessian = NULL, control = solver.control,
# LB = ipars[estidx, 5], UB = ipars[estidx, 6], arglist = arglist))
sol = solution$sol
hess = solution$hess
timer = Sys.time()-tic
convergence = sol$convergence
mpars[which(eidx[,(m+1)]==1, arr.ind = TRUE),m+1] = sol$pars
ipars[estidx, 1] = sol$pars
arglist$mpars = mpars
arglist$ipars = ipars
} else{
hess = NULL
timer = Sys.time()-tic
convergence = 0
sol = list()
sol$message = "all parameters fixed"
}
fit = list()
# check convergence else write message/return
if( convergence == 0 ){
fit = switch(model$modeldesc$distribution,
mvnorm = .dccmakefitmodel(garchmodel = "fdccnorm", f = normal.fdccLLH2,
arglist = arglist, timer = 0, message = sol$message, fname = "normal.fdccLLH2"))
#mvlaplace = .dccmakefitmodel(garchmodel = "fdcclaplace", f = laplace.fdccLLH2,
# arglist = arglist, timer = 0, message = sol$message, fname = "laplace.fdccLLH2"),
#mvt =.dccmakefitmodel(garchmodel = "fdccstudent", f = student.fdccLLH2,
# arglist = arglist, timer = 0, message = sol$message, fname = "student.fdccLLH2"))
fit$timer = Sys.time() - tic
} else{
fit$message = sol$message
fit$convergence = 1
}
# make model list to return some usefule information which
model$mpars = mpars
model$ipars = ipars
model$pars[,1] = ipars[,1]
model$midx = midx
model$eidx = eidx
model$umodel = umodel
fit$Qbar = Qbar
fit$realizedVol = realizedVol
ans = new("DCCfit",
mfit = fit,
model = model)
return(ans)
}
getfdccpars = function(pars, model){
modelinc = model$modelinc
pidx = model$pidx
sidx = model$sidx
if(modelinc[4]>0){
A = matrix(pars[pidx["fdcca",1]:pidx["fdcca",2],1], ncol = modelinc[4], byrow=FALSE)
A = apply(A, 2, function(x) sidx %*% x)
} else{
A = matrix(0)
}
if(modelinc[5]>0){
B = matrix(pars[pidx["fdccb",1]:pidx["fdccb",2],1], ncol = modelinc[5], byrow=FALSE)
B = apply(B, 2, function(x) sidx %*% x)
} else{
B = matrix(0)
}
# Because only the DCC (1,1) model is allowed, this is simplified:
BB = (B %*% t(B))
AA = (A %*% t(A))
C = matrix(1, NROW(BB), NCOL(BB)) - BB - AA
return(list(A = A, B = B, C = C))
}
.fdccfilter = function(spec, data, out.sample = 0, filter.control = list(n.old = NULL),
cluster = NULL, varcoef = NULL, realizedVol = NULL, ...)
{
tic = Sys.time()
model = spec@model
umodel = spec@umodel
n.old = filter.control$n.old
#-----------------------------------------------------------------------------------
# Data Extraction
m = dim(data)[2]
if( is.null( colnames(data) ) ) cnames = paste("Asset_", 1:m, sep = "") else cnames = colnames(data)
colnames(umodel$modelinc) = cnames
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]
if( (n-n.start) < 100)
stop("\ndccfilter-->error: function requires at least 100 data\n points to run\n")
data = xdata$data
index = xdata$index
period = xdata$period
# save the data to the model spec
model$modeldata$data = data
model$modeldata$index = index
model$modeldata$period = period
T = model$modeldata$T = n - n.start
model$modeldata$n.start = n.start
model$modeldata$asset.names = cnames
#-----------------------------------------------------------------------------------
# VAR model
if( spec@model$modelinc[1]>0 ){
tmp = mvmean.varfilter(model = model, data = data, varcoef = varcoef,
T = T, out.sample = out.sample)
model = tmp$model
zdata = tmp$zdata
mu = tmp$mu
p = tmp$p
N = tmp$N
} else{
zdata = data
ex = NULL
}
T = dim(zdata)[1] - out.sample
#-----------------------------------------------------------------------------------
if(is.null(filter.control$n.old)) n.old = T
#-----------------------------------------------------------------------------------
# Univariate GARCH filter
if(model$modelinc[1]>0){
for(i in 1:m){
if(sum(umodel$modelinc[1:6,i])>0)
stop("\nThe user supplied univariate spec object has a non-null mean specification but VAR already chosen for mean filtration!!!")
}
}
mspec = .makemultispec(umodel$modelinc, umodel$modeldesc$vmodel, umodel$modeldesc$vsubmodel,
umodel$modeldata$mexdata, umodel$modeldata$vexdata, umodel$start.pars,
umodel$fixed.pars, NULL)
filterlist = multifilter(multifitORspec = mspec, data = xts(zdata, index[1:nrow(zdata)]),
out.sample = out.sample, cluster = cluster, n.old = n.old,
realizedVol = realizedVol)
if(spec@model$modelinc[1]>0) model$mu = mu else model$mu = fitted(filterlist)
model$residuals = res = residuals(filterlist)
model$sigma = sig = sigma(filterlist)
stdresid = res/sig
if(is.null(filter.control$n.old)) dcc.old = dim(stdresid)[1] else dcc.old = n.old
#-----------------------------------------------------------------------------------
# Create the Full Model Pars and Indices
modelinc = model$modelinc
# create full par matrix
tmp = .fullinc_fdcc(modelinc, umodel)
midx = tmp$midx
fdccAidx = tmp$fdccAidx
fdccBidx = tmp$fdccBidx
mpars = midx*0
eidx = midx
# This is the estimated parameter index
unipars = sapply(filterlist@filter, FUN = function(x) x@filter$ipars[x@filter$ipars[,3]==1,1])
if(is.list(unipars)){
for(i in 1:length(unipars)){
uninames = names(unipars[[i]])
mpars[uninames, i] = unipars[[i]]
}
} else{
uninames = rownames(unipars)
mpars[uninames, 1:NCOL(unipars)] = unipars
}
# add include pars from DCC spec (includes the fixed pars)
mpars[which(midx[,m+1]==1, arr.ind = TRUE), m+1] = as.numeric( model$pars[model$pars[,3]==1,1] )
# DCC parameters
ipars = model$pars
estidx = as.logical( ipars[,3] )
npars = sum(estidx)
#-----------------------------------------------------------------------------------
Qbar = cov(stdresid[1:dcc.old, ])
H = sig^2
#-----------------------------------------------------------------------------------
mgarchenv = new.env(hash = TRUE)
arglist = list()
arglist$mgarchenv = mgarchenv
arglist$verbose = FALSE
arglist$cluster = cluster
arglist$filter.control = filter.control
arglist$cnames = cnames
arglist$m = m
arglist$T = T
arglist$n.old = n.old
arglist$dcc.old = dcc.old
arglist$data = zdata
arglist$index = index
arglist$realizedVol = realizedVol
arglist$model = model
arglist$filterlist = filterlist
arglist$umodel = umodel
arglist$midx = midx
arglist$eidx = eidx
arglist$mpars = mpars
arglist$ipars = ipars
arglist$estidx = estidx
arglist$dccN = npars
arglist$stdresid = stdresid
arglist$H = H
arglist$Qbar = Qbar
# FDCC specific
arglist$fdccAidx = fdccAidx
arglist$fdccBidx = fdccBidx
arglist$sidx = model$sidx
assign("rmgarch_llh", 1, envir = mgarchenv)
#-----------------------------------------------------------------------------------
filt = switch(model$modeldesc$distribution,
mvnorm = .dccmakefiltermodel(garchmodel = "fdccnorm", f = normalfilter.fdccLLH2,
arglist = arglist, timer = 0, message = 0, fname = "normalfilter.fdccLLH2"))
model$mpars = mpars
model$ipars = ipars
model$pars[,1] = ipars[,1]
model$midx = midx
model$eidx = eidx
model$umodel = umodel
filt$Qbar = Qbar
filt$realizedVol = realizedVol
filt$timer = Sys.time() - tic
ans = new("DCCfilter",
mfilter = filt,
model = model)
return(ans)
}
.fdccforecast = function(fit, n.ahead = 1, n.roll = 0,
external.forecasts = list(mregfor = NULL, vregfor = NULL), cluster = NULL, ...)
{
# checks for out.sample wrt n.roll
model = fit@model
modelinc = model$modelinc
ns = model$modeldata$n.start
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]
}
if( modelinc[1]>0 ){
if( modelinc[2] > 0 ){
if( is.null(external.forecasts$mregfor ) ){
warning("\nExternal Regressor Forecasts Matrix NULL...setting to zero...\n")
mregfor = matrix(0, ncol = modelinc[2], nrow = (n.roll + n.ahead) )
} else{
if( dim(mregfor)[2] != modelinc[2] ) 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.ahead!=1) stop("\ndccforecast-->error: only n.ahead==1 supported\n")
if( n.roll > ns ) stop("\ndccforecast-->error: n.roll greater than out.sample!", call. = FALSE)
#VARf = .varpredict(fit, n.ahead, n.roll, mregfor)$Mu
mu = varxforecast(X = fit@model$modeldata$data, Bcoef = model$varcoef, p = modelinc[1],
out.sample = ns, n.ahead = n.ahead, n.roll = n.roll, mregfor = mregfor)
} else{
mu = NULL
}
exf = external.forecasts
if( modelinc[1] > 0 ){
exf$mregfor = NULL
}
ans = .fdccforecastm(fit, n.ahead = n.ahead, n.roll = n.roll, external.forecasts = exf, cluster = cluster, ...)
if(modelinc[1]==0) mu = ans$mu
model$n.roll = n.roll
model$n.ahead = n.ahead
# 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) )
mforecast = list( H = ans$H, R = ans$R, Q = ans$Q, mu = mu )
ans = new("DCCforecast",
mforecast = mforecast,
model = model)
return( ans )
}
# 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
.fdccforecastm = function(fit, n.ahead = 1, n.roll = 10,
external.forecasts = list(mregfor = NULL, vregfor = NULL),
cluster = NULL, ...)
{
model = fit@model
modelinc = model$modelinc
umodel = model$umodel
m = dim(umodel$modelinc)[2]
ns = fit@model$modeldata$n.start
# first check whether the data needs to be filtered by VAR:
Data = fit@model$modeldata$data
if(modelinc[1]>0){
zdata = varxfilter(Data, p = model$modelinc[1], Bcoef = model$varcoef,
exogen = fit@model$modeldata$mexdata, postpad = c("constant"))$xresiduals
} else{
zdata = Data
}
fpars = lapply(1:m, FUN = function(i) fit@model$mpars[fit@model$midx[,i]==1,i])
mspec = .makemultispec(umodel$modelinc, umodel$modeldesc$vmodel, umodel$modeldesc$vsubmodel,
umodel$modeldata$mexdata, umodel$modeldata$vexdata, umodel$start.pars,
fpars, NULL)
filterlist = multifilter(multifitORspec = mspec, data = xts(zdata, fit@model$modeldata$index[1:nrow(zdata)]), out.sample = 0,
n.old = fit@model$modeldata$T, cluster = cluster, realizedVol = fit@mfit$realizedVol)
# all.equal(head(sigma(filterlist)), head(fit@model$sigma) )
n.roll = n.roll + 1
m = length(mspec@spec)
out.sample = fit@model$modeldata$n.start
mo = max(modelinc[4:5])
# we augmented n.roll before, now subtract
# forclist = multiforecast(multifitORspec = fitlist, n.ahead = n.ahead, n.roll = n.roll - 1)
forclist = multiforecast(multifitORspec = mspec, data = xts(zdata, fit@model$modeldata$index[1:nrow(zdata)]), n.ahead = n.ahead,
out.sample = ns, n.roll = n.roll - 1, external.forecasts = external.forecasts,
cluster = cluster, realizedVol = fit@mfit$realizedVol, ...)
# ToDo:
if(modelinc[1] == 0){
mu = array(NA, dim=c(n.ahead, m, n.roll))
f = lapply(forclist@forecast, function(x) fitted(x))
for(i in 1:n.roll) mu[,,i] = matrix( sapply( f, function(x) x[,i] ), ncol = m)
} else{
mu = NULL
}
sig = sigma(filterlist)
resid = residuals(filterlist)
stdresid = resid/sig
T = dim(fit@mfit$H)[3]
## 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)
f = lapply(forclist@forecast, function(x) sigma(x))
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( f, function(x) x[,i] ), ncol = m)
ans = .fdccforecastn(model, Qstart, Rstart, Hstart, xQbar, xstdresids, xfsig, n.ahead, mo)
Rtfor[[i]] = ans$Rtfor
Qtfor[[i]] = ans$Qtfor
Htfor[[i]] = ans$Htfor
Qstart = last( .abind_safe(Qstart, ans$Qtfor[, , 1]), mo )
Rstart = last( .abind_safe(Rstart, ans$Rtfor[, , 1]), mo )
Hstart = last( .abind_safe(Hstart, ans$Htfor[, , 1]), mo )
}
forc = list( H = Htfor, R = Rtfor, Q = Qtfor, mu = mu )
return(forc)
}
.fdccforecastn = function(model, Qstart, Rstart, Hstart, Qbar, stdresids, fsig, n.ahead, mo)
{
m = dim(Qbar)[1]
modelinc = model$modelinc
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]
tmp = getfdccpars(model$ipars, model)
A = tmp$A
B = tmp$B
BB = (B %*% t(B))
AA = (A %*% t(A))
C = matrix(1, NROW(BB), NCOL(BB)) - BB - AA
Q = C * Qbar
Qstar = diag( 1/sqrt( diag(Q) ) , m, m)
Rbar = Qstar %*% Q %*% t(Qstar)
for(i in 1:n.ahead){
if(i==1){
Qtfor[, , mo + i] = Q
if(modelinc[4]>0){
for(j in 1:modelinc[4]){
Qtfor[ , , mo + i] = Qtfor[ , , mo + i] + (A[,j]%*%t(A[,j])) * (stdresids[(mo + i - j), ] %*% t(stdresids[(mo + i - j), ]))
}
}
if(modelinc[5]>0){
for(j in 1:modelinc[5]){
Qtfor[ , , mo + i] = Qtfor[ , , mo + i] + (B[,j]%*%t(B[,j])) * Qtfor[ , , mo + i - j]
}
}
Qtmp = diag( 1/sqrt( diag(Qtfor[ , , mo + i]) ) , m, m)
Rtfor[ , , mo + i] = Qtmp %*% Qtfor[ , , mo + i] %*% t(Qtmp)
Dtmp = diag(fsig[i, ], m, m)
Htfor[ , , mo + i] = Dtmp %*% Rtfor[ , , mo + i] %*% Dtmp
} else{
Rtfor[ , , mo + i] = Rbar
if(modelinc[4]>0){
for(j in 1:modelinc[4]){
Rtfor[ , , mo + i] = Rtfor[ , , mo + i] + (A[,j]%*%t(A[,j])) * (Rtfor[ , , mo + i - j]-Rbar)
}
}
if(modelinc[5]>0){
for(j in 1:modelinc[5]){
Rtfor[ , , mo + i] = Rtfor[ , , mo + i] + (B[,j]%*%t(B[,j])) * (Rtfor[ , , mo + i - j]-Rbar)
}
}
Dtmp = diag(fsig[i, ], m, m)
Htfor[ , , mo + i] = Dtmp %*% Rtfor[ , , mo + i] %*% Dtmp
}
}
return( list( Rtfor = Rtfor[, , -(1:mo), drop = FALSE], Htfor = Htfor[, , -(1:mo), drop = FALSE],
Qtfor = Qtfor[, , -(1:mo), drop = FALSE] ) )
}
.fdccforecastn2 = function(model, Qstart, Rstart, Hstart, Qbar, stdresids, fsig, n.ahead, mo)
{
m = dim(Qbar)[1]
modelinc = model$modelinc
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]
tmp = getfdccpars(model$ipars, model)
A = tmp$A
B = tmp$B
BB = (B %*% t(B))
AA = (A %*% t(A))
C = matrix(1, NROW(BB), NCOL(BB)) - BB - AA
Q = C * Qbar
for(i in 1:n.ahead){
if(i==1){
Qtfor[, , mo + i] = Q
if(modelinc[4]>0){
for(j in 1:modelinc[4]){
Qtfor[ , , mo + i] = Qtfor[ , , mo + i] + (A[,j]%*%t(A[,j])) * (stdresids[(mo + i - j), ] %*% t(stdresids[(mo + i - j), ]))
}
}
if(modelinc[5]>0){
for(j in 1:modelinc[5]){
Qtfor[ , , mo + i] = Qtfor[ , , mo + i] + (B[,j]%*%t(B[,j])) * Qtfor[ , , mo + i - j]
}
}
Qtmp = diag( 1/sqrt( diag(Qtfor[ , , mo + i]) ) , m, m)
Rtfor[ , , mo + i] = Qtmp %*% Qtfor[ , , mo + i] %*% t(Qtmp)
Dtmp = diag(fsig[i, ], m, m)
Htfor[ , , mo + i] = Dtmp %*% Rtfor[ , , mo + i] %*% Dtmp
} else{
Qtfor[, , mo + i] = Q
if(modelinc[4]>0){
for(j in 1:modelinc[4]){
Qtfor[ , , mo + i] = Qtfor[ , , mo + i] + (A[,j]%*%t(A[,j])) * (Qtfor[ , , mo + i - j])
}
}
if(modelinc[5]>0){
for(j in 1:modelinc[5]){
Qtfor[ , , mo + i] = Qtfor[ , , mo + i] + (B[,j]%*%t(B[,j])) * (Qtfor[ , , mo + i - j])
}
}
Qtmp = diag( 1/sqrt( diag(Qtfor[ , , mo + i]) ) , m, m)
Rtfor[ , , mo + i] = Qtmp %*% Qtfor[ , , mo + i] %*% t(Qtmp)
Dtmp = diag(fsig[i, ], m, m)
Htfor[ , , mo + i] = Dtmp %*% Rtfor[ , , mo + i] %*% Dtmp
}
}
return( list( Rtfor = Rtfor[, , -(1:mo), drop = FALSE], Htfor = Htfor[, , -(1:mo), drop = FALSE],
Qtfor = Qtfor[, , -(1:mo), drop = FALSE] ) )
}
.fdccsim.fit = function(fitORspec, n.sim = 1000, n.start = 0, m.sim = 1,
startMethod = c("unconditional", "sample"), presigma = NULL,
preresiduals = NULL, prereturns = NULL, preQ = NULL, preZ = NULL,
Qbar = NULL, rseed = NULL, mexsimdata = NULL,
vexsimdata = NULL, cluster = NULL, VAR.fit = NULL, prerealized = NULL, ...)
{
fit = fitORspec
T = fit@model$modeldata$T
Data = fit@model$modeldata$data[1:T,]
m = dim(Data)[2]
mo = max(fit@model$modelinc[4:5])
mg = fit@model$maxgarchOrder
startMethod = startMethod[1]
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] )
}
model = fit@model
umodel = model$umodel
modelinc = model$modelinc
fpars = lapply(1:m, FUN = function(i) fit@model$mpars[fit@model$midx[,i]==1,i])
mspec = .makemultispec(umodel$modelinc, umodel$modeldesc$vmodel, umodel$modeldesc$vsubmodel,
umodel$modeldata$mexdata, umodel$modeldata$vexdata, umodel$start.pars,
fpars, NULL)
# Technically, we should pass an array of size (m, m, mo) for the
# sample but how many use a dccOrder > 1? (and it would mean modifying
# the C++ code to deal with arrays/lists for this input).
if(startMethod == "sample"){
if(is.null(preZ)){
preZ = matrix(tail(residuals(fit)/sigma(fit), mo), ncol = m)
} else{
preZ = matrix(tail(preZ, 1), ncol = m, nrow = mo, byrow = TRUE)
}
if(is.null(preQ)){
preQ = fit@mfit$Q[[length(fit@mfit$Q)]]
} else{
dcc.symcheck(preQ, m, d = NULL)
}
Rbar = preQ/(sqrt(diag(preQ)) %*% t(sqrt(diag(preQ))) )
} else{
if(is.null(preZ)){
preZ = matrix(0, ncol = m, nrow = mo)
} else{
preZ = matrix(tail(preZ, 1), ncol = m, nrow = mo, byrow = TRUE)
}
Rbar = cor(Data)
if(is.null(preQ)){
preQ = Rbar
} else{
dcc.symcheck(preQ, m, d = NULL)
Rbar = preQ/(sqrt(diag(preQ)) %*% t(sqrt(diag(preQ))) )
}
}
if(is.null(Qbar)){
Qbar = fit@mfit$Qbar
} else{
dcc.symcheck(Qbar, m, d = NULL)
}
uncv = sapply(mspec@spec, FUN = function(x) uncvariance(x))
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( dim(presigma)[1] != mg )
stop(paste("\ndccsim-->error: wrong row dimension for presigma (need ", mg, " rows.", sep = ""))
} else{
if(startMethod == "sample"){
mx = max(sapply(mspec@spec, FUN = function(x) x@model$maxOrder))
presigma = matrix(NA, ncol = m, nrow = mx)
tmp = last(fit@mfit$H, mx)
for(i in 1:mx) presigma[i,] = sqrt(diag(tmp[,,i]))
}
}
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( dim(preresiduals)[1] != mg )
stop(paste("\ndccsim-->error: wrong row dimension for preresiduals (need ", mg, " rows.", sep = ""))
} else{
if(startMethod == "sample"){
mx = max(sapply(mspec@spec, FUN = function(x) x@model$maxOrder))
preresiduals = matrix(NA, ncol = m, nrow = mx)
tmp = tail(fit@model$residuals, mx)
for(i in 1:mx) preresiduals[i,] = tmp[i,]
}
}
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.")
if( dim(prereturns)[1] != mg )
stop(paste("\ndccsim-->error: wrong row dimension for prereturns (need ", mg, " rows.", sep = ""))
} else{
if(startMethod == "sample"){
mx = max(sapply(mspec@spec, FUN = function(x) x@model$maxOrder))
prereturns = matrix(NA, ncol = m, nrow = mx)
tmp = tail(Data, mx)
for(i in 1:mx) prereturns[i,] = tmp[i,]
}
}
if(fit@model$umodel$modeldesc$vmodel[1]=="realGARCH"){
if( !is.null(prerealized) ){
if( !is.matrix(prerealized) )
stop("\ndccsim-->error: prerealized must be a matrix.")
if( dim(prerealized)[2] != m )
stop("\ndccsim-->error: wrong column dimension for prerealized.")
if( dim(prerealized)[1] != mg )
stop(paste("\ndccsim-->error: wrong row dimension for prerealized (need ", mg, " rows.", sep = ""))
} else{
if(startMethod == "sample"){
mx = max(sapply(mspec@spec, FUN = function(x) x@model$maxOrder))
prerealized = matrix(NA, ncol = m, nrow = mx)
tmp = tail(fit@mfit$realizedVol[1:T,], mx)
for(i in 1:mx) prerealized[i,] = tmp[i,]
}
}
} else{
mx = max(sapply(mspec@spec, FUN = function(x) x@model$maxOrder))
prerealized = matrix(NA, ncol = m, nrow = mx)
}
# switch distributions
if(fit@model$modeldesc$distribution == "mvnorm"){
if(length(rseed) == 1){
set.seed( rseed )
tmp = matrix(rnorm(m * (n.sim + n.start) * m.sim, 0, 1), ncol = m, nrow = n.sim+n.start)
z = array(NA, dim = c(n.sim + n.start + mo, m, m.sim))
for(i in 1:m.sim) z[,,i] = rbind(preZ, tmp)
} 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] = rbind(preZ, matrix(rnorm(m * (n.sim + n.start), 0, 1), nrow = n.sim + n.start, ncol = m))
}
}
} else if(fit@model$modeldesc$distribution == "mvlaplace"){
if(length(rseed) == 1){
set.seed( rseed )
tmp = matrix(rdist(distribution = "ged", n = m * (n.sim + n.start) * m.sim, 0, 1, shape = 1), ncol = m, nrow = n.sim+n.start)
z = array(NA, dim = c(n.sim + n.start + mo, m, m.sim))
for(i in 1:m.sim) z[,,i] = rbind(preZ, tmp)
} 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] = rbind(preZ, matrix(rdist(distribution = "ged", n = m * (n.sim + n.start), 0, 1, shape = 1), nrow = n.sim + n.start, ncol = m))
}
}
} else{
if(length(rseed) == 1){
set.seed( rseed )
tmp = matrix(rdist(distribution = "std", n = m * (n.sim + n.start) * m.sim, 0, 1, shape = rshape(fit)), ncol = m, nrow = n.sim+n.start)
z = array(NA, dim = c(n.sim + n.start + mo, m, m.sim))
for(i in 1:m.sim) z[,,i] = rbind(preZ, tmp)
} 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] = rbind(preZ, matrix(rdist(distribution = "std", n = m * (n.sim + n.start), 0, 1, shape = rshape(fit)), nrow = n.sim + n.start, ncol = m))
}
}
}
# ok now to expand rseed (but still allow for replication given initial seed)
if(length(rseed) == 1){
rseed = c(rseed, (1:m.sim)*(rseed+1))
}
tmp = getfdccpars(model$ipars, model)
A = tmp$A
B = tmp$B
BB = (B %*% t(B))
AA = (A %*% t(A))
C = matrix(1, NROW(BB), NCOL(BB)) - BB - AA
# model, Z, A, B, C, Qbar, preQ, Rbar, mo, n.sim, n.start, m, rseed
simRes = simX = simR = simQ = simH = simSeries = vector(mode = "list", length = m.sim)
if( !is.null(cluster) ){
simH = vector(mode = "list", length = m.sim)
simX = vector(mode = "list", length = m.sim)
clusterEvalQ(cluster, require(rmgarch))
clusterExport(cluster, c("model", "z", "A", "B", "C",
"preQ", "Rbar", "Qbar", "mo", "n.sim", "n.start", "m",
"rseed",".fdccsimf"), envir = environment())
mtmp = parLapply(cluster, as.list(1:m.sim), fun = function(j){
.fdccsimf(model, Z = z[,,j], A = A, B = B, C = C,
Qbar = Qbar, preQ = preQ, Rbar = Rbar, mo = 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))
#simZ = lapply(mtmp, FUN = function(x) x$Z)
simZ = vector(mode = "list", length = m)
for(i in 1:m) simZ[[i]] = sapply(mtmp, FUN = function(x) x$Z[,i])
clusterExport(cluster, c("fit", "n.sim", "n.start", "m.sim",
"startMethod", "simZ", "presigma", "preresiduals",
"prereturns", "mexsimdata", "vexsimdata","prerealized"),
envir = environment())
xtmp = parLapply(cluster, as.list(1:m), fun = function(j){
maxx = mspec@spec[[j]]@model$maxOrder;
htmp = ugarchpath(mspec@spec[[j]], n.sim = n.sim + n.start, n.start = 0, m.sim = m.sim,
custom.dist = list(name = "sample", distfit = matrix(simZ[[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) | model$modelinc[1]>0 ) NA else tail(prereturns[,j], maxx),
mexsimdata = if( model$modelinc[1]==0 ) mexsimdata[[j]] else NULL,
vexsimdata = vexsimdata[[j]], prerealized = tail(prerealized[,j], maxx));
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))
})
} 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){
.fdccsimf(model, Z = z[,,j], A = A, B = B, C = C,
Qbar = Qbar, preQ = preQ, Rbar = Rbar, mo = 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))
#simZ = lapply(mtmp, FUN = function(x) x$Z)
simZ = vector(mode = "list", length = m)
for(i in 1:m) simZ[[i]] = sapply(mtmp, FUN = function(x) x$Z[,i])
xtmp = lapply(as.list(1:m), FUN = function(j){
maxx = mspec@spec[[j]]@model$maxOrder;
htmp = ugarchpath(mspec@spec[[j]], n.sim = n.sim + n.start, n.start = 0, m.sim = m.sim,
custom.dist = list(name = "sample", distfit = matrix(simZ[[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) | model$modelinc[1]>0 ) NA else tail(prereturns[,j], maxx),
mexsimdata = if( model$modelinc[1]==0 ) mexsimdata[[j]] else NULL,
vexsimdata = vexsimdata[[j]], prerealized = tail(prerealized[,j], maxx));
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
}
if(model$modelinc[1]>0){
simxX = array(NA, dim = c(n.sim+n.start, m, m.sim))
for(i in 1:m.sim) simxX[,,i] = 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+n.start)
} else{
simxX = array(NA, dim = c(n.sim+n.start, m, m.sim))
for(i in 1:m.sim) simxX[,,i] = 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(tail(matrix(simxX[,,i], ncol = m), n.sim), nrow = n.sim)
}
if( model$modelinc[1]>0 ){
simRes = simX
simX = mvmean.varsim(model = model, Data = Data, res = simX,
mexsimdata = mexsimdata, prereturns = prereturns, m.sim = m.sim,
n.sim = n.sim, n.start = n.start, startMethod = startMethod,
cluster = cluster)
} 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$simZ = simZ
msim$simRes = simRes
#msim$Z = stdresid
msim$rseed = rseed
model$n.sim = n.sim
model$m.sim = m.sim
model$n.start = n.start
model$startMethod = startMethod[1]
ans = new("DCCsim",
msim = msim,
model = model)
return( ans )
}
.fdccsim.spec = function(fitORspec, n.sim = 1000, n.start = 0, m.sim = 1,
startMethod = c("unconditional", "sample"), presigma = NULL,
preresiduals = NULL, prereturns = NULL, preQ = NULL, preZ = NULL,
Qbar = NULL, rseed = NULL, mexsimdata = NULL,
vexsimdata = NULL, cluster = NULL, VAR.fit = NULL,
prerealized = NULL, ...)
{
spec = fitORspec
startMethod = startMethod[1]
# check that VAR.fit is included if used
if( spec@model$modelinc[1]>0 ){
if( is.null(VAR.fit) ) stop("\ndccsim-->error: VAR.fit must not be NULL for VAR method when calling dccsim using spec!", call. = FALSE)
}
model = spec@model
model$umodel = spec@umodel
m = dim(spec@umodel$modelinc)[2]
mo = max(spec@model$modelinc[4:5])
mg = spec@model$maxgarchOrder
model$modeldata$asset.names = paste("Asset", 1:m, sep = "")
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] )
}
if(is.null(preZ)){
preZ = matrix(0, ncol = m, nrow = mo)
} else{
preZ = matrix(tail(preZ, 1), ncol = m, nrow = mo, byrow = TRUE)
}
if(is.null(preQ)){
stop("\ndccsim-->error: preQ cannot be NULL when method uses spec!")
} else{
dcc.symcheck(preQ, m, d = NULL)
Rbar = preQ/(sqrt(diag(preQ)) %*% t(sqrt(diag(preQ))) )
}
if(is.null(Qbar)){
stop("\ndccsim-->error: Qbar cannot be NULL when method uses spec!")
} else{
dcc.symcheck(Qbar, m, d = NULL)
}
model = spec@model
umodel = spec@umodel
modelinc = model$modelinc
midx = .fullinc(modelinc, umodel)
mspec = .makemultispec(umodel$modelinc, umodel$modeldesc$vmodel, umodel$modeldesc$vsubmodel,
umodel$modeldata$mexdata, umodel$modeldata$vexdata, umodel$start.pars,
umodel$fixed.pars, NULL)
# uncvariance has uGARCHfit AND uGARCHspec dispatch methods
uncv = sapply(mspec@spec, FUN = function(x) uncvariance(x))
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( dim(presigma)[1] != mg )
stop(paste("\ndccsim-->error: wrong row dimension for presigma (need ", mg, " rows.", sep = ""))
}
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( dim(preresiduals)[1] != mg )
stop(paste("\ndccsim-->error: wrong row dimension for preresiduals (need ", mg, " rows.", sep = ""))
}
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.")
if( dim(prereturns)[1] != mg )
stop(paste("\ndccsim-->error: wrong row dimension for prereturns (need ", mg, " rows.", sep = ""))
}
if(spec@umodel$modeldesc$vmodel[1]=="realGARCH"){
if( !is.null(prerealized) ){
if( !is.matrix(prerealized) )
stop("\ndccsim-->error: prerealized must be a matrix.")
if( dim(prerealized)[2] != m )
stop("\ndccsim-->error: wrong column dimension for prerealized.")
if( dim(prerealized)[1] != mg )
stop(paste("\ndccsim-->error: wrong row dimension for prerealized (need ", mg, " rows.", sep = ""))
}
} else{
prerealized = matrix(NA, ncol = m, nrow = mg)
}
# we use the GED distribution with nu = 1 which corresponds to the Laplace case.
if(model$modeldesc$distribution == "mvnorm"){
if(length(rseed) == 1){
set.seed( rseed )
tmp = matrix(rnorm(m * (n.sim + n.start) * m.sim, 0, 1), ncol = m, nrow = n.sim+n.start)
z = array(NA, dim = c(n.sim + n.start + mo, m, m.sim))
for(i in 1:m.sim) z[,,i] = rbind(preZ, tmp)
} 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] = rbind(preZ, matrix(rnorm(m * (n.sim + n.start), 0, 1), nrow = n.sim + n.start, ncol = m))
}
}
} else if(model$modeldesc$distribution == "mvlaplace"){
if(length(rseed) == 1){
set.seed( rseed )
tmp = rdist(distribution = "ged", n = m * (n.sim + n.start) * m.sim, 0, 1, shape = 1)
z = array(NA, dim = c(n.sim + n.start + mo, m, m.sim))
for(i in 1:m.sim) z[,,i] = rbind(preZ, tmp)
} 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] = rbind(preZ, matrix(rdist(distribution = "ged", n = m * (n.sim + n.start), 0, 1, shape = 1), nrow = n.sim + n.start, ncol = m))
}
}
} else{
if(length(rseed) == 1){
set.seed( rseed )
tmp = rdist(distribution = "std", n = m * (n.sim + n.start) * m.sim, 0, 1, shape = model$mpars["mshape",1])
z = array(NA, dim = c(n.sim + n.start + mo, m, m.sim))
for(i in 1:m.sim) z[,,i] = rbind(preZ, tmp)
} 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] = rbind(preZ, matrix(rdist(distribution = "std", n = m * (n.sim + n.start), 0, 1, shape = model$mpars["mshape",1]), nrow = n.sim + n.start, ncol = m))
}
}
}
# ok now to expand rseed
if(length(rseed) == 1){
rseed = c(rseed, as.integer(runif(m.sim, 1, Sys.time())))
}
tmp = getfdccpars(model$ipars, model)
A = tmp$A
B = tmp$B
BB = (B %*% t(B))
AA = (A %*% t(A))
C = matrix(1, NROW(BB), NCOL(BB)) - BB - AA
simRes = simX = simR = simQ = simH = simSeries = vector(mode = "list", length = m.sim)
if( !is.null(cluster) ){
simH = vector(mode = "list", length = m.sim)
simX = vector(mode = "list", length = m.sim)
clusterEvalQ(cluster, require(rmgarch))
clusterExport(cluster, c("model", "z", "A", "B", "C", "preQ",
"Rbar", "Qbar", "mo", "n.sim", "n.start", "m",
"rseed",".fdccsimf"), envir = environment())
mtmp = parLapply(cluster, as.list(1:m.sim), fun = function(j){
.fdccsimf(model, Z = z[,,j], A = A, B = B, C = C,
Qbar = Qbar, preQ = preQ, Rbar = Rbar, mo = 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))
#simZ = lapply(mtmp, FUN = function(x) x$Z)
simZ = vector(mode = "list", length = m)
for(i in 1:m) simZ[[i]] = sapply(mtmp, FUN = function(x) x$Z[,i])
clusterExport(cluster, c("mspec", "n.sim", "n.start", "m.sim",
"startMethod", "simZ", "presigma", "preresiduals",
"prereturns", "mexsimdata", "vexsimdata","prerealized"),
envir = environment())
xtmp = parLapply(cluster, as.list(1:m), fun = function(j){
maxx = mspec@spec[[j]]@model$maxOrder;
htmp = ugarchpath(mspec@spec[[j]], n.sim = n.sim + n.start, n.start = 0, m.sim = m.sim,
custom.dist = list(name = "sample", distfit = matrix(simZ[[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) | model$modelinc[1]>0 ) NA else tail(prereturns[,j], maxx),
mexsimdata = if( model$modelinc[1]==0 ) mexsimdata[[j]] else NULL,
vexsimdata = vexsimdata[[j]], prerealized = tail(prerealized[,j], maxx))
h = matrix(tail(htmp@path$sigmaSim^2, n.sim), nrow = n.sim)
x = matrix(htmp@path$seriesSim, nrow = n.sim + n.start)
xres = matrix(htmp@path$residSim, nrow = n.sim + n.start)
return(list(h = h, x = x, xres = xres))
})
} 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){
.fdccsimf(model, Z = z[,,j], A = A, B = B, C = C,
Qbar = Qbar, preQ = preQ, Rbar = Rbar, mo = 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))
#simZ = lapply(mtmp, FUN = function(x) x$Z)
simZ = vector(mode = "list", length = m)
for(i in 1:m) simZ[[i]] = sapply(mtmp, FUN = function(x) x$Z[,i])
xtmp = lapply(as.list(1:m), FUN = function(j){
maxx = mspec@spec[[j]]@model$maxOrder;
htmp = ugarchpath(mspec@spec[[j]], n.sim = n.sim + n.start, n.start = 0, m.sim = m.sim,
custom.dist = list(name = "sample", distfit = matrix(simZ[[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) | model$modelinc[1]>0 ) NA else tail(prereturns[,j], maxx),
mexsimdata = if( model$modelinc[1]==0 ) mexsimdata[[j]] else NULL, vexsimdata = vexsimdata[[j]],
prerealized = tail(prerealized[,j], maxx));
h = matrix(tail(htmp@path$sigmaSim^2, n.sim), nrow = n.sim);
x = matrix(htmp@path$seriesSim, nrow = n.sim + n.start);
xres = matrix(htmp@path$residSim, nrow = n.sim + n.start);
return(list(h = h, x = x, xres = xres))
})
}
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
}
if(model$modelinc[1]>0){
simxX = array(NA, dim = c(n.sim+n.start, m, m.sim))
for(i in 1:m.sim) simxX[,,i] = 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+n.start)
} else{
simxX = array(NA, dim = c(n.sim+n.start, m, m.sim))
for(i in 1:m.sim) simxX[,,i] = 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(tail(matrix(simxX[,,i], ncol = m), n.sim), nrow = n.sim)
}
simxRes = array(NA, dim = c(n.sim, m, m.sim))
for(i in 1:n.sim) simxRes[i,,] = t(sapply(xtmp, FUN = function(x) as.numeric(x$xres[i,])))
simRes = vector(mode = "list", length = m.sim)
for(i in 1:m.sim) simRes[[i]] = matrix(simxRes[,,i], nrow = n.sim)
if( model$modelinc[1]>0 ){
model$varcoef = VAR.fit$Bcoef
Data = VAR.fit$xfitted
simRes = simX
simX = mvmean.varsim(model = model, Data = Data, res = simX,
mexsimdata = mexsimdata, prereturns = prereturns, m.sim = m.sim,
n.sim = n.sim, n.start = n.start, startMethod = startMethod,
cluster = cluster)
} 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$simZ = simZ
msim$simRes = simRes
#msim$Z = stdresid
msim$rseed = rseed
msim$model$Data = NULL
model$n.sim = n.sim
model$m.sim = m.sim
model$n.start = n.start
model$startMethod = "unconditional"
ans = new("DCCsim",
msim = msim,
model = model)
return( ans )
}
.fdccsimf = function(model, Z, A, B, C, Qbar, preQ, Rbar, mo, n.sim, n.start, m, rseed)
{
modelinc = model$modelinc
ipars = model$pars
idx = model$pidx
n = n.sim + n.start + mo
set.seed(rseed[1]+1)
# All distributions here make use of this since they are mixtures of normals
NZ = matrix(rnorm(m * (n.sim+n.start+mo)), nrow = n.sim + n.start + mo, ncol = m)
# SEXP model, SEXP A, SEXP B, SEXP C, SEXP Qbar, SEXP preQ, SEXP Rbar, SEXP Z, SEXP NZ
res = switch(model$modeldesc$distribution,
mvnorm = .Call("fdccsimmvn",
model = as.integer(modelinc),
A = A,
B = B,
C = C,
Qbar = as.matrix(Qbar),
preQ = as.matrix(preQ),
Rbar = as.matrix(Rbar),
Z = as.matrix(Z),
NZ = as.matrix(NZ),
PACKAGE = "rmgarch")
)
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[[1]][[i]]
}
ans = list( Q = Q, R = R, Z = res[[3]])
return( ans )
}
.fullinc_fdcc = function(modelinc, umodel){
m = dim(umodel$modelinc)[2]
# the par matrix will have the max parameter vis Asset (if an asset does not have
# the parameter it will be zero...for matrix estimation.
# the index places the solver estimates into their place in the matrix
###################################################################################
# Asset1 .... Asset2 .... AssetN ... ........ Joint
# mu
# ar
# ma
# arfima
# archm
# mxreg
# omega
# alpha
# beta
# gamma
# delta
# lambda
# eta1
# eta2
# vxreg
# skew
# shape
# ghlambda
# xi
# fdccC
# fdccA
# fdccB
# mshape
# mskew
vecmax = rep(0, 19)
names(vecmax) = rownames(umodel$modelinc[1:19,])
vecmax = apply(umodel$modelinc, 1, FUN = function(x) max(x) )
maxOrder = apply(umodel$modelinc, 2, FUN = function(x) max(c(x[2], x[3], x[8], x[9])))
sumv = 19 + sum(pmax(1, vecmax[c(2,3,6,8,9,10,11,12,13,15,16)])) - 11
tmpmat = matrix(0, ncol = m+1, nrow = sumv)
nx = 0
pnames = NULL
# mu
if(vecmax[1]>0){
tmpmat[1, 1:m] = umodel$modelinc[1,]
}
nx = nx + max(1, vecmax[1])
pnames = c(pnames, "mu")
# ar
if(vecmax[2]>0){
for(i in 1:vecmax[2]){
tmpmat[nx+i, 1:m] = as.integer( umodel$modelinc[2,] >= i)
pnames = c(pnames, paste("ar", i, sep = ""))
}
} else{
pnames = c(pnames, "ar")
}
nx = nx + max(1, vecmax[2])
# ma
if(vecmax[3]>0){
for(i in 1:vecmax[3]){
tmpmat[nx+i, 1:m] = as.integer( umodel$modelinc[3,] >= i)
pnames = c(pnames, paste("ma", i, sep = ""))
}
} else{
pnames = c(pnames, "ma")
}
nx = nx + max(1, vecmax[3])
# arfima
if(vecmax[4]>0){
tmpmat[nx+1, 1:m] = umodel$modelinc[4, ]
}
nx = nx + max(1, vecmax[4])
pnames = c(pnames, "arfima")
# archm no allowed
nx = nx + max(1, vecmax[5])
pnames = c(pnames, "archm")
# mxreg
if(vecmax[6]>0){
for(i in 1:vecmax[6]){
tmpmat[nx+i, 1:m] = as.integer(umodel$modelinc[6, ] >= i)
pnames = c(pnames, paste("mxreg", i, sep = ""))
}
} else{
pnames = c(pnames, "mxreg")
}
nx = nx + max(1, vecmax[6])
# omega
if(vecmax[7]>0){
tmpmat[nx+1, 1:m] = umodel$modelinc[7, ]
}
nx = nx + max(1, vecmax[7])
pnames = c(pnames, "omega")
# alpha
if(vecmax[8]>0){
for(i in 1:vecmax[8]){
tmpmat[nx+i, 1:m] = as.integer( umodel$modelinc[8, ] >= i)
pnames = c(pnames, paste("alpha", i, sep = ""))
}
} else{
pnames = c(pnames, "alpha")
}
nx = nx + max(1, vecmax[8])
# beta
if(vecmax[9]>0){
for(i in 1:vecmax[9]){
tmpmat[nx+i, 1:m] = as.integer( umodel$modelinc[9, ] >= i)
pnames = c(pnames, paste("beta", i, sep = ""))
}
} else{
pnames = c(pnames, "beta")
}
nx = nx + max(1, vecmax[9])
# gamma
if(vecmax[10]>0){
for(i in 1:vecmax[10]){
tmpmat[nx+i, 1:m] = as.integer( umodel$modelinc[10, ] >= i)
pnames = c(pnames, paste("gamma", i, sep = ""))
}
} else{
pnames = c(pnames, "gamma")
}
nx = nx + max(1, vecmax[10])
# eta1
if(vecmax[11]>0){
for(i in 1:vecmax[11]){
tmpmat[nx+i, 1:m] = as.integer( umodel$modelinc[11, ] >= i)
pnames = c(pnames, paste("eta1", i, sep = ""))
}
} else{
pnames = c(pnames, "eta1")
}
nx = nx + max(1, vecmax[11])
# eta2
if(vecmax[12]>0){
for(i in 1:vecmax[12]){
tmpmat[nx+i, 1:m] = as.integer( umodel$modelinc[12, ] >= i)
pnames = c(pnames, paste("eta2", i, sep = ""))
}
} else{
pnames = c(pnames, "eta2")
}
nx = nx + max(1, vecmax[12])
# delta
if(vecmax[13]>0){
tmpmat[nx+1, 1:m] = umodel$modelinc[13, ]
}
nx = nx + max(1, vecmax[13])
pnames = c(pnames, "delta")
# lambda
if(vecmax[14]>0){
tmpmat[nx+1, 1:m] = umodel$modelinc[14, ]
}
nx = nx + max(1, vecmax[14])
pnames = c(pnames, "lambda")
# vxreg
if(vecmax[15]>0){
for(i in 1:vecmax[15]){
tmpmat[nx+i, 1:m] = as.integer( umodel$modelinc[15, ] >= i)
pnames = c(pnames, paste("vxreg", i, sep = ""))
}
} else{
pnames = c(pnames, "vxreg")
}
nx = nx + max(1, vecmax[15])
# skew
if(vecmax[16]>0){
tmpmat[nx+1, 1:m] = umodel$modelinc[16, ]
}
nx = nx + max(1, vecmax[16])
pnames = c(pnames, "skew")
# shape
if(vecmax[17]>0){
tmpmat[nx+1, 1:m] = umodel$modelinc[17, ]
}
nx = nx + max(1, vecmax[17])
pnames = c(pnames, "shape")
# skew
if(vecmax[18]>0){
tmpmat[nx+1, 1:m] = umodel$modelinc[18, ]
}
nx = nx + max(1, vecmax[18])
pnames = c(pnames, "ghlambda")
if(vecmax[19]>0){
tmpmat[nx+1, 1:m] = umodel$modelinc[19, ]
}
nx = nx + max(1, vecmax[19])
pnames = c(pnames, "xi")
sumdcc = max(1, (modelinc[4]*modelinc[3]))+max(1, (modelinc[5]*modelinc[3])) +
max(1, modelinc[6]) + max(1, modelinc[6])
tmpmat = rbind(tmpmat, matrix(0, ncol = m+1, nrow = sumdcc))
# intercept
if(modelinc[4]>0){
xid = 0
for(i in 1:modelinc[4]){
for(j in 1:modelinc[3]){
tmpmat[(nx+xid+j), m+1] = 1
pnames = c(pnames, paste("fdcca[g=", j,",P=", i, "]", sep=""))
}
xid = xid + modelinc[3]
}
} else{
pnames = c(pnames, "fdcca")
}
fdccAidx = (nx+1):max( (nx+1), (nx+modelinc[3]*modelinc[4]) )
nx = nx + max(1, modelinc[4]*modelinc[3])
if(modelinc[5]>0){
xid = 0
for(i in 1:modelinc[5]){
for(j in 1:modelinc[3]){
tmpmat[(nx+xid+j), m+1] = 1
pnames = c(pnames, paste("fdccb[g=", j,",Q=", i, "]", sep=""))
}
xid = xid + modelinc[3]
}
} else{
pnames = c(pnames, "fdccb")
}
fdccBidx = (nx+1):max( (nx+1), (nx+modelinc[3]*modelinc[5]) )
nx = nx + max(1, modelinc[5]*modelinc[3])
# we allow for possibility of vector valued shape and skew parameters for future expansion
if(modelinc[6]>0){
for(i in 1:modelinc[6]){
tmpmat[nx+i, m+1] = 1
if(modelinc[6]>1) pnames = c(pnames, paste("mshape", i, sep = "")) else pnames = c(pnames, "mshape")
}
} else{
pnames = c(pnames, "mshape")
}
nx = nx + max(1, modelinc[6])
if(modelinc[7]>0){
for(i in 1:modelinc[7]){
tmpmat[nx+i, m+1] = 1
if(modelinc[7]>1) pnames = c(pnames, paste("mskew", i, sep = "")) else pnames = c(pnames, "mskew")
}
} else{
pnames = c(pnames, "mskew")
}
# NOW we have and indicator matrix
colnames(tmpmat) = c(paste("Asset", 1:m, sep = ""), "Joint")
rownames(tmpmat) = pnames
return(list(midx = tmpmat, fdccAidx = fdccAidx, fdccBidx = fdccBidx))
}
.fdcccon = function(pars, arglist){
ipars = arglist$ipars
estidx = arglist$estidx
idx = arglist$model$pidx
ipars[estidx, 1] = pars
modelinc = arglist$model$modelinc
if(modelinc[4]>0) fdccA = matrix(ipars[idx["fdcca", 1]:idx["fdcca", 2], 1], ncol = modelinc[4]) else fdccA = matrix(0, nrow = modelinc[3], ncol = modelinc[5])
if(modelinc[5]>0) fdccB = matrix(ipars[idx["fdccb", 1]:idx["fdccb", 2], 1], ncol = modelinc[5]) else fdccB = matrix(0, nrow = modelinc[3], ncol = modelinc[4])
if(modelinc[3]==1){
ans = fdccA + fdccB
} else{
# each pairwise combination (only allow DCC(1,1))
ans = apply(expand.grid(fdccA, fdccA, fdccB, fdccB), 1, function(x) x[1]*x[2] + x[3]*x[4])
}
return( as.numeric(ans) )
}
Try the rmgarch package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.