Package 'copula'

Description

The copula package provides (S4) classes of commonly used elliptical, (nested) Archimedean, extreme value and other copula families; methods for density, distribution, random number generation, and plots.

Fitting copula models and goodness-of-fit tests. Independence and serial (univariate and multivariate) independence tests, and other copula related tests.

Details

The DESCRIPTION file:

Index of help topics:

.pairsCond              Pairs Plot of a cu.u Object (Internal Use)
A..Z                    Sinc, Zolotarev's, and Other Mathematical
                        Utility Functions
An                      Nonparametric Rank-based Estimators of the
                        Pickands Dependence Function
Bernoulli               Compute Bernoulli Numbers
Copula                  Density, Evaluation, and Random Number
                        Generation for Copula Functions
Eulerian                Eulerian and Stirling Numbers of First and
                        Second Kind
K                       Kendall Distribution Function for Archimedean
                        Copulas
Mvdc                    Multivariate Distributions Constructed from
                        Copulas
RSpobs                  Pseudo-Observations of Radial and Uniform Part
                        of Elliptical and Archimedean Copulas
SMI.12                  SMI Data - 141 Days in Winter 2011/2012
Sibuya                  Sibuya Distribution - Sampling and
                        Probabilities
absdPsiMC               Absolute Value of Generator Derivatives via
                        Monte Carlo
acopula-class           Class "acopula" of Archimedean Copula Families
acopula-families        Specific Archimedean Copula Families ("acopula"
                        Objects)
allComp                 All Components of a (Inner or Outer) Nested
                        Archimedean Copula
archmCopula             Construction of Archimedean Copula Class Object
archmCopula-class       Class "archmCopula"
beta.                   Sample and Population Version of Blomqvist's
                        Beta for Archimedean Copulas
cCopula                 Conditional Distributions and Their Inverses
                        from Copulas
cloud2-methods          Cloud Plot Methods ('cloud2') in Package
                        'copula'
coeffG                  Coefficients of Polynomial used for Gumbel
                        Copula
contour-methods         Methods for Contour Plots in Package 'copula'
contourplot2-methods    Contour Plot Methods 'contourplot2' in Package
                        'copula'
copula-class            Mother Classes "Copula", etc of all Copulas in
                        the Package
copula-package          Multivariate Dependence Modeling with Copulas
corKendall              (Fast) Computation of Pairwise Kendall's Taus
dDiag                   Density of the Diagonal of (Nested) Archimedean
                        Copulas
describeCop             Copula (Short) Description as String
dnacopula               Density Evaluation for (Nested) Archimedean
                        Copulas
ebeta                   Various Estimators for (Nested) Archimedean
                        Copulas
ellipCopula             Construction of Elliptical Copula Class Objects
ellipCopula-class       Class "ellipCopula" of Elliptical Copulas
emde                    Minimum Distance Estimators for (Nested)
                        Archimedean Copulas
emle                    Maximum Likelihood Estimators for (Nested)
                        Archimedean Copulas
empCopula               The Empirical Copula
empCopula-class         Class "empCopula" of Empirical Copulas
enacopula               Estimation Procedures for (Nested) Archimedean
                        Copulas
evCopula                Construction of Extreme-Value Copula Objects
evCopula-class          Classes Representing Extreme-Value Copulas
evTestA                 Bivariate Test of Extreme-Value Dependence
                        Based on Pickands' Dependence Function
evTestC                 Large-sample Test of Multivariate Extreme-Value
                        Dependence
evTestK                 Bivariate Test of Extreme-Value Dependence
                        Based on Kendall's Distribution
exchEVTest              Test of Exchangeability for Certain Bivariate
                        Copulas
exchTest                Test of Exchangeability for a Bivariate Copula
fgmCopula               Construction of a fgmCopula Class Object
fgmCopula-class         Class "fgmCopula" - Multivariate Multiparameter
                        Farlie-Gumbel-Morgenstern Copulas
fhCopula                Construction of Fréchet-Hoeffding Bound Copula
                        Objects
fhCopula-class          Class "fhCopula" of Fréchet-Hoeffding Bound
                        Copulas
fitCopula               Fitting Copulas to Data - Copula Parameter
                        Estimation
fitCopula-class         Classes of Fitted Multivariate Models: Copula,
                        Mvdc
fitLambda               Non-parametric Estimators of the Matrix of
                        Tail-Dependence Coefficients
fitMvdc                 Estimation of Multivariate Models Defined via
                        Copulas
fixParam                Fix a Subset of a Copula Parameter Vector
gasoil                  Daily Crude Oil and Natural Gas Prices from
                        2003 to 2006
getAcop                 Get "acopula" Family Object by Name
getIniParam             Get Initial Parameter Estimate for Copula
getTheta                Get the Parameter(s) of a Copula
gnacopula               Goodness-of-fit Testing for (Nested)
                        Archimedean Copulas
gofBTstat               Various Goodness-of-fit Test Statistics
gofCopula               Goodness-of-fit Tests for Copulas
gofEVCopula             Goodness-of-fit Tests for Bivariate
                        Extreme-Value Copulas
gofTstat                Goodness-of-fit Test Statistics
htrafo                  GOF Testing Transformation of Hering and Hofert
iPsi                    Generator Functions for Archimedean and
                        Extreme-Value Copulas
indepCopula             Construction of Independence Copula Objects
indepCopula-class       Class "indepCopula"
indepTest               Test Independence of Continuous Random
                        Variables via Empirical Copula
initOpt                 Initial Interval or Value for Parameter
                        Estimation of Archimedean Copulas
interval                Construct Simple "interval" Object
interval-class          Class "interval" of Simple Intervals
khoudrajiCopula         Construction of copulas using Khoudraji's
                        device
khoudrajiCopula-class   Class '"khoudrajiCopula"' and its Subclasses
log1mexp                Compute f(a) = log(1 +/- exp(-a)) Numerically
                        Optimally
loss                    LOSS and ALAE Insurance Data
margCopula              Marginal copula of a Copula With Specified
                        Margins
mixCopula               Create Mixture of Copulas
mixCopula-class         Class '"mixCopula"' of Copula Mixtures
moCopula                The Marshall-Olkin Copula
moCopula-class          Class "moCopula" of Marshall-Olkin Copulas
multIndepTest           Independence Test Among Continuous Random
                        Vectors Based on the Empirical Copula Process
multSerialIndepTest     Serial Independence Test for Multivariate Time
                        Series via Empirical Copula
mvdc-class              Class "mvdc": Multivariate Distributions from
                        Copulas
nacFrail.time           Timing for Sampling Frailties of Nested
                        Archimedean Copulas
nacPairthetas           Pairwise Thetas of Nested Archimedean Copulas
nacopula-class          Class "nacopula" of Nested Archimedean Copulas
nesdepth                Nesting Depth of a Nested Archimedean Copula
                        ("nacopula")
onacopula               Constructing (Outer) Nested Archimedean Copulas
opower                  Outer Power Transformation of Archimedean
                        Copulas
p2P                     Tools to Work with Matrices
pacR                    Distribution of the Radial Part of an
                        Archimedean Copula
pairs2                  Scatter-Plot Matrix ('pairs') for Copula
                        Distributions with Nice Defaults
pairsRosenblatt         Plots for Graphical GOF Test via Pairwise
                        Rosenblatt Transforms
pairwiseCcop            Computations for Graphical GOF Test via
                        Pairwise Rosenblatt Transforms
persp-methods           Methods for Function 'persp' in Package
                        'copula'
plackettCopula          Construction of a Plackett Copula
plackettCopula-class    Class "plackettCopula" of Plackett Copulas
plot-methods            Methods for 'plot' in Package 'copula'
pnacopula               Evaluation of (Nested) Archimedean Copulas
pobs                    Pseudo-Observations
polylog                 Polylogarithm Li_s(z) and Debye Functions
polynEval               Evaluate Polynomials
printNacopula           Print Compact Overview of a Nested Archimedean
                        Copula ("nacopula")
