lrnv / copulas.jl

A fully `Distributions.jl`-compliant copula package

Copulas.jl

A fully `Distributions.jl`-compliant copula package

Copulas.jl brings most standard copula features into native Julia: random number generation, pdf and cdf, fitting, copula-based multivariate distributions through Sklar's theorem, etc. Since copulas are distribution functions, we fully comply with the Distributions.jl API. This allows interoperability with the broader ecosystem, 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 package R::copula. While well-maintained and regularly updated, R::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. One of the notable benefits of such a native implementation (among others) is the floating point type agnosticity, i.e. compatibility with BigFloat, DoubleFloats, MultiFloats, etc.

The package revolves around two main types:

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

Warning: This is fairly experimental work, use with caution.

Getting started

The package is registered in Julia's General registry so you may simply install the package by running :

] add Copulas

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 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)

The list of availiable copula models is very large, check it out on our documentation !

The general implementation philosophy is for the code to follow the mathematical boundaries of the implemented concepts. For example, this is the only implementation we know (in any language) that allows for all Archimedean copulas to be sampled: we use the Williamson transformation for non-standard generators, including user-provided black-box ones.

Feature comparison

There are competing packages in Julia, such as BivariateCopulas.jl which only deals with a few models in bivariate settings but has very nice graphs, or DatagenCopulaBased.jl, which only provides sampling and does not have exactly the same models as Copulas.jl. While not fully covering out both of these package's functionality (mostly because the three projects chose different implementation paths), Copulas.jl brings, as a key feature, the compliance with the broader ecosystem. The following table provides a feature comparison between the three:

	Copulas.jl	DatagenCopulaBased.jl	BivariateCopulas.jl
Distributions.jl's API	✔️	❌	✔️
Fitting	✔️	❌	❌
Plotting	❌	❌	✔️
Available copulas
- Classic Bivariate	✔️	✔️	✔️
- Classic Multivariate	✔️	✔️	❌
- Archimedeans	✔️ (All of them)	⚠️ Selected ones	⚠️Selected ones
- Obscure Bivariate	✔️	❌	❌
- Archimedean Chains	❌	✔️	❌

Since our primary target is maintainability and readability of the implementation, we did not consider the efficiency and the performance of the code yet. Proper benchmarks will come in the near future.

Contributions are welcome

If you want to contribute to the package, ask a question, found a bug or simply want to chat, do not hesitate to open an issue on this repo. General guidelines on collaborative practices (colprac) are available at https://github.com/SciML/ColPrac.

Citation

Do not hesitate to star this repository to show support. If you use this package in your researches, please cite it with the following bibtex code:

