##### https://github.com/cran/nacopula
Tip revision: 161411b
nacopula.R
``````## Copyright (C) 2010 Marius Hofert and Martin Maechler
##
## This program is free software; you can redistribute it and/or modify it under
## Foundation; either version 3 of the License, or (at your option) any later
## version.
##
## This program 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.
##
## You should have received a copy of the GNU General Public License along with
## this program; if not, see <http://www.gnu.org/licenses/>.

#### Implementation of function evaluations and random number generation for
#### nested Archimedean copulas

##' Returns the copula density at u
##'
##' @title Density of nested Archimedean copulas
##' @param x nacopula
##' @param u argument of the copula x (can be a matrix)
##' @param log if TRUE the log-density is evaluated
##' @param ... further arguments passed to the copula specific 'dacopula' function;
##'   typically  'n.MC' (Monte Carlo size), but potentially more
##' @author Marius Hofert and Martin Maechler
dnacopula <- function(x, u, log=FALSE, ...) {
stopifnot(is(x, "outer_nacopula"))
if(length(x@childCops))
stop("currently, only Archimedean copulas are provided")
if(!is.matrix(u)) u <- rbind(u)
if((d <- ncol(u)) < 2) stop("u should be at least bivariate")
x@copula@dacopula(u, x@copula@theta, log=log, ...)
}

##' Returns the copula density at u. Generic algorithm
##'
##' @title Density of nested Archimedean copulas (generic form)
##' @param x nacopula
##' @param u argument of the copula x
##' @param n.MC if > 0 a Monte Carlo approach is applied with sample size equal
##'        to n.MC; otherwise the exact formula is used
##' @param log if TRUE the log-density is evaluated
##' @author Marius Hofert
dnacopulag <- function(x, u, n.MC=0, log = FALSE) {
stopifnot(is(x, "outer_nacopula"), 0 <= u, u <= 1)
if(length(x@childCops))
stop("currently, only Archimedean copulas are provided")
if(!is.matrix(u)) u <- rbind(u)
if((d <- ncol(u)) < 2) stop("u should be at least bivariate")
acop <- x@copula
th <- acop@theta
res <- rep(NaN,n <- nrow(u)) # density not defined on the margins
n01 <- apply(u,1,function(x) all(0 < x, x < 1)) # indices for which density has to be evaluated
if(any(n01)) {
u. <- u[n01,, drop=FALSE]
psiI <- acop@psiInv(u.,th)
psiID <- acop@psiInvD1abs(u.,th)
res[n01] <- acop@psiDabs(rowSums(psiI),theta = th,degree = d, n.MC = n.MC, log = log)
res[n01] <- if(log) res[n01] + rowSums(log(psiID)) else res[n01] * apply(psiID,1,prod)
}
res
}

##' Returns the copula value at u
##'
##' @title Evaluation of nested Archimedean copula
##' @param x nacopula
##' @param u argument of the copula x (can be a matrix)
##' @return f_x(u)
##' @author Marius Hofert, Martin Maechler
pnacopula <- function(x,u) {
if(!is.matrix(u)) u <- rbind(u)
stopifnot(ncol(u) >= dim(x)) # will be larger for children
C <- x@copula
th <- C@theta
res <- C@psi(rowSums(## use u[,j, drop=FALSE] for the direct components 'comp':
cbind(C@psiInv(u[,x@comp, drop=FALSE], theta=th),
## and recurse down for the children:
C@psiInv(unlist(lapply(x@childCops, pnacopula, u=u)), theta=th))),
theta=th)
## if u is a vector, res contains a names attribute, which we have to remove
## otherwise all.equal() in nac-experi.R fails
names(res) <- NULL
res
}

