Version:	1.1-6
VersionNote:	Last CRAN: 1.1-5 on 2025-02-09, 1.1-4 on 2024-08-12, 1.1-3 on 2023-12-07
Date:	2025-03-05
Title:	Multivariate Dependence with Copulas
Depends:	R (≥ 3.5.0)
Imports:	stats, graphics, methods, stats4, Matrix (≥ 1.5-0), lattice, colorspace, gsl, ADGofTest, stabledist (≥ 0.6-4), mvtnorm, pcaPP, pspline, numDeriv
Suggests:	MASS, KernSmooth, sfsmisc, scatterplot3d, Rmpfr, bbmle, knitr, rmarkdown, animation, abind, crop, gridExtra, HAC, lcopula, mev, mvnormtest, parallel, partitions, polynom, qrng, randtoolbox, rugarch, Runuran, tseries, VGAM, VineCopula, zoo
SuggestsNote:	packages {animation, ..., zoo} (last lines above) are only used in vignettes, demos and a few tests.
VignetteBuilder:	knitr
Enhances:	nor1mix, copulaData
SystemRequirements:	pdfcrop (part of TexLive) is required to rebuild the vignettes.
Description:	Classes (S4) of commonly used elliptical, Archimedean, extreme-value and other copula families, as well as their rotations, mixtures and asymmetrizations. Nested Archimedean copulas, related tools and special functions. Methods for density, distribution, random number generation, bivariate dependence measures, Rosenblatt transform, Kendall distribution function, perspective and contour plots. Fitting of copula models with potentially partly fixed parameters, including standard errors. Serial independence tests, copula specification tests (independence, exchangeability, radial symmetry, extreme-value dependence, goodness-of-fit) and model selection based on cross-validation. Empirical copula, smoothed versions, and non-parametric estimators of the Pickands dependence function.
License:	GPL (≥ 3) | file LICENCE
Collate:	AllClass.R Classes.R AllGeneric.R Auxiliaries.R aux-acopula.R asymCopula.R mixCopula.R rotCopula.R Copula.R special-func.R amhCopula.R claytonCopula.R frankCopula.R cop_objects.R nacopula.R dC-dc.R amhExpr.R An.R archmCopula.R cCopula.R claytonExpr.R ellipCopula.R empCopula.R empPsi.R acR.R estimation.R evCopula.R evTests.R exchTests.R fgmCopula.R fitCopula.R fitLambda.R fitMvdc.R fixedPar.R frankExpr.R galambosCopula.R galambosExpr-math.R galambosExpr.R ggraph-tools.R pairsRosenblatt.R prob.R gofTrafos.R gofEVTests.R gofCopula.R graphics.R gumbelCopula.R gumbelExpr.R huslerReissCopula.R huslerReissExpr.R indepCopula.R fhCopula.R lowfhCopula.R upfhCopula.R indepTests.R joeCopula.R K.R logseries.R mvdc.R margCopula.R matrix_tools.R normalCopula.R obs.R opower.R plackettCopula.R plackettExpr.R moCopula.R rstable1.R safeUroot.R schlatherCopula.R stable.R timing.R tCopula.R tawnCopula.R tawnExpr.R tevCopula.R varianceReduction.R wrapper.R xvCopula.R zzz.R
Encoding:	UTF-8
URL:	https://copula.r-forge.r-project.org/, https://r-forge.r-project.org/projects/copula/, https://CRAN.r-project.org/package=copula
BugReports:	https://r-forge.r-project.org/tracker/?func=add&group_id=2140&atid=5417
NeedsCompilation:	yes
Packaged:	2025-03-18 09:43:59 UTC; maechler
Author:	Marius Hofert [aut], Ivan Kojadinovic [aut], Martin Maechler [aut, cre], Jun Yan [aut], Johanna G. Nešlehová [ctb] (evTestK()), Rebecca Morger [ctb] (fitCopula.ml(): code for free mixCopula weight parameters)
Maintainer:	Martin Maechler <maechler@stat.math.ethz.ch>
Repository:	CRAN
Date/Publication:	2025-03-18 13:30:06 UTC

Multivariate Dependence Modeling with Copulas

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

Pairs Plot of a cu.u Object (Internal Use)

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().

Nonparametric Rank-based Estimators of the Pickands Dependence Function

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))

Compute Bernoulli Numbers

Description

Compute the nth Bernoulli number, or generate all Bernoulli numbers up to the nth, 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

Density, Evaluation, and Random Number Generation for Copula Functions

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])

Kendall Distribution Function for Archimedean Copulas

Description

The Kendall distribution of an Archimedean copula is defined by

K(u) = P(C(U_1,U_2,\dots,U_d) \le u),

where u \in [0,1], and the d-dimensional (U_1,U_2,\dots,U_d) is distributed according to the copula C. Note that the random variable C(U_1,U_2,\dots,U_d) is known as “probability integral transform”. Its distribution function K is equal to the identity if d = 1, but is non-trivial for d \ge 2.

Kn() computes the empirical Kendall distribution function, pK() the distribution function (so K() itself), qK() the quantile function, dK() the density, and rK() random number generation from K() for an Archimedean copula.

Usage

Kn(u, x, method = c("GR", "GNZ")) # empirical Kendall distribution function
dK(u, copula, d, n.MC = 0, log.p = FALSE) # density
pK(u, copula, d, n.MC = 0, log.p = FALSE) # df
qK(p, copula, d, n.MC = 0, log.p = FALSE, # quantile function
   method = c("default", "simple", "sort", "discrete", "monoH.FC"),
   u.grid, xtraChecks = FALSE, ...)
rK(n, copula, d) # random number generation

Arguments

Details

For a completely monotone Archimedean generator \psi,

K(u)=\sum_{k=0}^{d-1} \frac{\psi^{(k)}(\psi^{-1}(u))}{k!} (-\psi^{-1}(u))^k,\ u\in[0,1];

see Barbe et al. (1996). The corresponding density is