prob                    Computing Probabilities of Hypercubes
qqplot2                 Q-Q Plot with Rugs and Pointwise Asymptotic
                        Confidence Intervals
rAntitheticVariates     Variance-Reduction Methods
rF01Frank               Sample Univariate Distributions Involved in
                        Nested Frank and Joe Copulas
rFFrank                 Sampling Distribution F for Frank and Joe
radSymTest              Test of Exchangeability for a Bivariate Copula
rdj                     Daily Returns of Three Stocks in the Dow Jones
retstable               Sampling Exponentially Tilted Stable
                        Distributions
rlog                    Sampling Logarithmic Distributions
rnacModel               Random nacopula Model
rnacopula               Sampling Nested Archimedean Copulas
rnchild                 Sampling Child 'nacopula's
rotCopula               Construction and Class of Rotated aka Reflected
                        Copulas
rstable1                Random numbers from (Skew) Stable Distributions
safeUroot               One-dimensional Root (Zero) Finding - Extra
                        "Safety" for Convenience
serialIndepTest         Serial Independence Test for Continuous Time
                        Series Via Empirical Copula
setTheta                Specify the Parameter(s) of a Copula
show-methods            Methods for 'show()' in Package 'copula'
splom2-methods          Methods for Scatter Plot Matrix 'splom2' in
                        Package 'copula'
tau                     Dependence Measures for Bivariate Copulas
tauAMH                  Ali-Mikhail-Haq ("AMH")'s and Joe's Kendall's
                        Tau
uranium                 Uranium Exploration Dataset of Cook & Johnson
                        (1986)
wireframe2-methods      Perspective Plots via 'wireframe2'
xvCopula                Model (copula) selection based on 'k'-fold
                        cross-validation

Further information is available in the following vignettes:

The copula package provides

• Classes (S4) of commonly used copulas including elliptical (normal and t; ellipCopula), Archimedean (Clayton, Gumbel, Frank, Joe, and Ali-Mikhail-Haq; ; archmCopula and acopula), extreme value (Gumbel, Husler-Reiss, Galambos, Tawn, and t-EV; evCopula), and other families (Plackett and Farlie-Gumbel-Morgenstern).

• Methods for density, distribution, random number generation (dCopula, pCopula and rCopula); bivariate dependence measures (rho, tau, etc), perspective and contour plots.

• Functions (and methods) for fitting copula models including variance estimates (fitCopula).

• Independence tests among random variables and vectors.

• Serial independence tests for univariate and multivariate continuous time series.

• Goodness-of-fit tests for copulas based on multipliers, and the parametric bootstrap, with several transformation options.

• Bivariate and multivariate tests of extreme-value dependence.

• Bivariate tests of exchangeability.

Now with former package nacopula for working with nested Archimedean copulas. Specifically,

• it provides procedures for computing function values and cube volumes (prob),

• characteristics such as Kendall's tau and tail dependence coefficients (via family objects, e.g., copGumbel),

• efficient sampling algorithms (rnacopula),

• various estimators and goodness-of-fit tests.

• The package also contains related univariate distributions and special functions such as the Sibuya distribution (Sibuya), the polylogarithm (polylog), Stirling and Eulerian numbers (Eulerian).

Further information is available in the following vignettes:

For a list of exported functions, use help(package = "copula").

References

Yan, J. (2007) Enjoy the Joy of Copulas: With a Package copula. Journal of Statistical Software 21(4), 1–21. https://www.jstatsoft.org/v21/i04/.

Kojadinovic, I. and Yan, J. (2010). Modeling Multivariate Distributions with Continuous Margins Using the copula R Package. Journal of Statistical Software 34(9), 1–20. https://www.jstatsoft.org/v34/i09/.

Hofert, M. and Mächler, M. (2011), Nested Archimedean Copulas Meet R: The nacopula Package., Journal of Statistical Software 39(9), 1–20. https://www.jstatsoft.org/v39/i09/.

Nelsen, R. B. (2006) An introduction to Copulas. Springer, New York.

See Also

The following CRAN packages currently use (‘depend on’) copula: CoClust, copulaedas, Depela, HAC, ipptoolbox, vines.

Examples

## Some of the more important functions (and their examples) are

example(fitCopula)## fitting Copulas
example(fitMvdc)  ## fitting multivariate distributions via Copulas
example(nacopula) ## nested Archimedean Copulas

## Independence Tests:  These also draw a 'Dependogram':
example(indepTest)       ## Testing for Independence
example(serialIndepTest) ## Testing for Serial Independence

Description

.pairsCond() is an internal function for plotting the pairwise Rosenblatt transforms, i.e., the pairwise conditional distributions, as returned by pairwiseCcop(), via the principal function pairsRosenblatt().

The intention is that pairsRosenblatt() be called, rather than this auxiliary function.

Usage

.pairsCond(gcu.u, panel = points, colList,
     col = par("col"), bg = par("bg"), labels, ...,
     text.panel = textPanel, label.pos = 0.5,
     cex.labels = NULL, font.labels = 1, gap = 0,
     axes = TRUE, panel.border = TRUE, key = TRUE,
     keyOpt = list(space= 2.5, width= 1.5, axis= TRUE,
                   rug.at= numeric(), title= NULL, line= 5),
     main = NULL, main.centered = FALSE,
     line.main = if(is.list(main)) 5/4*par("cex.main")* rev(seq_along(main)) else 2,
     sub = NULL, sub.centered = FALSE, line.sub = 4)

Arguments

Note

based on pairs.default() and filled.contour() from R-2.14.1 - used in Hofert and Maechler (2013)

Author(s)

Marius Hofert and Martin Maechler

See Also

pairsRosenblatt(), the prinicipal function, calling .pairsCond().

Description

Computes the absolute values of the $d$th generator derivative $\psi^{(d)}$ via Monte Carlo simulation.

Usage

absdPsiMC(t, family, theta, degree = 1, n.MC,
          method = c("log", "direct", "pois.direct", "pois"),
          log = FALSE, is.log.t = FALSE)

Arguments

Details

The absolute value of the $d$th derivative of the Laplace-Stieltjes transform $\psi=\mathcal{LS}[F]$ can be approximated via

$(-1)^d\psi^{(d)}(t)=\int_0^\infty x^d\exp(-tx)\,dF(x)\approx\frac{1}{N}\sum_{k=1}^NV_k^d\exp(-V_kt),\ t> 0,$

where $V_k\sim F,\ k\in\{1,\dots,N\}$. This approximation is used where $d=$degree and $N=$n.MC. Note that this is comparably fast even if t contains many evaluation points, since the random variates $V_k\sim F,\ k\in\{1,\dots,N\}$ only have to be generated once, not depending on t.

Value

numeric vector of the same length as t containing the absolute values of the generator derivatives.

References

Hofert, M., Mächler, M., and McNeil, A. J. (2013). Archimedean Copulas in High Dimensions: Estimators and Numerical Challenges Motivated by Financial Applications. Journal de la Société Française de Statistique 154(1), 25–63.

See Also

acopula-families.

Examples

t <- c(0:100,Inf)
set.seed(1)
(ps <- absdPsiMC(t, family="Gumbel", theta=2, degree=10, n.MC=10000, log=TRUE))
## Note: The absolute value of the derivative at 0 should be Inf for
## Gumbel, however, it is always finite for the Monte Carlo approximation
set.seed(1)
ps2 <- absdPsiMC(log(t), family="Gumbel", theta=2, degree=10,
                 n.MC=10000, log=TRUE, is.log.t = TRUE)
stopifnot(all.equal(ps[-1], ps2[-1], tolerance=1e-14))
## Now is there an advantage of using "is.log.t" ?
sapply(eval(formals(absdPsiMC)$method), function(MM)
       absdPsiMC(780, family="Gumbel", method = MM,
                 theta=2, degree=10, n.MC=10000, log=TRUE, is.log.t = TRUE))
## not really better, yet...

Description

This class "acopula" of Archimedean Copula Families is mainly used for providing objects of known Archimedean families with all related functions.

Objects from the Class

Objects can be created by calls of the form new("acopula", ...). For several well-known Archimedean copula families, the package copula already provides such family objects.

