R/drkdevinecop.R

In kdevine: Multivariate Kernel Density Estimation with Vine Copulas

#' Working with a \code{kdevinecop} object
#'
#' A vine copula density estimate (stored in a \code{kdevinecop} object)
#' can be evaluated on arbitrary points with \code{dkevinecop}. Furthermore,
#' you can simulate from the estimated density with \code{rkdevinecop}.
#'
#' @aliases dkevinecop rkdevinecop
#'
#' @param u \eqn{m x 2} matrix of evaluation points.
#' @param obj \code{kdevinecop} object.
#' @param stable logical; option for stabilizing the estimator: the estimated
#' pair copula density is cut off at \eqn{50}.
#'
#' @return A numeric vector of the density/cdf or a \eqn{n x 2} matrix of
#' simulated data.
#'
#' @author Thomas Nagler
#'
#' @seealso
#' \code{\link{kdevinecop}},
#' \code{\link[kdecopula:dkdecop]{dkdecop}},
#' \code{\link[kdecopula:rkdecop]{rkdecop}},
#' \code{\link[qrng:ghalton]{ghalton}}
#'
#' @references
#' Nagler, T., Czado, C. (2016) \cr
#' Evading the curse of dimensionality in nonparametric density estimation. \cr
#' Journal of Multivariate Analysis 151, 69-89 (doi:10.1016/j.jmva.2016.07.003)
#'
#' Dissmann, J., Brechmann, E. C., Czado, C., and Kurowicka, D. (2013). \cr
#' Selecting and estimating regular vine copulae and application to financial returns. \cr
#' Computational Statistics & Data Analysis, 59(0):52--69.
#'
#' @examples
#' data(wdbc, package = "kdecopula")                    # load data
#' u <- VineCopula::pobs(wdbc[, 5:7], ties = "average") # rank-transform
#'
#' fit <- kdevinecop(u)                # estimate density
#' dkdevinecop(c(0.1, 0.1, 0.1), fit)  # evaluate density estimate
#'
#' @importFrom kdecopula dkdecop hkdecop
#' @export
dkdevinecop <- function(u, obj, stable = FALSE) {
    stopifnot(is.numeric(u))
    stopifnot(ncol(u) == ncol(obj$matrix))
    if (any(u > 1) || any(u < 0))
        stop("Values have to be in the interval (0,1).")

    u <- as.matrix(u)
    if (ncol(u) == 1)
        u <- matrix(u, 1, nrow(u))
    n <- nrow(u)
    d <- ncol(u)
    if (ncol(u) != ncol(obj$matrix))
        stop("Dimensions of 'u' and 'matrix' do not match.")

    ##
    M <- as.matrix(obj$matrix)
    Mold <- M
    o <- diag(M)
    M <- reorderRVineMatrix(M)
    uev <- as.matrix(u[, o[length(o):1]])
    if (ncol(uev) == 1)
        uev <- matrix(uev, 1, nrow(uev))

    MaxMat <- createMaxMat(M)
    CondDistr <- neededCondDistr(M)

    ## initialize objects
    res <- as.list(numeric(d - 1))
    for (i in 1:(d - 1))
        res[[i]] <- as.list(numeric(d - i))
    V <- list()
    V$direct <- array(NA, dim = c(d, d, n))
    V$indirect <- array(NA, dim = c(d, d, n))
    V$direct[d, , ] <- t(uev[, d:1])

    val <- array(1, dim = c(d, d, n))

    for (k in d:2) {
        for (i in 1:(k - 1)) {
            # get names and match object
            name <- split_num(naming(Mold[c(i, k:d), i]))
            nums <- lapply(extract_nums(obj[[d - k + 1]]), split_num)
            match <- sapply(nums, function(x) all(name %in% x))
            if (any(match)) {
                ## get object and pseudo-observations
                res.ki <- obj[[d - k + 1]][[which(match)]]
                m <- MaxMat[k, i]
                zr1 <- V$direct[k, i, ]
                zr2 <- if (m == M[k, i]) {
                    V$direct[k, (d - m + 1), ]
                } else {
                    V$indirect[k, (d - m + 1), ]
                }
                ev <- cbind(zr2, zr1)
                cfit <- res.ki$c

                ## flip grid if arguments in object are in different order
                if (!all(name[1:2] == nums[[which(match)]][1:2]))
                    cfit$flip <- TRUE

                ## evaluate density and h-functions
                if (any(class(cfit) == "indep.copula")) {
                    val[k-1, i, ] <- rep(1, nrow(ev))
                } else {
                    val[k-1, i, ] <- dkdecop(ev,
                                             obj = cfit,
                                             stable = stable)
                }

                if (k > 2) {
                    if(CondDistr$direct[k - 1, i]) {
                        if (any(class(cfit) == "indep.copula")) {
                            V$direct[k - 1, i, ] <- ev[,2]
                        } else {
                            V$direct[k - 1, i, ] <- hkdecop(ev,
                                                            obj = cfit,
                                                            cond.var = 1L)
                        }
                    }

                    if(CondDistr$indirect[k - 1, i]) {
                        if (any(class(cfit) == "indep.copula")) {
                            V$indirect[k - 1, i, ] <- ev[,1]
                        } else {
                            V$indirect[k - 1, i, ] <- hkdecop(ev,
                                                              obj = cfit,
                                                              cond.var = 2L)
                        }
                    }

                }
            }
        }
    }
    out <- exp(apply(log(val), 3, sum))
    if (stable)
        out <- pmin(out, 10^(1 + d/2))
    out
}