\frac{(-1)^d\psi^{(d)}(\psi^{-1}(u))}{(d-1)!} (-(\psi^{-1})'(u))(\psi^{-1}(u))^{d-1}

Value

The empirical Kendall distribution function, density, distribution function, quantile function and random number generator.

Note

Currently, the "default" method of qK() is fast but not very accurate, see the ‘Examples’ for more accuracy (with more CPU effort).

References

Barbe, P., Genest, C., Ghoudi, K., and Rémillard, B. (1996), On Kendall's Process, Journal of Multivariate Analysis 58, 197–229.

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. doi:10.1016/j.jmva.2012.02.019

Genest, C. and Rivest, L.-P. (1993). Statistical inference procedures for bivariate Archimedean copulas. Journal of the American Statistical Association 88, 1034–1043.

Genest, C., G. Nešlehová, J., and Ziegel, J. (2011). Inference in multivariate Archimedean copula models. TEST 20, 223–256.

See Also

htrafo or emde (where K is used); splinefun(*, "monoHC") for that method.

Examples

tau <- 0.5
(theta <- copGumbel@iTau(tau)) # 2
d <- 20
(cop <- onacopulaL("Gumbel", list(theta,1:d)))

## Basic check of the empirical Kendall distribution function
set.seed(271)
n <- 1000
U <- rCopula(n, copula = cop)
X <- qnorm(U)
K.sample <- pCopula(U, copula = cop)
u <- seq(0, 1, length.out = 256)
edfK <- ecdf(K.sample)
plot(u, edfK(u), type = "l", ylim = 0:1,
     xlab = quote(italic(u)), ylab = quote(K[n](italic(u)))) # simulated
K.n <- Kn(u, x = X)
lines(u, K.n, col = "royalblue3") # Kn
## Difference at 0
edfK(0) # edf of K at 0
K.n[1] # K_n(0); this is > 0 since K.n is the edf of a discrete distribution
## => therefore, Kn(K.sample, x = X) is not uniform
plot(Kn(K.sample, x = X), ylim = 0:1)
## Note: Kn(0) -> 0 for n -> Inf

## Compute Kendall distribution function
u <- seq(0,1, length.out = 255)
Ku    <- pK(u, copula = cop@copula, d = d) # exact
Ku.MC <- pK(u, copula = cop@copula, d = d, n.MC = 1000) # via Monte Carlo
stopifnot(all.equal(log(Ku),
		    pK(u, copula = cop@copula, d = d, log.p=TRUE)))# rel.err 3.2e-16

## Build sample from K
set.seed(1)
n <- 200
W <- rK(n, copula = cop)

## Plot empirical distribution function based on W
## and the corresponding theoretical Kendall distribution function
## (exact and via Monte Carlo)
plot(ecdf(W), col = "blue", xlim = 0:1, verticals=TRUE,
     main = quote("Empirical"~ F[n](C(U)) ~
                     "and its Kendall distribution" ~ K(u)),
     do.points = FALSE, asp = 1)
abline(0,1, lty = 2); abline(h = 0:1, v = 0:1, lty = 3, col = "gray")
lines(u, Ku.MC, col = "red") # not quite monotone
lines(u, Ku, col = "black")  # strictly  monotone:
stopifnot(diff(Ku) >= 0)
legend(.25, .75, expression(F[n], K[MC](u), K(u)),
       col=c("blue" , "red", "black"), lty = 1, lwd = 1.5, bty = "n")

if(require("Rmpfr")) { # pK() now also works with high precision numbers:
 uM <- mpfr(0:255, 99)/256
 if(FALSE) {
   # not yet, now fails in  polyG() :
   KuM <- pK(uM, copula = cop@copula, d = d)
  ##  debug(copula:::.pK)
  debug(copula:::polyG)
 }
}# if( Rmpfr )


## Testing qK
pexpr <- quote( 0:63/63 );  p <- eval(pexpr)
d <- 10
cop <- onacopulaL("Gumbel", list(theta = 2, 1:d))
system.time(qK0 <- qK(p, copula = cop@copula, d = d)) # "default" - fast


system.time(qK1  <- qK(p, copula= cop@copula, d=d, method = "simple"))
system.time(qK1. <- qK(p, copula= cop@copula, d=d, method = "simple", tol = 1e-12))
system.time(qK2  <- qK(p, copula= cop@copula, d=d, method = "sort"))
system.time(qK2. <- qK(p, copula= cop@copula, d=d, method = "sort",   tol = 1e-12))
system.time(qK3  <- qK(p, copula= cop@copula, d=d, method = "discrete", u.grid = 0:1e4/1e4))
system.time(qK4  <- qK(p, copula= cop@copula, d=d, method = "monoH.FC",
                       u.grid = 0:5e2/5e2))
system.time(qK4. <- qK(p, copula= cop@copula, d=d, method = "monoH.FC",
                       u.grid = 0:5e2/5e2, tol = 1e-12))
system.time(qK5  <- qK(p, copula= cop@copula, d=d, method = "monoH.FC",
                       u.grid = 0:5e3/5e3))
system.time(qK5. <- qK(p, copula= cop@copula, d=d, method = "monoH.FC",
                       u.grid = 0:5e3/5e3, tol = 1e-12))
system.time(qK6  <- qK(p, copula= cop@copula, d=d, method = "monoH.FC",
                       u.grid = (0:5e3/5e3)^2))
system.time(qK6. <- qK(p, copula= cop@copula, d=d, method = "monoH.FC",
                       u.grid = (0:5e3/5e3)^2, tol = 1e-12))

## Visually they all coincide :
cols <- adjustcolor(c("gray50", "gray80", "light blue",
                      "royal blue", "purple3", "purple4", "purple"), 0.6)
matplot(p, cbind(qK0, qK1, qK2, qK3, qK4, qK5, qK6), type = "l", lwd = 2*7:1, lty = 1:7, col = cols,
        xlab = bquote(p == .(pexpr)), ylab = quote({K^{-1}}(u)),
        main = "qK(p, method = *)")
legend("topleft", col = cols, lwd = 2*7:1, lty = 1:7, bty = "n", inset = .03,
       legend= paste0("method= ",
             sQuote(c("default", "simple", "sort",
                      "discrete(1e4)", "monoH.FC(500)", "monoH.FC(5e3)", "monoH.FC(*^2)"))))

## See they *are* inverses  (but only approximately!):
eqInv <- function(qK) all.equal(p, pK(qK, cop@copula, d=d), tol=0)

eqInv(qK0 ) # "default"	       0.03  worst
eqInv(qK1 ) # "simple"	       0.0011 - best
eqInv(qK1.) # "simple", e-12   0.00000 (8.73 e-13) !
eqInv(qK2 ) # "sort"	       0.0013 (close)
eqInv(qK2.) # "sort", e-12     0.00000 (7.32 e-12)
eqInv(qK3 ) # "discrete"       0.0026
eqInv(qK4 ) # "monoH.FC(500)"  0.0095
eqInv(qK4.) # "m.H.FC(5c)e-12" 0.00963
eqInv(qK5 ) # "monoH.FC(5e3)"  0.001148
eqInv(qK5.) # "m.H.FC(5k)e-12" 0.000989
eqInv(qK6 ) # "monoH.FC(*^2)"  0.001111
eqInv(qK6.) # "m.H.FC(*^2)e-12"0.00000 (1.190 e-09)

## and ensure the differences are not too large
stopifnot(
 all.equal(qK0, qK1, tol = 1e-2) # !
 ,
 all.equal(qK1, qK2, tol = 1e-4)
 ,
 all.equal(qK2, qK3, tol = 1e-3)
 ,
 all.equal(qK3, qK4, tol = 1e-3)
 ,
 all.equal(qK4, qK0, tol = 1e-2) # !
)


stopifnot(all.equal(p, pK(qK0, cop@copula, d=d), tol = 0.04))

Multivariate Distributions Constructed from Copulas

Description

Density, distribution function, and random generator for a multivariate distribution via copula and parametric margins.

For likelihood and fitting these distributions to data, see fitMvdc.

Usage

mvdc(copula, margins, paramMargins, marginsIdentical = FALSE,
     check = TRUE, fixupNames = TRUE)
dMvdc(x, mvdc, log=FALSE)
pMvdc(x, mvdc)
rMvdc(n, mvdc)

Arguments

Details

The characters in argument margins are used to construct density, distribution, and quantile function names. For example, norm can be used to specify marginal distribution, because dnorm, pnorm, and qnorm are all available.

A user-defined distribution, for example, fancy, can be used as margin provided that dfancy, pfancy, and qfancy are available.

Each component list in argument paramMargins is a list with named components which are used to specify the parameters of the marginal distributions. For example, the list

    paramMargins = list(list(mean = 0, sd = 2), list(rate = 2))

can be used to specify that the first margin is normal with mean 0 and standard deviation 2, and the second margin is exponential with rate 2.

Value

mvdc() constructs an object of class "mvdc". dMvdc() gives the density, pMvdc() gives the cumulative distribution function, and rMvdc() generates random variates.

Note

mvdc(), fitMvdc, etc, are only for parametric margins. If you do not want to model all margins parametrically, use the standard copula approach, transforming the data by their empirical margins via pobs and modelling the copula alone, e.g., using fitCopula, i.e., conceptually, using

     fitCopula(.., pobs(x))

See Also

ellipCopula, archmCopula; the classes mvdc and copula.

Examples

## construct a bivariate distribution whose marginals
## are normal and exponential respectively, coupled
## together via a normal copula
mv.NE <- mvdc(normalCopula(0.75), c("norm", "exp"),
              list(list(mean = 0, sd =2), list(rate = 2)))
dim(mv.NE)
mv.NE  # using its print() / show() method

persp  (mv.NE, dMvdc, xlim = c(-4, 4), ylim=c(0, 2), main = "dMvdc(mv.NE)")
persp  (mv.NE, pMvdc, xlim = c(-4, 4), ylim=c(0, 2), main = "pMvdc(mv.NE)")
contour(mv.NE, dMvdc, xlim = c(-4, 4), ylim=c(0, 2))

# Generate (bivariate) random numbers from that, and visualize
x.samp <- rMvdc(250, mv.NE)
plot(x.samp)
summary(fx <- dMvdc(x.samp, mv.NE))
summary(Fx <- pMvdc(x.samp, mv.NE))
op <- par(mfcol=c(1,2))
pp <- persp(mv.NE, pMvdc, xlim = c(-5,5), ylim=c(0,2),
            main = "pMvdc(mv.NE)", ticktype="detail")
px <- copula:::perspMvdc(x.samp, FUN = F.n, xlim = c(-5, 5), ylim = c(0, 2),
                         main = "F.n(x.samp)", ticktype="detail")
par(op)
all.equal(px, pp)# about 5% difference

Pseudo-Observations of Radial and Uniform Part of Elliptical and Archimedean Copulas

Description

Given a matrix of iid multivariate data from a meta-elliptical or meta-Archimedean model, RSpobs() computes pseudo-observations of the radial part R and the vector \bm{S} which follows a uniform distribution on the unit sphere (for elliptical copulas) or the unit simplex (for Archimedean copulas). These quantities can be used for (graphical) goodness-of-fit tests, for example.

Usage

RSpobs(x, do.pobs = TRUE, method = c("ellip", "archm"), ...)

Arguments

Details

The construction of the pseudo-obersvations of the radial part and the uniform distribution on the unit sphere/simplex is described in Genest, Hofert, G. Nešlehová (2014).

Value

A list with components R (an n-vector containing the pseudo-observations of the radial part) and S (an (n, d)-matrix containing the pseudo-observations of the uniform distribution (on the unit sphere/simplex)).

References

Genest, C., Hofert, M., G. Nešlehová, J., (2014). Is the dependence Archimedean, elliptical, or what? To be submitted.

See Also

pobs() for computing the “classical” pseudo-observations.

Examples

set.seed(100)
n <- 250 # sample size
d <- 5 # dimension
nu <- 3 # degrees of freedom

## Build a mean vector and a dispersion matrix,
## and generate multivariate t_nu data:
mu <- rev(seq_len(d)) # d, d-1, .., 1
L <- diag(d) # identity in dim d
L[lower.tri(L)] <- 1:(d*(d-1)/2)/d # Cholesky factor (diagonal > 0)
Sigma <- crossprod(L) # pos.-def. dispersion matrix (*not* covariance of X)
X <- rep(mu, each=n) + mvtnorm::rmvt(n, sigma=Sigma, df=nu) # multiv. t_nu data
## note: this is *wrong*: mvtnorm::rmvt(n, mean=mu, sigma=Sigma, df=nu)

## compute pseudo-observations of the radial part and uniform distribution
## once for 1a), once for 1b) below
RS.t    <- RSpobs(X, method="ellip", qQg=function(p) qt(p, df=nu)) # 'correct'
RS.norm <- RSpobs(X, method="ellip", qQg=qnorm) # for testing 'wrong' distribution
stopifnot(length(RS.norm$R) == n, length(RS.t$R) == n,
          dim(RS.norm$S) == c(n,d), dim(RS.t$S) == c(n,d))

## 1) Graphically testing the radial part

## 1a) Q-Q plot of R against the correct quantiles
qqplot2(RS.t$R, qF=function(p) sqrt(d*qf(p, df1=d, df2=nu)),
        main.args=list(text= substitute(bold(italic(F[list(d.,nu.)](r^2/d.))~~"Q-Q Plot"),
                                        list(d.=d, nu.=nu)),
		       side=3, cex=1.3, line=1.1, xpd=NA))

## 1b) Q-Q plot of R against the quantiles of F_R for a multivariate normal
##     distribution
qqplot2(RS.norm$R, qF=function(p) sqrt(qchisq(p, df=d)),
        main.args=list(text= substitute(bold(italic(chi[D_]) ~~ "Q-Q Plot"), list(D_=d)),
               side=3, cex=1.3, line=1.1, xpd=NA))

## 2) Graphically testing the angular distribution

## auxiliary function
qqp <- function(k, Bmat) {
    d <- ncol(Bmat) + 1
    qqplot2(Bmat[,k],
            qF = function(p) qbeta(p, shape1=k/2, shape2=(d-k)/2),
            main.args=list(text= substitute(plain("Beta")(s1,s2) ~~ bold("Q-Q Plot"),
                                            list(s1 = k/2, s2 = (d-k)/2)),
  	                   side=3, cex=1.3, line=1.1, xpd=NA))
}

## 2a) Q-Q plot of the 'correct' angular distribution
##     (Bmat[,k] should follow a Beta(k/2, (d-k)/2) distribution)
Bmat.t <- gofBTstat(RS.t$S)
qqp(1, Bmat=Bmat.t) # k=1
qqp(3, Bmat=Bmat.t) # k=3

## 2b) Q-Q plot of the 'wrong' angular distribution
Bmat.norm <- gofBTstat(RS.norm$S)
qqp(1, Bmat=Bmat.norm) # k=1
qqp(3, Bmat=Bmat.norm) # k=3

## 3) Graphically check independence between radial part and B_1 and B_3

## 'correct' distributions (multivariate t)
plot(pobs(cbind(RS.t$R, Bmat.t[,1])), # k = 1
          xlab=quote(italic(R)), ylab=quote(italic(B)[1]),
          main=quote(bold("Rank plot between"~~italic(R)~~"and"~~italic(B)[1])))
plot(pobs(cbind(RS.t$R, Bmat.t[,3])), # k = 3
	  xlab=quote(italic(R)), ylab=quote(italic(B)[3]),
          main=quote(bold("Rank plot between"~~italic(R)~~"and"~~italic(B)[3])))

## 'wrong' distributions (multivariate normal)
plot(pobs(cbind(RS.norm$R, Bmat.norm[,1])), # k = 1
          xlab=quote(italic(R)), ylab=quote(italic(B)[1]),
          main=quote(bold("Rank plot between"~~italic(R)~~"and"~~italic(B)[1])))
plot(pobs(cbind(RS.norm$R, Bmat.norm[,3])), # k = 3
	  xlab=quote(italic(R)), ylab=quote(italic(B)[3]),
          main=quote(bold("Rank plot between"~~italic(R)~~"and"~~italic(B)[3])))

SMI Data – 141 Days in Winter 2011/2012

Description

SMI.12 contains the close prices of all 20 constituents of the Swiss Market Index (SMI) from 2011-09-09 to 2012-03-28.

Usage

data(SMI.12)

Format

SMI.12 is conceptually a multivariate time series, here simply stored as numeric matrix, where the rownames are dates (of week days).

The format is:

num [1:141, 1:20] 16.1 15.7 15.7 16.1 16.6 ... - attr(*, "dimnames")=List of 2 ..$ : chr [1:141] "2011-09-09" "2011-09-12" "2011-09-13" "2011-09-14" ... ..$ : chr [1:20] "ABBN" "ATLN" "ADEN" "CSGN" ...

... from 2011-09-09 to 2012-03-28

lSMI is the list of the original data (before NA “imputation”).

Source

The data was drawn from Yahoo! Finance.

Examples

data(SMI.12)
## maybe
head(SMI.12)

str(D.12 <- as.Date(rownames(SMI.12)))
summary(D.12)

matplot(D.12, SMI.12, type="l", log = "y",
        main = "The 20 SMI constituents  (2011-09 -- 2012-03)",
        xaxt="n", xlab = "2011  /  2012")
Axis(D, side=1)



if(FALSE) { ##--- This worked up to mid 2012, but no longer ---
 begSMI <- "2011-09-09"
 endSMI <- "2012-03-28"
 ##-- read *public* data ------------------------------
 stopifnot(require(zoo), # -> to access all the zoo methods
           require(tseries))
 symSMI <- c("ABBN.VX","ATLN.VX","ADEN.VX","CSGN.VX","GIVN.VX","HOLN.VX",
	     "BAER.VX","NESN.VX","NOVN.VX","CFR.VX", "ROG.VX", "SGSN.VX",
	     "UHR.VX", "SREN.VX","SCMN.VX","SYNN.VX","SYST.VX","RIGN.VX",
	     "UBSN.VX","ZURN.VX")
 lSMI <- sapply(symSMI, function(sym)
		get.hist.quote(instrument = sym, start= begSMI, end= endSMI,
			       quote = "Close", provider = "yahoo",
			       drop=TRUE))
 ## check if stock data have the same length for each company.
 sapply(lSMI, length)
 ## "concatenate" all:
 SMIo <- do.call(cbind, lSMI)
 ## and fill in the NAs :
 SMI.12 <- na.fill(SMIo, "extend")
 colnames(SMI.12) <- sub("\\.VX", "", colnames(SMI.12))
 SMI.12 <- as.matrix(SMI.12)
}##----       --- original download

zoo.there <- "package:zoo" %in% search()
if(zoo.there || require("zoo")) {
  stopifnot(identical(SMI.12,
     local({ S <- as.matrix(na.fill(do.call(cbind, lSMI), "extend"))
             colnames(S) <- sub("\\.VX", "", colnames(S)); S })))
  if(!zoo.there) detach("package:zoo")
}

Sibuya Distribution - Sampling and Probabilities

Description

The Sibuya distribution \mathrm{Sib}(\alpha) can be defined by its Laplace transform

1-(1-\exp(-t))^\alpha,\ t\in[0,\infty),

its distribution function

F(k)=1-(-1)^k{\alpha-1\choose k}=1-\frac{1}{kB(k,1-\alpha)},\ k\in\mathbf{N}

(where B denotes the beta function) or its probability mass function

p_k={\alpha\choose k}(-1)^{k-1},\ k\in\mathbf{N},

where \alpha\in(0,1].

pSibuya evaluates the distribution function.

dSibuya evaluates the probability mass function.

rSibuya generates random variates from \mathrm{Sib}(\alpha) with the algorithm described in Hofert (2011), Proposition 3.2.

dsumSibuya gives the probability mass function of the n-fold convolution of Sibuya variables, that is, the sum of n independent Sibuya random variables, S = \sum_{i=1}^n X_i, where X_i \sim \mathrm{Sib}(\alpha).

This probability mass function can be shown (see Hofert (2010, pp. 99)) to be

\sum_{j=1}^n{n\choose j}{j\alpha\choose k} (-1)^{k-j},\ k\in\{n,n+1,\dots\}.

Usage

rSibuya(n, alpha)
dSibuya(x, alpha, log=FALSE)
pSibuya(x, alpha, lower.tail=TRUE, log.p=FALSE)

dsumSibuya(x, n, alpha,
           method=c("log", "direct", "diff", "exp.log",
                    "Rmpfr", "Rmpfr0", "RmpfrM", "Rmpfr0M"),
           mpfr.ctrl = list(minPrec = 21, fac = 1.25, verbose=TRUE),
           log=FALSE)

Arguments

Details

The Sibuya distribution has no finite moments, that is, specifically infinite mean and variance.

For documentation and didactical purposes, rSibuyaR is a pure-R implementation of rSibuya, of course slower than rSibuya as the latter is implemented in C.

Note that the sum to evaluate for dsumSibuya is numerically highly challenging, even already for small \alpha values (for example, n \ge 10), and therefore should be used with care. It may require high-precision arithmetic which can be accessed with method="Rmpfr" (and the Rmpfr package).

Value

rSibuya:

A vector of positive integers of length n containing the generated random variates.

dSibuya, pSibuya:

a vector of probabilities of the same length as x.

dsumSibuya:

a vector of probabilities, positive if and only if x >= n and of the same length as x (or n if that is longer).

References

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

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

See Also

rFJoe and rF01Joe (where rSibuya is applied).

Examples

## Sample n random variates from a Sibuya(alpha) distribution and plot a
## histogram
n <- 1000
alpha <- .4
X <- rSibuya(n, alpha)
hist(log(X), prob=TRUE); lines(density(log(X)), col=2, lwd=2)

Eulerian and Stirling Numbers of First and Second Kind

Description

Compute Eulerian numbers and Stirling numbers of the first and second kind, possibly vectorized for all k “at once”.

Usage

Stirling1(n, k)
Stirling2(n, k, method = c("lookup.or.store", "direct"))
Eulerian (n, k, method = c("lookup.or.store", "direct"))

Stirling1.all(n)
Stirling2.all(n)
Eulerian.all (n)

Arguments

Details

Eulerian numbers:
A(n,k) = the number of permutations of 1,2,...,n with exactly k ascents (or exactly k descents).

Stirling numbers of the first kind:
s(n,k) = (-1)^{n-k} times the number of permutations of 1,2,...,n with exactly k cycles.

Stirling numbers of the second kind:
S^{(k)}_n is the number of ways of partitioning a set of n elements into k non-empty subsets.

Value

A(n,k), s(n,k) or S(n,k) = S^{(k)}_n, respectively.

Eulerian.all(n) is the same as sapply(0:(n-1), Eulerian, n=n) (for n > 0),
Stirling1.all(n) is the same as sapply(1:n, Stirling1, n=n), and
Stirling2.all(n) is the same as sapply(1:n, Stirling2, n=n), but more efficient.

Note

For typical double precision arithmetic,
Eulerian*(n, *) overflow (to Inf) for n \ge 172,
Stirling1*(n, *) overflow (to \pmInf) for n \ge 171, and
Stirling2*(n, *) overflow (to Inf) for n \ge 220.

References

Eulerians:

NIST Digital Library of Mathematical Functions, 26.14: https://dlmf.nist.gov/26.14

Stirling numbers:

Abramowitz and Stegun 24,1,4 (p. 824-5 ; Table 24.4, p.835); Closed Form : p.824 "C."

NIST Digital Library of Mathematical Functions, 26.8: https://dlmf.nist.gov/26.8

Examples

Stirling1(7,2)
Stirling2(7,3)

Stirling1.all(9)
Stirling2.all(9)

Absolute Value of Generator Derivatives via Monte Carlo

Description

Computes the absolute values of the dth 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 dth 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...

Distribution of the Radial Part of an Archimedean Copula

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])

Class "acopula" of Archimedean Copula Families

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")

All Components of a (Inner or Outer) Nested Archimedean Copula

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)

Construction of Archimedean Copula Class Object

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)

Class "archmCopula"

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.

Dependence Measures for Bivariate Copulas

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])

Sample and Population Version of Blomqvist's Beta for Archimedean Copulas

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")

Conditional Distributions and Their Inverses from Copulas

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

(Xtras <- copula:::doExtras()) # determine whether examples will be extra (long)

## 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 <- if(Xtras) 1000 else 250
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 <- if(Xtras) 6 else 3
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 <- if(Xtras) 1000 else 400
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

Cloud Plot Methods ('cloud2') in Package 'copula'

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)

Coefficients of Polynomial used for Gumbel Copula

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

% latex 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)

Methods for Contour Plots in Package 'copula'

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))

Contour Plot Methods 'contourplot2' in Package 'copula'

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"))

Specific Archimedean Copula Families ("acopula" Objects)

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)")

Mother Classes "Copula", etc of all Copulas in the Package

Description

A copula is a multivariate distribution with uniform margins. The virtual class "Copula" is the mother (or “super class”) of all copula classes in the package copula which encompasses classes of the former packages nacopula and copula.

The virtual class "parCopula" extends "Copula" and is the super class of all copulas that can be fully parametrized and hence fitted to data. For these, at least the dim() method must be well defined.

The virtual class "dimCopula" extends "Copula" and has an explicit slot dimension, with corresponding trivial dim() method.

The virtual class "copula" extends bot "dimCopula" and "parCopula" and is the mother of all copula classes from former package copula. It has set of slots for (the dimension and) parameter vector, see below.

The virtual class "Xcopula" contains a slot copula of class "parCopula".

The virtual class "xcopula" extends "parCopula" and "Xcopula"; an (“actual”) class example are the rotated copulas, rotCopula.

Objects from the Class

Objects are typically created by are by tCopula(), evCopula(), etc.

Slots

Class "dimCopula" and its subclasses, notably "copula", have a slot

dimension:

an "integer" (of length 1), the copula dimension d.

Class "copula" (and its subclasses) have additional slots

parameters:

numeric vector of parameter values, can be NA (i.e., NA_real_).

param.names:

"character" vector of parameter names (and hence of the same length as parameters).

param.lowbnd:

lower bounds for the parameters, of class "numeric".

param.upbnd:

upper bounds for the parameters, of class "numeric".

fullname:

deprecated; object of class "character", family names of the copula.

Warning

