Copulas
What is this package ?
Warning: This is fairly untested and experimental work and the API might change without notice.
This package brings most standard copula features into native Julia: random number generation, pdf and cdf, fitting, copulabased multivariate distributions through Sklar's theorem, etc., while fully complying with the 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
.
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 3variate Clayton 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 theArchimedeanCopula
abstract type, see below,WCopula
andMCopula
, which are FréchetHoeffding 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 generallyPatchworkCopula
)
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 and (log)density in any dimension, 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.  We plan on implementing the Williamson transformations so that
radialdist
can be automaticlaly deduced fromϕ
and vice versa, if you dont know much about your archimedean family
Dev Roadmap
Urgent things
 Add tests and documentation
Next
 Extensive documentation and tests for the current implementation.
 Implement archimedean density generally.
 Docs: show how to implement another archimedean.
 Give the user the choice of fitting method via
fit(dist,data; method="MLE")
orfit(dist,data; method="itau")
orfit(dist,data; method="irho")
.  Fitting a generic archimedean : should provide an empirical generator
 Make
Archimedean
more generic : inputing onlyradial_dist
or onlyphi
shoudl be enough to getpdf, cdf, rand, tau, rho, itau, irho, fit, radial_dist
, etc... Williamson dtransform and inverse dtransform should be implemented. The checking of nesting possibility should be done automatically with some rules (is phi_inv \circ phi complementely monotonous ? with obviously shortcut for interfamily nestings.)
Maybe later

Vines
? 
NestedArchimedean
and very easy implementation of new archimeean copulas via the radial dist or the phi/invphi + Williamson transform. 
BernsteinCopula
andBetaCopula
could also be implemented. 
PatchworkCopula
andCheckerboardCopula
: could be nice things to have :)  Goodness of fits tests ?
Contributions are welcome
Do not hesitate to open an issue to discuss :)