Slots

name:

A string (class "character") describing the copula family, for example, "AMH" (or simply "A"), "Clayton" ("C"), "Frank" ("F"), "Gumbel" ("G"), or "Joe" ("J").

theta:

Parameter value, a numeric, where NA means “unspecified”.

psi, iPsi:

The (Archimedean) generator $\psi$ (with $\psi$(t)=$\exp$(-t) being the generator of the independence copula) and its inverse (function). iPsi has an optional argument log which, if TRUE returns the logarithm of inverse of the generator.

absdPsi:

A function which computes the absolute value of the derivative of the generator $\psi$ for the given parameter theta and of the given degree (defaults to 1). Note that there is no informational loss by computing the absolute value since the derivatives alternate in sign (the generator derivative is simply (-1)^degree$*$absdPsi). The number n.MC denotes the sample size for a Monte Carlo evaluation approach. If n.MC is zero (the default), the generator derivatives are evaluated with their exact formulas. The optional parameter log (defaults to FALSE) indicates whether or not the logarithmic value is returned.

absdiPsi:

a function computing the absolute value of the derivative of the generator inverse (iPsi()) for the given parameter theta. The optional parameter log (defaults to FALSE) indicates whether the logarithm of the absolute value of the first derivative of iPsi() is returned.

dDiag:

a function computing the density of the diagonal of the Archimedean copula at u with parameter theta. The parameter log is as described before.

dacopula:

a function computing the density of the Archimedean copula at u with parameter theta. The meanings of the parameters n.MC and log are as described before.

score:

a function computing the derivative of the density with respect to the parameter $\theta$.

uscore:

a function computing the derivative of the density with respect to the each of the arguments.

paraInterval:

Either NULL or an object of class "interval", which is typically obtained from a call such as interval("[a,b)").

paraConstr:

A function of theta returning TRUE if and only if theta is a valid parameter value. Note that paraConstr is built automatically from the interval, whenever the paraInterval slot is valid. "interval".

nestConstr:

A function, which returns TRUE if and only if the two provided parameters theta0 and theta1 satisfy the sufficient nesting condition for this family.

V0:

A function which samples n random variates from the distribution $F$ with Laplace-Stieltjes transform $\psi$ and parameter theta.

dV0:

A function which computes either the probability mass function or the probability density function (depending on the Archimedean family) of the distribution function whose Laplace-Stieltjes transform equals the generator $\psi$ at the argument x (possibly a vector) for the given parameter theta. An optional argument log indicates whether the logarithm of the mass or density is computed (defaults to FALSE).

V01:

A function which gets a vector of realizations of V0, two parameters theta0 and theta1 which satisfy the sufficient nesting condition, and which returns a vector of the same length as V0 with random variates from the distribution function $F_{01}$ with Laplace-Stieltjes transform $\psi_{01}$ (see dV01) and parameters $\theta_0=\,$theta0, $\theta_1=\,$theta1.

dV01:

Similar to dV0 with the difference being that this function computes the probability mass or density function for the Laplace-Stieltjes transform

$\psi_{01}(t;V_0) = \exp(-V_0\psi_0^{-1}(\psi_1(t))),$

corresponding to the distribution function $F_{01}$.

Arguments are the evaluation point(s) x, the value(s) V0, and the parameters theta0 and theta1. As for dV0, the optional argument log can be specified (defaults to FALSE). Note that if x is a vector, V0 must either have length one (in which case V0 is the same for every component of x) or V0 must be of the same length as x (in which case the components of V0 correspond to the ones of x).

tau, iTau:

Compute Kendall's tau of the bivariate Archimedean copula with generator $\psi$ as a function of theta, respectively, theta as a function of Kendall's tau.

lambdaL, lambdaU, lambdaLInv, lambdaUInv:

Compute the lower (upper) tail-dependence coefficient of the bivariate Archimedean copula with generator $\psi$ as a function of theta, respectively, theta as a function of the lower (upper) tail-dependence coefficient.

For more details about Archimedean families, corresponding distributions and properties, see the references.

Methods

initialize

signature(.Object = "acopula"): is used to automatically construct the function slot paraConstr, when the paraInterval is provided (typically via interval()).

show

signature("acopula"): compact overview of the copula.

References

See those of the families, for example, copGumbel.

See Also

Specific provided copula family objects, for example, copAMH, copClayton, copFrank, copGumbel, copJoe.
To access these, you may also use getAcop.

A nested Archimedean copula without child copulas (see class "nacopula") is a proper Archimedean copula, and hence, onacopula() can be used to construct a specific parametrized Archimedean copula; see the example below.

Alternatively, setTheta also returns such a (parametrized) Archimedean copula.

Examples

## acopula class information
showClass("acopula")

## Information and structure of Clayton copulas
copClayton
str(copClayton)

## What are admissible parameters for Clayton copulas?
copClayton@paraInterval

## A Clayton "acopula" with Kendall's tau = 0.8 :
(cCl.2 <- setTheta(copClayton, iTau(copClayton, 0.8)))

## Can two Clayton copulas with parameters theta0 and theta1 be nested?
## Case 1: theta0 = 3, theta1 = 2
copClayton@nestConstr(theta0 = 3, theta1 = 2)
## -> FALSE as the sufficient nesting criterion is not fulfilled
## Case 2: theta0 = 2, theta1 = 3
copClayton@nestConstr(theta0 = 2, theta1 = 3) # TRUE

## For more examples, see  help("acopula-families")

Description

pacR() computes the distribution function $F_R$ of the radial part of an Archimedean copula, given by

$F_R(x)=1-\sum_{k=0}^{d-1} \frac{(-x)^k\psi^{(k)}(x)}{k!},\ x\in[0,\infty);$

The formula (in a slightly more general form) is given by McNeil and G. Nešlehová (2009).

qacR() computes the quantile function of $F_R$.

Usage

pacR(x, family, theta, d, lower.tail = TRUE, log.p = FALSE, ...)
qacR(p, family, theta, d, log.p = FALSE, interval,
     tol = .Machine$double.eps^0.25, maxiter = 1000, ...)

Arguments

Value

The distribution function of the radial part evaluated at x, or its inverse, the quantile at p.

References

McNeil, A. J., G. Nešlehová, J. (2009). Multivariate Archimedean copulas, $d$-monotone functions and $l_1$-norm symmetric distributions. The Annals of Statistics 37(5b), 3059–3097.

Examples

## setup
family <- "Gumbel"
tau <- 0.5
m <- 256
dmax <- 20
x <- seq(0, 20, length.out=m)

## compute and plot pacR() for various d's
y <- vapply(1:dmax, function(d)
            pacR(x, family=family, theta=iTau(archmCopula(family), tau), d=d),
            rep(NA_real_, m))
plot(x, y[,1], type="l", ylim=c(0,1),
     xlab = quote(italic(x)), ylab = quote(F[R](x)),
     main = substitute(italic(F[R](x))~~ "for" ~ d==1:.D, list(.D = dmax)))
for(k in 2:dmax) lines(x, y[,k])

Description

Given the nested Archimedean copula x, return an integer vector of the indices of all components of the corresponding outer_nacopula which are components of x, either direct components or components of possible child copulas. This is typically only used by programmers investigating the exact nesting structure.

For an outer_nacopula object x, allComp(x) must be the same as 1:dim(x), whereas its “inner” component copulas will each contain a subset of those indices only.

Usage

allComp(x)

Arguments

Value

An integer vector of indices $j$ of all components $u_j$ as described in the description above.

Examples

C3 <- onacopula("AMH", C(0.7135, 1, C(0.943, 2:3)))
 allComp(C3) # components are 1:3
 allComp(C3@childCops[[1]]) # for the child, only  (2, 3)

Description

Bivariate and multivariate versions of the nonparametric rank-based estimators of the Pickands dependence function $A$, studied in Genest and Segers (2009) and Gudendorf and Segers (2011).

Usage

An.biv(x, w, estimator = c("CFG", "Pickands"), corrected = TRUE,
       ties.method = eval(formals(rank)$ties.method))