This implementation is still at the experimental stage and is subject to change during the development.

Note

The "copula" class is extended by the evCopula, archmCopula, and ellipCopula classes. Instances of such copulas can be created via functions evCopula, archmCopula and ellipCopula.

"plackettCopula" and fgmCopula are special types of copulas which do not belong to either one of the three classes above.

See Also

Help for the (sub)classes archmCopula, ellipCopula, evCopula, and fgmCopula.

The Archimedean and nested Archimedean classes (from former package nacopula), with a more extensive list of slots (partly instead of methods), acopula, and nacopula.

Examples

hc <- evCopula("husler", 1.25)
dim(hc)
smoothScatter(u <- rCopula(2^11, hc))
lambda   (hc)
tau (hc)
rho(hc)
str(hc)

Internal Copula Functions

Description

Internal Copula functions

Details

These are not to be called by the user (or in some cases are just waiting for proper documentation to be written).

(Fast) Computation of Pairwise Kendall's Taus

Description

For a data matrix x, compute the Kendall's tau “correlation” matrix, i.e., all pairwise Kendall's taus between the columns of x.

By default and when x has no missing values (NAs), the fast O(n log(n)) algorithm of cor.fk() is used.

Usage

corKendall(x, checkNA = TRUE,
           use = if(checkNA && anyNA(x)) "pairwise" else "everything")

Arguments

Value

The p \times p matrix K of pairwise Kendall's taus, with K[i,j] := tau(x[,i], x[,j]).

See Also

cor.fk() from pcaPP (used by default when there are no missing values (NAs) in x).

etau() or fitCopula(*, method = "itau") make use of corKendall().

Examples

## If there are no NA's, corKendall() is faster than cor(*, "kendall")
## and gives the same :

system.time(C1 <- cor(swiss, method="kendall"))
system.time(C2 <- corKendall(swiss))
stopifnot(all.equal(C1, C2,  tol = 1e-5))

## In the case of missing values (NA), corKendall() reverts to
## cor(*, "kendall", use = "pairwise") {no longer very fast} :

swM <- swiss # shorter names and three missings:
colnames(swM) <- abbreviate(colnames(swiss), min=6)
swM[1,2] <- swM[7,3] <- swM[25,5] <- NA
(C3 <- corKendall(swM)) # now automatically uses the same as
stopifnot(identical(C3, cor(swM, method="kendall", use="pairwise")))
## and is quite close to the non-missing "truth":
stopifnot(all.equal(unname(C3), unname(C2), tol = 0.06)) # rel.diff.= 0.055

try(corKendall(swM, checkNA=FALSE)) # --> Error
## the error is really from  pcaPP::cor.fk(swM)

Density of the Diagonal of (Nested) Archimedean Copulas

Description

Evaluate the density of the diagonal of a d-dimensional (nested) Archimedean copula. Note that the diagonal of a copula is a cumulative distribution function. Currently, only Archimedean copulas are implemented.

Usage

dDiag(u, cop, log=FALSE)

Arguments

Value

A numeric vector containing the values of the density of the diagonal of the Archimedean copula at u.

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 class, dnacopula.

Examples

th. <- c(0.1, 0.2, 0.5, 0.8, 1.4, 2., 5.)
curve(dDiag(x, cop=onacopulaL("Clayton", list(th.[1], 1:3))), 0, 1,
      n=1000, ylab="dDiag(x, *)", main="Diagonal densities of Clayton")
abline(h=0, lty=3)
for(j in 2:length(th.))
  curve(dDiag(x, cop=onacopulaL("Clayton", list(th.[j], 1:3))), add=TRUE,
	     col=j, n=1000)
legend("topleft", do.call(expression, lapply(th., function(th)
                                 substitute(theta == TH, list(TH=th)))),
       lty = 1, col=seq_along(th.), bty="n")

Copula (Short) Description as String

Description

Describe a copula object, i.e., its basic properties as a string. This is a utility used when print()ing or plot()ting copulas, e.g., after a fitting.

Usage

describeCop(x, kind = c("short", "very short", "long"), prefix = "", ...)

Arguments

Value

a character string.

Methods

signature(x = "archmCopula", kind = "ANY")

..

signature(x = "copula", kind = "character")

..

signature(x = "copula", kind = "missing")

..

signature(x = "ellipCopula", kind = "character")

..

signature(x = "fgmCopula", kind = "ANY")

..

signature(x = "xcopula", kind = "ANY")

..

See Also

Copula class definition copula;

Examples

## FIXME

Density Evaluation for (Nested) Archimedean Copulas

Description

For a (nested) Archimedean copula (object of class nacopula) x, dCopula(u, x) (or also currently still dnacopula(x, u)) evaluates the density of x at the given vector or matrix u.

Usage

## S4 method for signature 'matrix,nacopula'
dCopula(u, copula, log=FALSE, ...)

## *Deprecated*:
dnacopula(x, u, log=FALSE, ...)

Arguments

Details

If it exists, the density of an Archimedean copula C with generator \psi at \bm{u}\in(0,1)^d is given by

c(\bm{u})=\psi^{(d)}(\psi^{-1}(u_1)+\dots+\psi^{-1}(u_d)) \prod_{j=1}^d(\psi^{-1}(u_j))^\prime = \frac{\psi^{(d)}(\psi^{-1}(u_1)+\dots+\psi^{-1}(u_d))}{ \prod_{j=1}^d\psi^\prime(\psi^{-1}(u_j))}.

Value

A numeric vector containing the values of the density of the Archimedean copula at u.

Note

dCopula(u, copula) is a generic function with methods for all our copula classes, see dCopula.

References

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.

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

For more details about the derivatives of an Archimedean generator, see, for example, absdPsi in class acopula.

Examples

## Construct a twenty-dimensional Gumbel copula with parameter chosen
## such that Kendall's tau of the bivariate margins is 0.25.
theta <- copJoe@iTau(.25)
C20 <- onacopula("J", C(theta, 1:20))

## Evaluate the copula density at the point u = (0.5,...,0.5)
u <- rep(0.5, 20)
dCopula(u, C20)

## the same with Monte Carlo based on 10000 simulated "frailties"
dCopula(u, C20, n.MC = 10000)

## Evaluate the exact log-density at several points
u <- matrix(runif(100), ncol=20)
dCopula(u, C20, log = TRUE)

## Back-compatibility check
stopifnot(identical( dCopula (u, C20), suppressWarnings(
                    dnacopula(C20, u))),
          identical( dCopula (u, C20, log = TRUE), suppressWarnings(
                    dnacopula(C20, u, log = TRUE))))

Construction of Elliptical Copula Class Objects

Description

Creating elliptical copula objects with corresponding dimension and parameters, including the dispersion structure P (pronounced “Rho”).

Usage

ellipCopula (family, param, dim = 2, dispstr = "ex", df = 4, ...)

normalCopula(param, dim = 2, dispstr = "ex")

     tCopula(param, dim = 2, dispstr = "ex", df = 4,
             df.fixed = FALSE, df.min = 0.01)

dispstrToep(perm = NULL, check = TRUE)

## S4 method for signature 'matrix,normalCopula'
pCopula(u, copula, algorithm=NULL, keepAttr=FALSE, ...)
## S4 method for signature 'matrix,tCopula'
pCopula(u, copula, algorithm=NULL, keepAttr=FALSE, ...)

Arguments

Value

For

ellipCopula(), normalCopula(), or tCopula():

an elliptical copula object of class "normalCopula" or "tCopula".

dispstrToep():

the character string "toep", optionally with attribute (see attributes, attr) "perm" with a permutation p of 1:d, such that (the column permuted) U_p, or in the data case the column-permuted matrix U[,p] has as dispersion matrix toeplitz(c(1, par)), with par the respective parameter vector of bivariate “correlations” \rho_{i,j}.

Note that the result of dispstrToep() is currently stored in the dispstr slot of the copula object.

pCopula(u, *):

the numerical vector of copula probabilites of length nrow(u).

Note

ellipCopula() is a wrapper for normalCopula() and tCopula().

The pCopula() methods for the normal- and t-copulas accept optional arguments to be passed to the underlying (numerical integration) algorithms from package mvtnorm's pmvnorm and pmvt, respectively, notably algorithm, see GenzBretz, or abseps which defaults to 0.001.

See Also

p2P(), and getSigma() for construction and extraction of the dispersion matrix P or Sigma matrix of (generalized) correlations.

archmCopula, fitCopula.

Examples

normalCopula(c(0.5, 0.6, 0.7), dim = 3, dispstr = "un")
t.cop <- tCopula(c(0.5, 0.3), dim = 3, dispstr = "toep",
                 df = 2, df.fixed = TRUE)
getSigma(t.cop) # P matrix (with diagonal = 1)
stopifnot(all.equal(toeplitz(c(1, .5, .3)), getSigma(t.cop)))

## dispersion "AR1" :
nC.7 <- normalCopula(0.8, dim = 7, dispstr = "ar1")
getSigma(nC.7)
stopifnot(all.equal(toeplitz(.8^(0:6)), getSigma(nC.7)))

## from the wrapper
norm.cop <- ellipCopula("normal", param = c(0.5, 0.6, 0.7),
                        dim = 3, dispstr = "un")
if(require("scatterplot3d") && dev.interactive(orNone=TRUE)) {
  ## 3d scatter plot of 1000 random observations
  scatterplot3d(rCopula(1000, norm.cop))
  scatterplot3d(rCopula(1000, t.cop))
}
set.seed(12)
uN <- rCopula(512, norm.cop)
set.seed(2); pN1 <- pCopula(uN, norm.cop)
set.seed(3); pN2 <- pCopula(uN, norm.cop)
stopifnot(identical(pN1, pN2)) # no longer random for dim = 3
(Xtras <- copula:::doExtras())
if(Xtras) { ## a bit more accurately:
  set.seed(4); pN1. <- pCopula(uN, norm.cop, abseps = 1e-9)
  set.seed(5); pN2. <- pCopula(uN, norm.cop, abseps = 1e-9)
  stopifnot(all.equal(pN1., pN2., 1e-5))# see 3.397e-6
  ## but increasing the required precision (e.g., abseps=1e-15) does *NOT* help
}

## For smaller copula dimension 'd', alternatives are available and
## non-random, see ?GenzBretz from package 'mvtnorm' :
has_mvtn <- "package:mvtnorm" %in% search() #% << (workaround ESS Rd render bug)
if(!has_mvtn)
  require("mvtnorm")# -> GenzBretz(), Miva(), and TVPACK() are available
## Note that Miwa() would become very slow for dimensions 5, 6, ..
set.seed(4); pN1.M <- pCopula(uN, norm.cop, algorithm = Miwa(steps = 512))
set.seed(5); pN2.M <- pCopula(uN, norm.cop, algorithm = Miwa(steps = 512))
stopifnot(all.equal(pN1.M, pN2.M, tol= 1e-15))# *no* randomness
set.seed(4); pN1.T <- pCopula(uN, norm.cop, algorithm = TVPACK(abseps = 1e-10))
set.seed(5); pN2.T <- pCopula(uN, norm.cop, algorithm = TVPACK(abseps = 1e-14))
stopifnot(all.equal(pN1.T, pN2.T, tol= 1e-15))# *no* randomness (but no effect of 'abseps')
if(!has_mvtn)
   detach("package:mvtnorm")# (revert)


## Versions with unspecified parameters:
tCopula()
allEQ <- function(u,v) all.equal(u, v, tolerance=0)
stopifnot(allEQ(ellipCopula("norm"), normalCopula()),
          allEQ(ellipCopula("t"), tCopula()))
tCopula(dim=3)
tCopula(dim=4, df.fixed=TRUE)
tCopula(dim=5, disp = "toep", df.fixed=TRUE)
normalCopula(dim=4, disp = "un")

## Toeplitz after *permutation* dispersions (new in copula 1.1-0) ---------
tpar <- c(7,5,3)/8 # *gives* pos.def.:
(ev <- eigen(toeplitz(c(1, tpar)), symmetric=TRUE, only.values=TRUE)$values)
stopifnot(ev > 0)
N4.   <- ellipCopula("normal", dim=4, param=tpar, dispstr = "toep") #"regular"
## reversed order is "the same" for toeplitz structure:
N4.pr <- ellipCopula("normal", dim=4, param=tpar, dispstr = dispstrToep(4:1))
N4.p1 <- ellipCopula("normal", dim=4, param=tpar, dispstr = dispstrToep(c(4,1:3)))
N4.p2 <- ellipCopula("normal", dim=4, param=tpar, dispstr = dispstrToep(c(4:3,1:2)))
N4.p3 <- ellipCopula("normal", dim=4, param=tpar, dispstr = dispstrToep(c(2,4,1,3)))

(pm <- attr(N4.p3@dispstr, "perm")) # (2 4 1 3)
ip <- c(3,1,4,2) # the *inverse* permutation of (2 4 1 3) = Matrix::invPerm(pm)
(Sp3 <- getSigma(N4.p3)) # <-- "permuted toeplitz"
Sp3[ip, ip] # re-ordered rows & columns => *is* toeplitz :
stopifnot(exprs = {
  ## permutations  pm  and  ip  are inverses:
  pm[ip] == 1:4
  ip[pm] == 1:4
  is.matrix(T4 <- toeplitz(c(1, tpar)))
  identical(getSigma(N4.),   T4)
  identical(getSigma(N4.pr), T4) # 4:1 and 1:4 is "the same" for Rho
  identical(Sp3[ip, ip]  ,   T4)
  identical(Sp3,      T4[pm,pm])
})
## Data generation -- NB: The U matrices are equal only "in distribution":
set.seed(7); U.p3 <- rCopula(1000, N4.p3)
set.seed(7); U.   <- rCopula(1000, N4.)
stopifnot(exprs = {
 all.equal(loglikCopula(tpar, u=U.p3,      copula= N4.p3),
           loglikCopula(tpar, u=U.p3[,ip], copula= N4.) -> LL3)
 all.equal(loglikCopula(tpar, u=U.,      copula= N4.),
           loglikCopula(tpar, u=U.[,pm], copula= N4.p3) -> LL.)
})
c(LL. , LL3)# similar but different
if(Xtras) {
  fm.  <- fitCopula(N4.  , U.  )
  fm.3 <- fitCopula(N4.p3, U.p3)
  summary(fm.3)
  stopifnot(all.equal(coef(fm.), coef(fm.3), tol = 0.01))# similar but different
}

Class "ellipCopula" of Elliptical Copulas

Description

Copulas generated from elliptical multivariate distributions, notably Normal- and t-copulas (of specific class "normalCopula" or "tCopula", respectively).

Objects from the Class

Objects are typically created by ellipCopula(), normalCopula(), or tCopula().

Slots

dispstr:

"character" string indicating how the dispersion matrix is parameterized; one of "ex", "ar1", "toep", or "un", see the dispstr argument of ellipCopula().

dimension:

Object of class "numeric", dimension of the copula.

parameters:

a numeric, (vector of) the parameter value(s).

param.names:

character vector with names for the parameters slot, of the same length.

param.lowbnd:

numeric vector of lower bounds for the parameters slot, of the same length.

param.upbnd:

upper bounds for parameters, analogous to parm.lowbnd.

fullname:

deprecated; object of class "character", family names of the copula.

Extends

Class "ellipCopula" extends class copula directly. Classes "normalCopula" and "tCopula" extend "ellipCopula" directly.

Methods

Many methods are available, notably dCopula, pCopula, and rCopula. Use, e.g., methods(class = "tCopula") to find others.

See Also

ellipCopula which also documents tCopula() and normalCopula(); copula-class.

Minimum Distance Estimators for (Nested) Archimedean Copulas

Description

Compute minimum distance estimators for (nested) Archimedean copulas.

Usage

emde(u, cop,
     method=c("mde.chisq.CvM", "mde.chisq.KS",
              "mde.gamma.CvM", "mde.gamma.KS"),
     interval=initOpt(cop@copula@name),
     include.K = FALSE, repara = TRUE, ...)

Arguments

Details