@article{Laverny2024,
    author = {Oskar Laverny and Santiago Jimenez}, 
    title = {Copulas.jl: A fully Distributions.jl-compliant copula package},
    journal = {Journal of Open Source Software},
    doi = {10.21105/joss.06189}, 
    url = {https://doi.org/10.21105/joss.06189}, 
    year = {2024}, 
    publisher = {The Open Journal}, 
    volume = {9}, 
    number = {94}, 
    pages = {6189}
}

copulas.jl's Issues

[New copulas] Vines Copulas

It would be nice to propose vines copulas.

Step One would be to find a right source to devise the code structure necessary to iplement such models.

Step Two would be to implement the class and the samplers.

Step Three would be to propose an estimation procedure, if possible generic enough :

• maybe one estimation procedure to estimate the structure itself, see Uytendale or some name like that%.
• and then another one, given the structure, to estimate each of the conditional copulas ? BTW, w need ot be able to devise conditional copulas for that to work correctly.

[Bug] Empirical Copula

Hi, I wanted to ask if the Copula empirica pdf and cdf is working correctly or not yet.

i used

emp = EmpiricalCopula(u) den = pdf(emp, [0.5, 0.5])

ERROR: MethodError: no method matching _logpdf(::EmpiricalCopula{2, Adjoint{Float32, Matrix{Float32}}}, ::Vector{Float64})

emp = EmpiricalCopula(u) den = cdf(emp, [0.5, 0.5])

ERROR: UndefVarError: mean not defined

[Contrib guideline] Add a sentence mentioning the colprac guide.

Link to JOSS review : openjournals/joss-reviews#6189

The contribution guidelines section is very short, but I think this is enough. Since you have the colprac badge, you could add a sentence there saying that more info can be found in the ColPrac guide.

[Fitting] Implement more fitting procedures

Some of the models are still lacking estimation procedures at all.

All models that have a density should be able to provide a MLE estimation procedure to perform IFM via SklarDist. For the moment it is not the case.

Some very used models could provide gradieents of the LLh directly to increase performance of the MLE.

just thoughts.

Archimedean sampling through McNeil & Neslehova's Radial/Simplex decompositions

Followup to the disucssion in #42

The goal is to switch for all / some / a few archimedean copulas (to be determined) from the frailty sampling that we currently have to the Williamson sampling from this article :

https://github.com/Article{mcneil2009multivariate,
title={Multivariate Archimedean copulas, d-monotone functions and ℓ 1-norm symmetric distributions},
author={McNeil, Alexander J and Ne{\v{s}}lehov{'a}, Johanna},
year={2009}
}

It requires the use of Williamson D-transform (defined in the article) that I implemented there : lrnv/WilliamsonTransforms.jl. The package is currently pending for registry approval, but is basically working. The method of samplign derscribed in the paper, together with this implementation, would allow sampling from ALL archimedean copulas, while the current method is restricted to those with completely monotonous generators.

With this generic fallback, we could allow instantiation and sampling of an archimedean copula from any given generator. On the other hand, some Williamson d-transform of some generators are known and could be used directly, to make computations more efficients.

Maybe we need a class for generators or something to cary with them their monotonicity property.

We definitely need bindings to make that possible, and maybe change the name of the radial dist functions htat are already there... Dunno yet.

PS : it will also allow to sample from negatively dependent claytons, which is not possible yet.

[Benchmark] Speed comparison with Matlab & R

Is there some benchmark? I wonder how much speed improvement when it's compared with copulapdf function in Matlab.

[New copulas] Extreme values copulas

As noted by @lidiamandre in #50 , implementing the option to use extreme values copulas would be neat :)

We may take a look at the R evd package to get insights.

Error on non-correlation matrices

Instead of making eveything a correlation matrix, let's error on them, would be easier.

[Student copula] fitting

See there : https://jds-online.org/journal/JDS/article/1266/info

to be able to fit student copula.

pdf/logpdf for Archimedean copula

Looking at R copula library, the calculation of the pdf/logpdf is quite a pain.

If you don't mind the dependency, could use ForwardDiff.jl to get the derivatives of the generator, then everything is pretty easy. The main question is how to get the higher order derivatives, I see some discussion e.g. https://discourse.julialang.org/t/automatic-differentiation-for-higher-order-derivatives/62337.
https://discourse.julialang.org/t/forwarddiff-3rd-or-even-higher-order-derivatives/8011/6
The easy option is to just write

∂ϕ(cop::CT, t::Real) where {CT<:Copulas.ArchimedeanCopula} = ForwardDiff.derivative(x -> Copulas.ϕ(cop,x),t) ∂²ϕ(cop::CT, t::Real) where {CT<:Copulas.ArchimedeanCopula} = ForwardDiff.derivative(x -> ∂ϕ(cop,x), t)
etc. There might be a smarter way to do this, but I don't see anything obvious. Since, you only need the d'th and first derivative of the generator for the pdf, it may be fine for now, but there's probably some way to get the first derivative along the way to the d'th derivative.

The R library also stops after d=6 for Archimedean, so could follow that convention, since at some point Archimedean in high dimension is silly.

Add sampling possibility for negatively dependent claytons

This will allow better coverage of the method, e.g. negatively correlated claytons.

Add generic coherence check for archimedeans

From McNeil & Neslehova 2009, we have that the williamson d-transfrom has the same distribution of \sum\phi(u_i) where \phi is the generator and u_i are samples from the copula.

Maybe we could use this property as a test: sample from Wd transfrom on one side, and from the copula on the other and use a KS test to ensure that the distributions matches ? We could make that generic and test a lot of different cases of archimedean copulas :)

[Docs] Detail a bit more the Turing example (add plot).

Link to JOSS review : openjournals/joss-reviews#6189

Luca wrote:

In the turing integration example, the code ends with a plot, would be good to include the resulting image in the documentation.

If possible, some visual representation would also be nice to have for the example of the README

So first adding the image in the documentation should be easy, the only think is to tell Documenter.jl to actually compile this code (which takes a while, this is why i removed it the other time).

But then plotting a copula with Plots.jl should also be a possibility, a nice thing to have.

[Paper] Unicode rendering in code blocks is broken.

Link to JOSS review : openjournals/joss-reviews#6189

[Bug] Kendall's tau for TCopula & non-archimedeans more generally

Hello, when I try to calculate the kendal tau of any copula, it appears that it is not defined
`julia> Σ = [1.0 0.8; 0.8 1.0]
2×2 Matrix{Float64}:
1.0 0.8
0.8 1.0

julia> Cop = GaussianCopula(Σ)
GaussianCopula{2, Matrix{Float64}}(
Σ: [1.0 0.8; 0.8 1.0]
)

julia> τ(Cop)
ERROR: UndefVarError: τ not defined
Stacktrace:
[1] top-level scope
@ REPL[7]:1`

[API choices] `fit_mle` for insteances, `_tau`

function Distributions.fit(D::SklarDist, x)
    # The first thing to do is to fit the marginals : 
    @assert length(D) == size(x, 1)
    # One could put fit but EM works for with fit_mle only
    m = Tuple(Distributions.fit_mle(D.m[i], x[i, :]) for i in axes(x, 1))
    u = pseudos(x)

    # τ⁻¹(FrankCopula, τ) works but not τ⁻¹(FrankCopula{3, Float64}, τ)
    # I found https://discourse.julialang.org/t/extract-type-name-only-from-parametric-type/14188/20
    # to just extract the typename.
    # Again defining fit(C::ArchimedeanCopula, u) = fit(::Type{CT},u) where {CT <: ArchimedeanCopula} would simplify the writting
    C = Distributions.fit(Base.typename(typeof(D.C)).wrapper, u)
    return SklarDist(C, m)
end

D̂mix = fit(D, simu)

Note that I had to define the weighted fit_mle(LogNormal, x, w) (only fit_mle(LogNormal, x) is in Distributions.jl).
I believe this does the job

using StatsBase
function Distributions.fit_mle(::Type{<:LogNormal}, x::AbstractArray{T}, w::AbstractArray) where T<:Real
    lx = log.(x)
    μ, σ = mean_and_std(lx, weights(w))
    LogNormal(μ, σ)
end

(Note that StochasticEM() version only use unweighted version of fit_mle.)

Hence, IMO, it seems that for compatibility between EM.jl and Copula.jl, it is Copula.jl that needs to simply add Distributions.fit(D::SklarDist, x) as shown here. I believe it should not have any consequences on your existing functions.

Originally posted by @dmetivie in dmetivie/ExpectationMaximization.jl#8 (comment)

Difficulty generating independent copula samples and also evaluating density and cumulative functions for Fréchet–Hoeffding copula bounds

Hello, nice to talk to you.

I have a problem when I want to generate random samples and also to calculate the "pdf" and the "cdf" of the independent copula and the Fréchet–Hoeffding bounds. I know that densities for the bounds do not exist, however the cumulative function should be able to be calculated. In the first code I show a small example where it works well for the Gaussian copula. It involves generating a 3D graph of the density and cumulative function of the copula. To change between density and cumulative function it is only necessary to comment and uncomment

####################################################################

using Copulas, Distributions, Plotly

###Define a valid correlation matrix for the bivariate Gaussian copula
Σ = [1.0 0.5; 0.5 1.0] # example

###Create the bivariate Gaussian copula with the correlation matrix Σ
G = GaussianCopula(Σ)
#G = IndependentCopula(2)
#G = MCopula(2)
###Generate a grid of two-dimensional values (u, v) in [0, 1]x[0, 1]
u = range(0, stop=1, length=100)
v = range(0, stop=1, length=100)
grid = [(ui, vi) for ui in u, vi in v]
###Calculate the cumulative function or density of the copula in the grid matrix
density = zeros(length(grid))
###acumulate = zeros(length(grid))
for (i, (ui, vi)) in enumerate(grid)
###acumulate[i] = cdf(G, [ui, vi])
density[i] = pdf(G, [ui, vi])
end

####Reformat the density into a two-dimensional array
z_data = reshape(density, 100, 100)

layout = Layout(
title="Gaussian Copula",
autosize=false,
scene_camera_eye=attr(x=1.87, y=0.88, z=-0.64),
width=500, height=500,
margin=attr(l=65, r=50, b=65, t=90)
)

plot(surface(
z=z_data,
contours_z=attr(
show=true,
usecolormap=true,
highlightcolor="limegreen",
project_z=true
)
), layout)

###############################################################################

The second code is a little more basic, it simply consists of simulating random samples of the independent copula and the Fréchet–Hoeffding bounds. Now, as you can notice, when simulating the Gaussian copula I have no problems, however for the independent copula I get the error
" MethodError: no method matching Random.Sampler(::Type{Random.TaskLocalRNG}, :Random.SamplerType{IndependentCopula{2}}, ::Val{1})"

##########################

using Copulas, Distributions

###Gaussian copula simulation
Σ = [1.0 0.5; 0.5 1.0]
G = GaussianCopula(Σ)
simu1 = rand(G,1)

###independent copula simulation
I = IndependentCopula(2)
simu2 = rand(I,1) #Does not generate samples and gives the mentioned error

###Mcopula simulation
M = MCopula(2)
simu3 = rand(M,1)

###WCopula simulation
W = Copulas.WCopula(2)
simu4 = rand(W,1)

###########################

The code does not generate samples of the independent copula, nor does it generate samples of "WCopula(2)" since an error of the type "UndefVarError: WCopula not defined" appears. I could solve it by changing "WCopula(2) for Copulas.WCopula(2) )" . However, the error of independent copulation continues. For "MCopula()" and the other elliptic and Archimedean copulas it works well.

Thank you very much for your attention and I look forward to your prompt response.

Elliptical CDF

Gaussian and student have no cdf.

from discourse : https://discourse.julialang.org/t/cdf-of-a-copula-from-copulas-jl/85786

I fixed Copulas.jl (in the latest master commit), so that it provides a more accurate error: we rely on Distributions.jl's implementation of the multivariate Gaussian and Student distributions to handle the complexity behind the Gaussian and Student copulas, and Distributions.jl has no multivariate Gaussian cdf.
Here is a Github issue for this : JuliaStats/Distributions.jl#260
Also, seems relevant there : JuliaStats/StatsFuns.jl#114
But it does not look like it is going to be fixed soon. If you find somewhere an integration routine that computes the cdf of a multivariate Gaussian, I am willing to include it as a temporary workaround though, until those issues are fixed in StatsFun.jl and Distributions.jl

Move archimedean generators to their own structs.

The idea is to have

abstract type Generator end
max_monotony(Generator) = ...
\phi(Generator) = ....
\phi_inv(Generator) = ... # with an adequate default method
williamson_dist(Generator,d) = ... # with an adequate default method. 

then

struct WilliamsonGenerator <: Generator
    X::TX # the random variable 
    d::Int # the dimension in which the transformation has been done. 
end
struct ClaytonGenerator <: Generator
    theta::T
end

with the appropriate methods overloads.

And then we can construct copulas as

# for archimedeans: 
ClaytonCopula{d,T} = ArchimedeanCopula{d,ClatyonGenerator{T}}
# or even 
struct LiouvilleCopula{d,V,G}
    alpha::V # vector of dirichlet parameters
    phi::G # archimedean generator
end

it is still unclear how the other methods (fitting, tau, invtau, rho, invrho) will be implemented for specific archimedeans..

It will also easily allow for nested archimedeans, by defining e.g.

can_nest(inner::Generator, outer::Generator) = ...
struct NestedLiouvilleCopula
    phi::Generator
    content::Vector{Union{NestedLiouvilleCopula,LiouvilleCopula,Distributions.Uniform}}
end

that recursively checks possibility by using the nesting conditions through e.g.

can_nest(inner::ClaytonGenerator, outer::ClaytonGenerator) = inner.theta < outer.theta
can_nest(inner::ClaytonGenerator,outer::GumbelGenerator) = # more complicated but doable

maybe even a generic can_nest could be implemented ? I dunno if it is even possible.

So maybe ArchimedeanCopula<:LiouvilleCopula should be true ?

[Internals] move to FrailtyGenerator

This will allow to hide the ugly WilliamsonFromFrailty, with the same kind of interface as WilliamsonGenerator.

Moreover, the \phi function can be directly taken as the laplace transfrom of the random variable, I think that Distributions.jl has laplace transfrom of variables right ? Since laplace transfrom of positive random variables are always completely monotoonous, that would make a lot of sense !

problem with Copula Frank samples

Hi, I had a very strange problem, when I want to sample the two-dimensional Frank copula I get infinite values.

F = FrankCopula(2,98.11)
r = rand(F,1)
Just with one value, it is giving infinity

I also tried with parameters like 86.44 and 75.

TagBot trigger issue

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

[Docs] Convert docstring references to the right format

The right format is :

* [bibtexkey](@cite) My super reference (Year), jounral XXX...

Lack of documentation on PlackettCopula

@Santymax98 I did not check on the pull-request, but you did not write the documentation of the new copula type as the other ones (granted, some of them are missing).

Would you like to do it yourself ? Copy the format of the documentation from some other copula model, and include in it your references.

Or do you prefer me to do it ?

example on front page gives error

(using Julia 1.8.5)

Not sure what's causing this, but this is the output from the fit() command in the example:

 D̂ = fit(SklarDist{FrankCopula,Tuple{Gamma,Normal,LogNormal}}, simu)
┌ Warning: GSL Error 1 (input domain error): domain error at debye.c:202
└ @ GSL.C C:\Users\ethan\.julia\packages\GSL\T2BEO\src\error_handling.jl:23
ERROR: DomainError with [unknown value]:
domain error at debye.c:202
Stacktrace:
  [1] custom_error_handler(reason::String, file::String, line::Int32, errno::Int32)
    @ GSL.C C:\Users\ethan\.julia\packages\GSL\T2BEO\src\error_handling.jl:26
  [2] custom_error_handler(reason::Ptr{UInt8}, file::Ptr{UInt8}, line::Int32, errno::Int32)
    @ GSL.C C:\Users\ethan\.julia\packages\GSL\T2BEO\src\error_handling.jl:9
  [3] sf_debye_1(x::Float64)
    @ GSL.C C:\Users\ethan\.julia\packages\GSL\T2BEO\src\gen\direct_wrappers\gsl_sf_debye_h.jl:35
  [4] (::Copulas.var"#19#20"{Float64})(x::Float64)
    @ Copulas C:\Users\ethan\.julia\packages\Copulas\NGFd0\src\ArchimedeanCopulas\FrankCopula.jl:39
  [5] (::Roots.FnWrapper)(x::Float64)
    @ Roots C:\Users\ethan\.julia\packages\Roots\qr9ed\src\alternative_interfaces.jl:98
  [6] (::Roots.Callable_Function{Val{1}, Val{false}, Roots.FnWrapper, Nothing})(x::Float64)
    @ Roots C:\Users\ethan\.julia\packages\Roots\qr9ed\src\functions.jl:28
  [7] update_state(::Roots.Secant, F::Roots.Callable_Function{Val{1}, Val{false}, Roots.FnWrapper, Nothing}, o::Roots.UnivariateZeroState{Float64, Float64}, options::Roots.UnivariateZeroOptions{Float64, Float64, Float64, Float64}, l::Roots.NullTracks)
    @ Roots C:\Users\ethan\.julia\packages\Roots\qr9ed\src\DerivativeFree\secant.jl:39
  [8] update_state(::Roots.Secant, F::Roots.Callable_Function{Val{1}, Val{false}, Roots.FnWrapper, Nothing}, o::Roots.UnivariateZeroState{Float64, Float64}, options::Roots.UnivariateZeroOptions{Float64, Float64, Float64, Float64})
    @ Roots C:\Users\ethan\.julia\packages\Roots\qr9ed\src\DerivativeFree\secant.jl:29
  [9] solve!(𝐙 ::Roots.ZeroProblemIterator{Roots.Secant, Roots.AlefeldPotraShi, Roots.Callable_Function{Val{1}, Val{false}, Roots.FnWrapper, Nothing}, Roots.UnivariateZeroState{Float64, Float64}, Roots.UnivariateZeroOptions{Float64, Float64, Float64, Float64}, Roots.NullTracks}; verbose::Bool)
    @ Roots C:\Users\ethan\.julia\packages\Roots\qr9ed\src\hybrid.jl:54
 [10] solve(𝑭𝑿  ::Roots.ZeroProblem{Roots.FnWrapper, Float64}, M::Roots.Order0, p::Nothing; verbose::Bool, kwargs::Base.Pairs{Symbol, Roots.NullTracks, Tuple{Symbol}, NamedTuple{(:tracks,), Tuple{Roots.NullTracks}}})
    @ Roots C:\Users\ethan\.julia\packages\Roots\qr9ed\src\find_zero.jl:493
 [11] find_zero(f::Roots.FnWrapper, x0::Float64, M::Roots.Order0, p′::Nothing; p::Nothing, verbose::Bool,
 tracks::Roots.NullTracks, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ Roots C:\Users\ethan\.julia\packages\Roots\qr9ed\src\find_zero.jl:221
 [12] find_zero (repeats 2 times)
    @ C:\Users\ethan\.julia\packages\Roots\qr9ed\src\find_zero.jl:210 [inlined]
 [13] derivative_free(f::Function, x0::Float64; order::Int64, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ Roots C:\Users\ethan\.julia\packages\Roots\qr9ed\src\alternative_interfaces.jl:243
 [14] derivative_free(f::Function, x0::Float64)
    @ Roots C:\Users\ethan\.julia\packages\Roots\qr9ed\src\alternative_interfaces.jl:222
 [15] #fzero#133
    @ C:\Users\ethan\.julia\packages\Roots\qr9ed\src\alternative_interfaces.jl:157 [inlined]
 [16] fzero
    @ C:\Users\ethan\.julia\packages\Roots\qr9ed\src\alternative_interfaces.jl:154 [inlined]
 [17] τ⁻¹
    @ C:\Users\ethan\.julia\packages\Copulas\NGFd0\src\ArchimedeanCopulas\FrankCopula.jl:39 [inlined]
 [18] fit(#unused#::Type{FrankCopula}, u::Transpose{Float64, Matrix{Float64}})
    @ Copulas C:\Users\ethan\.julia\packages\Copulas\NGFd0\src\ArchimedeanCopula.jl:72
 [19] fit(#unused#::Type{SklarDist{FrankCopula, Tuple{Gamma, Normal, LogNormal}}}, x::Matrix{Float64})
    @ Copulas C:\Users\ethan\.julia\packages\Copulas\NGFd0\src\SklarDist.jl:36
 [20] top-level scope
    @ REPL[271]:1

What could be causing this? I don't understand the error message. Any help appreciated, and thanks for making this awesome functionality available.

Doubt about the Williamson transform

Dear Oskar

I hope everething ok.

I wanted to ask you something regarding your implementation of the Williamson transform. Is it possible to obtain the accumulated function "quantile()" for the Radial function of the Archimedean family provided by the implementation?

[Bug] RafteryCopula broken tests

@Santymax98 There are broken tests in Rafterycopula, as when we merged your PR we forgot to run them...

On the pdf here

On the cdf there

This is problematic, I had the pdf retuninrg negative values on some test cases. Would you mind checking again the coirrectness of your code ?

Cheers

[Bug] Brake ties at random

The pseudo function should break ties at random (in curent version it does not), while still using an efficient O(n logn) algorithm. See https://github.com/oheil/NormalizeQuantiles.jl/blob/master/src/NormalizeQuantiles.jl for an implementation of sampleranks that might work for us.

Design automatic soundness test for all copulas.

The idea is that a few basic checks should exists and be the same for every copulas. In these tests, we might include :

• Marginal uniformity, to ensure that the model is indeed a copula. This can include verifying that the cdf has uniform margins AND that marginals samples are uniformely distributed (so a check on the cdf and a check on the samples). We alreayd have this one but we could generlize it a bit more and make it automatic

• check mandatory values of pdf&cdf on the frontier of the hypercube

• check that cdf is in [0,1] and that pdf is positive on some random points

• check numerically that the pdf is indeed the derivative of the cdf indeed

• check that tau \circ tau_inv is the identity on the domain

• check that rho \circ rho_inv is identitdy in the domain

• check that fitting the model works under a sklardist construct

• Archimedean have specific tests themselves.

• Generators could be checked numerically to be indeed monotone enough (compute theoretically the derivatives up to d-2 for exemple and check signs).

Some of these checks need to filter the models on the ones that implements a certain method... this looks complicated to do and i do not really want to do it by hand, but it is the right thing to do i think. Other ideas on stuff that could be automaticcally checked ?

Bug in simulation algorithm

Simple duplication of the issue noted at org-arl/AlphaStableDistributions.jl#46

Maybe we could simply depend on AlphaStableDistributions.jl again when they fix their stuff?

add Elliptical CDF using MvNormalCDF.jl or something similar

Copulas.jl does not compute cdf values for elliptical copulas. Per the discussion here ( https://discourse.julialang.org/t/cdf-of-a-copula-from-copulas-jl/85786/6 ), using MvNormalCDF.jl (or something similar) to calculate the CDF may be one way to go. Thanks for @PharmCat 's volunteer.

[Paper] JOSS review checklist

This started as a copy of Luca's and Ander's review but now is our checklist of things to do before publishing the paper.

Link to JOSS review : openjournals/joss-reviews#6189

Paper

• #115
• #120
• #121
• #120
• #121
• #116
• #135
• #139
• #140
• #141

Bugs

• #11
• #50
• #105
• #80
• #98

Documentation

• #88
• #118
• #119
• #145
• #143
• #142
• #117

Features

• #150
• #151
• #152

TODO at the end of the review

• #129
• #163

[Paper] Add reference to `DataGenCopulaBased.jl` and `BivariateCopulas.jl`, both have zenodo DOIs.

Link to JOSS review : openjournals/joss-reviews#6189

Add a precompile workflow

We should use a precompile workflow that exercise some of the copulas and most of the usefull functions -- pdf, cdf, fitting, etc.. Maybe based on the test workflows ?

[Interface] Clarify properties of several copulas

The following properties of the copulas:

• kendall tau
• spearman rho
• tail dependence coefficients

should be clarified and included into the general API. it is weird that we dont have them yet for all copulas. For the moment, only several of them implements kendall tau and spearman rho, and the tail coefficients are lacking.

[New copulas] Fréchet-Hoeffding mixture family

So that W and M are in the same parametric family, ref AnderGray/ProbabilityBoundsAnalysis.jl#42

Inverse Gaussian Copula

Hi!!,

In the Gaussian inverse copula when the parameter is zero it should return the independent copula (special case) as the verification is below the requirement that the parameter must be positive it is removing the opportunity to return the independent copula when the parameter is equal to zero

[Bug] Frank Copula when parameter is negative

I was using this package for the first time as I need to simulate from some copulas in Julia rather than in R. The Frank copula is defined for all the real values but 0 but here it doesn't allow me to generate when my parameter is negative (and hence when I have negative dependence in the data), could this be fixed please, many thanks :)

A good addition to this package would be to have more Extreme Value copulas, such as the ones we can generate from the evd package in R

Accept interval parameters

Every copula family should accept an interval parameter, and return for e.g. Distributions.cdf an interval. Dunno if this will work out of the box but should be tested. Ref AnderGray/ProbabilityBoundsAnalysis.jl#42

Some Questions

If some of the questions seem ill-informed please let me know. I'm exploring the territory as I experiment with EconomicScenarioGenerators.jl (specifically, working on integrating Copulas.jl in JuliaActuary/EconomicScenarioGenerators.jl#39

• when sampling from the copulas with rand do all the copulas return a value in [0,1] representing the cumulative probability to that point?
• Do you know of more convenient ways to construct a symmetric covariance matrix than to use the whole thing?

Gumbel with parameter > 10

Hello, I noticed that the Gumbel copula does not allow obtaining samples larger than copulas with a parameter greater than 10, since the \alpha - stable distribution has a limitation of 0.1 to 2 and I think it should be 0 < \alpha < 1

Originally posted by @Santymax98 in #66 (comment)

H-volume

As a follow up to AnderGray/ProbabilityBoundsAnalysis.jl#42

Regarding the H-volume, yes, I think this should go to Copulas.jl : i'll add a generic there. I also have a R implementation in https://github.com/lrnv/cort which is based on the follinwg paper :

Cherubini U, Romagnoli S (2009). “Computing the Volume of n -Dimensional Copulas.” Applied Mathematical Finance, 16(4), 307--314. ISSN 1350-486X, 1466-4313.

https://github.com/lrnv/cort/blob/b5d0d7a0d7d214ab0fd7a167dc99fa2d4cbe8cac/R/generics.R#L121-L127

[Paper] Typos to fix + Eventually automatic proof reading ?

Link to JOSS review : openjournals/joss-reviews#6189

I spotted a few typos, you may consider using a tool like to check more systematically. This could also be setup on CI (actually the julia repo does it) but it might be an overkill.

• line 40 of the paper must -> most
• line 41 of the paper complience -> compliance
• line 46 of the paper specify -> specifies
• line 55 folloiwng -> following (I also think consists of is the grammatically correct form, instead of is consisting of but I am not native)
• line 63 exploit -> exploits
• Documentation navbar Exemples -> Examples
• here in the "Unfinished work" box, wtill -> still

[Documentation] Automatized proof reading the docs & docstrings ?

We could use something like that to proof-read the docs again : https://forem.julialang.org/danielvandh/efficiently-fixing-typos-in-documentation-1j88

MvNormalCDF 0.3 compatible

If it possible, please check MvNormalCDF.jl v0.3.0 compatible (on dev brunch).

PharmCat/MvNormalCDF.jl#9

Planning changes in types to make compatible with ForwardDiff.jl

[Paper] Add comparison table with BivariateCopulas.jl and DatagenCopulaBased.jl

Link to JOSS review : openjournals/joss-reviews#6189

Lucas wrote:

At line 40 a pretty strong statement is made

Copulas.jl is clearly the must fully featured, and brings, as a key feature, the complience with the broader ecosystem.

It is not clear to me from the paper alone that Copulas.jl is the most feature complete of the three and would be good to give some more arguments in the paper.

Ideally, a table comparing features of the three would be nice to have.

I agree that a simple table comparing key functionalities would be usefull. We have to decide what to include in that table.

[Paper] Publish a new version before JOSS paper with all changes from the review.

Do not forget to do that when reviewers comment are all dealt with

Fix broken kendall tau on GumbelBarnettCopula and InvGaussianCopula

@Santymax98 This is yours to do. the gsl_... functions you wrote do not have the right interface for Julia, where the funtions are called GSL.integration_qags() intead, but the parameters are not the same as what you have ? where is this gsl_integration_qags() call coming from ?

Fix type parameters of each copula consistently

For the moment type parameters of copulas are whatever.

It would be very nice to have at least, on every copula :

• its dimension d
• the types of values of generates (Float64, BigFloat, etc... as T

The mothertype could become Copula{d,T} instead of Copula{d} that is is right now.

[Docs] Find new applications examples for the docs

Link to JOSS review : openjournals/joss-reviews#6189

Luca wrote:

Currently the documentation feels to have more theory than code, it would be nice to expand the Examples section with a few more examples of how to use the code, ideally within some applications of how to solve some "real-world" problems. E.g. applying Copulas.jl to some popular datasets and concrete example problems.

Thus we need to find a use-case ona popular dataset or concrete problem to add it to the docs.

@Santymax98 Maybe you have in your work one simple use-case that could be shown ?