An(x, w, ties.method = eval(formals(rank)$ties.method))

Arguments

Details

More details can be found in the references.

Value

An.biv() returns a vector containing the values of the estimated Pickands dependence function at the points in w (and is the same as former Anfun()).

The function An computes simultaneously the three corrected multivariate estimators studied in Gudendorf and Segers (2011) at the points in w and retuns a list whose components are

References

C. Genest and J. Segers (2009). Rank-based inference for bivariate extreme-value copulas. Annals of Statistics 37, 2990–3022.

G. Gudendorf and J. Segers (2011). Nonparametric estimation of multivariate extreme-value copulas. Journal of Statistical Planning and Inference 142, 3073–3085.

See Also

evCopula, A, and evTestA. Further, evTestC, evTestK, exchEVTest, and gofEVCopula.

Examples

## True Pickands dependence functions
curve(A(gumbelCopula(4   ), x), 0, 1)
curve(A(gumbelCopula(2   ), x), add=TRUE, col=2)
curve(A(gumbelCopula(1.33), x), add=TRUE, col=3)

## CFG estimator
curve(An.biv(rCopula(1000, gumbelCopula(4   )), x), lty=2, add=TRUE)
curve(An.biv(rCopula(1000, gumbelCopula(2   )), x), lty=2, add=TRUE, col=2)
curve(An.biv(rCopula(1000, gumbelCopula(1.33)), x), lty=2, add=TRUE, col=3)

## Pickands estimator
curve(An.biv(rCopula(1000, gumbelCopula(4   )), x, estimator="Pickands"),
      lty=3, add=TRUE)
curve(An.biv(rCopula(1000, gumbelCopula(2   )), x, estimator="Pickands"),
      lty=3, add=TRUE, col=2)
curve(An.biv(rCopula(1000, gumbelCopula(1.33)), x, estimator="Pickands"),
      lty=3, add=TRUE, col=3)

legend("bottomleft",  paste0("Gumbel(", format(c(4, 2, 1.33)),")"),
       lwd=1, col=1:3, bty="n")
legend("bottomright", c("true", "CFG est.", "Pickands est."), lty=1:3, bty="n")

## Relationship between An.biv and An
u <- c(runif(100),0,1) # include 0 and 1
x <- rCopula(1000, gumbelCopula(4))
r <- An(x, cbind(1-u, u))
all.equal(r$CFG, An.biv(x, u))
all.equal(r$P, An.biv(x, u, estimator="Pickands"))

## A trivariate example
x <- rCopula(1000, gumbelCopula(4, dim = 3))
u <- matrix(runif(300), 100, 3)
w <- u / apply(u, 1, sum)
r <- An(x, w)

## Endpoint corrections are applied
An(x, cbind(1, 0, 0))
An(x, cbind(0, 1, 0))
An(x, cbind(0, 0, 1))

Description

Constructs an Archimedean copula class object with its corresponding parameter and dimension.

Usage

archmCopula(family, param = NA_real_, dim = 2, ...)

claytonCopula(param = NA_real_, dim = 2,
          use.indepC = c("message", "TRUE", "FALSE"))
frankCopula(param = NA_real_, dim = 2,
          use.indepC = c("message", "TRUE", "FALSE"))
gumbelCopula(param = NA_real_, dim = 2,
          use.indepC = c("message", "TRUE", "FALSE"))
amhCopula(param = NA_real_, dim = 2,
          use.indepC = c("message", "TRUE", "FALSE"))
joeCopula(param = NA_real_, dim = 2,
          use.indepC = c("message", "TRUE", "FALSE"))

Arguments

Details

archmCopula() is a wrapper for claytonCopula(), frankCopula(), gumbelCopula(), amhCopula() and joeCopula.

For the mathematical definitions of the respective Archimedean families, see copClayton.

For $d = 2$, i.e. dim = 2, the AMH, Clayton and Frank copulas allow to model negative Kendall's tau (tau) behavior via negative $\theta$, for AMH and Clayton $-1 \le \theta$, and for Frank $-\infty < \theta$. The respective boundary cases are

AMH:

$\tau(\theta = -1) = -0.1817258$,

Clayton:

$\tau(\theta = -1) = -1$,

Frank:

$\tau(\theta = -Inf) = -1$ (as limit).

For the Ali-Mikhail-Haq copula family ("amhCopula"), only the bivariate case is available; however copAMH has no such limitation.

Note that in all cases except for Frank and AMH, and $d = 2$ and $theta < 0$, the densities (dCopula() methods) are evaluated via the dacopula slot of the corresponding acopula-classed Archimedean copulas, implemented in a numeric stable way without any restriction on the dimension $d$.
Similarly, the (cumulative) distribution function (“"the copula"” $C()$) is evaluated via the corresponding acopula-classed Archimedean copula's functions in the psi and iPsi slots.

Value

An Archimedean copula object of class "claytonCopula", "frankCopula", "gumbelCopula", "amhCopula", or "joeCopula".

References

R.B. Nelsen (2006), An introduction to Copulas, Springer, New York.

See Also

acopula-classed Archimedean copulas, such as copClayton, copGumbel, etc, notably for mathematical definitions including the meaning of param.

fitCopula for fitting such copulas to data.

ellipCopula, evCopula.

Examples

clayton.cop <- claytonCopula(2, dim = 3)
## scatterplot3d(rCopula(1000, clayton.cop))

## negative param (= theta) is allowed for dim = 2 :
tau(claytonCopula(-0.5)) ## = -1/3
tauClayton <- Vectorize(function(theta) tau(claytonCopula(theta, dim=2)))
plot(tauClayton, -1, 10, xlab=quote(theta), ylim = c(-1,1), n=1025)
abline(h=-1:1,v=0, col="#11111150", lty=2); axis(1, at=-1)

tauFrank <- Vectorize(function(theta) tau(frankCopula(theta, dim=2)))
plot(tauFrank, -40, 50, xlab=quote(theta), ylim = c(-1,1), n=1025)
abline(h=-1:1,v=0, col="#11111150", lty=2)

## tauAMH() is function in our package
iTau(amhCopula(), -1) # -1 with a range warning
iTau(amhCopula(), (5 - 8*log(2)) / 3) # -1 with a range warning

ic <- frankCopula(0) # independence copula (with a "message")
stopifnot(identical(ic,
   frankCopula(0, use.indepC = "TRUE")))# indep.copula  withOUT message
(fC <- frankCopula(0, use.indepC = "FALSE"))
## a Frank copula which corresponds to the indep.copula (but is not)

frankCopula(dim = 3)# with NA parameters
frank.cop <- frankCopula(3)# dim=2
persp(frank.cop, dCopula)

gumbel.cop <- archmCopula("gumbel", 5)
stopifnot(identical(gumbel.cop, gumbelCopula(5)))
contour(gumbel.cop, dCopula)

amh.cop <- amhCopula(0.5)
u. <- as.matrix(expand.grid(u=(0:10)/10, v=(0:10)/10, KEEP.OUT.ATTRS=FALSE))
du <- dCopula(u., amh.cop)
stopifnot(is.finite(du) | apply(u. == 0, 1,any)| apply(u. == 1, 1,any))

## A 7-dim Frank copula
frank.cop <- frankCopula(3, dim = 7)
x <- rCopula(5, frank.cop)
## dCopula now *does* work:
dCopula(x, frank.cop)

## A 7-dim Gumbel copula
gumbelC.7 <- gumbelCopula(2, dim = 7)
dCopula(x, gumbelC.7)

## A 12-dim Joe copula
joe.cop <- joeCopula(iTau(joeCopula(), 0.5), dim = 12)
x <- rCopula(5, joe.cop)
dCopula(x, joe.cop)

Description

Archimedean copula class.

Objects from the Class

Created by calls of the form new("archmCopula", ...) or rather typically by archmCopula(). Implemented families are Clayton, Gumbel, Frank, Joe, and Ali-Mikhail-Haq.

Slots

exprdist:

Object of class "expression": expressions of the cdf and pdf of the copula. These expressions are used in function pCopula and dCopula.

dimension, parameters, etc:

all inherited from the super class copula.

Methods

dCopula

signature(copula = "claytonCopula"): ...

pCopula

signature(copula = "claytonCopula"): ...

rCopula

signature(copula = "claytonCopula"): ...

dCopula

signature(copula = "frankCopula"): ...

pCopula

signature(copula = "frankCopula"): ...

rCopula

signature(copula = "frankCopula"): ...

dCopula

signature(copula = "gumbelCopula"): ...

pCopula

signature(copula = "gumbelCopula"): ...

rCopula

signature(copula = "gumbelCopula"): ...

dCopula

signature(copula = "amhCopula"): ...

pCopula

signature(copula = "amhCopula"): ...

rCopula

signature(copula = "amhCopula"): ...

dCopula

signature(copula = "joeCopula"): ...

pCopula

signature(copula = "joeCopula"): ...

rCopula

signature(copula = "joeCopula"): ...

Extends

Class "archmCopula" extends class "copula" directly. Class "claytonCopula", "frankCopula", "gumbelCopula", "amhCopula" and "joeCopula" extends class "archmCopula" directly.

Note

"gumbelCopula" is also of class "evCopula".

See Also

archmCopula, for constructing such copula objects; copula-class.

Description

These functions compute Kendall's tau, Spearman's rho, and the tail dependence index for bivariate copulas. iTau and iRho, sometimes called “calibration” functions are the inverses: they determine (“calibrate”) the copula parameter (which must be one-dimensional!) given the value of Kendall's tau or Spearman's rho.

Usage

tau (copula, ...)
rho (copula, ...)
lambda(copula, ...)
iTau (copula, tau, ...)
iRho (copula, rho, ...)

Arguments

Details

The calibration functions iTau() and iRho() in fact return a moment estimate of the parameter for one-parameter copulas.

When there are no closed-form expressions for Kendall's tau or Spearman's rho, the calibration functions use numerical approximation techniques (see the last reference). For closed-form expressions, see Frees and Valdez (1998). For the t copula, the calibration function based on Spearman's rho uses the corresponding expression for the normal copula as an approximation.

References

E.W. Frees and E.A. Valdez (1998) Understanding relationships using copulas. North American Actuarial Journal 2, 1–25.

Iwan Kojadinovic and Jun Yan (2010) Comparison of three semiparametric methods for estimating dependence parameters in copula models. Insurance: Mathematics and Economics 47, 52–63. doi:10.1016/j.insmatheco.2010.03.008

See Also

The acopula class objects have slots, tau, lambdaL, and lambdaU providing functions for tau(), and the two tail indices lambda(), and slot iTau for iTau(), see the examples and copGumbel, etc.

Examples

gumbel.cop <- gumbelCopula(3)
tau(gumbel.cop)
rho(gumbel.cop)
lambda(gumbel.cop)
iTau(joeCopula(), 0.5)
stopifnot(all.equal(tau(gumbel.cop), copGumbel@tau(3)),

          all.equal(lambda(gumbel.cop),
                    c(copGumbel@lambdaL(3), copGumbel@lambdaU(3)),
                    check.attributes=FALSE),

          all.equal(iTau (gumbel.cop, 0.681),
                    copGumbel@iTau(0.681))
)

## let us compute the sample versions
x <- rCopula(200, gumbel.cop)
cor(x, method = "kendall")
cor(x, method = "spearman")
## compare with the true parameter value 3
iTau(gumbel.cop, cor(x, method="kendall" )[1,2])
iRho(gumbel.cop, cor(x, method="spearman")[1,2])

Description

Compute the $n$th Bernoulli number, or generate all Bernoulli numbers up to the $n$th, using diverse methods, that is, algorithms.

NOTE the current default methods will be changed – to get better accuracy!

Usage

Bernoulli    (n, method = c("sumBin", "sumRamanujan", "asymptotic"),
              verbose = FALSE)
Bernoulli.all(n, method = c("A-T", "sumBin", "sumRamanujan", "asymptotic"),
              precBits = NULL, verbose = getOption("verbose"))

Arguments

Value

Bernoulli():

a number

Bernoulli.all():

a numeric vector of length n, containing B(n)

References

Kaneko, Masanobu (2000) The Akiyama-Tanigawa algorithm for Bernoulli numbers; Journal of Integer Sequences 3, article 00.2.9

See Also

Eulerian, Stirling1, etc.

Examples

## The example for the paper
MASS::fractions(Bernoulli.all(8, verbose=TRUE))

B10 <- Bernoulli.all(10)
MASS::fractions(B10)

system.time(B50  <- Bernoulli.all(50))#  {does not cache} -- still "no time"
system.time(B100 <- Bernoulli.all(100))# still less than a milli second

## Using Bernoulli() is not much slower, but hopefully *more* accurate!
## Check first - TODO
system.time(B.1c <- Bernoulli(100))# caches ..
system.time(B1c. <- Bernoulli(100))# ==> now much faster
stopifnot(identical(B.1c, B1c.))

if(FALSE)## reset the cache:
assign("Bern.tab", list(), envir = copula:::.nacopEnv)

## More experiments in the source of the copula package ../tests/Stirling-etc.R

Description

Compute the population (beta.()) and sample (betan()) version of Blomqvist's beta for an Archimedean copula.

See the reference below for definitions and formulas.

Usage

beta.(cop, theta, d, scaling=FALSE)
betan(u, scaling=FALSE)

Arguments

Value

beta.:

a number, being the population version of Blomqvist's beta for the corresponding Archimedean copula;

betan:

a number, being the sample version of Blomqvist's beta for the given data.

References

Schmid and Schmidt (2007), Nonparametric inference on multivariate versions of Blomqvist's beta and related measures of tail dependence, Metrika 66, 323–354.

See Also

acopula

Examples

beta.(copGumbel, 2.5, d = 5)

d.set <- c(2:6, 8, 10, 15, 20, 30)
cols <- adjustcolor(colorRampPalette(c("red", "orange", "blue"),
                                     space = "Lab")(length(d.set)), 0.8)
## AMH:
for(i in seq_along(d.set))
   curve(Vectorize(beta.,"theta")(copAMH, x, d = d.set[i]), 0, .999999,
         main = "Blomqvist's beta(.) for  AMH",
         xlab = quote(theta), ylab = quote(beta(theta, AMH)),
         add = (i > 1), lwd=2, col=cols[i])
mtext("NB:  d=2 and d=3 are the same")
legend("topleft", paste("d =",d.set), bty="n", lwd=2, col=cols)

## Gumbel:
for(i in seq_along(d.set))
   curve(Vectorize(beta.,"theta")(copGumbel, x, d = d.set[i]), 1, 10,
         main = "Blomqvist's beta(.) for  Gumbel",
         xlab = quote(theta), ylab = quote(beta(theta, Gumbel)),
         add=(i > 1), lwd=2, col=cols[i])
legend("bottomright", paste("d =",d.set), bty="n", lwd=2, col=cols)

## Clayton:
for(i in seq_along(d.set))
   curve(Vectorize(beta.,"theta")(copClayton, x, d = d.set[i]), 1e-5, 10,
         main = "Blomqvist's beta(.) for  Clayton",
         xlab = quote(theta), ylab = quote(beta(theta, Gumbel)),
         add=(i > 1), lwd=2, col=cols[i])
legend("bottomright", paste("d =",d.set), bty="n", lwd=2, col=cols)

## Joe:
for(i in seq_along(d.set))
   curve(Vectorize(beta.,"theta")(copJoe, x, d = d.set[i]), 1, 10,
         main = "Blomqvist's beta(.) for  Joe",
         xlab = quote(theta), ylab = quote(beta(theta, Gumbel)),
         add=(i > 1), lwd=2, col=cols[i])
legend("bottomright", paste("d =",d.set), bty="n", lwd=2, col=cols)