First, htrafo is applied to map the n\times d-matrix of given realizations to a n\times d-matrix or n\times (d-1)-matrix, depending on whether the last component is included (include.K=TRUE) or not. Second, using either the sum of squares of the standard normal quantile function (method="mde.chisq.CvM" and method="mde.chisq.KS") or the sum of negative logarithms (method="mde.gamma.CvM" and method="mde.gamma.KS"), a map to a chi-square or an Erlang distribution is applied, respectively. Finally, a Cramér-von Mises (method="mde.chisq.CvM" and method="mde.gamma.CvM") or Kolmogorov-Smirnov (method="mde.chisq.KS" and method="mde.gamma.KS") distance is applied. This is repeated in an optimization until the copula parameter is found such that this distance is minimized.

Note that the same transformations as described above are applied for goodness-of-fit testing; see the ‘See Also’ section).

Value

list as returned by optimize, including the minimum distance estimator.

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.

Hering, C. and Hofert, M. (2014), Goodness-of-fit tests for Archimedean copulas in high dimensions, Innovations in Quantitative Risk Management.

See Also

enacopula (wrapper for different estimators), gofCopula (wrapper for different goodness-of-fit tests), htrafo (transformation to a multivariate uniform distribution), and K (Kendall distribution function).

Examples

tau <- 0.25
(theta <- copGumbel@iTau(tau)) # 4/3
d <- 20
(cop <- onacopulaL("Gumbel", list(theta,1:d)))

set.seed(1)
n <- 200
U <- rnacopula(n, cop)

(meths <- eval(formals(emde)$method)) # "mde.chisq.CvM", ...
fun <- function(meth, u, cop, theta){
	run.time <- system.time(val <- emde(u, cop=cop, method=meth)$minimum)
	list(value=val, error=val-theta, utime.ms=1000*run.time[[1]])
}
(res <- sapply(meths, fun, u=U, cop=cop, theta=theta))

Maximum Likelihood Estimators for (Nested) Archimedean Copulas

Description

Compute (simulated) maximum likelihood estimators for (nested) Archimedean copulas.

Usage

emle(u, cop, n.MC=0, optimizer="optimize", method,
     interval=initOpt(cop@copula@name),
     start=list(theta=initOpt(cop@copula@name, interval=FALSE, u=u)),
     ...)
.emle(u, cop, n.MC=0,
      interval=initOpt(cop@copula@name), ...)

Arguments

Details

Exact formulas for the generator derivatives were derived in Hofert et al. (2012). Based on these formulas one can compute the (log-)densities of the Archimedean copulas. Note that for some densities, the formulas are numerically highly non-trivial to compute and considerable efforts were put in to make the computations numerically feasible even in large dimensions (see the source code of the Gumbel copula, for example). Both MLE and SMLE showed good performance in the simulation study conducted by Hofert et al. (2013) including the challenging 100-dimensional case. Alternative estimators (see also enacopula) often used because of their numerical feasibility, might break down in much smaller dimensions.

Note: SMLE for Clayton currently faces serious numerical issues and is due to further research. This is only interesting from a theoretical point of view, since the exact derivatives are known and numerically non-critical to evaluate.

Value

emle

an R object of class "mle2" (and thus useful for obtaining confidence intervals) with the (simulated) maximum likelihood estimator.

.emle

list as returned by optimize() including the maximum likelihood estimator (does not confidence intervals but is typically faster).

References

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.

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

mle2 from package bbmle and mle from stats4 on which mle2 is modeled. enacopula (wrapper for different estimators). demo(opC-demo) and demo(GIG-demo) for examples of two-parameter families.

Examples

tau <- 0.25
(theta <- copGumbel@iTau(tau)) # 4/3
d <-  20
(cop <- onacopulaL("Gumbel", list(theta,1:d)))

set.seed(1)
n <- 200
U <- rnacopula(n,cop)

## Estimation
system.time(efm <- emle(U, cop))
summary(efm) # using bblme's 'mle2' method

## Profile likelihood plot [using S4 methods from bbmle/stats4] :
pfm <- profile(efm)
ci  <- confint(pfm, level=0.95)
ci
stopifnot(ci[1] <= theta, theta <= ci[2])
plot(pfm)               # |z| against theta, |z| = sqrt(deviance)
plot(pfm, absVal=FALSE, #  z  against theta
     show.points=TRUE) # showing how it's interpolated
## and show the true theta:
abline(v=theta, col="lightgray", lwd=2, lty=2)
axis(1, pos = 0, at=theta, label=quote(theta[0]))

## Plot of the log-likelihood, MLE  and  conf.int.:
logL <- function(x) -efm@minuslogl(x)
       # == -sum(copGumbel@dacopula(U, theta=x, log=TRUE))
logL. <- Vectorize(logL)
I <- c(cop@copula@iTau(0.1), cop@copula@iTau(0.4))
curve(logL., from=I[1], to=I[2], xlab=quote(theta),
      ylab="log-likelihood",
      main="log-likelihood for Gumbel")
abline(v = c(theta, efm@coef), col="magenta", lwd=2, lty=2)
axis(1, at=c(theta, efm@coef), padj = c(-0.5, -0.8), hadj = -0.2,
     col.axis="magenta", label= expression(theta[0], hat(theta)[n]))
abline(v=ci, col="gray30", lwd=2, lty=3)
text(ci[2], extendrange(par("usr")[3:4], f= -.04)[1],
     "95% conf. int.", col="gray30", adj = -0.1)

The Empirical Copula

Description

Computes the empirical copula (according to a provided method) and auxiliary tools.

Usage

empCopula(X, smoothing = c("none", "beta", "checkerboard",
                           "schaake.shuffle"), offset = 0,
          ties.method = c("max", "average", "first", "last", "random", "min"))
C.n(u, X, smoothing = c("none", "beta", "checkerboard"), offset = 0,
    ties.method = c("max", "average", "first", "last", "random", "min"))
dCn(u, U, j.ind = 1:d, b = 1/sqrt(nrow(U)), ...)
F.n(x, X, offset = 0, smoothing = c("none", "beta", "checkerboard"))
Cn(x, w) ## <-- deprecated!  use  C.n(w, x) instead!
toEmpMargins(U, x, ...)

Arguments

Details

Given pseudo-observations from a distribution with continuous margins and copula C, the empirical copula is the (default) empirical distribution function of these pseudo-observations. It is thus a natural nonparametric estimator of C. The function C.n() computes the empirical copula or two alternative smoothed versions of it: the empirical beta copula or the empirical checkerboard copula; see Eqs. (2.1) and (4.1) in Segers, Sibuya and Tsukahara (2017), and the references therein. empCopula() is the constructor of an object of class empCopula.

The function dCn() approximates first-order partial derivatives of the unknown copula using the empirical copula.

The function F.n() computes the empirical distribution function of a multivariate sample. Note that C.n(u, X, smoothing="none", *) simply calls F.n(u, pobs(X), *) after checking u.

There are several asymptotically equivalent definitions of the empirical copula. C.n(, smoothing = "none") is simply defined as the empirical distribution function computed from the pseudo-observations, that is,

C_n(\bm{u})=\frac{1}{n}\sum_{i=1}^n\mathbf{1}_{\{\hat{\bm{U}}_i\le\bm{u}\}},

where \hat{\bm{U}}_i, i\in\{1,\dots,n\}, denote the pseudo-observations and n the sample size. Internally, C.n(,smoothing = "none") is just a wrapper for F.n() and is expected to be fed with the pseudo-observations.

The approximation for the jth partial derivative of the unknown copula C is implemented as, for example, in Rémillard and Scaillet (2009), and given by

\hat{\dot{C}}_{jn}(\bm{u})=\frac{C_n(u_1,..,u_{j-1},min(u_j+b,1),u_{j+1},..,u_d)-C_n(u_1,..,u_{j-1},max(u_j-b,0),u_{j+1},..,u_d)}{2b},

where b denotes the bandwidth and C_n the empirical copula.

Value

empCopula() is the constructor for objects of class empCopula.

C.n() returns the empirical copula of the pseudo-observations X evaluated at u (or a smoothed version of it).

dCn() returns a vector (length(j.ind) is 1) or a matrix (with number of columns equal to length(j.ind)), containing the approximated first-order partial derivatives of the unknown copula at u with respect to the arguments in j.ind.

F.n() returns the empirical distribution function of X evaluated at x if smoothing = "none", the empirical beta copula evaluated at x if smoothing = "beta" and the empirical checkerboard copula evaluated at x if smoothing = "checkerboard".

toEmpMargins() transforms the copula sample U to the empirical margins based on the sample x.

Note

The first version of our empirical copula implementation, Cn(), had its two arguments reversed compared to C.n(), and is deprecated now. You must swap its arguments to transform to new code.

The use of the two smoothed versions assumes implicitly no ties in the component samples of the data.

References

Rüschendorf, L. (1976). Asymptotic distributions of multivariate rank order statistics, Annals of Statistics 4, 912–923.

Deheuvels, P. (1979). La fonction de dépendance empirique et ses propriétés: un test non paramétrique d'indépendance, Acad. Roy. Belg. Bull. Cl. Sci., 5th Ser. 65, 274–292.

Deheuvels, P. (1981). A non parametric test for independence, Publ. Inst. Statist. Univ. Paris 26, 29–50.

Clark, M., Gangopadhyay, S., Hay, L., Rajagopalan, B. and Wilby, R. (2004). The Schaake Shuffle: A Method for Reconstructing Space-Time Variability in Forecasted Precipitation and Temperature Fields. Journal of Hydrometeorology, pages 243-262.

Rémillard, B. and Scaillet, O. (2009). Testing for equality between two copulas. Journal of Multivariate Analysis, 100(3), pages 377-386.

Segers, J., Sibuya, M. and Tsukahara, H. (2017). The Empirical Beta Copula. Journal of Multivariate Analysis, 155, pages 35–51, https://arxiv.org/abs/1607.04430.

Kiriliouk, A., Segers, J. and Tsukahara, H. (2020). Resampling Procedures with Empirical Beta Copulas. https://arxiv.org/abs/1905.12466.

See Also

pobs() for computing pseudo-observations.

Examples

## Generate data X (from a meta-Gumbel model with N(0,1) margins)
n <- 100
d <- 3
family <- "Gumbel"
theta <- 2
cop <- onacopulaL(family, list(theta=theta, 1:d))
set.seed(1)
X <- qnorm(rCopula(n, cop)) # meta-Gumbel data with N(0,1) margins

## Evaluate empirical copula
u <- matrix(runif(n*d), n, d) # random points were to evaluate the empirical copula
ec <- C.n(u, X = X)

## Compare the empirical copula with the true copula
pc <- pCopula(u, copula = cop)
mean(abs(pc - ec)) # ~= 0.012 -- increase n to decrease this error

## The same for the two smoothed versions
beta <- C.n(u, X, smoothing = "beta")
mean(abs(pc - beta))
check <- C.n(u, X, smoothing = "checkerboard")
mean(abs(pc - check))

## Compare the empirical copula with F.n(pobs())
U <- pobs(X) # pseudo-observations
stopifnot(identical(ec, F.n(u, X = pobs(U)))) # even identical

## Compare the empirical copula based on U at U with the Kendall distribution
## Note: Theoretically, C(U) ~ K, so K(C_n(U, U = U)) should approximately be U(0,1)
plot(ecdf(pK(C.n(U, X), cop = cop@copula, d = d)), asp = 1, xaxs="i", yaxs="i")
segments(0,0, 1,1, col=adjustcolor("blue",1/3), lwd=5, lty = 2)
abline(v=0:1, col="gray70", lty = 2)

## Compare the empirical copula and the true copula on the diagonal
C.n.diag <- function(u) C.n(do.call(cbind, rep(list(u), d)), X = X) # diagonal of C_n
C.diag <- function(u) pCopula(do.call(cbind, rep(list(u), d)), cop) # diagonal of C
curve(C.n.diag, from = 0, to = 1, # empirical copula diagonal
      main = paste("True vs empirical diagonal of a", family, "copula"),
      xlab = "u", ylab = quote("True C(u,..,u) and empirical"~C[n](u,..,u)))
curve(C.diag, lty = 2, add = TRUE) # add true copula diagonal
legend("bottomright", lty = 2:1, bty = "n", inset = 0.02,
       legend = expression(C, C[n]))

## Approximate partial derivatives w.r.t. the 2nd and 3rd component
j.ind <- 2:3 # indices w.r.t. which the partial derivatives are computed
## Partial derivatives based on the empirical copula and the true copula
der23 <- dCn(u, U = pobs(U), j.ind = j.ind)
der23. <- copula:::dCdu(archmCopula(family, param=theta, dim=d), u=u)[,j.ind]
## Approximation error
summary(as.vector(abs(der23-der23.)))

## For an example of using F.n(), see help(mvdc)% ./Mvdc.Rd

## Generate a bivariate empirical copula object (various smoothing methods)
n <- 10 # sample size
d <- 2 # dimension
set.seed(271)
X <- rCopula(n, copula = claytonCopula(3, dim = d))
ecop.orig  <- empCopula(X) # smoothing = "none"
ecop.beta  <- empCopula(X, smoothing = "beta")
ecop.check <- empCopula(X, smoothing = "checkerboard")

## Sample from these (smoothed) empirical copulas
m <- 50
U.orig  <-  rCopula(m, copula = ecop.orig)
U.beta  <-  rCopula(m, copula = ecop.beta)
U.check <-  rCopula(m, copula = ecop.check)

## Plot
wireframe2(ecop.orig,  FUN = pCopula, draw.4.pCoplines = FALSE)
wireframe2(ecop.beta,  FUN = pCopula)
wireframe2(ecop.check, FUN = pCopula)
## Density (only exists when smoothing = "beta")
wireframe2(ecop.beta,  FUN = dCopula)

## Transform a copula sample to empirical margins
set.seed(271)
X <- qexp(rCopula(1000, copula = claytonCopula(2))) # multivariate distribution
U <- rCopula(917, copula = gumbelCopula(2)) # new copula sample
X. <- toEmpMargins(U, x = X) # tranform U to the empirical margins of X
plot(X.) # Gumbel sample with empirical margins of X

Class "empCopula" of Empirical Copulas

Description

Empirical Copula class.

Objects from the Class

Created by calls of the form new("empCopula", ...) or rather typically by empCopula() based on a matrix X of pseudo-observations. Smoothing options are available, see there.

Slots

X:

matrix of pseudo-observations based on which the empirical copula is constructed.

smoothing:

character string determining the smoothing method.

offset:

numeric giving the shift in the normalizing factor for computing the empirical copula.

ties.method:

a string indicating rank()'s ties method for computing the empirical copula.

See Also

The class constructor are empCopula(), also for examples.

Estimation Procedures for (Nested) Archimedean Copulas

Description

A set of ten different estimators, currently for one-parameter Archimedean copulas, of possibly quite high dimensions.

Usage

enacopula(u, cop,
          method = c("mle", "smle", "dmle",
                     "mde.chisq.CvM", "mde.chisq.KS",
                     "mde.gamma.CvM", "mde.gamma.KS",
                     "tau.tau.mean", "tau.theta.mean", "beta"),
          n.MC = if (method == "smle") 10000 else 0,
          interval = initOpt(cop@copula@name),
          xargs = list(), ...)

Arguments

Details

enacopula serves as a wrapper for the different implemented estimators and provides a uniform framework to utilize them. For more information, see the single estimators as given in the section ‘See Also’.

Note that Hofert, Mächler, and McNeil (2013) compared these estimators. Their findings include a rather poor performance and numerically challenging problems of some of these estimators. In particular, the estimators obtained by method="mde.gamma.CvM", method="mde.gamma.KS", method="tau.theta.mean", and method="beta" should be used with care (or not at all). Overall, MLE performed best (by far).

Value

the estimated parameter, \hat{\theta}, that is, currently a number as only one-parameter Archimedean copulas are considered.

References

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.

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

emle which returns an object of "mle" providing useful methods not available for other estimators. demo(opC-demo) and vignette("GIG", package="copula") for examples of two-parameter families. edmle for the diagonal maximum likelihood estimator. emde for the minimum distance estimators. etau for the estimators based on Kendall's tau. ebeta for the estimator based on Blomqvist's beta.

Examples

