Nothing
#################################################################################
##
## R package rgarch by Alexios Ghalanos Copyright (C) 2008, 2009, 2010, 2011
## This file is part of the R package rgarch.
##
## The R package rgarch is free software: you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation, either version 3 of the License, or
## (at your option) any later version.
##
## The R package rgarch is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
## GNU General Public License for more details.
##
#################################################################################
##########################################################################################
# Some functions imported directly from QRMlib and mvtnorm and adjusted for use in this package
##########################################################################################
#------------------------------------------------------------------------------------
# Copula densities take the form:
# dMultivariate/Product(dUnivariate)
# log(dMultivariate) - Sum(log(dUnivariate))
.Spearman = function(data)
{
Rho = cor(apply(data, 2, rank))
return( Rho )
}
.Pconstruct <- function(theta){
n <- length(theta)
d <- (1 + sqrt(1+8*n))/2
A <- matrix(0,nrow=d,ncol=d)
A[lower.tri(A)] <- theta
diag(A) <- rep(1,d)
Q <- A %*% t(A)
P <- cov2cor(Q)
P
}
.Pdeconstruct <- function(P){
A <- t(chol(P))
Adiag <- diag(diag(A))
Astar <- solve(Adiag) %*% A
Astar[lower.tri(Astar)]
}
.vechR = function(Rho)
{
Rho[lower.tri(Rho)]
}
.ivechR = function(theta)
{
n = length(theta)
n = (1 + sqrt(1+8*n))/2
Rho = matrix(NA, n, n)
Rho[lower.tri(Rho)] = theta
Rho[upper.tri(Rho)] = theta
diag(Rho) = 1
Rho
}
.vech2 = function(Rho)
{
A = t(chol(Rho))
Adiag = diag(diag(A))
Astar = solve(Adiag) %*% A
return( Astar[lower.tri(Astar)] )
}
.ivech2 = function(theta)
{
n = length(theta)
d = (1 + sqrt(1+8*n))/2
A = matrix(0, d, d)
A[lower.tri(A)] = theta
diag(A) = rep(1,d)
Q = A %*% t(A)
Rho = cov2cor(Q)
return( Rho )
}
# about 3 times faster than "cor" function
.Kendall = function(data)
{
n <- dim(data)[1]
d <- dim(data)[2]
Rho <- matrix(0, nrow = d, ncol = d)
gr = matrix(combn(1:d, 2))
z = apply(gr, 2, FUN = function(x) Kendall(data[,x[1]], data[, x[2]])$tau)
Rho[lower.tri(Rho)] = z
Rho = Rho + t(Rho)
diag(Rho) = 1
nms <- dimnames(data)[[2]]
dimnames(Rho) <- list(nms, nms)
return( Rho )
}
.fit.kendall = function(Udata)
{
Rho = .Kendall(Udata)
Rbar = sin(pi * Rho/2)
n = dim(Udata)[2]
dimn = n * (n - 1) / 2
vRbar = Rbar[lower.tri(Rbar)]
X = diag(dimn)
estimate = lm(vRbar ~ X - 1)$coef
A = matrix(0, ncol = n, nrow = n)
A[lower.tri(A)] = as.numeric(estimate)
A = A + t(A)
diag(A) = 1
return( A )
}
########################################################################
# From corpcor package:
# Method by Higham 1988
.makeposdef = function(m)
{
d = dim(m)[1]
es = eigen(m)
esv = es$values
delta = 2*d*max(abs(esv))*.Machine$double.eps
tau = pmax(0, delta - esv)
dm = es$vectors %*% diag(tau, d) %*% t(es$vectors)
return( m + dm )
}
#------------------------------------------------------------------------------------
#------------------------------------------------------------------------------------
#######################################################################################
# copula likelihood functions (cpp coded)
#------------------------------------------------------------------------------------
.tvtcopulafn = function(pars, data, type = "LLH")
{
nu = pars[1]
alpha = pars[2]
beta = pars[3]
# this is the slowest part of the routine (student quantile)
# qt is vectorized but no discernable difference in using apply instead
Qdata = apply(data, 2, FUN = function(x) qt(p = x, df = nu))
n = dim(data)[1]
m = dim(data)[2]
Qbar = cov(Qdata)
Rbar = cor(Qdata)
ans = try(.Call(name = "dccCopulaStudent", Qbar = Qbar, U = Qdata, Rbar = Rbar, dcca = alpha,
dccb = beta, tnu = nu, dccorder = c(1,1,1), dccsum = alpha + beta), silent = TRUE)
if(inherits(ans, "try-error")){
ret = Inf
} else {
ret = switch(type, LLH = -ans[[3]], ALL = ans)
}
return(ret)
}
.tcopulafn = function(pars, data, Rbar, type = "LLH")
{
nu = pars[1]
Qdata = apply(data, 2, FUN = function(x) qt(p = x ,df = nu))
#ans = try(.Call(name = "staticCopulaStudent", U = Qdata, Rbar = Rbar, tnu = nu), silent = TRUE)
ans = dcopula.student(U = Qdata, Corr = Rbar, df = nu, logvalue = TRUE)
ret = switch(type,
LLH = -ans,
ALL = list(LLH = ans, Rbar = Rbar))
return(ret)
}
.tvgausscopulafn = function(pars, data, type = "LLH")
{
alpha = pars[2]
beta = pars[3]
Qdata = apply(data, 2, FUN = function(x) qnorm(p = x))
n = dim(data)[1]
m = dim(data)[2]
Qbar = cov(Qdata)
Rbar = cor(Qdata)
ans = try(.Call(name = "dccCopulaNormal", Qbar = Qbar, U = Qdata, Rbar = Rbar, dcca = alpha,
dccb = beta, dccorder = c(1,1,1), dccsum = alpha + beta), silent = TRUE)
if(inherits(ans, "try-error")){
ret = Inf
} else {
ret = switch(type, LLH = -ans[[3]], ALL = ans)
}
return(ret)
}
.gausscopulafn = function(pars, data, Rbar, type = "LLH")
{
Qdata = apply(data, 2, FUN = function(x) qnorm(p = x))
ans = dcopula.gauss(U = Qdata, Sigma = Rbar, logvalue = TRUE)
ret = switch(type,
LLH = -ans,
ALL = list(LLH = ans, Rbar = Rbar))
return(ret)
}
#------------------------------------------------------------------------------------
#####################################################################################
# Auxiliary Functions
#------------------------------------------------------------------------------------
.cor2cov = function(corr, sigmas)
{
V = (sigmas)%o%(sigmas)*corr
V
}
.cgarchcvar = function(f, pars, garchfit, model, env, fname)
{
garchNames = vector(mode = "character")
assetNames = get("cnames", envir = env)
spec = get("spec", envir = env)
trace = get("trace", envir = env)
dccNames = spec@mspec$optimization.model$modelnames
m = length(spec@uspec@spec)
for(i in 1:m){
cf = coef(garchfit@fit[[i]])
cfnames = names(cf)
garchNames = c(garchNames, paste(assetNames[i], cfnames, sep = ".") )
}
allnames = c(garchNames, dccNames)
# we need to generate the hessian and scores of the copula-garch model
# by partitioning the matrix (as in the DCC model).
g.pars = coef(garchfit)
g.n = length(unlist(g.pars))
#m = dim(g.pars)[2]
c.pars = pars
c.n = length(c.pars)
n.pars = g.n + c.n
A = matrix(0, n.pars, n.pars)
tmp = garchfit@fit[[1]]@fit$scores
N = dim(tmp)[1]
# partition the first g.n block into diagonal sublock holding the standard
# errors of each GARCH model
tidx = 1
# parallel
for(i in 1:m){
cv = garchfit@fit[[i]]@fit$cvar
widx = dim(cv)[1]
A[(tidx:(tidx + widx - 1)),(tidx:(tidx + widx - 1))] = solve(cv)
tidx = tidx + widx
}
numd = .cgarchderiv(f, g.pars, c.pars, garchfit, model, N, env, fname)
H = numd$H
likplus = numd$likplus
likminus = numd$likminus
hstep = numd$hstep
pars = numd$pars
A[(n.pars - c.n + 1):n.pars,] = H
jointscores = zeros(N, n.pars)
tidx = 1
for(i in 1:m){
cv = garchfit@fit[[i]]@fit$cvar
widx = dim(cv)[1]
jointscores[,(tidx:(tidx + widx - 1))] = garchfit@fit[[i]]@fit$scores
tidx = tidx + widx
}
sctemp = matrix(likplus[,(n.pars - c.n + 1):n.pars] - likminus[,(n.pars - c.n + 1):n.pars])
copulascores = matrix(NA, ncol = dim(sctemp)[2], nrow = dim(sctemp)[1])
sdtemp = 2*repmat(t(hstep[(n.pars - c.n + 1):n.pars]), N, 1)
for(i in 1:dim(sctemp)[2])
{
copulascores[, i] = sctemp[,i]/sdtemp[,i]
}
jointscores[,(n.pars - c.n + 1):n.pars] = copulascores
B = cov(jointscores)
A = A/N
copulavar = (solve(A) %*% B %*% solve(A))/N
se.coef = sqrt(diag(abs(copulavar)))
tval = as.numeric( pars/se.coef )
pval = 2* ( 1 - pnorm( abs( tval ) ) )
matcoef = matrix(NA, nrow = length(pars), ncol = 4)
matcoef[, 1] = pars
matcoef[, 2] = se.coef
matcoef[, 3] = tval
matcoef[, 4] = pval
dimnames(matcoef) = list(allnames, c(" Estimate",
" Std. Error", " t value", "Pr(>|t|)"))
sol = list(cvar = copulavar, scores = jointscores, A = A, n.pars = n.pars, pars = pars, matcoef = matcoef)
return(sol)
}
.cgarchderiv = function(f, g.pars, c.pars, garchfit, model, N, env, fname)
{
parallel = get("parallel", env)
parallel.control = get("parallel.control", env)
if(is.matrix(g.pars)) ng = length(g.pars) else ng = length(unlist(g.pars))
if(is.matrix(g.pars)) mg = dim(g.pars)[2] else mg = length((g.pars))
if(is.matrix(g.pars)) gpars = lapply(1:mg, FUN = function(i) g.pars[,i]) else gpars = g.pars
nc = length(c.pars)
ngidx = sapply(gpars, FUN = function(x) length(x))
ngidx = cbind(cumsum(c(1,ngidx[-length(ngidx)])), cumsum(ngidx))
ngidx = rbind(ngidx, matrix(c(ng+1, ng+nc), ncol = 2))
n = ng + nc
pars = c(unlist(g.pars), c.pars)
ures = get("ures", envir = env)
trace = get("trace", envir = env)
transformation = get("transformation", envir = env)
# Compute the stepsize (h)
fx = .cgarchfullLLH(f, pars, ngidx, garchfit, model, env)$llh
.eps = .Machine$double.eps
hstep = .eps ^ (1/3) * pmax( abs( pars ), 1 )
hpars = pars + hstep
hstep = hpars - pars
ee = Matrix(diag(hstep), sparse = TRUE)
g = vector(mode = "numeric", length = n)
# parallel
g = vector(mode = "numeric", length = n)
if( parallel ){
if( parallel.control$pkg == "multicore" ){
tmp = mclapply(1:n, FUN = function(i){
cat(paste("Evaluating StepValue ",i," out of ",n,"\n",sep=""))
tmpg = .cgarchfullLLH(f, as.numeric(pars) + ee[,i], ngidx, garchfit, model, env)$llh
return( tmpg )
}, mc.cores = parallel.control$cores)
g = as.numeric( unlist(tmp) )
} else{
sfInit(parallel = TRUE, cpus = parallel.control$cores)
sfExport("f", "pars", "ee", "ngidx", "garchfit", "model", "n", "env", local = TRUE)
tmp = sfLapply(as.list(1:n), fun = function(i){
cat(paste("Evaluating StepValue ",i," out of ",n,"\n",sep=""))
tmpg = .cgarchfullLLH(f, as.numeric(pars) + ee[,i], ngidx, garchfit, model, env)$llh
return( tmpg )
})
# sfStop() we'll stop it further down
g = as.numeric( unlist(tmp) )
}
} else{
tmp = lapply(as.list(1:n), FUN = function(i){
cat(paste("Evaluating StepValue ",i," out of ",n,"\n",sep=""))
tmpg = .cgarchfullLLH(f, as.numeric(pars) + ee[,i], ngidx, garchfit, model, env)$llh
return( tmpg )
})
g = as.numeric( unlist(tmp) )
}
H = hstep %*% t(hstep)
idx = 1
# Compute "double" forward and backward steps
# parallel
Nx = n*nc
if( parallel ){
if( parallel.control$pkg == "multicore" ){
tmp = mclapply(1:n, FUN = function(i){
Htmp = H
for(j in (n - Nx + 1):n){
if(i <= j){
Htmp[i, j] = (.cgarchfullLLH(f, as.numeric(pars) + ee[,i] + ee[,j], ngidx, garchfit, model, env)$llh -g[i] - g[j] + fx)/Htmp[i,j]
Htmp[j, i] = Htmp[i, j]
}
}
return(Htmp)
})
} else{
sfInit(parallel = TRUE, cpus = parallel.control$cores)
sfExport("f", "pars", "ee", "ngidx", "garchfit", "model", "fx", "n", "Nx", "env", local = TRUE)
tmp = sfLapply(as.list(1:n), fun = function(i){
Htmp = H
for(j in (n - Nx + 1):n){
if(i <= j){
Htmp[i, j] = (.cgarchfullLLH(f, as.numeric(pars) + ee[,i] + ee[,j], ngidx, garchfit, model, env)$llh -g[i] - g[j] + fx)/Htmp[i,j]
Htmp[j, i] = Htmp[i, j]
}
}
return(Htmp)
})
sfStop()
}
for(i in 1:n){
for(j in (n - Nx + 1):n){
if(i <= j){
H[i, j] = tmp[[i]][i, j]
H[j, i] = tmp[[i]][j, i]
}
}
}
} else{
tmp = lapply(as.list(1:n), FUN = function(i){
Htmp = H
for(j in (n - Nx + 1):n){
if(i <= j){
Htmp[i, j] = (.cgarchfullLLH(f, as.numeric(pars) + ee[,i] + ee[,j], ngidx, garchfit, model, env)$llh -g[i] - g[j] + fx)/Htmp[i,j]
Htmp[j, i] = Htmp[i, j]
}
}
return(Htmp)
})
for(i in 1:n){
for(j in (n - Nx + 1):n){
if(i <= j){
H[i, j] = tmp[[i]][i, j]
H[j, i] = H[i, j]
}
}
}
}
newH = H[(n - nc + 1):n, ]
hstep = pmax( abs( pars/2 ), 1e-2 ) * .eps^(1/3)
hplus = pars + hstep
hminus = pars - hstep
likplus = zeros(N, n)
likminus = zeros(N, n)
# parallel
for(i in (n - nc + 1):n){
hparsplus = pars
hparsplus[i] = hplus[i]
LHT = .cgarchfullLLH(f, hparsplus, ngidx, garchfit, model, env)$llhvec
likplus[,i] = LHT
hparsminus = pars
hparsminus[i] = hminus[i]
LHT = .cgarchfullLLH(f, hparsminus, ngidx, garchfit, model, env)$llhvec
likminus[,i] = LHT
}
sol = list(H = newH, likplus = likplus, likminus = likminus, hstep = hstep, pars = pars)
return(sol)
}
.cgarchfullLLH = function(f, pars, ngidx, garchfit, model, env)
{
# we need to simulate with different garch parameters
# gpars is passed as a list of pertrubed coefficients
# lapply(1:n, FUN = function(i) gpars[,i])
ng = dim(ngidx)[1]
cpars = pars[ngidx[ng,1]:ngidx[ng,2]]
gpars = pars[ngidx[1,1]:ngidx[(ng-1),2]]
m = ng - 1
n = length(gpars)
gpars = .relist(gpars, ngidx)
mflt = vector(mode = "list", length = m)
for(i in 1:m){
pr = gpars[[i]]
names(pr) = garchfit@fit[[i]]@model$modelnames
pr = as.list(pr)
spec = getspec(garchfit@fit[[i]])
spec@optimization.model$fixed.pars = pr
mflt[[i]] = ugarchfilter(spec = spec, data = garchfit@fit[[i]]@fit$data)
}
garch.llhvec = sapply(mflt, FUN = function(x) x@filter$log.likelihoods)
zres = sapply(mflt, FUN = function(x) x@filter$z)
ures = matrix(NA, ncol = 5, nrow = dim(zres)[1])
transformation = model$transformation
spd.control = model$spd.control
ans = switch(transformation,
parametric = .pparametric.filter(mflt, zres),
empirical = .pempirical(zres),
spd = .pspd(zres, spd.control))
if(transformation == "spd") {
ures = ans$ures
sfit = ans$sfit
} else{
ures = ans
sfit = NULL
}
# make tail adjustments in order to avoid problem with the quantile functions in
# optimization
if(any(ures > 0.99999)){
xn = which(ures > 0.99999)
ures[xn] = 0.99999
}
if(any(ures < eps)){
xn = which(ures < (1.5*eps))
ures[xn] = eps
}
cfilter = f(pars = cpars, data = ures, returnType = "ALL", garchenv = env)
c.llhvec = cfilter$lik
full.llhvec = apply(cbind(garch.llhvec, c.llhvec), 1, "sum")
sol = list()
sol$llhvec = full.llhvec
sol$llh = sum(full.llhvec)
return( sol )
}
.relist = function(gpars, idx)
{
n = dim(idx)[1]
idx = idx[-n,]
pars = lapply(1:(n-1), FUN = function(i) gpars[idx[i,1]:idx[i,2]])
return( pars )
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.