## Frank:
for(i in seq_along(d.set))
   curve(Vectorize(beta.,"theta")(copFrank, x, d = d.set[i]), 1e-5, 50,
         main = "Blomqvist's beta(.) for  Frank",
         xlab = quote(theta), ylab = quote(beta(theta, Gumbel)),
         add=(i > 1), lwd=2, col=cols[i])
legend("bottomright", paste("d =",d.set), bty="n", lwd=2, col=cols)

## Shows the numeric problems:
curve(Vectorize(beta.,"theta")(copFrank, x, d = 29), 35, 42, col="violet")

Description

Compute the conditional distribution function $C(u_d\,|\,u_1,\dots, u_{d-1})$ of $u_d$ given $u_1,\dots,u_{d-1}$.

Usage

cCopula(u, copula, indices = 1:dim(copula), inverse = FALSE,
        log = FALSE, drop = FALSE, ...)

## Deprecated (use cCopula() instead):
rtrafo(u, copula, indices = 1:dim(copula), inverse = FALSE, log = FALSE)
cacopula(u, cop, n.MC = 0, log = FALSE)

Arguments

Details

By default and if fed with a sample of the corresponding copula, cCopula() computes the Rosenblatt transform; see Rosenblatt (1952). The involved high-order derivatives for Archimedean copulas were derived in Hofert et al. (2012).

Sampling, that is, random number generation, can be achieved by using inverse=TRUE. In this case, the inverse Rosenblatt transformation is used, which, for sampling purposes, is also known as conditional distribution method. Note that, for Archimedean copulas not being Clayton, this can be slow as it involves numerical root finding in each (but the first) component.

Value

An $(n, k)$-matrix (unless n == 1 and drop is true, where a $k$-vector is returned) where $k$ is the length of indices. This matrix contains the conditional copula function values $C_{j|1,\dots,j-1}(u_j\,|\,u_1,\dots,u_{j-1})$ or, if inverse = TRUE, their inverses $C^-_{j|1,\dots,j-1}(u_j\,|\,u_1,\dots,u_{j-1})$ for all $j$ in indices.

Note

For some (but not all) families, this function also makes sense on the boundaries (if the corresponding limits can be computed).

References

Genest, C., Rémillard, B., and Beaudoin, D. (2009). Goodness-of-fit tests for copulas: A review and a power study. Insurance: Mathematics and Economics 44, 199–213.

Rosenblatt, M. (1952). Remarks on a Multivariate Transformation, The Annals of Mathematical Statistics 23, 3, 470–472.

Hofert, M., Mächler, M., and McNeil, A. J. (2012). Likelihood inference for Archimedean copulas in high dimensions under known margins. Journal of Multivariate Analysis 110, 133–150.

See Also

htrafo; acopula-families.

Examples

## 1) Sampling from a conditional distribution of a Clayton copula given u_1

## Define the copula
tau <- 0.5
theta <- iTau(claytonCopula(), tau = tau)
d <- 2
cc <- claytonCopula(theta, dim = d)
n <- 1000
set.seed(271)

## A small u_1
u1 <- 0.05
U <- cCopula(cbind(u1, runif(n)), copula = cc, inverse = TRUE)
plot(U[,2], ylab = quote(U[2]))

## A large u_1
u1 <- 0.95
U <- cCopula(cbind(u1, runif(n)), copula = cc, inverse = TRUE)
plot(U[,2], ylab = quote(U[2]))


## 2) Sample via conditional distribution method and then apply the
##    Rosenblatt transform
##    Note: We choose the numerically more involved (and thus slower)
##          Gumbel case here

## Define the copula
tau <- 0.5
theta <- iTau(gumbelCopula(), tau = tau)
d <- 5
gc <- gumbelCopula(theta, dim = d)
n <- 200
set.seed(271)
U. <- matrix(runif(n*d), ncol = d) # U(0,1)^d


## Transform to Gumbel sample via conditional distribution method
U <- cCopula(U., copula = gc, inverse = TRUE) # slow for ACs except Clayton
splom2(U) # scatter-plot matrix copula sample

## Rosenblatt transform back to U(0,1)^d (as a check)
U. <- cCopula(U, copula = gc)
splom2(U.) # U(0,1)^d again


## 3) cCopula() for elliptical copulas

tau <- 0.5
theta <- iTau(claytonCopula(), tau = tau)
d <- 5
cc <- claytonCopula(theta, dim = d)
set.seed(271)
n <- 1000
U <- rCopula(n, copula = cc)
X <- qnorm(U) # X now follows a meta-Clayton model with N(0,1) marginals
U <- pobs(X) # build pseudo-observations

fN <- fitCopula(normalCopula(dim = d), data = U) # fit a Gauss copula
U.RN <- cCopula(U, copula = fN@copula)
splom2(U.RN, cex = 0.2) # visible but not so clearly

f.t <- fitCopula(tCopula(dim = d), U)
U.Rt <- cCopula(U, copula = f.t@copula) # transform with a fitted t copula
splom2(U.Rt, cex = 0.2) # still visible but not so clear

## Inverse (and check consistency)
U.N <- cCopula(U.RN, copula = fN @copula, inverse = TRUE)
U.t <- cCopula(U.Rt, copula = f.t@copula, inverse = TRUE)

tol <- 1e-14
stopifnot(
    all.equal(U, U.N),
    all.equal(U, U.t),
    all.equal(log(U.RN),
              cCopula(U, copula = fN @copula, log = TRUE), tolerance = tol),
    all.equal(log(U.Rt),
              cCopula(U, copula = f.t@copula, log = TRUE), tolerance = tol)
)

## 4) cCopula() for a more sophisticated mixture copula (bivariate case only!)

tau <- 0.5
cc <- claytonCopula(iTau(claytonCopula(), tau = tau)) # first mixture component
tc <- tCopula(iTau(tCopula(), tau = tau), df = 3) # t_3 copula
tc90 <- rotCopula(tc, flip = c(TRUE, FALSE)) # t copula rotated by 90 degrees
wts <- c(1/2, 1/2) # mixture weights
mc <- mixCopula(list(cc, tc90), w = wts) # mixture copula with one copula rotated

set.seed(271)
U <- rCopula(n, copula = mc)
U. <- cCopula(U, copula = mc) # Rosenblatt transform back to U(0,1)^2 (as a check)
plot(U., xlab = quote(U*"'"[1]), ylab = quote(U*"'"[2])) # check for uniformity

Description

Function and Methods cloud2() to draw (lattice) cloud plots of two-dimensional distributions from package copula.

Usage

## S4 method for signature 'matrix'
cloud2(x,
      xlim = range(x[,1], finite = TRUE),
      ylim = range(x[,2], finite = TRUE),
      zlim = range(x[,3], finite = TRUE),
      xlab = NULL, ylab = NULL, zlab = NULL,
      scales = list(arrows = FALSE, col = "black"),
      par.settings = standard.theme(color = FALSE), ...)
## S4 method for signature 'data.frame'
cloud2(x,
      xlim = range(x[,1], finite = TRUE),
      ylim = range(x[,2], finite = TRUE),
      zlim = range(x[,3], finite = TRUE),
      xlab = NULL, ylab = NULL, zlab = NULL,
      scales = list(arrows = FALSE, col = "black"),
      par.settings = standard.theme(color = FALSE), ...)
## S4 method for signature 'Copula'
cloud2(x, n,
      xlim = 0:1, ylim = 0:1, zlim = 0:1,
      xlab = quote(U[1]), ylab = quote(U[2]), zlab = quote(U[3]), ...)
## S4 method for signature 'mvdc'
cloud2(x, n,
      xlim = NULL, ylim = NULL, zlim = NULL,
      xlab = quote(X[1]), ylab = quote(X[2]), zlab = quote(X[3]), ...)

Arguments

Value

An object of class “trellis” as returned by cloud().

Methods

Cloud plots for objects of class "matrix" , "data.frame", "Copula" or "mvdc".

See Also

The lattice-based splom2-methods for data, and wireframe2-methods and contourplot2-methods for functions.

Examples

## For 'matrix' objects
cop <- gumbelCopula(2, dim = 3)
n <- 1000
set.seed(271)
U <- rCopula(n, copula = cop)
cloud2(U, xlab = quote(U[1]), ylab = quote(U[2]), zlab = quote(U[3]))

