#' Kernel estimation of vine copula densities
#'
#' The function estimates a vine copula density using kernel estimators for the
#' pair copulas (based on the \link{kdecopula} package).
#'
#' @param data (\eqn{n x d}) matrix of copula data (have to lie in \eqn{[0,1^d]}).
#' @param matrix R-Vine matrix (\eqn{n x d}) specifying the structure of the vine;
#' if \code{NA} (default) the structure selection heuristic of Dissman et al.
#' (2013) is applied.
#' @param method see \code{\link[kdecopula:kdecop]{kdecop}}.
#' @param renorm.iter see \code{\link[kdecopula:kdecop]{kdecop}}.
#' @param mult see \code{\link[kdecopula:kdecop]{kdecop}}.
#' @param test.level significance level for independence test. If you provide a
#' number in \eqn{[0, 1]}, an independence test
#' (\code{\link[VineCopula:BiCopIndTest]{BiCopIndTest}}) will be performed for
#' each pair; if the null hypothesis of independence cannot be rejected, the
#' independence copula will be set for this pair. If \code{test.level = NA}
#' (default), no independence test will be performed.
#' @param trunc.level integer; the truncation level. All pair copulas in trees
#' above the truncation level will be set to independence.
#' @param treecrit criterion for structure selection; defaults to \code{"tau"}.
#' @param cores integer; if \code{cores > 1}, estimation will be parallized
#' within each tree (using \code{\link[foreach]{foreach}}).
#' @param info logical; if \code{TRUE}, additional information about the
#' estimate will be gathered (see \code{\link[kdecopula:kdecop]{kdecop}}).
#'
#' @return An object of class \code{kdevinecop}. That is, a list containing
#' \item{T1, T2, ...}{lists of the estimted pair copulas in each tree,}
#' \item{matrix}{the structure matrix of the vine,}
#' \item{info}{additional information about the fit (if \code{info = TRUE}).}
#'
#' @references
#' Nagler, T., Czado, C. (2016) \cr Evading the curse of
#' dimensionality in nonparametric density estimation with simplified vine
#' copulas. \cr \emph{Journal of Multivariate Analysis 151, 69-89
#' (doi:10.1016/j.jmva.2016.07.003)}
#'
#' Nagler, T., Schellhase, C. and Czado, C. (2017) \cr Nonparametric
#' estimation of simplified vine copula models: comparison of methods
#' arXiv:1701.00845
#'
#' 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.
#'
#' @seealso
#' \code{\link{dkdevinecop}},
#' \code{\link[kdecopula:kdecop]{kdecop}},
#' \code{\link[VineCopula:BiCopIndTest]{BiCopIndTest}},
#' \code{\link[foreach]{foreach}}
#'
#' @examples
#' data(wdbc, package = "kdecopula")
#' # rank-transform to copula data (margins are uniform)
#' u <- VineCopula::pobs(wdbc[, 5:7], ties = "average")
#' \dontshow{u <- u[1:30, ]}
#' fit <- kdevinecop(u) # estimate density
#' dkdevinecop(c(0.1, 0.1, 0.1), fit) # evaluate density estimate
#' contour(fit) # contour matrix (Gaussian scale)
#' pairs(rkdevinecop(500, fit)) # plot simulated data
#'
#' @importFrom kdecopula kdecop hkdecop
#' @importFrom VineCopula BiCopIndTest RVineMatrix TauMatrix
#' @importFrom foreach foreach %dopar%
#' @importFrom parallel detectCores makeCluster stopCluster
#' @importFrom doParallel registerDoParallel
#' @export
kdevinecop <- function(data, matrix = NA, method = "TLL2", renorm.iter = 3L,
mult = 1, test.level = NA, trunc.level = NA,
treecrit = "tau", cores = 1, info = FALSE) {
## adjust input
if (is.null(info))
info <- FALSE
if (is.null(matrix))
matrix <- NA
if (is.null(method))
method <- "TLL2"
if (is.null(mult))
mult <- 1
if (is.null(test.level))
test.level <- 1
if (is.na(test.level))
test.level <- 1
if (is.null(trunc.level))
trunc.level <- ncol(data)
if (is.na(trunc.level))
trunc.level <- ncol(data)
if (is.null(treecrit))
treecrit <- "tau"
if (is.na(treecrit))
treecrit <- "tau"
if (is.null(cores))
cores <- 1
data <- as.matrix(data)
data <- pobs(data, ties.method = "first")
matrix <- as.matrix(matrix)
## sanity checks
d <- ncol(data)
n <- nrow(data)
if (any(data >= 1) || any(data < 0))
stop("Data has be in the interval [0,1].")
if (n < 2)
stop("Number of observations has to be at least 2.")
if (d < 2)
stop("Dimension has to be at least 2.")
if (!(treecrit %in% c("tau", "AIC", "cAIC", "hoeffd")))
stop("'treecrit' not available; please choose either 'tau', 'AIC', 'cAIC', or 'hoeffd'")
## call structure selection routine if no matrix given
if (any(is.na(matrix)) & d > 2) {
return(structure_select2(data,
type = 0,
method = method,
mult = mult,
info = info,
struct.crit = treecrit,
test.level = test.level,
trunc.level = trunc.level,
renorm.iter = renorm.iter,
cores = cores))
} else if (any(is.na(matrix)) & d == 2) {
# set standard matrix for d = 2 (there is only on possible structure)
matrix <- matrix(c(2, 1, 0, 1), 2, 2)
} else if (nrow(matrix) == 1 && matrix == 1) { # select C-Vine
return(structure_select2(data,
type = 1,
method = method,
mult = mult,
info = info,
struct.crit = treecrit,
test.level = test.level,
trunc.level = trunc.level,
renorm.iter = renorm.iter,
cores = cores))
}
## check matrix if provided
if (nrow(matrix) != ncol(matrix))
stop("Structure matrix has to be quadratic.")
if (ncol(data) != ncol(matrix))
stop("Dimensions of data and matrix don't match.")
if (max(matrix) > nrow(matrix))
stop("Error in the structure matrix.")
## register parallel backend
if (cores != 1 | is.na(cores)) {
if (is.na(cores))
cores <- max(1, detectCores() - 1)
if (cores > 1) {
cl <- makeCluster(cores)
registerDoParallel(cl)
on.exit(try(stopCluster(), silent = TRUE))
}
}
##
M <- ToLowerTri(matrix)
Mold <- M
o <- diag(M)
M <- reorderRVineMatrix(M)
data <- data[, o[length(o):1]]
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))
llikv <- array(0, dim = c(d, d, n))
llik <- matrix(0, d, d)
effp <- matrix(0, d, d)
AIC <- matrix(0, d, d)
cAIC <- matrix(0, d, d)
BIC <- matrix(0, d, d)
nms <- matrix("", d - 1, d - 1)
V <- list()
V$direct <- array(NA, dim = c(d, d, n))
V$indirect <- array(NA, dim = c(d, d, n))
V$direct[d, , ] <- t(data[, d:1])
## For independence pair-copulas
indepinfo <- list(effp = 0,
likvalues = rep(1, n),
loglik = 0,
effp = 0,
AIC = 0,
cAIC = 0,
BIC = 0)
for (k in d:2) {
doEst <- function(i) {
if (k > i) {
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), ]
}
samples <- cbind(zr2, zr1)
indep <- ifelse(test.level < 1,
BiCopIndTest(zr2, zr1)$p.value >= test.level,
FALSE)
if (trunc.level <= (d-k))
indep <- TRUE
if (indep) {
cfit <- list()
if (info)
cfit$info <- indepinfo
class(cfit) <- c("kdecopula", "indep.copula")
} else {
cfit <- kdecop(samples,
mult = mult,
method = method,
renorm.iter = renorm.iter,
info = info)
}
hfit <- list()
direct <- indirect <- NULL
if (CondDistr$direct[k - 1, i]) {
direct <- hkdecop(samples,
obj = cfit,
cond.var = 1L)
}
if (CondDistr$indirect[k - 1, i]) {
indirect <- hkdecop(samples,
obj = cfit,
cond.var = 2L)
}
names <- naming(Mold[c(i, k:d), i])
res.ki <- list(c = cfit, name = names)
return(list(direct = direct,
indirect = indirect,
res.ki = res.ki))
} else {
return(NULL)
}
}
res.k <- if (cores > 1) {
foreach(i = 1:(k-1),
.export = c("naming"),
.packages = c("kdevine", "kdecopula")) %dopar% doEst(i)
} else lapply(1:(k-1), doEst)
## clean up and finalize
for (i in 1:(d-1)) {
nums <- Mold[c(i, k:d), i]
#nums[1:2] <- nums[2:1]
name <- naming(nums)
if (any(extract_nums(res.k) == name)) {
ind <- which(extract_nums(res.k) == name)
res.ki <- res.k[[ind]]
res[[d + 1 - k]][[i]] <- res.ki$res.ki
if (!is.null(res.ki$direct)) {
V$direct[k - 1, i, ] <- res.ki$direct
}
if (!is.null(res.ki$indirect)) {
V$indirect[k - 1, i, ] <- res.ki$indirect
}
if (info) {
cfit <- res.ki$res.ki$c
llikv[k, i, ] <- log(cfit$info$likvalues)
llik[k, i] <- cfit$info$loglik
effp[k, i] <- cfit$info$effp
AIC[k, i] <- -2 * cfit$info$loglik + 2 * effp[k, i]
cAIC[k, i] <-
AIC[k, i] + (2 * effp[k, i] * (effp[k, i] + 1)) /
(n - effp[k, i] - 1)
BIC[k, i] <- -2 * cfit$info$loglik + log(n) * effp[k, i]
}
}
} # end i = 1:(d-1)
} # end k = d:2
## finalize results
res[[d]] <- data
res[[d + 1]] <- Mold
res[[d + 2]] <- list(llikv = apply(llikv, 3, sum),
loglik = sum(llik),
pair.loglik = llik,
effp = sum(effp),
pair.effp = effp,
AIC = sum(AIC),
pair.AIC = AIC,
cAIC = sum(AIC) +
(2 * sum(effp) * (sum(effp) + 1)) /
(n - sum(effp) - 1),
pair.cAIC = cAIC,
BIC = sum(BIC),
pair.BIC = BIC)
names(res) <- vapply(1:(d - 1), function(x) paste("T", x, sep = ""), "")
names(res)[d:(d + 2)] <- c("data", "matrix", "info")
if (!info)
res[[d + 2]] <- NULL
## return results
class(res) <- "kdevinecop"
res
}
## required functions from VineCopula package ----------------------------------
normalizeRVineMatrix <- function(RVM) {
oldOrder <- diag(RVM$Matrix)
Matrix <- reorderRVineMatrix(RVM$Matrix)
return(RVineMatrix(Matrix,
RVM$family,
RVM$par,
RVM$par2,
names = rev(RVM$names[oldOrder])))
}
reorderRVineMatrix <- function(Matrix) {
oldOrder <- diag(Matrix)
O <- apply(t(1:nrow(Matrix)), 2, "==", Matrix)
for (i in 1:nrow(Matrix)) {
Matrix[O[, oldOrder[i]]] <- nrow(Matrix) - i + 1
}
return(Matrix)
}
createMaxMat <- function(Matrix) {
if (dim(Matrix)[1] != dim(Matrix)[2])
stop("Structure matrix has to be quadratic.")
MaxMat <- reorderRVineMatrix(Matrix)
n <- nrow(MaxMat)
for (j in 1:(n - 1)) {
for (i in (n - 1):j) {
MaxMat[i, j] <- max(MaxMat[i:(i + 1), j])
}
}
tMaxMat <- MaxMat
tMaxMat[is.na(tMaxMat)] <- 0
oldSort <- diag(Matrix)
oldSort <- oldSort[n:1]
for (i in 1:n) {
MaxMat[tMaxMat == i] <- oldSort[i]
}
return(MaxMat)
}
neededCondDistr <- function(Vine) {
if (dim(Vine)[1] != dim(Vine)[2])
stop("Structure matrix has to be quadratic.")
Vine <- reorderRVineMatrix(Vine)
MaxMat <- createMaxMat(Vine)
d <- nrow(Vine)
M <- list()
M$direct <- matrix(FALSE, d, d)
M$indirect <- matrix(FALSE, d, d)
M$direct[2:d, 1] <- TRUE
for (i in 2:(d - 1)) {
v <- d - i + 1
bw <- as.matrix(MaxMat[i:d, 1:(i - 1)]) == v
direct <- Vine[i:d, 1:(i - 1)] == v
M$indirect[i:d, i] <- apply(as.matrix(bw & (!direct)), 1, any)
M$direct[i:d, i] <- TRUE
M$direct[i, i] <- any(as.matrix(bw)[1, ] & as.matrix(direct)[1, ])
}
return(M)
}
ToLowerTri <- function(x) {
## only change matrix if not already lower triagonal
if(all(x[lower.tri(x)] == 0)) {
x[nrow(x):1, ncol(x):1]
} else {
x
}
}
