General Discussion
An important parametric class of copulas is the class of Archimedean copulas. To define Archimedean copulas, we must take a look at their generators, which are unrelated to spherical generators, and must be $d$-monotone functions.
Generators and d-monotony
Archimedean generators can be defined as follows:
Definition (Archimedean generator): A $d$-Archimedean generator is a $d$-monotone function
\[\phi :\mathbb R_+ \to [0,1]\]
such that $\phi(0) = 1$ and $\phi(+\infty) = 0$.
where the notion of $d$-monotone function can be defined as follows:
Definition (d-monotony [20]): A function $\phi$ is said to be $d$-monotone if it has $d-2$ derivatives which satisfy
\[(-1)^k \phi^{(k)} \ge 0 \;\forall k \in \{1,..,d-2\},\]
and if $(-1)^{d-2}\phi^{(d-2)}$ is a non-increasing and convex function.
A function that is $d$-monotone for all $d$ is called completely monotone.
In this package, there is an abstract class Generator
that contains those generators. Many Archimedean generators are already implemented for you ! See the list of implemented archimedean generator to get an overview.
If you do not find the one you need, you may define it yourself by subtyping Generator
. The API does not ask for much information, which is really convenient. Only the two following methods are required:
- The
ϕ(G::MyGenerator,t)
function returns the value of the archimedean generator itself. - The
max_monotony(G::MyGenerator)
returns its maximum monotony, that is the greater integer $d$ making the generator $d$-monotonous.
Thus, a new generator implementation may simply look like:
struct MyGenerator{T} <: Generator
θ::T
end
ϕ(G::MyGenerator,t) = exp(-G.θ * t) # can you recognise this one ?
max_monotony(G::MyGenerator) = Inf
These two functions are enough to sample the corresponding Archimedean copula (see how in the Inverse Williamson $d$-transforms section of the documentation). However, if you know a bit more about your generator, implementing a few more simple methods can largely fasten the algorithms. You'll find more details on these methods in the Generator
docstring.
For example, Here is a graph of a few Clayton Generators:
using Copulas: ϕ,ClaytonGenerator,IndependentGenerator
using Plots
plot( x -> ϕ(ClaytonGenerator(-0.5),x), xlims=(0,5), label="ClaytonGenerator(-0.5)")
plot!(x -> ϕ(IndependentGenerator(),x), label="IndependentGenerator()")
plot!(x -> ϕ(ClaytonGenerator(0.5),x), label="ClaytonGenerator(0.5)")
plot!(x -> ϕ(ClaytonGenerator(1),x), label="ClaytonGenerator(1)")
plot!(x -> ϕ(ClaytonGenerator(5),x), label="ClaytonGenerator(5)")
And the corresponding inverse functions:
using Copulas: ϕ⁻¹,ClaytonGenerator,IndependentGenerator
using Plots
plot( x -> ϕ⁻¹(ClaytonGenerator(-0.5),x), xlims=(0,1), ylims=(0,5), label="ClaytonGenerator(-0.5)")
plot!(x -> ϕ⁻¹(IndependentGenerator(),x), label="IndependentGenerator()")
plot!(x -> ϕ⁻¹(ClaytonGenerator(0.5),x), label="ClaytonGenerator(0.5)")
plot!(x -> ϕ⁻¹(ClaytonGenerator(1),x), label="ClaytonGenerator(1)")
plot!(x -> ϕ⁻¹(ClaytonGenerator(5),x), label="ClaytonGenerator(5)")
Copulas.Generator
— TypeGenerator
Abstract type. Implements the API for archimedean generators.
An Archimedean generator is simply a function $\phi :\mathbb R_+ \to [0,1]$ such that $\phi(0) = 1$ and $\phi(+\infty) = 0$.
To generate an archimedean copula in dimension $d$, the function also needs to be $d$-monotone, that is :
- $\phi$ is $d-2$ times derivable.
- $(-1)^k \phi^{(k)} \ge 0 \;\forall k \in \{1,..,d-2\},$ and if $(-1)^{d-2}\phi^{(d-2)}$ is a non-increasing and convex function.
The access to the function $\phi$ itself is done through the interface:
ϕ(G::Generator, t)
We do not check algorithmically that the proposed generators are d-monotonous. Instead, it is up to the person implementing the generator to tell the interface how big can $d$ be through the function
max_monotony(G::MyGenerator) = # some integer, the maximum d so that the generator is d-monotonous.
More methods can be implemented for performance, althouhg there are implement defaults in the package :
ϕ⁻¹( G::Generator, x)
gives the inverse function of the generator.ϕ⁽¹⁾(G::Generator, t)
gives the first derivative.ϕ⁽ᵏ⁾(G::Generator, k, t)
gives the kth derivative.williamson_dist(G::Generator, d)
gives the Wiliamson d-transform of the generator, see WilliamsonTransforms.jl.
References:
- [20] McNeil, A. J., & Nešlehová, J. (2009). Multivariate Archimedean copulas, d-monotone functions and ℓ 1-norm symmetric distributions.
Note that the rate at which these functions are reaching 0 (and their inverse reaching infinity on the left boundary) can vary a lot from one to the other. Note also that the difference between each of them is easier to grasp on the inverse plot.
Williamson d-transform
An easy way to construct new $d$-monotonous generators is the use of the Williamson $d$-transform.
Definition (Williamson d-transformation): For a univariate non-negative random variable $X$, with cumulative distribution function $F$ and an integer $d\ge 2$, the Williamson-d-transform of $X$ is the real function supported on $[0,\infty[$ given by:
\[\phi(t) = 𝒲_{d}(X)(t)\]
\[=\int_{t}^{\infty} \left(1 - \frac{t}{x}\right)^{d-1} dF(x)\]
\[= \mathbb E\left( (1 - \frac{t}{X})^{d-1}_+\right) \mathbb 1_{t > 0} + \left(1 - F(0)\right)\mathbb 1_{t <0}\]
In this package, we implemented it through the WilliamsonGenerator
class. It can be used as follows:
WilliamsonGenerator(X::UnivariateRandomVariable, d)
.
This function computes the Williamson d-transform of the provided random variable $X$ using the WilliamsonTransforms.jl
package. See [20, 21] for the literature.
The $d$-transform of a positive random variable is $d$-monotonous but not $k$-monotonous for any $k > d$. Its max monotony is therefore $d$. This has a few implications, one of the biggest one is that the $d$-variate Archimedean copula that corresponds has no density.
More genrally, if you want your Archimedean copula to have a density, you have to use a generator that is more-monotonous that the dimension of your model.
Copulas.WilliamsonGenerator
— TypeWilliamsonGenerator{TX}
i𝒲{TX}
Fields:
X::TX
– a random variable that represents its Williamson d-transformd::Int
– the dimension of the transformation.
Constructor
WilliamsonGenerator(X::Distributions.UnivariateDistribution, d)
i𝒲(X::Distributions.UnivariateDistribution,d)
The WilliamsonGenerator
(alias i𝒲
) allows to construct a d-monotonous archimedean generator from a positive random variable X::Distributions.UnivariateDistribution
. The transformation, which is called the inverse Williamson transformation, is implemented in WilliamsonTransforms.jl.
For a univariate non-negative random variable $X$, with cumulative distribution function $F$ and an integer $d\ge 2$, the Williamson-d-transform of $X$ is the real function supported on $[0,\infty[$ given by:
\[\phi(t) = 𝒲_{d}(X)(t) = \int_{t}^{\infty} \left(1 - \frac{t}{x}\right)^{d-1} dF(x) = \mathbb E\left( (1 - \frac{t}{X})^{d-1}_+\right) \mathbb 1_{t > 0} + \left(1 - F(0)\right)\mathbb 1_{t <0}\]
This function has several properties:
- We have that $\phi(0) = 1$ and $\phi(Inf) = 0$
- $\phi$ is $d-2$ times derivable, and the signs of its derivatives alternates : $\forall k \in 0,...,d-2, (-1)^k \phi^{(k)} \ge 0$.
- $\phi^{(d-2)}$ is convex.
These properties makes this function what is called a d-monotone archimedean generator, able to generate archimedean copulas in dimensions up to $d$. Our implementation provides this through the Generator
interface: the function $\phi$ can be accessed by
G = WilliamsonGenerator(X, d)
ϕ(G,t)
Note that you'll always have:
max_monotony(WilliamsonGenerator(X,d)) === d
References:
Inverse Williamson d-transform
The Williamson d-transform is a bijective transformation[1] from the set of positive random variables to the set of generators. It therefore has an inverse transformation (called, surprisingly, the inverse Williamson $d$-transform) that construct the positive random variable R from a generator $\phi$.
This transformation is implemented through one method in the Generator interface that is worth talking a bit about : williamson_dist(G::Generator, d)
. This function computes the inverse Williamson d-transform of the d-monotone archimedean generator ϕ, still using the WilliamsonTransforms.jl
package. See [20, 21].
To put it in a nutshell, for $\phi$ a $d$-monotone archimedean generator, the inverse Williamson-d-transform of $\\phi$ is the cumulative distribution function $F$ of a non-negative random variable $R$, defined by :
\[F(x) = 𝒲_{d}^{-1}(\phi)(x) = 1 - \frac{(-x)^{d-1} \phi_+^{(d-1)}(x)}{k!} - \sum_{k=0}^{d-2} \frac{(-x)^k \phi^{(k)}(x)}{k!}\]
The WilliamsonTransforms.jl
package implements this transformation (and its inverse, the Williamson d-transfrom) in all generality. It returns this cumulative distribution function in the form of the corresponding random variable <:Distributions.ContinuousUnivariateDistribution
from Distributions.jl
. You may then compute :
- The cdf via
Distributions.cdf
- The pdf via
Distributions.pdf
and the logpdf viaDistributions.logpdf
- Samples from the distribution via
rand(X,n)
.
As an example of a generator produced by the Williamson transformation and its inverse, we propose to construct a generator from a LogNormal distribution:
using Distributions
using Copulas: i𝒲, ϕ⁻¹, IndependentGenerator
using Plots
G = i𝒲(LogNormal(), 2)
plot(x -> ϕ⁻¹(G,x), xlims=(0.1,0.9), label="G")
plot!(x -> ϕ⁻¹(IndependentGenerator(),x), label="Independence")
The i𝒲
alias stands for WiliamsonGenerator
. To stress the generality of the approach, remark that any positive distribution is allowed, including discrete ones:
using Distributions
using Copulas: i𝒲, ϕ⁻¹
using Plots
G1 = i𝒲(Binomial(10,0.3), 2)
G2 = i𝒲(Binomial(10,0.3), 3)
plot(x -> ϕ⁻¹(G1,x), xlims=(0.1,0.9), label="G1")
plot!(x -> ϕ⁻¹(G2,x), xlims=(0.1,0.9), label="G2")
As obvious from the definition of the Williamson transform, using a discrete distribution produces piecewise-linear generators, where the number of pieces is dependent on the order of the transformation.
Archimedean copulas
Let's first define formally archimedean copulas:
Definition (Archimedean copula): If $\phi$ is a $d$-monotonous Archimedean generator, then the function
\[C(\bm u) = \phi\left(\sum\limits_{i=1}^d \phi^{-1}(u_i)\right)\]
is a copula.
There are a few archimedean generators that are worth noting since they correspond to known archimedean copulas families:
IndependentGenerator
: $\phi(t) =e^{-t} \text{ generates } \Pi$.ClaytonGenerator
: $\phi_{\theta}(t) = \left(1+t\theta\right)^{-\theta^{-1}}$ generates the $\mathrm{Clayton}(\theta)$ copula.GumbelGenerator
: $\phi_{\theta}(t) = \exp\{-t^{\theta^{-1}}\}$ generates the $\mathrm{Gumbel}(\theta)$ copula.FrankGenerator
: $\phi_{\theta}(t) = -\theta^{-1}\ln\left(1+e^{-t-\theta}-e^{-t}\right)$ generates the $\mathrm{Franck}(\theta)$ copula.
There are a lot of others implemented in the package, see our large list of implemented archimedean generator.
Archimedean copulas have a nice decomposition, called the Radial-simplex decomposition:
Property (Radial-simplex decomposition [20, 22]): A $d$-variate random vector $\bm U$ following an Archimedean copula with generator $\phi$ can be decomposed into
\[\bm U = \phi.(\bm S R),\]
where $\bm S$ is uniform on the $d$-variate simplex and $R$ is a non-negative random variable, independent form $\bm S$, defined as the inverse Williamson $d$-transform of $\phi$.
This is why williamson_dist(G::Generator,d)
is such an important function in the API: it allows to generator the radial part and sample the Archimedean copula. You may call this function directly to see what distribution will be used:
using Copulas: williamson_dist, FrankCopula
williamson_dist(FrankCopula(3,7))
Copulas.WilliamsonFromFrailty{Copulas.Logarithmic{Float64}, 3}(
frailty_dist: Copulas.Logarithmic{Float64}(α=0.9990881180344455, h=-7.0)
)
For the Frank Copula, as for many classic copulas, the distribution used is known. We pull some of them from Distributions.jl
but implement a few more, as this Logarithmic one. Another useful example are negatively-dependent Clayton copulas:
using Copulas: williamson_dist, ClaytonCopula
williamson_dist(ClaytonCopula(3,-0.2))
Copulas.ClaytonWilliamsonDistribution{Float64, Int64}(θ=-0.2, d=3)
for which the corresponding distribution is known but has no particular name, thus we implemented it under the ClaytonWilliamsonDistribution
name.
It is well-known that completely monotone generators are Laplace transforms of non-negative random variables. This gives rise to another decomposition:
Property (Frailty decomposition [23]: When $\phi$ is completely monotone, it is the Laplace transform of a non-negative random variable $W$ such that
\[\bm U = \phi(\bm Y / W),\]
where $\bm Y$ is a vector of independent and identically distributed (i.i.d.) exponential distributions.
The link between the distribution of $R$ and the distribution of $W$ can be explicited. We exploit this link and provide the WilliamsonFromFrailty()
constructor that construct the distribution of $R$ from the distribution of $W$ and returns the corresponding WilliamsonGenerator
from the frailty distribution itself. The corresponding ϕ is simply the laplace transform of $W$. This is another potential way of constructing new archimedean copulas !
We use this fraily approach for several generators, since sometimes it is faster, including e.g. the Clayton one with positive dependence:
using Copulas: williamson_dist, ClaytonCopula
williamson_dist(ClaytonCopula(3,10))
Copulas.WilliamsonFromFrailty{Distributions.Gamma{Float64}, 3}(
frailty_dist: Distributions.Gamma{Float64}(α=0.1, θ=10.0)
)
Copulas.ArchimedeanCopula
— TypeArchimedeanCopula{d, TG}
Fields: - G::TG : the generator <: Generator.
Constructor:
ArchimedeanCopula(d::Int,G::Generator)
For some Archimedean Generator
G::Generator
and some dimenson d
, this class models the archimedean copula which has this generator. The constructor checks for validity by ensuring that max_monotony(G) ≥ d
. The $d$-variate archimedean copula with generator $\phi$ writes:
\[C(\mathbf u) = \phi^{-1}\left(\sum_{i=1}^d \phi(u_i)\right)\]
The default sampling method is the Radial-simplex decomposition using the Williamson transformation of $\phi$.
There exists several known parametric generators that are implement in the package. For every NamedGenerator <: Generator
implemented in the package, we provide a type alias `NamedCopula{d,...} = ArchimedeanCopula{d,NamedGenerator{...}}
to be able to manipulate the classic archimedean copulas without too much hassle for known and usefull special cases.
A generic archimdean copula can be constructed as follows:
struct MyGenerator <: Generator end
ϕ(G::MyGenerator,t) = exp(-t) # your archimedean generator, can be any d-monotonous function.
max_monotony(G::MyGenerator) = Inf # could depend on generators parameters.
C = ArchimedeanCopula(d,MyGenerator())
The obtained model can be used as follows:
spl = rand(C,1000) # sampling
cdf(C,spl) # cdf
pdf(C,spl) # pdf
loglikelihood(C,spl) # llh
Bonus: If you know the Williamson d-transform of your generator and not your generator itself, you may take a look at WilliamsonGenerator
that implements them. If you rather know the frailty distribution, take a look at WilliamsonFromFrailty
.
References:
<!– TODO: Make a few graphs of bivariate archimedeans pdfs and cdfs. And provide a few more standard tools for these copulas ? –>
- [20]
- A. J. McNeil and J. Nešlehová. Multivariate Archimedean Copulas, d -Monotone Functions and L1 -Norm Symmetric Distributions. The Annals of Statistics 37, 3059–3097 (2009).
- [21]
- R. E. Williamson. On multiply monotone functions and their laplace transforms (Mathematics Division, Office of Scientific Research, US Air Force, 1955).
- [22]
- A. J. McNeil. Sampling Nested Archimedean Copulas. Journal of Statistical Computation and Simulation 78, 567–581 (2008).
- [23]
- M. Hofert, M. Mächler and A. J. McNeil. Archimedean Copulas in High Dimensions: Estimators and Numerical Challenges Motivated by Financial Applications. Journal de la Société Française de Statistique 154, 25–63 (2013).
- 1This bijection is to be taken carefuly: the bijection is between random variables with unit scales and generators with common value at 1, sicne on both rescaling does not change the underlying copula.