tau <- 0.25
(theta <- copGumbel@iTau(tau)) # 4/3
d <- 12
(cop <- onacopulaL("Gumbel", list(theta,1:d)))

set.seed(1)
n <- 100
U <- rnacopula(n, cop)

meths <- eval(formals(enacopula)$method)

fun <- function(meth, u, cop, theta) {
	run.time <- system.time(val <- enacopula(u, cop=cop, method=meth))
	list(value=val, error=val-theta, utime.ms=1000*run.time[[1]])
}
t(res <- sapply(meths, fun, u=U, cop=cop, theta=theta))

Various Estimators for (Nested) Archimedean Copulas

Description

Various Estimators for (Nested) Archimedean Copulas, namely,

ebeta

Method-of-moments-like estimator based on (a multivariate version of) Blomqvist'sbeta.

edmle

Maximum likelihood estimator based on the diagonal of a (nested) Archimedean copula.

etau

Method-of-moments-like estimators based on (bivariate) Kendall's tau.

Usage

ebeta(u, cop, interval = initOpt(cop@copula@name), ...)
edmle(u, cop, interval = initOpt(cop@copula@name), warn=TRUE, ...)
 etau(u, cop, method = c("tau.mean", "theta.mean"), warn=TRUE, ...)

Arguments

Details

For ebeta, the parameter is estimated with a method-of-moments-like procedure such that the population version of the multivariate Blomqvist's beta matches its sample version.

Note that the copula diagonal is a distribution function and the maximum of all components of a random vector following the copula is distributed according to this distribution function. For edmle, the parameter is estimated via maximum-likelihood estimation based on the diagonal.

For etau, corKendall(u, ...) is used and if there are no NAs in u, by default (if no additional arguments are provided), corKendall() calls the O(n log(n)) fast cor.fk() from package pcaPP instead of the O(n^2) cor(*, method="kendall"). Conversely, when u has NAs, by default, corKendall(u, ...) will use cor(u, method="kendall", use = "pairwise") such that etau(u, *) will work.
Furthermore, method="tau.mean" means that the average of sample versions of Kendall's tau are computed first and then the parameter is determined such that the population version of Kendall's tau matches this average (if possible); the method="theta.mean" stands for first computing all pairwise Kendall's tau estimators and then returning the mean of these estimators.

For more details, see Hofert et al. (2013).

Note that these estimators should be used with care; see the performance results in Hofert et al. (2013). In particular, etau should be used with the (default) method "tau.mean" since "theta.mean" is both slower and more prone to errors.

Value

ebeta

the return value of safeUroot (that is, typically almost the same as the value of uniroot) giving the Blomqvist beta estimator.

edmle

list as returned by optimize, including the diagonal maximum likelihood estimator.

etau

method-of-moments-like estimator based on Kendall's tau for the chosen method.

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

corKendall().

The more sophisticated estimators emle (Maximum Likelihood) and emde (Minimum Distance). enacopula (wrapper for different estimators).

Examples

tau <- 0.25
(theta <- copGumbel@iTau(tau)) # 4/3 = 1.333..
d <- 20
(cop <- onacopulaL("Gumbel", list(theta,1:d)))

set.seed(1)
n <- 200
U <- rnacopula(n, cop)

system.time(theta.hat.beta <- ebeta(U, cop=cop))
theta.hat.beta$root

system.time(theta.hat.dmle <- edmle(U, cop=cop))
theta.hat.dmle$minimum

system.time(theta.hat.etau <- etau(U, cop=cop, method="tau.mean"))
theta.hat.etau

system.time(theta.hat.etau. <- etau(U, cop=cop, method="theta.mean"))
theta.hat.etau.

## etau()  in the case of missing values (NA's)
## ------                 ---------------------
##' @title add Missing Values completely at random
##' @param x  matrix or vector
##' @param prob desired probability ("fraction") of missing values (\code{\link{NA}}s).
##' @return x[] with some (100*prob percent) entries replaced by \code{\link{NA}}s.
addNAs <- function(x, prob) {
    np <- length(x)
    x[sample.int(np, prob*np)] <- NA
    x
}

## UM[] := U[] with 5% missing values
set.seed(7); UM <- addNAs(U, p = 0.05)
mean(is.na(UM)) # 0.05
## This error if x has NA's was happening for  etau(UM, cop=cop)
## before copula version 0.999-17 (June 2017) :
try(eM <- etau(UM, cop=cop, use = "everything"))
        #  --> Error ... NA/NaN/Inf in foreign function call
## The new default:
eM0 <- etau(UM, cop=cop, use = "pairwise")
eM1 <- etau(UM, cop=cop, use = "complete")
##  use = "complete" is really equivalent to dropping all obs. with with missing values:
stopifnot(all.equal(eM1, etau(na.omit(UM), cop=cop), tol = 1e-15))
## but  use = "pairwise" ---> cor(*, use = "pairwise") is much better:
rbind(etau.U = theta.hat.etau, etau.UM.pairwise = eM0, etau.UM.complete = eM1)

Construction of Extreme-Value Copula Objects

Description

Constructs an extreme-value copula class object with its corresponding parameter.

Usage

evCopula(family, param, dim = 2, ...)
galambosCopula(param)
huslerReissCopula(param)
tawnCopula(param)
tevCopula(param, df = 4, df.fixed = FALSE)

Arguments

Value

An object of class "gumbelCopula", "galambosCopula", "huslerReissCopula", "tawnCopula", or "tevCopula".

Note

The Gumbel copula is both an Archimedean and an extreme-value copula, with principal documentation on gumbelCopula (or archmCopula).

See Also

ellipCopula, archmCopula, gofEVCopula, An.

Examples

## Gumbel is both
stopifnot(identical(   evCopula("gumbel"), gumbelCopula()),
          identical(archmCopula("gumbel"), gumbelCopula()))

## For a given degree of dependence these copulas are strikingly similar :

tau <- 1/3

gumbel.cop      <- gumbelCopula     (iTau(gumbelCopula(),      tau))
galambos.cop    <- galambosCopula   (iTau(galambosCopula(),    tau))
huslerReiss.cop <- huslerReissCopula(iTau(huslerReissCopula(), tau))
tawn.cop        <- tawnCopula       (iTau(tawnCopula(),        tau))
tev.cop         <- tevCopula        (iTau(tevCopula(),         tau))

curve(A(gumbel.cop, x), 0, 1, ylab = "A(<cop>( iTau(<cop>(), tau)), x)",
      main = paste("A(x) for five Extreme Value cop. w/  tau =", format(tau)))
curve(A(galambos.cop, x), lty=2, add=TRUE)
curve(A(huslerReiss.cop, x), lty=3, add=TRUE)
curve(A(tawn.cop, x), lty=4, add=TRUE)
curve(A(tev.cop, x), lty=5, col=2, add=TRUE)# very close to Gumbel

## And look at the differences
curve(A(gumbel.cop, x) - A(tawn.cop, x), ylim = c(-1,1)*0.005,
      ylab = '', main = "A(<Gumbel>, x) - A(<EV-Cop.>, x)")
abline(h=0, lty=2)
curve(A(gumbel.cop, x) - A(galambos.cop, x), add=TRUE, col=2)
curve(A(gumbel.cop, x) - A(huslerReiss.cop, x), add=TRUE, col=3)
curve(A(gumbel.cop, x) - A(tev.cop, x), add=TRUE, col=4, lwd=2)


## the t-EV-copula has always positive tau :
curve(vapply(x, function(x) tau(tevCopula(x)), 0.), -1, 1,
      n=257, ylim=0:1, xlab=quote(rho),ylab=quote(tau),
      main= quote(tau( tevCopula(rho) )), col = 2, lwd = 2)
rect(-1,0,1,1, lty = 2, border = adjustcolor("black", 0.5))

Classes Representing Extreme-Value Copulas

Description

Class evCopula is the virtual (mother) class of all extreme-value copulas. There currently are five subclasses, "galambosCopula", "huslerReissCopula", "tawnCopula", "tevCopula", and "gumbelCopula", the latter of which is also an Archimedean copula, see the page for class "archmCopula".

Objects from the Class

evCopula is a virtual class: No objects may be created from it. Objects of class "galambosCopula" etc, can be created by calls of the form new("galambosCopula", ...), but typically rather by galambosCopula(), etc, see there.

Slots

All slots are inherited from the mother class "copula", see there.

Methods

dCopula

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

pCopula

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

rCopula

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

dCopula

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

pCopula

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

rCopula

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

Extends

Class "evCopula" extends class "copula" directly. Classes "galambosCopula", "huslerReissCopula", "tawnCopula", and "tevCopula" extend class "evCopula" directly.

Note

Objects of class "gumbelCopula" are also of class "archmCopula".

See Also

evCopula, evTestC, evTestK, gofEVCopula, copula-class.

Bivariate Test of Extreme-Value Dependence Based on Pickands' Dependence Function

Description

Test of bivariate extreme-value dependence based on the process comparing the empirical copula with a natural nonparametric estimator of the unknown copula derived under extreme-value dependence. The test statistics are defined in the third reference. Approximate p-values for the test statistics are obtained by means of a multiplier technique.

Usage

evTestA(x, N = 1000, derivatives = c("An","Cn"),
        ties.method = eval(formals(rank)$ties.method),
        trace.lev = 0, report.err = FALSE)

Arguments

Details

More details are available in the third reference. See also Genest and Segers (2009) and Remillard and Scaillet (2009).

Value

An object of class htest which is a list, some of the components of which are

Note

This test was derived under the assumption of continuous margins, which implies that ties occur with probability zero. The presence of ties in the data might substantially affect the approximate p-value.

References

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

Rémillard, B. and Scaillet, O. (2009). Testing for equality between two copulas. Journal of Multivariate Analysis, 100(3), pages 377-386.

Kojadinovic, I. and Yan, J. (2010). Nonparametric rank-based tests of bivariate extreme-value dependence. Journal of Multivariate Analysis 101, 2234–2249.

See Also

evTestK, evTestC, evCopula, gofEVCopula, An.

Examples

## Do these data come from an extreme-value copula?
set.seed(63)
uG <- rCopula(100, gumbelCopula (3))
uC <- rCopula(100, claytonCopula(3))
## these two take 21 sec on nb-mm4 (Intel Core i7-5600U @ 2.60GHz):
evTestA(uG)
evTestA(uC) # significant even though Clayton is *NOT* an extreme value copula

## These are fast:
evTestA(uG, derivatives = "Cn")
evTestA(uC, derivatives = "Cn") # small p-value even though Clayton is *NOT* an EV copula.

Large-sample Test of Multivariate Extreme-Value Dependence

Description

Test of multivariate extreme-value dependence based on the empirical copula and max-stability. The test statistics are defined in the second reference. Approximate p-values for the test statistics are obtained by means of a multiplier technique.

Usage

evTestC(x, N = 1000)

Arguments

Details

More details are available in the second reference. See also Remillard and Scaillet (2009).

Value

An object of class htest which is a list, some of the components of which are

Note

This test was derived under the assumption of continuous margins, which implies that ties occur with probability zero. The presence of ties in the data might substantially affect the approximate p-value.

References

Rémillard, B. and Scaillet, O. (2009). Testing for equality between two copulas. Journal of Multivariate Analysis, 100(3), pages 377-386.

Kojadinovic, I., Segers, J., and Yan, J. (2011). Large-sample tests of extreme-value dependence for multivariate copulas. The Canadian Journal of Statistics 39, 4, pages 703-720.

See Also

evTestK, evTestA, evCopula, gofEVCopula, An.

Examples

## Do these data come from an extreme-value copula?
evTestC(rCopula(200, gumbelCopula(3)))
evTestC(rCopula(200, claytonCopula(3)))

## Three-dimensional examples
evTestC(rCopula(200, gumbelCopula(3, dim=3)))
evTestC(rCopula(200, claytonCopula(3, dim=3)))

Bivariate Test of Extreme-Value Dependence Based on Kendall's Distribution

Description

Test of extreme-value dependence based on the bivariate probability integral transformation. The test statistic is defined in Ben Ghorbal, G. Nešlehová, and Genest (2009).

Usage

evTestK(x, method = c("fsample","asymptotic","jackknife"), ties = NA, N = 100)

Arguments

Details

The code for this test was generously provided by Johanna G. Nešlehová. More details are available in Appendix B of Ben Ghorbal, G. Nešlehová and Genest (2009).

Value

An object of class htest which is a list, some of the components of which are

References

Ghorbal, M. B., Genest, C., and G. Nešlehová, J. (2009) On the test of Ghoudi, Khoudraji, and Rivest for extreme-value dependence. The Canadian Journal of Statistics 37, 1–9.

Kojadinovic, I. (2017). Some copula inference procedures adapted to the presence of ties. Computational Statistics and Data Analysis 112, 24–41, https://arxiv.org/abs/1609.05519.

See Also

evTestC, evTestA, evCopula, gofEVCopula, An.

Examples

set.seed(321)
## Do the data come from an extreme-value copula?
evTestK(Ug <- rCopula(200, gumbelCopula(3))) # not significant => yes, EV
    dim(Uc <- rCopula(200, claytonCopula(3)))
## Clayton:                       tests are highly significant => no, not EV
(K1 <- evTestK(Uc))
(K2 <- evTestK(Uc, method = "asymptotic"))

system.time(print(K3 <- evTestK(Uc, method = "jackknife")))
## slower: 1.06 sec (2015 intel core i7)

Test of Exchangeability for Certain Bivariate Copulas

Description

Test for assessing the exchangeability of the underlying bivariate copula when it is either extreme-value or left-tail decreasing. The test uses the nonparametric estimators of the Pickands dependence function studied by Genest and Segers (2009).

The test statistic is defined in the second reference. An approximate p-value for the test statistic is obtained by means of a multiplier technique if there are no ties in the component series of the bivariate data, or by means of an appropriate bootstrap otherwise.

Usage

exchEVTest(x, N = 1000,  estimator = c("CFG", "Pickands"),
           ties = NA, ties.method = eval(formals(rank)$ties.method),
           m = 100, derivatives = c("Cn", "An"))

Arguments

Details

More details are available in the references.

Value

An object of class htest which is a list, some of the components of which are

References

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

Kojadinovic, I. and Yan, J. (2012) A nonparametric test of exchangeability for extreme-value and left-tail decreasing bivariate copulas. The Scandinavian Journal of Statistics 39:3, 480–496.

Kojadinovic, I. (2017). Some copula inference procedures adapted to the presence of ties. Computational Statistics and Data Analysis 112, 24–41, https://arxiv.org/abs/1609.05519.

See Also

exchTest, radSymTest, gofCopula.

Examples

## Data from an exchangeable left-tail decreasing copulas
exchEVTest(rCopula(200,  gumbelCopula(3)))
exchEVTest(rCopula(200, claytonCopula(3)))

## An asymmetric Khoudraji-Gumbel-Hougaard copula
kc <- khoudrajiCopula(copula1 = indepCopula(),
                      copula2 = gumbelCopula(4),
                      shapes = c(0.4, 0.95))
exchEVTest(rCopula(200, kc))

Test of Exchangeability for a Bivariate Copula

Description

Test for assessing the exchangeability of the underlying bivariate copula based on the empirical copula. The test statistics are defined in the first two references. Approximate p-values for the test statistics are obtained by means of a multiplier technique if there are no ties in the component series of the bivariate data, or by means of an appropriate bootstrap otherwise.

Usage

exchTest(x, N = 1000, ties = NA,
         ties.method = eval(formals(rank)$ties.method), m = 0)

Arguments

Details

More details are available in the references.

Value

An object of class htest which is a list, some of the components of which are

References

Genest, C., G. Nešlehová, J. and Quessy, J.-F. (2012). Tests of symmetry for bivariate copulas. Annals of the Institute of Statistical Mathematics 64, 811–834.

Kojadinovic, I. and Yan, J. (2012). A nonparametric test of exchangeability for extreme-value and left-tail decreasing bivariate copulas. The Scandinavian Journal of Statistics 39:3, 480–496.

Kojadinovic, I. (2017). Some copula inference procedures adapted to the presence of ties. Computational Statistics and Data Analysis 112, 24–41, https://arxiv.org/abs/1609.05519.

