R/copula-fn.R

In rgarch: Flexible GARCH modelling in R

#################################################################################
##
##   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 )
}