## For 'Copula' objects
set.seed(271)
cloud2(cop, n = n) # same as above

## For 'mvdc' objects
mvNN <- mvdc(cop, c("norm", "norm", "exp"),
             list(list(mean = 0, sd = 1), list(mean = 1), list(rate = 2)))
cloud2(mvNN, n = n)

Description

Compute the coefficients $a_{d,k}(\theta)$ involved in the generator (psi) derivatives and the copula density of Gumbel copulas.

For non-small dimensions $d$, these are numerically challenging to compute accurately.

Usage

coeffG(d, alpha,
       method = c("sort", "horner", "direct", "dsumSibuya",
                  paste("dsSib", eval(formals(dsumSibuya)$method), sep = ".")),
       log = FALSE, verbose = FALSE)

Arguments

Value

a numeric vector of length d, of values

$a_k(\theta, d) = (-1)^{d-k}\sum_{j=k}^d \alpha^j * s(d,j) * S(j,k), k \in \{1,\ldots,d\}.$

Note

There are still known numerical problems (with non-"Rmpfr" methods; and those are slow), e.g., for d=100, alpha=0.8 and $sign(s(n,k)) = (-1)^{n-k}$.

As a consequence, the methods and its defaults may change in the future, and so the exact implementation of coeffG() is still considered somewhat experimental.

Examples

a.k  <- coeffG(16, 0.55)
plot(a.k, xlab = quote(k), ylab = quote(a[k]),
     main = "coeffG(16, 0.55)", log = "y", type = "o", col = 2)
a.kH <- coeffG(16, 0.55, method = "horner")
stopifnot(all.equal(a.k, a.kH, tol = 1e-11))# 1.10e-13 (64-bit Lnx, nb-mm4)

Description

Methods for function contour to draw contour lines aka a level plot for objects from package copula.

Usage

## S4 method for signature 'Copula'
contour(x, FUN,
                   n.grid = 26, delta = 0,
                   xlab = quote(u[1]), ylab = quote(u[2]),
                   box01 = TRUE, ...)
## S4 method for signature 'mvdc'
contour(x, FUN, xlim, ylim, n.grid = 26,
                   xlab = quote(x[1]), ylab = quote(x[2]),
		   box01 = FALSE, ...)

Arguments

Methods

Contour lines are drawn for "Copula" or "mvdc" objects, see x in the Arguments section.

See Also

The persp-methods for “perspective” aka “3D” plots.

Examples

contour(frankCopula(-0.8), dCopula)
contour(frankCopula(-0.8), dCopula, delta=1e-6)
contour(frankCopula(-1.2), pCopula)
contour(claytonCopula(2), pCopula)

## the Gumbel copula density is "extreme"
## --> use fine grid (and enough levels):
r <- contour(gumbelCopula(3), dCopula, n=200, nlevels=100)
range(r$z)# [0, 125.912]
## Now superimpose contours of three resolutions:
contour(r, levels = seq(1, max(r$z), by=2), lwd=1.5)
contour(r, levels = (1:13)/2, add=TRUE, col=adjustcolor(1,3/4), lty=2)
contour(r, levels = (1:13)/4, add=TRUE, col=adjustcolor(2,1/2),
        lty=3, lwd=3/4)

x <- mvdc(gumbelCopula(3), c("norm", "norm"),
          list(list(mean = 0, sd =1), list(mean = 1)))
contour(x, dMvdc, xlim=c(-2, 2), ylim=c(-1, 3))
contour(x, pMvdc, xlim=c(-2, 2), ylim=c(-1, 3))

Description

Methods for contourplot2(), a version of contourplot() from lattice, to draw contour plots of two dimensional distributions from package copula.

Usage

## NB: The 'matrix' and 'data.frame' methods are identical - documenting the former
## S4 method for signature 'matrix'
contourplot2(x, aspect = 1,
      xlim = range(x[,1], finite = TRUE),
      ylim = range(x[,2], finite = TRUE),
      xlab = NULL, ylab = NULL,
      cuts = 16, labels = !region, pretty = !labels,
      scales = list(alternating = c(1,1), tck = c(1,0)),
      region = TRUE, ...,
      col.regions = gray(seq(0.4, 1, length.out = max(100, 4*cuts))))

## S4 method for signature 'Copula'
contourplot2(x, FUN, n.grid = 26, delta = 0,
      xlim = 0:1, ylim = 0:1,
      xlab = quote(u[1]), ylab = quote(u[2]), ...)
## S4 method for signature 'mvdc'
contourplot2(x, FUN, n.grid = 26, xlim, ylim,
      xlab = quote(x[1]), ylab = quote(x[2]), ...)

Arguments

Value

An object of class “trellis” as returned by contourplot().

Methods

Contourplot plots for objects of class "matrix" , "data.frame", "Copula" or "mvdc".

See Also

The contour-methods for drawing perspective plots via base graphics.

The lattice-based wireframe2-methods to get perspective plots for functions, and cloud2-methods and splom2-methods for data.

Examples

## For 'matrix' objects
## The Frechet--Hoeffding bounds W and M
n.grid <- 26
u <- seq(0, 1, length.out = n.grid)
grid <- expand.grid("u[1]" = u, "u[2]" = u)
W <- function(u) pmax(0, rowSums(u)-1) # lower bound W
M <- function(u) apply(u, 1, min) # upper bound M
x.W <- cbind(grid, "W(u[1],u[2])" = W(grid)) # evaluate W on 'grid'
x.M <- cbind(grid, "M(u[1],u[2])" = M(grid)) # evaluate M on 'grid'
contourplot2(x.W) # contour plot of W
contourplot2(x.M) # contour plot of M

## For 'Copula' objects
cop <- frankCopula(-4)
contourplot2(cop, pCopula) # the copula
contourplot2(cop, pCopula, xlab = "x", ylab = "y") # adjusting the labels
contourplot2(cop, dCopula) # the density

## For 'mvdc' objects
mvNN <- mvdc(gumbelCopula(3), c("norm", "norm"),
             list(list(mean = 0, sd = 1), list(mean = 1)))
xl <- c(-2, 2)
yl <- c(-1, 3)
contourplot2(mvNN, FUN = dMvdc, xlim = xl, ylim = yl, contour = FALSE)
contourplot2(mvNN, FUN = dMvdc, xlim = xl, ylim = yl)
contourplot2(mvNN, FUN = dMvdc, xlim = xl, ylim = yl, region = FALSE, labels = FALSE)
contourplot2(mvNN, FUN = dMvdc, xlim = xl, ylim = yl, region = FALSE)
contourplot2(mvNN, FUN = dMvdc, xlim = xl, ylim = yl,
             col.regions = colorRampPalette(c("royalblue3", "maroon3"), space="Lab"))

Description

Specific Archimedean families ("acopula" objects) implemented in the package copula.

These families are “classical” as from p. 116 of Nelsen (2007). More specifially, see Table 1 of Hofert (2011).

Usage

copAMH
copClayton
copFrank
copGumbel
copJoe

Details

All these are objects of the formal class "acopula".

copAMH:

Archimedean family of Ali-Mikhail-Haq with parametric generator

$\psi(t)=(1-\theta)/(\exp(t)-\theta),\ t\in[0,\infty],$

with $\theta\in[0,1)$. The range of admissible Kendall's tau is $[0,1/3)$.

Note that the lower and upper tail-dependence coefficients are both zero, that is, this copula family does not allow for tail dependence.

copClayton:

Archimedean family of Clayton with parametric generator

$\psi(t)=(1+t)^{-1/\theta},\ t\in[0,\infty],$

with $\theta\in(0,\infty)$. The range of admissible Kendall's tau, as well as that of the lower tail-dependence coefficient, is (0,1). For dimension $d = 2$, $\theta\in(-1,\infty)$ is admissible where negative $\theta$ allow negative Kendall's taus. Note that this copula does not allow for upper tail dependence.

copFrank:

Archimedean family of Frank with parametric generator