##' Compute the probability P[l < U <= u]  where U ~ copula x.
##'
##' @title Computing probabilities under nested Archimedean dependence
##' @param x outer nested archimedean copula
##' @param l d-vector of lower "integration" limits
##' @param u d-vector of upper "integration" limits
##' @return the probability that a random vector following the given copula
##'  falls in the hypercube with lower and upper corner l and u, respectively.
##' @author Marius Hofert, Martin Maechler
setGeneric("prob", function(x, l, u) standardGeneric("prob"))

setMethod("prob", signature(x ="outer_nacopula"),
function(x, l,u) {
d <- dim(x)
stopifnot(is.numeric(l), is.numeric(u),
length(u) == d, d == length(l),
0 <= l, l <= u, u <= 1)
if(d > 30)
stop("prob() for copula dimensions > 30 are not supported (yet)")
D <- 2^d
m <- 0:(D - 1)
## digitsBase() from package 'sfsmisc' {slightly simplified} :
## Purpose: Use binary representation of 0:N
## Author: Martin Maechler, Date:  Wed Dec  4 14:10:27 1991
II <- matrix(0, nrow = D, ncol = d)
for (i in d:1L) {
II[,i] <- m %% 2L + 1L
if (i > 1) m <- m %/% 2L
}
## Sign: the ("u","u",...,"u") case has +1; = c(2,2,...,2)
Sign <- c(1,-1)[1L + (- rowSums(II)) %% 2]
U <- array(cbind(l,u)[cbind(c(col(II)), c(II))], dim = dim(II))
sum(Sign * pnacopula(x, U))
})

##' Returns (n x d)-matrix of random variates
##'
##' @title Random number generation for nested Archimedean copulas
##' @param n number of vectors of random variates to generate
##' @param x outer_nacopula
##' @return matrix of random variates
##' @author Marius Hofert, Martin Maechler
rnacopula <- function(n, x, ...)
{
Cp <- x@copula			# outer copula
theta <- Cp@theta			# theta for outer copula
V0 <- Cp@V0(n,theta)		# generate V0's
childL <- lapply(x@childCops, rnchild, # <-- start recursion
theta0=theta,V0=V0,...)
dns <- length(x@comp)	 # dimension of the non-sectorial part
r <- matrix(rexp(n*dns), n, dns) # generate the non-sectorial part
## put pieces together
mat <- Cp@psi(r/V0, theta=theta)	# transform
mat <- cbind(mat, do.call(cbind,lapply(childL, `[[`, "U")))
## get correct sorting order:
j <- c(x@comp, unlist(lapply(childL, `[[`, "indCol")))
## extra check
stopifnot(length(j) == ncol(mat))
m <- mat[,order(j), drop=FALSE] # permute data and return
## extra checks
stopifnot(length(dm <- dim(m)) == 2, dm == dim(mat))
m
}

##' Returns a list with an (n x d)-matrix of random variates and a vector of
##' indices.
##'
##' @title Random number generation for children of nested Archimedean copulas
##' @param x nacopula
##' @param n number of vectors of random variates to generate
##' @param theta0 parameter theta0
##' @param V0 vector of V0's
##' @return list(U = matrix(*,n,d), indCol = vector of length d)
##' @author Marius Hofert, Martin Maechler
rnchild <- function(x, theta0, V0,...)
{
n <- length(V0)
Cp <- x@copula # inner copula
## Consistency checks -- for now {comment later} :
stopifnot(is(Cp, "acopula"), is.numeric(n), n == as.integer(n),
is.numeric(V0), length(V0) == n, is.numeric(theta0))
theta1 <- Cp@theta # theta_1 for inner copula
## generate V01's (only for one sector since the
## recursion in rnacopula() takes care of all sectors):
V01 <- Cp@V01(V0, theta0,theta1,...)
childL <- lapply(x@childCops, rnchild, # <-- recursion
theta0=theta1, V0=V01,...)
dns <- length(x@comp)	# dimension of the non-sectorial part
r <- matrix(rexp(n*dns), n, dns) # generate the non-sectorial part
## put pieces together: first own comp.s, then the children's :
mat <- Cp@psi(r/V01, theta1) # transform
if(length(childL) && length(U <- lapply(childL, `[[`, "U")))
mat <- cbind(mat, do.call(cbind, U))
## get correct sorting order:
j <- c(x@comp, unlist(lapply(childL, `[[`, "indCol")))
list(U = mat, indCol = j) # get list and return
}

