drkdevinecop.R
``````#' 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
#'
#' @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
#'
#' @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]
}
``````