https://github.com/cran/nacopula
Raw File
Tip revision: 5bc804b03d066ef4a2ab9cf3af6f4f2df5bcda4e authored by Martin Maechler on 23 September 2011, 00:00 UTC
version 0.7-9-1
Tip revision: 5bc804b
nacopula.R
## Copyright (C) 2010 Marius Hofert and Martin Maechler
##
## This program 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.
##
## 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)
}

back to top