See Also

radSymTest, exchEVTest, gofCopula.

Examples

## Data from an exchangeable copulas
exchTest(rCopula(200,  gumbelCopula(3)))
exchTest(rCopula(200, claytonCopula(3)))

## An asymmetric Khoudraji-Clayton copula
kc <- khoudrajiCopula(copula1 = indepCopula(),
                      copula2 = claytonCopula(6),
                      shapes = c(0.4, 0.95))
exchTest(rCopula(200, kc))

Construction of a fgmCopula Class Object

Description

Constructs a multivariate multiparameter Farlie-Gumbel-Morgenstern copula class object with its corresponding parameters and dimension.

Usage

fgmCopula(param, dim = 2)

Arguments

Value

A Farlie-Gumbel-Morgenstern copula object of class "fgmCopula".

Note

Note that a d-dimensional fgmCopula must have npar= 2^d - d - 1 parameters. The verification of the validity of the parameter values is of high complexity and may not work for high dimensional copulas.

The random number generation needs to be properly tested, especially for dimensions higher than 2.

References

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

See Also

Copula, copula-class, fitCopula.

Examples

## length(param) = #{parameters}  for  d-dimensional  FGM copula:
d <- 2:10; rbind(d, npar = 2^d - d - 1)
## d       2    3    4    5    6    7    8    9   10
## npar    1    4   11   26   57  120  247  502 1013

## a bivariate example
fgm.cop <- fgmCopula(1)
x <- rCopula(1000, fgm.cop)
cor(x, method = "kendall")
tau(fgm.cop)
cor(x, method = "spearman")
rho(fgm.cop)
persp  (fgm.cop, dCopula)
contour(fgm.cop, dCopula)

## a trivariate example with wrong parameter values
try(
 fgm2.cop <- fgmCopula(c(1,1,1,1), dim = 3)
) # Error: "Bad vector of parameters"

## a trivariate example with satisfactory parameter values
fgm2.cop <- fgmCopula(c(.2,-.2,-.4,.6), dim = 3)
fgm2.cop

Class "fgmCopula" - Multivariate Multiparameter Farlie-Gumbel-Morgenstern Copulas

Description

The class of multivariate multiparameter Farlie-Gumbel-Morgenstern copulas are typically created via fgmCopula(..).

Objects from the Class

Objects are typically created by fgmCopula(..), or more low-level by (careful) calls to new("fgmCopula", ..).

Slots

exprdist:

Object of class "expression", expressions for the cdf and pdf of the copula. These expressions are used in function pCopula() and dCopula().

subsets.char:

Object of class "character", containing the subsets of integers used for naming the parameters.

dimension:

Object of class "numeric", the dimension of the copula.

parameters:

Object of class "numeric", parameter values.

param.names:

Object of class "character", parameter names.

param.lowbnd:

Object of class "numeric", parameter lower bound.

param.upbnd:

Object of class "numeric", parameter upper bound.

fullname:

Object of class "character", family names of the copula (deprecated).

Methods

dCopula

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

pCopula

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

rCopula

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

Extends

Class "fgmCopula" extends class "copula" directly.

Note

The verification of the validity of the parameter values is of high complexity and may not work for high dimensional copulas.

The random number generation needs to be properly tested, especially for dimensions higher than 2.

References

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

See Also

copula-class; to create such objects, use fgmCopula(); see there, also for examples.

Construction of Fréchet-Hoeffding Bound Copula Objects

Description

Constructs the Fréchet-Hoeffding lower and upper bound copulas aka W and M.

Usage

fhCopula(family = c("upper", "lower"), dim = 2L)

lowfhCopula(dim = 2L)
 upfhCopula(dim = 2L)

Arguments

Value

A copula object of class "lowfhCopula" or "upfhCopula".

Note

fhCopula() is a wrapper for lowfhCopula() and upfhCopula().

The dCopula() method will simply return an error for these copulas (as their density does not exist). Also, since the Fréchet-Hoeffding bound copulas are not parametric, certain methods available for parametric copulas are not available.

Examples

## Lower bound  W : -------------------------

try(W <- lowfhCopula(dim = 3)) # lower bound is *not* a copula for dim > 2
W <- lowfhCopula()
wireframe2(W, FUN = pCopula)
plot(W, n=100) # perfect anti-correlation ( rho = tau = -1 )

## Upper bound  M : -------------------------

wireframe2(upfhCopula(dim = 2), pCopula)
M <- upfhCopula(dim = 3)
set.seed(271)
splom2(M, n = 100) # "random" data: all perfectly correlated

Class "fhCopula" of Fréchet-Hoeffding Bound Copulas

Description

Fréchet-Hoeffding bound copula class.

Objects from the Class

Created by calls of the form new("fhCopula", ...) or rather typically by fhCopula(), lowfhCopula(), or upfhCopula(). Actual (sub) classes are the lower and upper Fréchet-Hoeffding bound copulas lowfhCopula() (W), and upfhCopula() (M).

Slots

dimension:

inherited from super class "dimCopula".

exprdist:

an expression of length two, named "cdf" with the R expression of the CDF, and "pdf" which is empty as the PDF does not exist (everywhere).

See Also

ellipCopula, archmCopula, evCopula.

The class constructors are fhCopula(), lowfhCopula(), and upfhCopula(). See there for examples.

Fitting Copulas to Data – Copula Parameter Estimation

Description

Parameter estimation of copulas, i.e., fitting of a copula model to multivariate (possibly “pseudo”) observations.

Usage

loglikCopula(param = getTheta(copula), u, copula,
             error = c("-Inf", "warn-Inf", "let-it-be"))

loglikCopulaMany(pList, u, copula)

## Generic [and "parCopula" method] : %- ../R/fitCopula.R
fitCopula(copula, data, ...)
## S4 method for signature 'parCopula'
fitCopula(copula, data,
          method = c("mpl", "ml", "itau", "irho", "itau.mpl"),
          posDef = is(copula, "ellipCopula"),
          start = NULL, lower = NULL, upper = NULL,
          optim.method = optimMeth(copula, method, dim = d),
          optim.control = list(maxit=1000),
          estimate.variance = NA, hideWarnings = FALSE, ...)

optimMeth(copula, method, dim)

Arguments

Details

The only difference between "mpl" and "ml" is in the variance-covariance estimate, not in the parameter (\theta) estimates.

If method "mpl" in fitCopula() is used and if start is not assigned a value, estimates obtained from method "itau" are used as initial values in the optimization. Standard errors are computed as explained in Genest, Ghoudi and Rivest (1995); see also Kojadinovic and Yan (2010, Section 3). Their estimation requires the computation of certain partial derivatives of the (log) density. These have been implemented for six copula families thus far: the Clayton, Gumbel-Hougaard, Frank, Plackett, normal and t copula families. For other families, numerical differentiation based on grad() from package numDeriv is used (and a warning message is displayed).

In the multiparameter elliptical case and when the estimation is based on Kendall's tau or Spearman's rho, the estimated correlation matrix may not always be positive-definite. In that case, nearPD(*, corr=TRUE) (from Matrix) is applied to get a proper correlation matrix.

For normal and t copulas, fitCopula(, method = "mpl") and fitCopula(, method = "ml") maximize the log-likelihood based on mvtnorm's dmvnorm() and dmvt(), respectively. The latter two functions set the respective densities to zero if the correlation matrices of the corresponding distributions are not positive definite. As such, the estimated correlation matrices will be positive definite.

If methods "itau" or "irho" are used in fitCopula(), an estimate of the asymptotic variance (if available for the copula under consideration) will be correctly computed only if the argument data consists of pseudo-observations (see pobs()).

Consider the t copula with df.fixed=FALSE (see ellipCopula()). In this case, the methods "itau" and "irho" cannot be used in fitCopula() as they cannot estimate the degrees of freedom parameter df. For the methods "mpl" and "itau.mpl" the asymptotic variance cannot be (fully) estimated (yet). For the methods "ml" and "mpl", when start is not specified, the starting value for df is set to copula@df, typically 4.

To implement the Inference Functions for Margins (IFM) method (see, e.g., Joe 2005), set method="ml" and note that data need to be parametric pseudo-observations obtained from fitted parametric marginal distribution functions. The returned large-sample variance will then underestimate the true variance (as the procedure cannot take into account the (unknown) estimation error for the margins).

The fitting procedures based on optim() generate warnings because invalid parameter values are tried during the optimization process. When the number of parameters is one and the parameter space is bounded, using optim.method="Brent" is likely to give less warnings. Furthermore, from experience, optim.method="Nelder-Mead" is sometimes a more robust alternative to optim.method="BFGS" or "L-BFGS-B".

There are methods for vcov(), coef(), logLik(), and nobs().

Value

loglikCopula() returns the copula log-likelihood evaluated at the parameter (vector) param given the data u.

loglikCopulaMany() returns a numeric vector of such log-likelihoods; it assumes consistent parameter values, corresponding to loglikCopula()'s error = "let-it-be", for speed.

The return value of fitCopula() is an object of class "fitCopula" (inheriting from hidden class "fittedMV"), containing (among others!) the slots

estimate

The parameter estimates.

var.est

The large-sample (i.e., asymptotic) variance estimate of the parameter estimator unless estimate.variance=FALSE where it is matrix(numeric(), 0,0) (to be distinguishable from cases when the covariance estimates failed partially).

copula

The fitted copula object.

The summary() method for "fitCopula" objects returns an S3 “class” "summary.fitCopula", which is simply a list with components method, loglik and convergence, all three from the corresponding slots of the "fitCopula" objects, and coefficients (a matrix of estimated coefficients, standard errors, t values and p-values).

References

Genest, C. (1987). Frank's family of bivariate distributions. Biometrika 74, 549–555.

Genest, C. and Rivest, L.-P. (1993). Statistical inference procedures for bivariate Archimedean copulas. Journal of the American Statistical Association 88, 1034–1043.

Rousseeuw, P. and Molenberghs, G. (1993). Transformation of nonpositive semidefinite correlation matrices. Communications in Statistics: Theory and Methods 22, 965–984.

Genest, C., Ghoudi, K., and Rivest, L.-P. (1995). A semiparametric estimation procedure of dependence parameters in multivariate families of distributions. Biometrika 82, 543–552.

Joe, H. (2005). Asymptotic efficiency of the two-stage estimation method for copula-based models. Journal of Multivariate Analysis 94, 401–419.

Zeevi, Assaf and Mashal, Roy (2002) Beyond Correlation: Extreme Co-Movements between Financial Assets. doi:10.2139/ssrn.317122

Demarta, S. and McNeil, A. J. (2005). The t copula and related copulas. International Statistical Review 73, 111–129.

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.

Kojadinovic, I. and Yan, J. (2010). Comparison of three semiparametric methods for estimating dependence parameters in copula models. Insurance: Mathematics and Economics 47, 52–63.

See Also

Copula, fitMvdc for fitting multivariate distributions including the margins, gofCopula for goodness-of-fit tests.

For maximum likelihood of (nested) Archimedean copulas, see emle, etc.

Examples

(Xtras <- copula:::doExtras()) # determine whether examples will be extra (long)
n <- if(Xtras) 200 else 64 # sample size

## A Gumbel copula
set.seed(7) # for reproducibility
gumbel.cop <- gumbelCopula(3, dim=2)
x <- rCopula(n, gumbel.cop) # "true" observations (simulated)
u <- pobs(x)                # pseudo-observations
## Inverting Kendall's tau
fit.tau <- fitCopula(gumbelCopula(), u, method="itau")
fit.tau
confint(fit.tau) # works fine !
confint(fit.tau, level = 0.98)
summary(fit.tau) # a bit more, notably "Std. Error"s
coef(fit.tau)# named vector
coef(fit.tau, SE = TRUE)# matrix

## Inverting Spearman's rho
fit.rho <- fitCopula(gumbelCopula(), u, method="irho")
summary(fit.rho)
## Maximum pseudo-likelihood
fit.mpl <- fitCopula(gumbelCopula(), u, method="mpl")
fit.mpl
## Maximum likelihood -- use 'x', not 'u' ! --
fit.ml <- fitCopula(gumbelCopula(), x, method="ml")
summary(fit.ml) # now prints a bit more than simple 'fit.ml'
## ... and what's the log likelihood (in two different ways):
(ll. <- logLik(fit.ml))
stopifnot(all.equal(as.numeric(ll.),
            loglikCopula(coef(fit.ml), u=x, copula=gumbel.cop)))

## A Gauss/normal copula

## With multiple/*un*constrained parameters
set.seed(6) # for reproducibility
normal.cop <- normalCopula(c(0.6, 0.36, 0.6), dim=3, dispstr="un")
x <- rCopula(n, normal.cop) # "true" observations (simulated)
u <- pobs(x)                # pseudo-observations
## Inverting Kendall's tau
fit.tau <- fitCopula(normalCopula(dim=3, dispstr="un"), u, method="itau")
fit.tau

if(Xtras) withAutoprint({ ## needs time
  ## Inverting Spearman's rho
  fit.rho <- fitCopula(normalCopula(dim=3, dispstr="un"), u, method="irho")
  fit.rho
  ## Maximum pseudo-likelihood
  fit.mpl <- fitCopula(normalCopula(dim=3, dispstr="un"), u, method="mpl")
  summary(fit.mpl)
  coef(fit.mpl) # named vector
  coef(fit.mpl, SE = TRUE) # the matrix, with SE
}) #
## Maximum likelihood (use 'x', not 'u' !)
fit.ml <- fitCopula(normalCopula(dim=3, dispstr="un"), x, method="ml", traceOpt=TRUE)
summary(fit.ml)
confint(fit.ml)
confint(fit.ml, level = 0.999) # clearly non-0

## Fix some of the parameters
param <- c(.6, .3, NA_real_)
attr(param, "fixed") <- c(TRUE, FALSE, FALSE)
ncp <- normalCopula(param = param, dim = 3, dispstr = "un")
fixedParam(ncp) <- c(TRUE, TRUE, FALSE)
## 'traceOpt = 5': showing every 5-th log likelihood evaluation:
summary(Fxf.mpl <- fitCopula(ncp, u, method = "mpl", traceOpt = 5))
Fxf.mpl@copula # reminding of the fixed param. values

## With dispstr = "toep" :
normal.cop.toep <- normalCopula(c(0, 0), dim=3, dispstr="toep")
## Inverting Kendall's tau
fit.tau <- fitCopula(normalCopula(dim=3, dispstr="toep"), u, method="itau")
fit.tau
## Inverting Spearman's rho
fit.rho <- fitCopula(normalCopula(dim=3, dispstr="toep"), u, method="irho")
summary(fit.rho)
## Maximum pseudo-likelihood
fit.mpl <- fitCopula(normalCopula(dim=3, dispstr="toep"), u, method="mpl")
fit.mpl
## Maximum likelihood (use 'x', not 'u' !)
fit.ml <- fitCopula(normalCopula(dim=3, dispstr="toep"), x, method="ml")
summary(fit.ml)

## With dispstr = "ar1"
normal.cop.ar1 <- normalCopula(c(0), dim=3, dispstr="ar1")
## Inverting Kendall's tau
summary(fit.tau <- fitCopula(normalCopula(dim=3, dispstr="ar1"), u, method="itau"))
## Inverting Spearman's rho
summary(fit.rho <- fitCopula(normalCopula(dim=3, dispstr="ar1"), u, method="irho"))
## Maximum pseudo-likelihood
summary(fit.mpl <- fitCopula(normalCopula(dim=3, dispstr="ar1"), u, method="mpl"))
## Maximum likelihood (use 'x', not 'u' !)
fit.ml <- fitCopula(normalCopula(dim=3, dispstr="ar1"), x, method="ml")
summary(fit.ml)

## A t copula with variable df (df.fixed=FALSE)
(tCop <- tCopula(c(0.2,0.4,0.6), dim=3, dispstr="un", df=5))
set.seed(101)
x <- rCopula(n, tCop) # "true" observations (simulated)
## Maximum likelihood (start = (rho[1:3], df))
summary(tc.ml <- fitCopula(tCopula(dim=3, dispstr="un"), x, method="ml",
                           start = c(0,0,0, 10)))