#' @param n integer; number of observations.
#' @param U (optional) \eqn{n x d} matrix of independent uniform random
#'  variables.
#' @param quasi logical; the default (\code{FALSE}) returns pseudo-random
#' numbers, use \code{TRUE} for quasi-random numbers (generalized Halton, see
#' \code{\link[qrng:ghalton]{ghalton}}).
#'
#' @rdname dkdevinecop
#' @importFrom kdecopula hkdecop
#' @importFrom stats runif
#' @export
rkdevinecop <- function(n, obj, U = NULL, quasi = FALSE) {
    n <- round(n)
    stopifnot(n > 0)
    stopifnot(is.logical(quasi))

    ## get structurte matrix and helper matrices,
    ## convert all to uper triangular form for simpler indexing
    M <- obj$matrix
    o <- diag(M)
    d <- length(o)
    Mold <- M[d:1, d:1]
    M <- reorderRVineMatrix(M)
    MaxMat <- createMaxMat(M)[d:1, d:1]
    CondDistr <- lapply(neededCondDistr(M), function(m) m[d:1, d:1])
    M <- M[d:1, d:1]

    ## initialize V matrices with independent uniform data
    ## (see Scherer, Mai (2014). "Simulating Copulas", Chapter 4)
    stopifnot(is.logical(quasi))
    if (!quasi) {
        # simulate independent uniform random variables
        if (is.null(U)) {
            U <- matrix(runif(n * d), ncol = d)
        } else {
            U <- matrix(U, ncol = d)
            U <- U[, rev(o), drop = FALSE]
        }
    } else {
        # generate quasi random numbers
        U <- ghalton(n, d = d)
    }

    Vdirect <- Vindirect <- array(dim = c(d, d, n))
    for (i in 1:d) {
        Vdirect[i, i, ] <- U[, i]
    }
    Vindirect[1, 1, ] <- Vdirect[1, 1, ]

    ## simulation algorithm
    for (i in 2:d) {  # loop through variables
        for (k in (i - 1):1) {  # loop over pairs involving variable i
            # find the pair copula object
            name <- split_num(naming(Mold[c(i, k:1), i]))
            nums <- lapply(extract_nums(obj[[k]]), split_num)
            match <- sapply(nums, function(x) all(name %in% x))
            res.ki <- obj[[k]][[which(match)]]

            # find the corresponding pseudo-data
            mm <- MaxMat[k, i]
            if (mm == M[k, i]) {
                zz <- Vdirect[k, mm, ]
            } else {
                zz <- Vindirect[k, mm, ]
            }

            # flip grid if arguments in object are in different order
            if (!all(name[1:2] == nums[[which(match)]][1:2]))
                res.ki$c$flip <- TRUE

            # invert h-function to create pseudo-data for next step
            Vdirect[k, i, ] <- hkdecop(cbind(zz, Vdirect[k + 1, i, ]),
                                       res.ki$c,
                                       cond.var = 1,
                                       inverse = TRUE)
            if ((i < d) & CondDistr$indirect[k + 1, i]) {
                Vindirect[k + 1, i, ] <- hkdecop(cbind(zz, Vdirect[k, i, ]),
                                                 res.ki$c,
                                                 cond.var = 2)
            }
        }
    }

    ## return result in initial ordering of the variables
    out <- t(Vdirect[1, , ])
    out[, sort(o[length(o):1], index.return = TRUE)$ix, drop = FALSE]
}