Copulas.jl 1.0: from a Distributions.jl experiment to a stable dependence-modeling toolkit

After more than four years of developpement, Copulas.jl has reached version 1.0 ! I wanted to come back to the developpement journey, talk a bit about the functionalities but most importantly about why: why tag the 1.0 version now. This is why :

Copulas.jl now has a clear public contract that can be maintained under semantic versioning.

Copulas.jl started from a simple design idea: copulas are distribution functions, so in Julia they should behave like distributions. That idea is still in the middle of the scene, while the surrounding implementation and API have expanded substantially. This post is therefore both a 1.0 announcement and a summary of the package’s evolution.

What can you do with it?

The basic usage has barely changed. Pick a copula, add margins, and you have a distribution:

using Copulas, Distributions

C = ClaytonCopula(3, 2.0)
D = SklarDist(C, (Gamma(2, 3), Beta(1, 4), Normal()))
X = rand(D, 1_000) # One observation per column

You can now extract marginals, condition on observations, and apply Rosenblatt transforms:

C23 = subsetdims(C, (2, 3))
Cgiven = condition(C, 2, 0.4)

S = rosenblatt(D, X)
X2 = inverse_rosenblatt(D, S)

The conditional object is itself a distribution, so you can use rand, cdf, pdf, etc. on it too.

Fitting has received quite a bit of work. CopulaModel stores a fitted model and exposes the usual statistical summaries through StatsBase. For example, to fit the copula to uniform observations:

using StatsBase

U = rand(C, 1_000)
M = fit(CopulaModel, ClaytonCopula, U; method=:mle)
Cfit = fitted_distribution(M)

I = infer(M)
confint(I; level=0.95)

# H₀: the observations follow a Clayton copula
test = GOFCopulaTest(M; N=999)
pvalue(test)

Here fit estimates the parameters, and infer computes uncertainty afterwards. GOFCopulaTest checks the fit using a parametric bootstrap, refitting the model in each of the 999 replicates. Maximum pseudo-likelihood, rank-based estimators, bootstrap and jackknife procedures are also available. You can compare a set of candidate families using information criteria, and test independence, exchangeability, radial symmetry, extreme-value dependence or goodness of fit.

There are plotting recipes, integrations with ExpectationMaximization.jl and PartitionedDistributions.jl, etc, all thanks to the Distributions.jl's API.

Four years later

In the first announcement, I had a few families, sampling, CDFs and some fitting methods. Densities for generic Archimedean copulas were still on the todo list. So were nested copulas, Bernstein and Beta copulas, and checkerboard constructions. These are now implemented :)

A lot of my time went into Archimedean generators, and particularly the Williamson transform. For a positive random variable $R$ with distribution function $F_R$, its Williamson transform of order $d>1$ is

$$ \mathcal{W}_d(R)(t) = \mathbb{E}\left[\left(1-\frac{t}{R}\right)_+^{d-1}\right] = \int_t^\infty \left(1-\frac{t}{r}\right)^{d-1},\mathrm{d}F_R(r), \qquad t\geq 0, $$

where $x_+=\max(x,0)$. This turns a radial distribution into an Archimedean generator; the inverse recovers the radial distribution from a suitable generator. Both are implemented, both for integer or real orders:

using Distributions
using Copulas: 𝒲, 𝒲₋₁

R = LogNormal() # Any positive distribution.
G = 𝒲(R, 3.5) # Williamson generator of order d
R_bis = 𝒲₋₁(G, 3.5) # Recover the law of R. 
R2 = 𝒲₋₁(G, 2) # Another distribution. 

This gives us a common mathematical basis for implementing Archimedean and Liouville copulas generically, across all admissible generators and parameters. Archimax and nested Archimedean constructions can reuse that machinery too, with their own validity conditions.

The extreme-value models also moved beyond the bivariate case. They are now built around stable tail dependence functions, with multivariate models and empirical estimation.

One important goals from the begining of the project was the compatibility with other floating point types than Float64. It took until now to make sampling and conditioning reliably preserve Float32 and BigFloat on supported paths, because this required heavy work under the hood. It is not very visible, but makes using the package much more pleasant on the long run :)

Why tag 1.0 now?

I wanted to settle which parts of the package users can depend on. Until recently, that boundary was too implicit: a function used internally could easily look like an interface intended for downstream code. The source and documentation have been reviewed with that question in mind.

This is what I am committing to for 1.x: keeping the documented interface compatible while continuing to add models and improve the algorithms. There is still work to do, of course. The package now has tests and regression benchmarks to help with that.

Thanks to everyone who contributed code, examples, bug reports and reviews, in particular @Santymax98 but other contributors too.

Suggestions and contributions are, as always, welcome :)

Documentation · GitHub · Discourse thread