## Maximum pseudo-likelihood (the asymptotic variance cannot be estimated)
u <- pobs(x)          # pseudo-observations
tc.mpl <- fitCopula(tCopula(dim=3, dispstr="un"),
                     u, method="mpl", estimate.variance=FALSE,
                     start= c(0,0,0, 10))
summary(tc.mpl)

Classes of Fitted Multivariate Models: Copula, Mvdc

Description

Classes and summary methods related to copula model fitting.

Objects from the Class

Objects can be created by calls to fitCopula or fitMvdc, respectively or to their summary methods.

Slots

The “mother class”, "fittedMV" has the slots

estimate:

numeric, the estimated parameters.

var.est:

numeric, variance matrix estimate of the parameter estimator. See note below.

loglik:

numeric, log likelihood evaluated at the maximizer.

nsample:

numeric, integer representing the sample size.

method:

character, method of estimation.

fitting.stats:

a list, currently containing the numeric convergence code from optim, the counts, message, and all the control arguments explicitly passed to optim(). Since copula version 1.0-1 also keeps information about parameter transformations, currently needed only for mixCopula fits with free weights.

In addition, the "fitCopula" class has a slot

copula:

the fitted copula, of class "copula".

whereas the "fitMvdc" has

mvdc:

the fitted distribution, of class "mvdc".

Extends

Classes "fitCopula" and "fitMvdc" extend class "fittedMV", directly.

Methods

summary

signature(object = "fitMvdc"): ...

summary

signature(object = "fitCopula"): ...

Further, there are S3 methods (class "fittedMV") for coef(), vcov() and logLik(), see fitMvdc.

References

Genest, C., Ghoudi, K., and Rivest, L.-P. (1995). A semiparametric estimation procedure of dependence parameters in multivariate families of distributions. Biometrika 82, 543–552.

Non-parametric Estimators of the Matrix of Tail-Dependence Coefficients

Description

Computing non-parametric estimators of the (matrix of) tail-dependence coefficients.

Usage

fitLambda(u, method = c("Schmid.Schmidt", "Schmidt.Stadtmueller", "t"),
          p = 1/sqrt(nrow(u)), lower.tail = TRUE, verbose = FALSE, ...)

Arguments

Details

As seen in the examples, be careful using nonparametric estimators, they might not perform too well (depending on p and in general). After all, the notion of tail dependence is a limit, finite sample sizes may not be able to capture limits well.

Value

The matrix of pairwise coefficients of tail dependence; for method = "t" a list with components Lambda, the matrix of pairwise estimated correlation coefficients P and the matrix of pairwise estimated degrees of freedom Nu.

References

Jaworski, P., Durante, F., Härdle, W. K., Rychlik, T. (2010). Copula Theory and Its Applications Springer, Lecture Notes in Statistics – Proceedings.

Schmid, F., Schmidt, R. (2007). Multivariate conditional versions of Spearman's rho and related measures of tail dependence. Journal of Multivariate Analysis 98, 1123–1140. doi:10.1016/j.jmva.2006.05.005

Schmidt, R., Stadtmueller, U. (2006). Non-parametric Estimation of Tail Dependence. Scandinavian Journal of Statistics 33(2), 307–335. doi:10.1111/j.1467-9469.2005.00483.x

See Also

lambda() computes the true (lower and upper) tail coefficients for a given copula.

Examples

n <- 10000 # sample size
p <- 0.01 # cut-off

## Bivariate case
d <- 2
cop <- claytonCopula(2, dim = d)
set.seed(271)
U <- rCopula(n, copula = cop) # generate observations (unrealistic)
(lam.true <- lambda(cop)) # true tail-dependence coefficients lambda
(lam.C <- c(lower = fitLambda(U, p = p)[2,1],
            upper = fitLambda(U, p = p, lower.tail = FALSE)[2,1])) # estimate lambdas
## => pretty good
U. <- pobs(U) # pseudo-observations (realistic)
(lam.C. <- c(lower = fitLambda(U., p = p)[2,1],
             upper = fitLambda(U., p = p, lower.tail = FALSE)[2,1])) # estimate lambdas
## => The pseudo-observations do have an effect...

## Higher-dimensional case
d <- 5
cop <- claytonCopula(2, dim = d)
set.seed(271)
U <- rCopula(n, copula = cop) # generate observations (unrealistic)
(lam.true <- lambda(cop)) # true tail-dependence coefficients lambda
(Lam.C <- list(lower = fitLambda(U, p = p),
               upper = fitLambda(U, p = p, lower.tail = FALSE))) # estimate Lambdas
## => Not too good
U. <- pobs(U) # pseudo-observations (realistic)
(Lam.C. <- list(lower = fitLambda(U., p = p),
                upper = fitLambda(U., p = p, lower.tail = FALSE))) # estimate Lambdas
## => Performance not too great here in either case

Estimation of Multivariate Models Defined via Copulas

Description

Fitting copula-based multivariate distributions ("mvdc") to multivariate data, estimating both the marginal and the copula parameters.

If you assume (non parametric) margins, in other words, take the empirical distributions for all margins, you can use fitCopula(*, pobs(x)) instead.

Usage

loglikMvdc(param, x, mvdc)
fitMvdc(data, mvdc, start, optim.control = list(), method = "BFGS",
        lower = -Inf, upper = Inf,
        estimate.variance = fit$convergence == 0, hideWarnings = TRUE)

## S3 method for class 'fittedMV'
coef(object, SE = FALSE, orig = TRUE, ...)
## S3 method for class 'fittedMV'
logLik(object, ...)
## S3 method for class 'fittedMV'
vcov(object, orig = TRUE, ...)

Arguments

Value

The return value loglikMvdc() is the log likelihood evaluated for the given value of param.

The return value of fitMvdc() is an object of class "fitMvdc" (see there), containing slots (among others!):

The summary() method for "fitMvdc" objects returns a S3 “class” "summary.fitMvdc", simply a list with components method, loglik, and convergence, all three from corresponding slots of the "fitMvdc" objects, and

Note

User-defined marginal distributions can be used as long as the "{dpq}" functions are defined. See vignette("AR_Clayton", package="copula").

When covariates are available for marginal distributions or for the copula, one can construct the loglikelihood function and feed it to "optim" to estimate all the parameters.

Finally, note that some of the fitting functions generate error messages because invalid parameter values are tried during the optimization process (see optim). This should be rarer since 2013, notably for likelihood based methods (as the likelihood is now rather set to -Inf than giving an error).

Previously, loglikMvdc() had an argument hideWarnings; nowadays, do use suppressWarnings(..) if you are sure you do not want to see them.

See Also

mvdc and mvdc; further, Copula, fitCopula, gofCopula.

For fitting univariate marginals, fitdistr().

Examples

G3 <- gumbelCopula(3, dim=2)
gMvd2  <- mvdc(G3, c("exp","exp"),
               param = list(list(rate=2), list(rate=4)))
## with identical margins:
gMvd.I <- mvdc(G3, "exp",
               param = list(rate=3), marginsIdentical=TRUE)

(Xtras <- copula:::doExtras()) # determine whether examples will be extra (long)
n <- if(Xtras) 10000 else 200 # sample size (realistic vs short for example)

set.seed(11)
x <- rMvdc(n, gMvd2)
## Default:     hideWarnings = FALSE .. i.e. show warnings here
fit2 <- fitMvdc(x, gMvd2, start = c(1,1, 2))
fit2
confint(fit2)
summary(fit2) # slightly more
## The estimated, asymptotic var-cov matrix [was used for confint()]:
vcov(fit2)

## with even more output for the "identical margin" case:
fitI <- fitMvdc(x, gMvd.I, start = c(3, 2),
                optim.control=list(trace= TRUE, REPORT= 2))
summary(fitI)
coef(fitI, SE = TRUE)
stopifnot(is.matrix(coef(fitI, SE = TRUE)),
          is.matrix(print( confint(fitI) )) )

## a wrong starting value can already be *the* problem:
f2 <- try(fitMvdc(x, gMvd.I, start = c(1, 1),
           optim.control=list(trace= TRUE, REPORT= 2)))
##--> Error in optim( ... ) : non-finite finite-difference value [2]

##==> "Solution" :  Using a more robust (but slower) optim() method:
fI.2 <- fitMvdc(x, gMvd.I, start = c(1, 1), method = "Nelder",
                optim.control=list(trace= TRUE))
fI.2

Fix a Subset of a Copula Parameter Vector

Description

It is sometimes useful to keep fixed some components of a copula parameter vector whereas the others are “free” and will be estimated, e.g., by fitCopula.

The first two functions set or modify the “fixedness”, whereas isFree(), isFreeP() and nParam() are utilities enquiring about the “fixedness” of the parameters (of a copula).

Usage

fixParam(param, fixed = TRUE)
fixedParam(copula) <- value

isFreeP(param)
## S4 method for signature 'copula'
isFree(copula)
## and specific '*Copula' methods
## S4 method for signature 'copula'
nParam(copula, freeOnly = FALSE)
## and specific '*Copula' methods

Arguments

Value

fixParam(param) returns a numeric vector with attribute "fixed"(a logical, either TRUE or vector of the same length as param) to indicate which components of param are to be held fixed or not.

fixedParam<-, a generic function, returns a "copula" object with a partly fixed parameter (slot), i.e., corresponding to fixParam() above.

See Also

fitCopula for fitting; t-copulas, tCopula(*, df.fixed=TRUE) now uses parameter fixing for "df".

setTheta() for setting or changing the non-fixed parameter values.

Examples

nc1 <- normalCopula(dim = 3, fixParam(c(.6,.3,.2), c(TRUE, FALSE,FALSE)),
                    dispstr = "un")
nc1
nc13 <- nc12 <- nc1
fixedParam(nc12) <- c(TRUE, TRUE, FALSE) ; nc12
fixedParam(nc13) <- c(TRUE, FALSE, TRUE) ; nc13
set.seed(17); x <- rCopula(100, nc1)
summary(f.13 <- fitCopula(nc13, x, method = "ml"))
f.13@copula  # 'rho.2' is estimated; the others kept fixed

## Setting to 'FALSE' (removes the "fixed" parts completely):
nc0 <- nc13; fixedParam(nc0) <- FALSE
nc0
stopifnot(is.null(attributes(nc0@parameters)))

Daily Crude Oil and Natural Gas Prices from 2003 to 2006

Description

Three years of daily prices (from July 2003 to July 2006) of crude oil and natural gas. These data should be very close to those analysed in Grégoire, Genest and Gendron (2008).

Usage

data(gasoil, package="copula")

Format

A data frame of 762 daily prices from 2003 to 2006.

date

date (of class Date).

oil

daily price of crude oil

gas

daily price of natural gas

References

Grégoire, V., Genest, C., and Gendron, M. (2008) Using copulas to model price dependence in energy markets. Energy Risk 5(5), 58–64.

Examples

data(gasoil)
## Log Scaled  Oil & Gas Prices :
lattice :: xyplot(oil + gas ~ date, data = gasoil, auto.key=TRUE,
                  type = c("l","r"),
                  scales=list(y = list(log = TRUE), equispaced.log = FALSE))

Generator Functions for Archimedean and Extreme-Value Copulas

Description

Methods to evaluate the generator function, the inverse generator function, and derivatives of the inverse of the generator function for Archimedean copulas. For extreme-value copulas, the “Pickands dependence function” plays the role of a generator function.

Usage

psi(copula, s)
iPsi(copula, u, ...)
diPsi(copula, u, degree=1, log=FALSE, ...)

A(copula, w)
dAdu(copula, w)

Arguments

Details

psi() and iPsi() are, respectively, the generator function \psi() and its inverse \psi^{(-1)} for an Archimedean copula, see pnacopula for definition and more details.

diPsi() computes (currently only the first two) derivatives of iPsi() (= \psi^{(-1)}).

A(), the “Pickands dependence function”, can be seen as the generator function of an extreme-value copula. For instance, in the bivariate case, we have the following result (see, e.g., Gudendorf and Segers 2009):

A bivariate copula C is an extreme-value copula if and only if

C(u,v) = (uv)^{A(\log(v) / \log(uv))}, \qquad (u,v) \in (0,1]^2 \setminus \{(1,1)\},

where A: [0, 1] \to [1/2, 1] is convex and satisfies \max(t, 1-t) \le A(t) \le 1 for all t \in [0, 1].
In the d-variate case, there is a similar characterization, except that this time, the Pickands dependence function A is defined on the d-dimensional unit simplex.

dAdu() returns a data.frame containing the 1st and 2nd derivative of A().

References

Gudendorf, G. and Segers, J. (2010). Extreme-value copulas. In Copula theory and its applications, Jaworski, P., Durante, F., Härdle, W. and Rychlik, W., Eds. Springer-Verlag, Lecture Notes in Statistics, 127–146, doi:10.1007/978-3-642-12465-5; preprint at https://arxiv.org/abs/0911.1015.

See Also

Nonparametric estimators for A() are available, see An.

Examples

## List the available methods (and their definitions):
showMethods("A")
showMethods("iPsi", incl=TRUE)

Get "acopula" Family Object by Name

Description

