R/copula-likelihoods.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.
##
#################################################################################


copula.tvnormalLLH = function(pars, data, returnType = "llh", garchenv)
{
	fixed = get("fixed", garchenv)
	fixid = get("fixid", garchenv)
	mspec = get("mspec", garchenv)
	npar = get("npar", garchenv)
	fit.control = get("fit.control", garchenv)
	m = dim(data)[2]
	if(!is.null(fixed)){
		pall = vector(mode = "numeric", length = npar)
		pall[fixid] = fixed
		pall[-fixid] = pars
		pars = pall
	}
	T = dim(data)[1]
	omodel = get("omodel", garchenv)
	mo = omodel$maxOrder
	N = c(mo, T)
	dccOrder = mspec$dccOrder
	pos = omodel$pos.matrix
	dcca = dccb = skew = shape = 0
	if(pos[1,3] == 1) dcca = pars[pos[1,1]:pos[1,2]]
	if(pos[2,3] == 1) dccb = pars[pos[2,1]:pos[2,2]]
	if(pos[3,3] == 1) skew = pars[pos[3,1]:pos[3,2]]
	if(pos[4,3] == 1) shape = pars[pos[4,1]:pos[4,2]]
	
	#if( shape < 2.001 ) shape = 2.001
	
	Z = qnorm(data)
	dnormZ = dnorm(Z, log = TRUE)
	
	#xx = 0
	# for(i in 1:4) xx = xx + lgamma(0.5*(shape+1))  - log(sqrt(pi*shape)) - lgamma(shape/2) - 0.5*(shape+1)*sum( log(1+( (Z[2,i]^2)/shape) ) )
	# sum(dtZ[2,])
	Z = rbind( matrix(0, nrow = mo, ncol = m), Z )
	dnormZ = rbind( matrix(0, nrow = mo, ncol = m), dnormZ )
	
	Qbar = cov(Z)
	
	dccsum = sum(dcca) + sum(dccb)
	dccax = dcca
	dccbx = dccb
	dccorderx = c(dccOrder, mo)
	dccsumx = sum(dcca) + sum(dccb)
	
	if( fit.control$stationarity ){
		if( !is.na(dccsumx) && dccsumx >= 1 ) return(llh = get(".llh", garchenv) + 0.1*(abs(get(".llh", garchenv))))
	}
	
	res = .Call(name = "copulaNormalC2", Qbar = Qbar, Z = Z, dnormZ = dnormZ, 
			dcca = as.double(dccax), dccb = as.double(dccbx), 
			dccorder = as.numeric(dccorderx), dccsum = as.double(dccsumx), PACKAGE = "rgarch")
	
	llh = -1 * res[[3]]
	
	if(is.finite(llh) | !is.na(llh) | !is.nan(llh)){
		assign(".llh", llh, envir = garchenv)
	} else {
		llh = (get(".llh", garchenv) + 0.1*(abs(get(".llh", garchenv))))
	}
	
	ans = switch(returnType,
			lik = res[[2]][(mo+1):(mo+T)],
			llh = llh,
			ALL = list(lik = res[[2]][(mo+1):(mo+T)], llh = res[[3]], Rt = res[[4]][(mo+1):(mo+T)], Qt = res[[1]][(mo+1):(mo+T)]))
	return(ans)
}

copula.normalLLH = function(pars, data, returnType = "llh", garchenv)
{
	# .ivech(theta)
	fixed = get("fixed", garchenv)
	fixid = get("fixid", garchenv)
	mspec = get("mspec", garchenv)
	npar = get("npar", garchenv)
	method = mspec$method
	m = dim(data)[2]
	if(!is.null(fixed)){
		pall = vector(mode = "numeric", length = npar)
		pall[fixid] = fixed
		pall[-fixid] = pars
		pars = pall
	}
	T = dim(data)[1]
	omodel = get("omodel", garchenv)
	mo = omodel$maxOrder
	N = c(mo, T)
	dccOrder = mspec$dccOrder
	pos = omodel$pos.matrix
	
	Z = qnorm(data)
	dnormZ = dnorm(Z, log = TRUE)
	Z = rbind( matrix(0, nrow = mo, ncol = m), Z )
	dnormZ = rbind( matrix(0, nrow = mo, ncol = m), dnormZ )
	
	if(pos[5,3] == 1){
		Rbar = .Pconstruct(pars[pos[5,1]:pos[5,2]])
	} else{
		Rbar = cor(Z)
	}
	diag(Rbar) = 1
	
	eigR = eigen(Rbar)$values
	#print(eigR)
	if(any(is.complex(eigR)) || min(eigR) < 0){
		warning("\nNon p.s.d. covariance matrix...adjusting")
		Rbar = .makeposdef(Rbar)
	}
	
	res = .Call(name = "copulaNormalC1", Rbar = Rbar, Z = Z, dnormZ = dnormZ, PACKAGE = "rgarch")
	
	llh = -1 * res[[2]]
	
	if(is.finite(llh) | !is.na(llh) | !is.nan(llh)){
		assign(".llh", llh, envir = garchenv)
	} else {
		llh = (get(".llh", garchenv) + 0.1*(abs(get(".llh", garchenv))))
	}
	
	ans = switch(returnType,
			lik = res[[1]][(mo+1):(mo+T)],
			llh = llh,
			ALL = list(lik = res[[1]][(mo+1):(mo+T)], llh = res[[2]], Rt = Rbar))
	return(ans)
}

