Nothing
############ plot of the theoretical and empirical lambda-function
#' Lambda-Function (Plot) for Bivariate Copula Data
#'
#' This function plots/returns the lambda-function of given bivariate copula
#' data.
#'
#' If the family and parameter specification is stored in a [BiCop()]
#' object `obj`, the alternative versions \cr
#' \preformatted{BiCopLambda(obj, PLOT = TRUE, ...)}
#' and
#' \preformatted{BiCopLambda((u1, u2, obj, PLOT = TRUE, ...)}
#' can be used.
#'
#' @param u1,u2 Data vectors of equal length with values in \eqn{[0,1]} (default:
#' `u1` and `u2 = NULL`).
#' @param family An integer defining the bivariate copula family or indicating
#' the empirical lambda-function: \cr
#' `"emp"` = empirical lambda-function
#' (default) \cr
#' `1` = Gaussian copula; the theoretical lambda-function is simulated
#' (no closed formula available) \cr
#' `2` = Student-t copula; the theoretical lambda-function is simulated
#' (no closed formula available) \cr
#' `3` = Clayton copula \cr
#' `4` = Gumbel copula \cr
#' `5` = Frank copula \cr
#' `6` = Joe copula \cr
#' `7` = BB1 copula \cr
#' `8` = BB6 copula \cr
#' `9` = BB7 copula \cr
#' `10` = BB8
#' copula
#' @param par Copula parameter; if the empirical lambda-function is chosen,
#' `par = NULL` or `0` (default).
#' @param par2 Second copula parameter for t-, BB1, BB6, BB7 and BB8 copulas
#' (default: `par2 = 0`).
#' @param PLOT Logical; whether the results are plotted. If `PLOT =
#' FALSE`, the values \cr `empLambda` and/or `theoLambda` are
#' returned (see below; default: `PLOT = TRUE`).
#' @param ... Additional plot arguments.
#' @param obj `BiCop` object containing the family and parameter
#' specification.
#'
#' @return \item{empLambda}{If the empirical lambda-function is chosen and
#' `PLOT = FALSE`, a vector of the empirical lambda's is returned.}
#' \item{theoLambda}{If the theoretical lambda-function is chosen and
#' `PLOT = FALSE`, a vector of the theoretical lambda's is returned.}
#' @note The \eqn{\lambda}-function is characteristic for each bivariate copula
#' family and defined by Kendall's distribution function \eqn{K}:
#' \deqn{ \lambda(v,\theta) := v - K(v,\theta) } with
#' \deqn{K(v,\theta) := P(C_{\theta}(U_1,U_2) \leq v),\ \ v\in [0,1]. }{
#' K(v,\theta) := P(C_{\theta}(U_1,U_2) <= v), v \in [0,1]. }
#' For Archimedean copulas one has the following closed form expression
#' in terms of the generator function \eqn{\varphi}{\phi} of the copula
#' \eqn{C_{\theta}}:
#' \deqn{\lambda(v,\theta) = \frac{\varphi(v)}{\varphi '(v)}, }{
#' \lambda(v,\theta) = \phi(v) / \phi'(v), }
#' where \eqn{\varphi '}{\phi'} is the derivative of \eqn{\varphi}{\phi}. For
#' more details see Genest and Rivest (1993) or Schepsmeier (2010).
#'
#' For the bivariate Gaussian and Student-t copula no closed form expression for
#' the theoretical \eqn{\lambda}-function exists. Therefore it is simulated based
#' on samples of size 1000. For all other implemented copula families there are
#' closed form expressions available.
#'
#' The plot of the theoretical \eqn{\lambda}-function also shows the limits of
#' the \eqn{\lambda}-function corresponding to Kendall's tau \eqn{=0} and
#' Kendall's tau \eqn{=1} (\eqn{\lambda=0}).
#'
#' For rotated bivariate copulas one has to transform the input arguments
#' `u1` and/or `u2`. In particular, for copulas rotated by 90 degrees
#' `u1` has to be set to `1-u1`, for 270 degrees `u2` to
#' `1-u2` and for survival copulas `u1` and `u2` to `1-u1`
#' and `1-u2`, respectively. Then \eqn{\lambda}-functions for the
#' corresponding non-rotated copula families can be considered.
#'
#' @author Ulf Schepsmeier
#'
#' @seealso [BiCopMetaContour()], [BiCopKPlot()],
#' [BiCopChiPlot()], [BiCop()]
#'
#' @references Genest, C. and L.-P. Rivest (1993). Statistical inference
#' procedures for bivariate Archimedean copulas. Journal of the American
#' Statistical Association, 88 (423), 1034-1043.
#'
#' Schepsmeier, U. (2010). Maximum likelihood estimation of C-vine pair-copula
#' constructions based on bivariate copulas from different families. Diploma
#' thesis, Technische Universitaet Muenchen.\cr
#' <https://mediatum.ub.tum.de/?id=1079296>.
#'
#' @examples
#' \dontshow{set.seed(123)}
#' # simulate from Clayton copula
#' cop <- BiCop(3, tau = 0.5)
#' dat <- BiCopSim(1000, cop)
#'
#' # create lambda-function plots
#' op <- par(mfrow = c(1, 3))
#' BiCopLambda(dat[, 1], dat[, 2]) # empirical lambda-function
#' BiCopLambda(cop) # theoretical lambda-function
#' BiCopLambda(dat[, 1], dat[, 2], cop) # both
#' par(op)
#'
#'
BiCopLambda <- function(u1 = NULL, u2 = NULL, family = "emp", par = 0, par2 = 0,
PLOT = TRUE, obj = NULL, ...) {
## extract family and parameters if BiCop object is provided
if (!is.null(obj)) {
stopifnot(inherits(obj, "BiCop"))
family <- obj$family
par <- obj$par
par2 <- obj$par2
}
if (inherits(u1, "BiCop")) {
# for short hand usage extract from u1
if (inherits(u2, "logical"))
PLOT <- u2
obj <- u1
family <- obj$family
par <- obj$par
par2 <- obj$par2
u1 <- NULL
}
if (inherits(family, "BiCop")) {
# for short hand usage extract from family
if (inherits(par, "logical"))
PLOT <- par
obj <- family
family <- obj$family
par <- obj$par
par2 <- obj$par2
}
if ((is.null(u1) == TRUE) && (is.null(u2) == TRUE) && (family == 0 || par == 0))
stop("Either 'u1' and 'u2' have to be set for the emp.
lambda-function or 'family' and 'par' for the theo. lambda-function.")
if (length(u1) != length(u2))
stop("Lengths of 'u1' and 'u2' do not match.")
if (!(family %in% c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, "emp")))
stop("Copula family not implemented.")
if ((family %in% c(2, 7, 8, 9, 10)) && (par2 == 0))
stop("For t-, BB1 and BB7 copulas, 'par2' must be set.")
if (any(is.na(u1 + u2))) {
# send warning message
warning("Some of the data are NA. ",
"Only complete observations are used.")
# remove NAs
na.ind <- which(is.na(u1 + u2))
u1 <- u1[-na.ind]
u2 <- u2[-na.ind]
}
if ((is.null(u1) == FALSE) && (any(u1 > 1) || any(u1 < 0)))
stop("Data has be in the interval [0,1].")
if ((is.null(u2) == FALSE) && (any(u2 > 1) || any(u2 < 0)))
stop("Data has be in the interval [0,1].")
if (PLOT != TRUE && PLOT != FALSE)
stop("The parameter 'PLOT' has to be set to 'TRUE' or 'FALSE'.")
## check for parameter consistency
if (family != "emp")
BiCopCheck(family, par, par2)
## create grids
if (!is.null(u1))
v <- seq(0.001, 1, length.out = length(u1)) else v <- seq(0.001, 1, 0.001)
v1 <- v
theoLambda <- rep(0, length(v))
lambdaFull <- rep(0, length(v))
## calculate theoretical lambda
main <- ""
for (i in 2:length(v)) {
if (family == 3) {
theoLambda[i] <- -v1[i] * (1 - v1[i]^(par))/par
main <- "Clayton copula"
lambdaFull[i] <- (v1[i] * log(v1[i], exp(1)))
} else if (family == 4) {
theoLambda[i] <- v1[i] * log(v1[i])/(par)
main <- "Gumbel copula"
lambdaFull[i] <- (v1[i] * log(v1[i], exp(1)))
} else if (family == 5) {
theoLambda[i] <- -log((1 - exp(-par))/(1 - exp(-par * v1[i]))) * (1 - exp(-par * v1[i]))/(par * exp(-par * v1[i]))
main <- "Frank copula"
lambdaFull[i] <- (-v1[i] * log(1/v1[i], exp(1)))
} else if (family == 6) {
theoLambda[i] <- (log(1 - (1 - v[i])^par) * (1 - (1 - v[i])^par))/(par * (1 - v[i])^(par - 1))
main <- "Joe copula"
lambdaFull[i] <- (v1[i] * log(v1[i], exp(1)))
} else if (family == 7) {
theta <- par
delta <- par2
theoLambda[i] <- -1/(theta * delta) * (v[i]^(-theta) - 1)/(v[i]^(-1 - theta))
main <- "BB1 copula"
lambdaFull[i] <- (v1[i] * log(v1[i], exp(1)))
} else if (family == 8) {
theta <- par
delta <- par2
theoLambda[i] <- -log(-(1 - v[i])^theta + 1) * (1 - v[i] - (1 - v[i])^(-theta) + (1 - v[i])^(-theta) * v[i])/(delta * theta)
main <- "BB6 copula"
lambdaFull[i] <- (v1[i] * log(v1[i], exp(1)))
} else if (family == 9) {
theta <- par
delta <- par2
theoLambda[i] <- -1/(theta * delta) * ((1 - (1 - v[i])^theta)^(-delta) - 1)/((1 - v[i])^(theta - 1) * (1 - (1 - v[i])^theta)^(-delta - 1))
main <- "BB7 copula"
lambdaFull[i] <- (v1[i] * log(v1[i], exp(1)))
} else if (family == 10) {
theta <- par
delta <- par2
theoLambda[i] <- -log(((1 - v[i] * delta)^theta - 1)/((1 - delta)^theta - 1)) * (1 - v[i] * delta - (1 - v[i] * delta)^(-theta) + (1 - v[i] * delta)^(-theta) * v[i] * delta)/(theta * delta)
main <- "BB8 copula"
lambdaFull[i] <- (v1[i] * log(v1[i], exp(1)))
}
}
if (all(is.null(u1) == FALSE))
len <- length(u1) else len <- 1000
if (family == 1) {
theoLambda <- gtLambda(1, par, len = len)
main <- "Gaussian copula"
lambdaFull <- gtLambda(1, 0, len = len)
} else if (family == 2) {
theoLambda <- gtLambda(2, c(par, par2), len = len)
main <- "t-copula"
lambdaFull <- gtLambda(2, c(0, par2), len = len)
}
## calculate empircal lambda
if (family == "emp" || is.null(u1) == FALSE) {
n <- length(v)
nn <- length(u1)
empLambda <- rep(0, n)
K <- rep(0, n)
V <- rep(0, nn)
# Berechnung der konkordanten Paare
V[1:nn] <- mapply(function(x, y) length(which(x > u1 & y > u2)), u1, u2)
# Berechnung des empirischen K's
V1 <- V/(n - 1)
V1 <- V1
K <- sapply(v1, function(x) (1/nn) * length(which(V1[1:nn] <= x)))
# Berechnung der emp. lambdas
empLambda <- v1 - K
}
## create plot
if (PLOT) {
if (family == "emp") {
# only empirical one
if (min(empLambda[2:(nn - 1)]) < (-0.4)) {
plot(v1,
empLambda,
type = "l",
ylab = expression(lambda(v)),
xlab = "v",
ylim = c(-0.6, 0),
...)
} else {
plot(v1,
empLambda,
type = "l",
ylab = expression(lambda(v)),
xlab = "v",
ylim = c(-0.4, 0),
...)
}
} else if (family != "emp" && is.null(u1) == TRUE) {
# only theoretical one
if (min(theoLambda[2:999]) < (-0.4)) {
plot(v1,
theoLambda,
type = "l",
ylab = expression(lambda(v)),
xlab = "v",
ylim = c(-0.6, 0),
main = main,
...)
} else {
plot(v1,
theoLambda,
type = "l",
ylab = expression(lambda(v)),
xlab = "v",
ylim = c(-0.4, 0),
main = main,
...)
}
lines(v1, lambdaFull, type = "l", lty = 2)
lines(c(0, 1), c(0, 0), type = "l", lty = 2)
} else {
# both
minLambda <- min(min(theoLambda[2:(nn - 1)]), min(empLambda))
if (minLambda < (-0.4)) {
plot(v1,
empLambda,
type = "l",
ylab = expression(lambda(v)),
xlab = "v",
ylim = c(-0.6, 0),
...)
} else {
plot(v1,
empLambda,
type = "l",
ylab = expression(lambda(v)),
xlab = "v",
ylim = c(-0.4, 0),
...)
}
lines(v1, theoLambda, xlab = "v", type = "l", lwd = 2, col = "grey")
lines(v1, lambdaFull, type = "l", lty = 2)
lines(c(0, 1), c(0, 0), type = "l", lty = 2)
}
} else {
out <- list()
if (family == "emp")
out$empLambda <- empLambda else if (family != "emp" && is.null(u1) == TRUE) {
out$theoLambda <- theoLambda
} else {
out$empLambda <- empLambda
out$theoLambda <- theoLambda
}
return(out)
}
}
###### lambda-function for Gaussian- and t-copula #
# Input:
# copula Copula family (1='N',2='t')
# param Parameter
# Output:
# lambda lambda-function #
gtLambda <- function(copula, param, len = 10000) {
v <- seq(0.001, 1, length.out = len)
v1 <- v
n <- length(v)
nn <- len
lambda <- rep(0, n)
K <- rep(0, n)
V <- rep(0, nn)
mu <- c(0, 0)
if (param[1] == 0) {
for (i in 1:n) {
lambda[i] <- (v1[i] * log(v1[i], exp(1)))
}
} else {
if (copula == 1) {
rho <- param
# sigma=matrix(c(1,rho,rho,1), c(2,2)) Z=rmvnorm(nn, mu, sigma) u1=pnorm(Z[,1])
# u2=pnorm(Z[,2])
uu1 <- runif(nn)
vv2 <- runif(nn)
uu2 <- .C("Hinv1",
as.integer(1),
as.integer(nn),
as.double(vv2),
as.double(uu1),
as.double(rho),
as.double(0),
as.double(rep(0, nn)),
PACKAGE = "VineCopula")[[7]]
# x1=qnorm(U1) x2=qnorm(U2) C=dmvnorm(cbind(x1,x2), mu, sigma)
} else if (copula == 2) {
rho <- param[1]
nu <- param[2]
# sigma=matrix(c(1,rho,rho,1), c(2,2)) Z=rmvt(nn, sigma, nu) u1=pt(Z[,1], df=nu)
# u2=pt(Z[,2], df=nu)
uu1 <- runif(nn)
vv2 <- runif(nn)
uu2 <- .C("Hinv1",
as.integer(2),
as.integer(nn),
as.double(vv2),
as.double(uu1),
as.double(rho),
as.double(nu),
as.double(rep(0, nn)),
PACKAGE = "VineCopula")[[7]]
}
# Berechnung der konkordanten Paare
V[1:nn] <- mapply(function(x, y) length(which(x > uu1 & y > uu2)), uu1, uu2)
# Berechnung des empirischen K's
V1 <- V/(n - 1)
V1 <- V1
K <- sapply(v1, function(x) (1/nn) * length(which(V1[1:nn] <= x)))
# Berechnung der emp. lambdas
lambda <- v1 - K
# lambda=smooth(lambda)
}
return(lambda)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.