Announce

I am proud to annonce the publication and the registration of my new Julia package, Copulas.jl.

As it’s name suggests, Copulas.jl is a package that implements methods and tools to work with an arround copulas in the Julia programming language. Quoting the github page :

This package brings most standard [copula](https://en.wikipedia.org/wiki/Copula_(probability_theory)) features into native Julia: random number generation, pdf and cdf, fitting, copula-based multivariate distributions through Sklar's theorem, etc., while fully complying with the [`Distributions.jl`](https://github.com/JuliaStats/Distributions.jl) API (after all, copulas are distributions functions) in order to provide interoperability with other packages based on this API such as, e.g., [`Turing.jl`](https://github.com/TuringLang/Turing.jl).

Usually, people that use and work with copulas turn to R, because of the amazing R package copula. While it is still well maintained and regularly updated, the R package copula is a mixture of obscure, heavily optimized C code and more standard R code, which makes it a complicated code base for readability, extensibility, reliability and maintenance.

This is an attempt to provide a very light, fast, reliable and maintainable copula implementation in native Julia (which means, in particular, floating point type agnostic, i.e. compatibility with BigFloat, DoubleFloats, MultiFloats and other kind of numbers). The two most important exported types are:

• Copula: an abstract mother type for all the copulas in the package.
• SklarDist: allows construction of a multivariate distribution by specifying the copula and the marginals through Sklar’s theorem.

What is already implemented

The API contains random number generation, cdf and pdf evaluation, and the fit function from Distributions.jl. A typical use case might look like this:

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

# This generates a (3,1000)-sized dataset from the multivariate distribution D
simu = rand(D,1000)

# While the following estimates the parameters of the model from a dataset: 
D̂ = fit(SklarDist{FrankCopula,Tuple{Gamma,Normal,LogNormal}}, simu)
# Increase the number of observations to get a beter fit (or not?)  

Available copula families are:

• GaussianCopula,
• TCopula,
• ArchimedeanCopula (for any generator),
• ClaytonCopula,FrankCopula, AMHCopula, JoeCopula, GumbelCopula as example of the ArchimedeanCopula abstract type, see below,
• WCopula and MCopula, which are Fréchet-Hoeffding bounds,
• EmpiricalCopula to follow closely a given dataset.

The next ones to be implemented will probably be:

• Nested archimedeans (general, with the possibility to nest any family with any family, assuming it is possible, with parameter checks.)
• Bernstein copula and more general Beta copula as smoothing of the Empirical copula.
• CheckerboardCopula (and more generally PatchworkCopula)

Adding a new ArchimedeanCopula is very easy. The Clayton implementation is as short as:

struct ClaytonCopula{d,T} <: ArchimedeanCopula{d}
    θ::T
end
ClaytonCopula(d,θ)            = ClaytonCopula{d,typeof(θ)}(θ)     # Constructor
ϕ(C::ClaytonCopula, t)        = (1+sign(C.θ)*t)^(-1/C.θ)          # Generator
ϕ⁻¹(C::ClaytonCopula,t)       = sign(C.θ)*(t^(-C.θ)-1)            # Inverse Generator
τ(C::ClaytonCopula)           = C.θ/(C.θ+2)                       # θ -> τ
τ⁻¹(::Type{ClaytonCopula},τ)  = 2τ/(1-τ)                          # τ -> θ
radial_dist(C::ClaytonCopula) = Distributions.Gamma(1/C.θ,1)      # Radial distribution

The Archimedean API is modular:

• To sample an archimedean, only radial_dist and ϕ are needed.
• To evaluate the cdf, only ϕ and ϕ⁻¹ are needed
• Currently, to fit the copula τ⁻¹ is needed as we use the inverse tau moment method. But we plan on also implementing inverse rho and MLE (density needed).
• Note that the generator ϕ follows the convention ϕ(0)=1, while others (e.g., https://en.wikipedia.org/wiki/Copula_(probability_theory)#Archimedean_copulas) use ϕ⁻¹ as the generator.
• The density and thus log-likelyhood of all archimedean copulas is upcoming soon.

Next step

The main thing to do on this package is to add test and documentations to the code that is already done. Suggestions and contributions are, of course, welcomed !