if(FALSE) { # evaluate the following into your R session if you need debugging:
trace(rnacopula, browser, exit=browser, signature=signature(x ="outer_nacopula"))

debug(rnchild)
}

##' @title Constructor for outer_nacopula
##' @param family either character string (short or longer form of
##'	 copula family name) or an "acopula" family object
##' @param nacStructure a "formula" of the form C(th, comp, list(C(..), C(..)))
##' @return a valid outer_nacopula object
##' @author Martin Maechler
onacopula <- function(family, nacStructure) {
## , envir = ... , enclos=parent.frame() or
## , envir=environment()
##
### FIXME: base this on  onacopulaL() -- replacing nacStructure by nacList
### ----- : (1) replacing  C() with list() *AND* by
###         (2) wrapping the 3rd argument with list(.) if it's not already
### Use a *recursive* function like this one (but add the "list(.)" wrapping:
###    repC <- function(e) { sC <- as.symbol("C"); if(identical(e[[1]], sC)) e[[1]] <- as.symbol("list"); if(length(e) == 4) e[[4]] <- repC(e[[4]]); e}
nacl <- substitute(nacStructure)
stopifnot(identical(nacl[[1]], as.symbol("C")))
COP <- getAcop(family)
nacl[[1]] <- as.symbol("oC")
### does not work ..>>>>>>>>>>> we should use onacopulaL() inside functions! <<<<
## needed, e.g., when onacopula() is called from inside a user function:
## for(j in 2:3) if(is.language(nacl[[j]]))
##     nacl[[j]] <- eval(nacl[[j]], parent.frame())
## pframe <- parent.frame()
mkC <- function(cClass, a,b,c) {
if(missing(b) || length(b) == 0) b <- integer()
if(missing(c) || length(c) == 0) c <- list()
else if(length(c) == 1 && !is.list(c)) c <- list(c)
### does not work ...
## else if(!is.list(c)) {
##     c <- eval(c, pframe) ; stopifnot(is.list(c))
## }
## if(!is.numeric(a)) {
##     a <- eval(a, pframe) ; stopifnot(is.numeric(a), length(a) == 1)
## }
## if(!is.numeric(b)) {
##     b <- eval(b, pframe) ; stopifnot(is.numeric(b))
## }
else stopifnot(is.list(c))
stopifnot(is.numeric(a) || is.na(a), length(a) == 1, is.numeric(b))
if(any(sapply(c, class) != "nacopula"))
stop("third entry of 'nacStructure' must be NULL or list of 'C(..)' terms")
new(cClass, copula = setTheta(COP, a),
comp = as.integer(b), childCops = c)
}
C <- function(a,b,c) mkC("nacopula", a,b,c)
oC <- function(a,b,c) mkC("outer_nacopula", a,b,c)
eval(nacl)
}

##' @title Constructor for outer_nacopula - with list() input and *recursively*
##' @param family
##' @param nacList
##' @return
##' @author Martin Maechler
onacopulaL <- function(family, nacList) {
COP <- getAcop(family)
mkC <- function(cClass, abc) {
stopifnot(is.list(abc), 2 <= (nL <- length(abc)))
if(nL == 2) abc[[3]] <- list() else
if(nL > 3) stop("nacLists must be of length 3 (or 2)")
new(cClass, copula = setTheta(COP, abc[[1]]),
comp = as.integer(abc[[2]]),
childCops = lapply(abc[[3]], function(.) mkC("nacopula", .)))
}
mkC("outer_nacopula", nacList)
}

``````