Multivariate random vectors

This section gives some general definitions and tools about dependence structures, multivariate random vectors and copulas. Along this journey through the mathematical theory of copulas, we link to the rest of the documentation for more specific and detailed arguments on particular points, or simply to the technical documentation of the actual implementation. The interested theoretical reader can take a look at the standard books on the subject [1–4] or more recently [5–8].

We start here by defining a few concepts about dependence structures and copulas. Consider a real valued random vector $\bm X = \left(X_1,...,X_d\right): \Omega \to \mathbb R^d$. The random variables $X_1,...,X_d$ are called the marginals of the random vector $\bm X$.

using Distributions
X₁ = Normal()       # A standard gaussian random variable
X₂ = Gamma(2,3)     # A Gamma random variable
X₃ = Pareto(1)      # A Pareto random variable with no variance.
X₄ = LogNormal(0,1) # A Lognormal random variable

The probability distribution of the random vector $\bm X$ can be characterized by its distribution function $F$:

\[\begin{align*} F(\bm x) &= \mathbb P\left(\bm X \le \bm x\right)\\ &= \mathbb P\left(\forall i \in \{1,...,d\},\; X_i \le x_i\right). \end{align*}\]

For a function $F$ to be the distribution function of some random vector, it should be $d$-increasing, right-continuous and left-limited. For $i \in \{1,...,d\}$, the random variables $X_1,...,X_d$, called the marginals of the random vector, also have distribution functions denoted $F_1,...,F_d$ and defined by :

\[F_i(x_i) = F(+\infty,...,+\infty,x_i,+\infty,...,+\infty).\]

Note that the range $\mathrm{Ran}(F)$ of a distribution function $F$, univariate or multivariate, is always contained in $[0,1]$. When the random vector or random variable is absolutely continuous with respect to (w.r.t.) the Lebesgue measure restricted to its domain, the range is exactly $[0,1]$. When the distribution is discrete with $n$ atoms, the range is a finite set of $n+1$ values in $[0,1]$.

Copulas and Sklar's Theorem

There is a fundamental functional link between the function $F$ and its marginals $F_1,...,F_d$. This link is expressed by the mean of copulas.

Definition (Copula) : A copula, usually denoted $C$, is the distribution function of a random vector with marginals that are all uniform on $[0,1]$, i.e.

\[C_i(u) = u\mathbb 1_{u \in [0,1]} \text{ for all }i \in 1,...,d.\]

You may define a copula object in Julia by simply calling its constructor:

using Copulas
d = 4 # The dimension of the model
θ = 7 # Parameter
C = ClaytonCopula(4,7) # A 4-dimensional clayton copula with parameter θ = 7.

ClaytonCopula{4, Int64}(
G: Copulas.ClaytonGenerator{Int64}(7)
)

This object is a random vector, and behaves exactly as you would expect a random vector from Distributions.jl to behave: you may sample it with rand(C,100), compute its pdf or cdf with pdf(C,x) and cdf(C,x), etc:

u = rand(C,10)

4×10 Matrix{Float64}:
 0.541689  0.904653  0.838058  0.991346  …  0.592455  0.315444  0.265217
 0.458433  0.706722  0.878939  0.807879     0.784658  0.26961   0.310556
 0.425349  0.540633  0.567612  0.822973     0.658182  0.248751  0.404481
 0.374389  0.598137  0.795738  0.881928     0.781849  0.175403  0.327209

cdf(C,u)

10-element Vector{Float64}:
 0.3463681797781641
 0.5040580359968292
 0.5564436823278354
 0.7325242970891881
 0.4369133697107717
 0.2840082050036614
 0.2533658925860069
 0.5498855064373659
 0.1718841000477894
 0.2477081195204077

One of the reasons that makes copulas so useful is discovered by Sklar [9] in 1959:

Theorem (Sklar): For every random vector $\bm X$, there exists a copula $C$ such that

\[\forall \bm x\in \mathbb R^d, F(\bm x) = C(F_{1}(x_{1}),...,F_{d}(x_{d})).\]

The copula $C$ is uniquely determined on $\mathrm{Ran}(F_{1}) \times ... \times \mathrm{Ran}(F_{d})$, where $\mathrm{Ran}(F_i)$ denotes the range of the function $F_i$. In particular, if all marginals are absolutely continuous, $C$ is unique.

This result allows to decompose the distribution of $\bm X$ into several components: the marginal distributions on one side, and the copula on the other side, which governs the dependence structure between the marginals. This object is central in our work, and therefore deserves a moment of attention.

Example (Independence): The function

\[\Pi : \bm x \mapsto \prod_{i=1}^d x_i = \bm x^{\bm 1}\]

is a copula, corresponding to independent random vectors.

The independence copula can be constructed using the IndependentCopula(d) syntax as follows:

Π = IndependentCopula(d) # A 4-variate independence structure.

We leverage the Sklar theorem to construct multivariate random vectors from a copula-marginals specification. This can be used as follows:

MyDistribution = SklarDist(Π, (X₁,X₂,X₃,X₄))
MyOtherDistribution = SklarDist(C, (X₁,X₂,X₃,X₄))

And the API is still the same:

rand(MyDistribution,10)

4×10 Matrix{Float64}:
 0.506068   1.1736   0.115507  0.761589  …  1.69936   0.832913  -0.92773
 3.82966    2.21575  5.34875   4.2604       1.70107   0.720446  16.6845
 1.29203   11.5124   2.76127   1.70539      1.47626   1.03674    1.3449
 0.621665   1.78511  0.184542  0.471222     0.210294  1.78814    2.83079

rand(MyOtherDistribution,10)

4×10 Matrix{Float64}:
  2.60396  -0.533472   1.27053  -0.349767  …  -0.270593  -1.24575   2.2978
 10.2072    3.55393   10.4297    3.47268       3.63315    2.11319   9.32348
  5.26268   1.42884    4.54443   2.2984        1.55669    1.14145   6.32431
  1.96467   0.786505   2.05202   0.809932      0.639576   0.279512  6.11862

On the other hand, the pseudo() function computes ranks, effectively using Sklar's theorem the other way around (from the marginal space to the unit hypercube).

Copulas are bounded functions with values in [0,1] since they correspond to probabilities. But their range can be bounded more precisely:

Property (Fréchet-Hoeffding bounds [10]): For all $\bm x \in [0,1]^d$, every copula $C$ satisfies :

\[\langle \bm 1, \bm x - 1 + d^{-1}\rangle_{+} \le C(\bm x) \le \min \bm x,\]

where $y_{+} = \max(0,y)$.

The function $M : \bm x \mapsto \min\bm x$, called the upper Fréchet-Hoeffding bound, is a copula. The function $W : \bm x \mapsto \langle \bm 1, \bm x - 1 + d^{-1}\rangle_{+}$, called the lower Fréchet-Hoeffding bound, is on the other hand a copula only when $d=2$. These two copulas can be constructed through MCopula(d) and WCopula(2).

The upper Fréchet-Hoeffding bound corresponds to the case of comonotone random vector: a random vector $\bm X$ is said to be comonotone, i.e., to have copula $M$, when each of its marginals can be written as a non-decreasing transformation of the same random variable (say with $\mathcal U\left([0,1]\right)$ distribution). This is a simple but important dependence structure. See e.g.,[11, 12] on this particular copula. Note that the implementation of their sampler was straightforward due to their particular shapes:

rand(MCopula(2),10) # sampled values are all equal, this is comonotony

2×10 Matrix{Float64}:
 0.435569  0.959905  0.429309  0.438606  …  0.77662  0.911544  0.272435
 0.435569  0.959905  0.429309  0.438606     0.77662  0.911544  0.272435

u = rand(WCopula(2),10)
sum(u, dims=1) # sum is always equal to one, this is anticomonotony

1×10 Matrix{Float64}:
 1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0

Since copulas are distribution functions, like distribution functions of real-valued random variables and random vectors, there exists classical and useful parametric families of copulas. This is mostly the content of this package, and we refer to the rest of the documentation for more details on the models and their implementations.

Fitting copulas and compound distributions.

Distributions.jl's API contains a fit function for random vectors and random variables. We propose an implementation of it for copulas and multivariate compound distributions (composed of a copula and some given marginals). It can be used as follows:

using Copulas, Distributions, Random
# Construct a given model:
X₁ = Gamma(2,3)
X₂ = Pareto()
X₃ = LogNormal(0,1)
C = ClaytonCopula(3,0.7) # A 3-variate Clayton Copula with θ = 0.7
D = SklarDist(C,(X₁,X₂,X₃)) # The final distribution

simu = rand(D,1000) # Generate a dataset

# You may estimate a copula using the `fit` function:
D̂ = fit(SklarDist{ClaytonCopula,Tuple{Gamma,Normal,LogNormal}}, simu)

SklarDist{ClaytonCopula{3, Float64}, Tuple{Distributions.Gamma{Float64}, Distributions.Normal{Float64}, Distributions.LogNormal{Float64}}}(
C: ClaytonCopula{3, Float64}(
G: Copulas.ClaytonGenerator{Float64}(0.741823930441656)
)

m: (Distributions.Gamma{Float64}(α=1.875593711812958, θ=3.3061357979386847), Distributions.Normal{Float64}(μ=7.972312768030163, σ=40.972726287187704), Distributions.LogNormal{Float64}(μ=-0.04584673741832515, σ=0.9867461240185885))
)

We see on the output that the parameters were correctly estimated from this sample. More details on the estimator, including, e.g., standard errors, may be obtained with more complicated estimation routines. For a Bayesian approach using Turing.jl, see this example.

Going further

There are a lot of available copula families in the package, that can be regrouped into a few classes:

• Elliptical Copulas, including the Gaussian and Student cases.
• Archimedean Copulas, leveraging their Archimedean Generators
• Fréchet-Hoeffding bounds,
• Other Copulas

Each of these classes more-or-less correspond to an abstract type in our type hierarchy, and to a section of this documentation.