copula.tvstudentLLH = function(pars, data, returnType = "llh", garchenv)
{
	fixed = get("fixed", garchenv)
	fixid = get("fixid", garchenv)
	mspec = get("mspec", garchenv)
	npar = get("npar", garchenv)
	fit.control = get("fit.control", garchenv)
	
	m = dim(data)[2]
	if(!is.null(fixed)){
		pall = vector(mode = "numeric", length = npar)
		pall[fixid] = fixed
		pall[-fixid] = pars
		pars = pall
	}
	T = dim(data)[1]
	omodel = get("omodel", garchenv)
	mo = omodel$maxOrder
	N = c(mo, T)
	dccOrder = mspec$dccOrder
	pos = omodel$pos.matrix
	dcca = dccb = skew = shape = 0
	if(pos[1,3] == 1) dcca = pars[pos[1,1]:pos[1,2]]
	if(pos[2,3] == 1) dccb = pars[pos[2,1]:pos[2,2]]
	if(pos[3,3] == 1) skew = pars[pos[3,1]:pos[3,2]]
	if(pos[4,3] == 1) shape = pars[pos[4,1]:pos[4,2]]
	
	if( shape < 2.001 ) shape = 2.001
	
	Z = qt(data, df = shape)
	dtZ = dt(Z, df = shape, log = TRUE)
	
	#xx = 0
	# for(i in 1:4) xx = xx + lgamma(0.5*(shape+1))  - log(sqrt(pi*shape)) - lgamma(shape/2) - 0.5*(shape+1)*sum( log(1+( (Z[2,i]^2)/shape) ) )
	# sum(dtZ[2,])
	Z = rbind( matrix(0, nrow = mo, ncol = m), Z )
	dtZ = rbind( matrix(0, nrow = mo, ncol = m), dtZ )
	
	Qbar = cov(Z)
	
	dccsum = sum(dcca) + sum(dccb)
	dccax = dcca
	dccbx = dccb
	dccorderx = c(dccOrder, mo)
	dccsumx = sum(dcca) + sum(dccb)
	
	if( fit.control$stationarity ){
		if( !is.na(dccsumx) && dccsumx >= 1 ) return(llh = get(".llh", garchenv) + 0.1*(abs(get(".llh", garchenv))))
	}
	
	res = .Call(name = "copulaStudentC2", Qbar = Qbar, Z = Z, dtZ = dtZ, 
			dcca = as.double(dccax), dccb = as.double(dccbx), nu = as.double(shape), 
			dccorder = as.numeric(dccorderx), dccsum = as.double(dccsumx), PACKAGE = "rgarch")
	
	llh = -1 * res[[3]]
	
	if(is.finite(llh) | !is.na(llh) | !is.nan(llh)){
		assign(".llh", llh, envir = garchenv)
	} else {
		llh = (get(".llh", garchenv) + 0.1*(abs(get(".llh", garchenv))))
	}
	
	ans = switch(returnType,
			lik = res[[2]][(mo+1):(mo+T)],
			llh = llh,
			ALL = list(lik = res[[2]][(mo+1):(mo+T)], llh = res[[3]], Rt = res[[4]][(mo+1):(mo+T)], Qt = res[[1]][(mo+1):(mo+T)]))
	return(ans)
}

copula.studentLLH = function(pars, data, returnType = "llh", garchenv)
{
	# .ivech(theta)
	fixed = get("fixed", garchenv)
	fixid = get("fixid", garchenv)
	mspec = get("mspec", garchenv)
	npar = get("npar", garchenv)
	method = mspec$method
	m = dim(data)[2]
	if(!is.null(fixed)){
		pall = vector(mode = "numeric", length = npar)
		pall[fixid] = fixed
		pall[-fixid] = pars
		pars = pall
	}
	T = dim(data)[1]
	omodel = get("omodel", garchenv)
	mo = omodel$maxOrder
	N = c(mo, T)
	dccOrder = mspec$dccOrder
	pos = omodel$pos.matrix
	
	shape = 0
	if(pos[4,3] == 1) shape = pars[pos[4,1]:pos[4,2]]

	
	if( shape < 2.001 ) shape = 2.001
	
	Z = qt(data, df = shape)
	dtZ = dt(Z, df = shape, log = TRUE)
	Z = rbind( matrix(0, nrow = mo, ncol = m), Z )
	dtZ = rbind( matrix(0, nrow = mo, ncol = m), dtZ )
	
	if(pos[5,3] == 1){
		Rbar = .Pconstruct(pars[pos[5,1]:pos[5,2]])
	} else{
		Rtau = .Kendall(Z)
		Rbar = sin(pi * Rtau/2)
	}
	diag(Rbar) = 1
	
	eigR = eigen(Rbar)$values
	#print(eigR)
	if(any(is.complex(eigR)) || min(eigR) < 0){
		warning("\nNon p.s.d. covariance matrix...adjusting")
		Rbar = .makeposdef(Rbar)
	}
	
	res = .Call(name = "copulaStudentC1", Rbar = Rbar, Z = Z, dtZ = dtZ, nu = as.double(shape), PACKAGE = "rgarch")
	
	llh = -1 * res[[2]]
	
	if(is.finite(llh) | !is.na(llh) | !is.nan(llh)){
		assign(".llh", llh, envir = garchenv)
	} else {
		llh = (get(".llh", garchenv) + 0.1*(abs(get(".llh", garchenv))))
	}
	
	ans = switch(returnType,
			lik = res[[1]][(mo+1):(mo+T)],
			llh = llh,
			ALL = list(lik = res[[1]][(mo+1):(mo+T)], llh = res[[2]], Rt = Rbar))
	return(ans)
}