$-\log(1-(1-e^{-\theta})\exp(-t))/\theta,\ t\in[0,\infty]$

with $\theta\in(0,\infty)$. The range of admissible Kendall's tau is (0,1). Note that this copula family does not allow for tail dependence.

copGumbel:

Archimedean family of Gumbel with parametric generator

$\exp(-t^{1/\theta}),\ t\in[0,\infty]$

with $\theta\in[1,\infty)$. The range of admissible Kendall's tau, as well as that of the upper tail-dependence coefficient, is [0,1). Note that this copula does not allow for lower tail dependence.

copJoe:

Archimedean family of Joe with parametric generator

$1-(1-\exp(-t))^{1/\theta},\ t\in[0,\infty]$

with $\theta\in[1,\infty)$. The range of admissible Kendall's tau, as well as that of the upper tail-dependence coefficient, is [0,1). Note that this copula does not allow for lower tail dependence.

Note that staying within one of these Archimedean families, all of them can be nested if two (generic) generator parameters $\theta_0$, $\theta_1$ satisfy $\theta_0\le\theta_1$.

Value

A "acopula" object.

References

Nelsen, R. B. (2007). An Introduction to Copulas (2nd ed.). Springer.

Hofert, M. (2010). Sampling Nested Archimedean Copulas with Applications to CDO Pricing. Suedwestdeutscher Verlag fuer Hochschulschriften AG & Co. KG.

Hofert, M. (2011). Efficiently sampling nested Archimedean copulas. Computational Statistics & Data Analysis 55, 57–70.

Hofert, M. and Mächler, M. (2011). Nested Archimedean Copulas Meet R: The nacopula Package. Journal of Statistical Software 39(9), 1–20. https://www.jstatsoft.org/v39/i09/.

See Also

The class definition, "acopula". onacopula and setTheta for such Archimedean copulas with specific parameters.
getAcop accesses these families “programmatically”.

Examples

## Print a copAMH object and its structure
copAMH
str(copAMH)

## Show admissible parameters for a Clayton copula
copClayton@paraInterval

## Generate random variates from a Log(p) distribution via V0 of Frank
p <- 1/2
copFrank@V0(100, -log(1-p))

## Plot the upper tail-dependence coefficient as a function in the
## parameter for Gumbel's family
curve(copGumbel@lambdaU(x), xlim = c(1, 10), ylim = c(0,1), col = 4)

## Plot Kendall's tau as a function in the parameter for Joe's family
curve(copJoe@tau(x), xlim = c(1, 10), ylim = c(0,1), col = 4)

## ------- Plot psi() and tau() - and properties of all families ----

## The copula families currently provided:
(famNms <- ls("package:copula", patt="^cop[A-Z]"))

op <- par(mfrow= c(length(famNms), 2),
          mar = .6+ c(2,1.4,1,1), mgp = c(1.1, 0.4, 0))
for(nm in famNms) { Cf <- get(nm)
   thet <- Cf@iTau(0.3)
   curve(Cf@psi(x, theta = thet), 0, 5,
         xlab = quote(x), ylab="", ylim=0:1, col = 2,
         main = substitute(list(NAM ~~~ psi(x, theta == TH), tau == 0.3),
                           list(NAM=Cf@name, TH=thet)))
   I <- Cf@paraInterval
   Iu <- pmin(10, I[2])
   curve(Cf@tau(x), I[1], Iu, col = 3,
         xlab = bquote(theta %in% .(format(I))), ylab = "",
         main = substitute(NAM ~~ tau(theta), list(NAM=Cf@name)))
}
par(op)

## Construct a bivariate Clayton copula with parameter theta
theta <- 2
C2 <- onacopula("Clayton", C(theta, 1:2))
C2@copula # is an "acopula" with specific parameter theta

curve(C2@copula@psi(x, C2@copula@theta),
      main = quote("Generator" ~~ psi ~~ " of Clayton A.copula"),
      xlab = quote(theta1), ylab = quote(psi(theta1)),
      xlim = c(0,5), ylim = c(0,1), col = 4)

## What is the corresponding Kendall's tau?
C2@copula@tau(theta) # 0.5

## What are the corresponding tail-dependence coefficients?
C2@copula@lambdaL(theta)
C2@copula@lambdaU(theta)

## Generate n pairs of random variates from this copula
U <- rnacopula(n = 1000, C2)
## and plot the generated pairs of random variates
plot(U, asp=1, main = "n = 1000 from  Clayton(theta = 2)")

Description

Density (dCopula), distribution function (pCopula), and random generation (rCopula) for a copula object.

Usage

dCopula(u, copula, log=FALSE, ...)
pCopula(u, copula, ...)
rCopula(n, copula, ...)

Arguments

Details

The density (dCopula) and distribution function (pCopula) methods for Archimedean copulas now use the corresponding function slots of the Archimedean copula objects, such as copClayton, copGumbel, etc.

If an $u_j$ is outside $(0,1)$ we declare the density to be zero, and this is true even when another $u_k, k \ne j$ is NA or NaN; see also the “outside” example.

The distribution function of a t copula uses pmvt from package mvtnorm; similarly, the density (dCopula) calls dmvt from CRAN package mvtnorm. The normalCopula methods use dmvnorm and pmvnorm from the same package.

The random number generator for an Archimedean copula uses the conditional approach for the bivariate case and the Marshall-Olkin (1988) approach for dimension greater than 2.

Value

dCopula() gives the density, pCopula() gives the distribution function, and rCopula() generates random variates.

References

Frees, E. W. and Valdez, E. A. (1998). Understanding relationships using copulas. North American Actuarial Journal 2, 1–25.

Genest, C. and Favre, A.-C. (2007). Everything you always wanted to know about copula modeling but were afraid to ask. Journal of Hydrologic Engineering 12, 347–368.

Joe, H. (1997). Multivariate Models and Dependence Concepts. Chapman and Hall, London.

Marshall, A. W. and Olkin, I. (1988) Families of multivariate distributions. Journal of the American Statistical Association 83, 834–841.

Nelsen, R. B. (2006) An introduction to Copulas. Springer, New York.

See Also

the copula and acopula classes, the acopula families, acopula-families. Constructor functions such as ellipCopula, archmCopula, fgmCopula.

Examples

norm.cop <- normalCopula(0.5)
norm.cop
## one d-vector =^= 1-row matrix, works too :
dCopula(c(0.5, 0.5), norm.cop)
pCopula(c(0.5, 0.5), norm.cop)

u <- rCopula(100, norm.cop)
plot(u)
dCopula(u, norm.cop)
pCopula(u, norm.cop)
persp  (norm.cop, dCopula)
contour(norm.cop, pCopula)

## a 3-dimensional normal copula
u <- rCopula(1000, normalCopula(0.5, dim = 3))
if(require(scatterplot3d))
  scatterplot3d(u)

## a 3-dimensional clayton copula
cl3 <- claytonCopula(2, dim = 3)
v <- rCopula(1000, cl3)
pairs(v)
if(require(scatterplot3d))
  scatterplot3d(v)

## Compare with the "nacopula" version :
fu1 <- dCopula(v, cl3)
fu2 <- copClayton@dacopula(v, theta = 2)
Fu1 <- pCopula(v, cl3)
Fu2 <- pCopula(v, onacopula("Clayton", C(2.0, 1:3)))
## The density and cumulative values are the same:
stopifnot(all.equal(fu1, fu2, tolerance= 1e-14),
          all.equal(Fu1, Fu2, tolerance= 1e-15))

## NA and "outside" u[]
u <- v[1:12,]
## replace some by values outside (0,1) and some by NA/NaN
u[1, 2:3] <- c(1.5, NaN); u[2, 1] <- 2; u[3, 1:2] <- c(NA, -1)
u[cbind(4:9, 1:3)] <- c(NA, NaN)
f <- dCopula(u, cl3)
cbind(u, f) # note: f(.) == 0 at [1] and [3] inspite of NaN/NA
stopifnot(f[1:3] == 0, is.na(f[4:9]), 0 < f[10:12])