Get one of our "acopula" family objects (see acopula-families by name.

Named strings for “translation” between different names and forms of Archimedean copulas.

Usage

getAcop (family, check = TRUE)
getAname(family, objName = FALSE)

.ac.shortNames
.ac.longNames
.ac.objNames
.ac.classNames

Arguments

Value

getAcop() returns an "acopula" family object, typically one of one of our predefined ones.

getAname() returns a character string, the name of an "acopula" family object.

.as.longnames etc are named string constants, useful in programming for all our (five) standard Archimedean families.

See Also

Our predefined acopula-families; the class definition "acopula".

Examples

getAcop("Gumbel")

## different ways of getting the same "acopula" family object:
stopifnot(## Joe (three ways):
          identical(getAcop("J"), getAcop("Joe")),
          identical(getAcop("J"), copJoe),
          ## Frank (yet another two different ways):
          identical(getAcop(frankCopula()), copFrank),
          identical(getAcop("frankCopula"), copFrank))

stopifnot(
  identical(getAname(claytonCopula()), getAname("C")),
  identical(getAname(copClayton), "Clayton"), identical(getAname("J"), "Joe"),
  identical(getAname(amhCopula(), TRUE), "copAMH"),
  identical(getAname(joeCopula(), TRUE), "copJoe")
)

.ac.shortNames
.ac.longNames
.ac.objNames
.ac.classNames

Get Initial Parameter Estimate for Copula

Description

A (S4) generic function and methods providing a typically cheap method to get valid parameters for the copula, given the data. This is used, e.g., in fitCopula() when start is not specified.

Usage

getIniParam(copula, data, default, named = TRUE, ...)

Arguments

Value

a numeric vector of correct length, say param, which should e.g., “work” in loglikCopula(param, u = data, copula).

Methods

signature(copula = "parCopula")

Close to a default method (as class "parCopula" contains most copulas), currently mainly trying to use a version of fitCopula(*, method = "itau") (itself based on moment matching iTau()).

signature(copula = "mixCopula")

a relatively simple method, for the copula parameters, trying getInitParam(cop[[k]]) for each component, and using equal weights w[k].

See Also

getTheta() gets such a vector from a copula object; fitCopula, loglikCopula.

Examples

 # TODO !

Get the Parameter(s) of a Copula

Description

Get the parameter (vector) \theta (theta) of a copula, see setTheta for more background.

Usage

getTheta(copula, freeOnly = TRUE, attr = FALSE, named = attr)

Arguments

Value

parameter vector of the copula, a numeric vector, possibly with names and other attributes (depending on the attr and named arguments).

Methods

signature(copula = "parCopula")

a default method, returning numeric(0).

signature(copula = "copula")

..

signature(copula = "acopula")

..

signature(copula = "khoudrajiCopula")

..

signature(copula = "mixCopula")

..

signature(copula = "rotCopula")

..

signature(copula = "Xcopula")

..

See Also

setTheta, its inverse.

Examples

showMethods("getTheta")

getTheta(setTheta(copClayton, 0.5)) # is 0.5

Computations for Graphical GOF Test via Pairwise Rosenblatt Transforms

Description

Tools for computing a graphical goodness-of-fit (GOF) test based on pairwise Rosenblatt transformed data.

pairwiseCcop()

computes a (n,d,d)-array which contains pairwise Rosenblatt-transformed data.

pairwiseIndepTest()

takes such an array as input and computes a (d,d)-matrix of test results from pairwise tests of independence (as by indepTest()).

pviTest()

can be used to extract the matrix of p-values from the return matrix of pairwiseIndepTest().

gpviTest()

takes such a matrix of p-values and computes a global p-value with the method provided.

Usage

pairwiseCcop(u, copula, ...)
pairwiseIndepTest(cu.u, N=256,
        iTest = indepTestSim(n, p=2, m=2, N=N, verbose = idT.verbose, ...),
        verbose=TRUE, idT.verbose = verbose, ...)

 pviTest(piTest)
gpviTest(pvalues, method=p.adjust.methods, globalFun=min)

Arguments

Value

pairwiseCcop

(n,d,d)-array cu.u with cu.u[i,j] containing C(u_i\,|\,u_j) for i\neq j and u_i for i=j.

pairwiseIndepTest

(d,d)-matrix of lists with test results as returned by indepTest(). The test results correspond to pairwise tests of independence as conducted by indepTest().

pviTest

(d,d)-matrix of p-values.

gpviTest

global p-values for the specified methods.

Note

If u are distributed according to or “perfectly” sampled from a copula, p-values on GOF tests for that copula should be uniformly distributed in [0,1].

References

Hofert and Mächler (2014), see pairsRosenblatt.

See Also

pairsRosenblatt for where these tools are used, including demo(gof_graph) for examples.

Examples

## demo(gof_graph)

Goodness-of-fit Testing for (Nested) Archimedean Copulas

Description

gnacopula() conducts a goodness-of-fit test for the given (H_0-)copula cop based on the (copula-)data u.

NOTE: gnacopula() is deprecated, call gofCopula() instead.

Usage

gnacopula(u, cop, n.bootstrap,
          estim.method = eval(formals(enacopula)$method),
          include.K=TRUE, n.MC=0, trafo=c("Hering.Hofert", "Rosenblatt"),
          method=eval(formals(gofTstat)$method), verbose=TRUE, ...)

Arguments

Details

The function gnacopula() performs a parametric bootstrap for the goodness-of-fit test specified by trafo and method. The transformation given by trafo specifies the multivariate transformation which is first applied to the (copula-) data u (typically, the pseudo-observations are used); see htrafo() or cCopula() for more details. The argument method specifies the particular goodness-of-fit test carried out, which is either the Anderson-Darling test for the univariate standard uniform distribution (for method="AnChisq" or method="AnGamma") in a one-dimensional setup or the tests described in Genest et al. (2009) for the multivariate standard uniform distribution directly in a multivariate setup. As estimation method, the method provided by estim.method is used.

Note that a finite-sample correction is made when computing p-values; see gofCopula() for details.

A word of warning: Do work carefully with the variety of different goodness-of-fit tests that can be performed with gnacopula(). For example, among the possible estimation methods at hand, only MLE is known to be consistent (under conditions to be verified). Furthermore, for the tests based on the Anderson-Darling test statistic, it is theoretically not clear whether the parametric bootstrap converges. Consequently, the results obtained should be treated with care. Moreover, several estimation methods are known to be prone to numerical errors (see Hofert et al. (2013)) and are thus not recommended to be used in the parametric bootstrap. A warning is given if gnacopula() is called with a method not being MLE.

Value

gnacopula returns an R object of class "htest". This object contains a list with the bootstrap results including the components

p.value:

the bootstrapped p-value;

statistic:

the value of the test statistic computed for the data u;

data.name:

the name of u;

method:

a character describing the goodness-of-fit test applied;

estimator:

the estimator computed for the data u;

bootStats:

a list with component estimator containing the estimators for all bootstrap replications and component statistic containing the values of the test statistic for each bootstrap replication.

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.

Hering, C. and Hofert, M. (2011), Goodness-of-fit tests for Archimedean copulas in large dimensions, submitted.

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

gofTstat() for the implemented test statistis, htrafo() and cCopula() involved and K() for the Kendall distribution function.

gofCopula() for other (parametric bootstrap) based goodness-of-fit tests.

Goodness-of-fit Tests for Copulas

Description

The goodness-of-fit tests are based, by default, on the empirical process comparing the empirical copula with a parametric estimate of the copula derived under the null hypothesis, the default test statistic, "Sn", being the Cramer-von Mises functional S_n defined in Equation (2) of Genest, Remillard and Beaudoin (2009). In that case, approximate p-values for the test statistic can be obtained either using a parametric bootstrap (see references two and three) or by means of a faster multiplier approach (see references four and five).

Alternative test statistics can be used, in particular if a parametric bootstrap is employed.

The prinicipal function is gofCopula() which, depending on simulation either calls gofPB() or gofMB().

Usage

## Generic [and "rotCopula" method] ------ Main function ------
gofCopula(copula, x, ...)
## S4 method for signature 'copula'
gofCopula(copula, x, N = 1000,
          method = c("Sn", "SnB", "SnC", "Rn"),
          estim.method = c("mpl", "ml", "itau", "irho", "itau.mpl"),
          simulation = c("pb", "mult"), test.method = c("family", "single"),
          verbose = interactive(), ties = NA,
          ties.method = c("max", "average", "first", "last", "random", "min"),
          fit.ties.meth = eval(formals(rank)$ties.method), ...)

## (Deprecated) internal 'helper' functions : ---
gofPB(copula, x, N, method = c("Sn", "SnB", "SnC"),
      estim.method = c("mpl", "ml", "itau", "irho", "itau.mpl"),
      trafo.method = if(method == "Sn") "none" else c("cCopula", "htrafo"),
      trafoArgs = list(), test.method = c("family", "single"),
      verbose = interactive(), useR = FALSE, ties = NA,
      ties.method = c("max", "average", "first", "last", "random", "min"),
      fit.ties.meth = eval(formals(rank)$ties.method), ...)

gofMB(copula, x, N, method = c("Sn", "Rn"),
      estim.method = c("mpl", "ml", "itau", "irho"),
      test.method = c("family", "single"), verbose = interactive(),
      useR = FALSE, m = 1/2, zeta.m = 0, b = 1/sqrt(nrow(x)),
      ties.method = c("max", "average", "first", "last", "random", "min"),
      fit.ties.meth = eval(formals(rank)$ties.method), ...)

Arguments

Details

If the parametric bootstrap is used, the dependence parameters of the hypothesized copula family can be estimated by any estimation method available for the family, up to a few exceptions. If the multiplier is used, any of the rank-based methods can be used in the bivariate case, but only maximum pseudo-likelihood estimation can be used in the multivariate (multiparameter) case.

The price to pay for the higher computational efficiency of the multiplier is more programming work as certain partial derivatives need to be computed for each hypothesized parametric copula family. When estimation is based on maximization of the pseudo-likelihood, these have been implemented for six copula families thus far: the Clayton, Gumbel-Hougaard, Frank, Plackett, normal and t copula families. For other families, numerical differentiation based on grad() from package numDeriv is used (and a warning message is displayed).

Although the empirical processes involved in the multiplier and the parametric bootstrap-based test are asymptotically equivalent under the null, the finite-sample behavior of the two tests might differ significantly.

Both for the parametric bootstrap and the multiplier, the approximate p-value is computed as

(0.5 +\sum_{b=1}^N\mathbf{1}_{\{T_b\ge T\}})/(N+1),

where T and T_b denote the test statistic and the bootstrapped test statistc, respectively. This ensures that the approximate p-value is a number strictly between 0 and 1, which is sometimes necessary for further treatments. See Pesarin (2001) for more details.

For the normal and t copulas, several dependence structures can be hypothesized: "ex" for exchangeable, "ar1" for AR(1), "toep" for Toeplitz, and "un" for unstructured (see ellipCopula()). For the t copula, "df.fixed" has to be set to TRUE, which implies that the degrees of freedom are not considered as a parameter to be estimated.

The former argument print.every is deprecated and not supported anymore; use verbose instead.

Value

An object of class htest which is a list, some of the components of which are

Note

These tests were theoretically studied and implemented under the assumption of continuous margins, which implies that ties in the component samples occur with probability zero. The presence of ties in the data might substantially affect the approximate p-values. Through argument ties, the user can however select a version of the parametric bootstrap adapted to the presence of ties. No such adaption exists for the multiplier for the moment.

References

Genest, C., Huang, W., and Dufour, J.-M. (2013). A regularized goodness-of-fit test for copulas. Journal de la Société française de statistique 154, 64–77.

Genest, C. and Rémillard, B. (2008). Validity of the parametric bootstrap for goodness-of-fit testing in semiparametric models. Annales de l'Institut Henri Poincare: Probabilites et Statistiques 44, 1096–1127.

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–214.

Kojadinovic, I., Yan, J., and Holmes M. (2011). Fast large-sample goodness-of-fit tests for copulas. Statistica Sinica 21, 841–871.

Kojadinovic, I. and Yan, J. (2011). A goodness-of-fit test for multivariate multiparameter copulas based on multiplier central limit theorems. Statistics and Computing 21, 17–30.

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/.

Kojadinovic, I. (2017). Some copula inference procedures adapted to the presence of ties. Computational Statistics and Data Analysis 112, 24–41, https://arxiv.org/abs/1609.05519.

Pesarin, F. (2001). Multivariate Permutation Tests: With Applications in Biostatistics. Wiley.

See Also

fitCopula() for the underlying estimation procedure and gofTstat() for details on *some* of the available test statistics.

Examples

## The following example is available in batch through
## demo(gofCopula)

n <- 200; N <- 1000 # realistic (but too large for interactive use)
n <-  60; N <-  200 # (time (and tree !) saving ...)

## A two-dimensional data example ----------------------------------
set.seed(271)
x <- rCopula(n, claytonCopula(3))

## Does the Gumbel family seem to be a good choice (statistic "Sn")?
gofCopula(gumbelCopula(), x, N=N)
## With "SnC", really s..l..o..w.. --- with "SnB", *EVEN* slower
gofCopula(gumbelCopula(), x, N=N, method = "SnC", trafo.method = "cCopula")
## What about the Clayton family?
gofCopula(claytonCopula(), x, N=N)

## Similar with a different estimation method
gofCopula(gumbelCopula (), x, N=N, estim.method="itau")
gofCopula(claytonCopula(), x, N=N, estim.method="itau")


## A three-dimensional example  ------------------------------------
x <- rCopula(n, tCopula(c(0.5, 0.6, 0.7), dim = 3, dispstr = "un"))

## Does the Gumbel family seem to be a good choice?
g.copula <- gumbelCopula(dim = 3)
gofCopula(g.copula, x, N=N)
## What about the t copula?
t.copula <- tCopula(dim = 3, dispstr = "un", df.fixed = TRUE)
if(FALSE) ## this is *VERY* slow currently
  gofCopula(t.copula, x, N=N)

## The same with a different estimation method
gofCopula(g.copula, x, N=N, estim.method="itau")
if(FALSE) # still really slow
  gofCopula(t.copula, x, N=N, estim.method="itau")

## The same using the multiplier approach
gofCopula(g.copula, x, N=N, simulation="mult")
gofCopula(t.copula, x, N=N, simulation="mult")
if(FALSE) # no yet possible
    gofCopula(t.copula, x, N=N, simulation="mult", estim.method="itau")

Goodness-of-fit Tests for Bivariate Extreme-Value Copulas

Description

Goodness-of-fit tests for extreme-value copulas based on the empirical process comparing one of the two nonparameteric rank-based estimator of the Pickands dependence function studied in Genest and Segers (2009) with a parametric estimate of the Pickands dependence function derived under the null hypothesis. The test statistic is the Cramer-von Mises functional Sn defined in Equation (5) of Genest, Kojadinovic, G. Nešlehová, and Yan (2010). Approximate p-values for the test statistic are obtained using a parametric bootstrap.

Usage

gofEVCopula(copula, x, N = 1000,
            method = c("mpl", "ml", "itau", "irho"),
            estimator = c("CFG", "Pickands"), m = 1000,
            verbose = interactive(),
            ties.method = c("max", "average", "first", "last", "random", "min"),
            fit.ties.meth = eval(formals(rank)$ties.method), ...)

Arguments

Details

More details can be found in the second reference.

The former argument print.every is deprecated and not supported anymore; use verbose instead.

Value

An object of class htest which is a list, some of the components of which are

Note

For a given degree of dependence, the most popular bivariate extreme-value copulas are strikingly similar.

References

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

Genest, C. Kojadinovic, I., G. Nešlehová, J., and Yan, J. (2011). A goodness-of-fit test for bivariate extreme-value copulas. Bernoulli 17(1), 253–275.

See Also

evCopula, evTestC, evTestA, evTestK, gofCopula, An.

Examples

n <- 100; N <- 1000 # realistic (but too large currently for CRAN checks)
n <-  60; N <-  200 # (time (and tree !) saving ...)
x <- rCopula(n, claytonCopula(3))


## Does the Gumbel family seem to be a good choice?
gofEVCopula(gumbelCopula(), x, N=N)


## The same with different (and cheaper) estimation methods:
gofEVCopula(gumbelCopula(), x, N=N, method="itau")
gofEVCopula(gumbelCopula(), x, N=N, method="irho")


## The same with different extreme-value copulas
gofEVCopula(galambosCopula(), x, N=N)
gofEVCopula(galambosCopula(), x, N=N, method="itau")
gofEVCopula(galambosCopula(), x, N=N, method="irho")

gofEVCopula(huslerReissCopula(), x, N=N)
gofEVCopula(huslerReissCopula(), x, N=N, method="itau")
gofEVCopula(huslerReissCopula(), x, N=N, method="irho")

gofEVCopula(tevCopula(df.fixed=TRUE), x, N=N)
gofEVCopula(tevCopula(df.fixed=TRUE), x, N=N, method="itau")
gofEVCopula(tevCopula(df.fixed=TRUE), x, N=N, method="irho")

Various Goodness-of-fit Test Statistics

Description

gofBTstat() computes supposedly Beta distributed test statistics for checking uniformity of u on the unit sphere.

Usage

gofBTstat(u)

Arguments

Value

An (n,d-1)-matrix where the (i,k)th entry is

B_{ik}=\frac{\sum_{j=1}^k u_{ij}^2}{\sum_{j=1}^d u_{ij}^2}.

References

Li, R.-Z., Fang, K.-T., and Zhu, L.-X. (1997). Some Q-Q probability plots to test spherical and elliptical symmetry. Journal of Computational and Graphical Statistics 6(4), 435–450.

Examples

## generate data on the unit sphere
n <- 360
d <- 5
set.seed(1)
x <- matrix(rnorm(n*d), ncol=d)
U <- x/sqrt(rowSums(x^2))

## compute the test statistics B_k, k in {1,..,d-1}
Bmat <- gofBTstat(U)

## (graphically) check if Bmat[,k] follows a Beta(k/2, (d-k)/2) distribution
qqp <- function(k, Bmat) {
    d <- ncol(Bmat)+1L
    tit <- substitute(plain("Beta")(.A.,.B.)~~ bold("Q-Q Plot"),
                      list(.A. = k/2, .B. = (d-k)/2))
    qqplot2(Bmat[,k], qF=function(p) qbeta(p, shape1=k/2, shape2=(d-k)/2),
            main.args=list(text=tit, side=3, cex=1.3, line=1.1, xpd=NA))
}
qqp(1, Bmat=Bmat) # k=1
qqp(3, Bmat=Bmat) # k=3

Goodness-of-fit Test Statistics

Description

gofTstat() computes various goodness-of-fit test statistics typically used in gofCopula(*, simulation = "pb"). gofT2stat() computes the two-sample goodness of fit test statistic of Rémillard and Scaillet (2009).

Usage

gofTstat(u, method = c("Sn", "SnB", "SnC", "AnChisq", "AnGamma"),
         useR = FALSE, ...)
gofT2stat(u1, u2, useR = FALSE)

Arguments