Fitting compound distributions

Through the SklarDist interface, it is possible to fit distributions constructed from a copula and marginals:

using Copulas
using Distributions

# Let's sample some datas:
X₁ = LogNormal()
X₂ = Pareto()
X₃ = Gamma()
X₄ = Normal()
C = SurvivalCopula(FrankCopula(4,7),(2,4))
D = SklarDist(C,(X₁,X₂,X₃,X₄))
data = rand(D,1000)

4×1000 Matrix{Float64}:
  0.401445  1.21509    7.81948  0.270553  …  1.64186   10.8205    0.294878
  4.60104   3.07032    1.26845  4.02309      3.69228    1.32904   5.32843
  0.52025   0.261444   2.83385  0.257615     0.020718   1.12437   0.31552
 -0.133422  0.598604  -1.26938  0.122044     1.25908   -0.618249  1.43523

The fit function uses a type as its first argument that describes the structure of the model :

MyCop = SurvivalCopula{4,ClaytonCopula}
MyMargs = Tuple{LogNormal,Pareto,Gamma,Normal}
MyD = SklarDist{MyCop, MyMargs}
fitted_model = fit(MyD, data; copula_method=:itau)

SklarDist{SurvivalCopula{4, ClaytonCopula{4, Float64}}, Tuple{Distributions.LogNormal{Float64}, Distributions.Pareto{Float64}, Distributions.Gamma{Float64}, Distributions.Normal{Float64}}}(
C: SurvivalCopula(ClaytonCopula{4, Float64}(0.3363188691771293,), (1, 2, 3, 4))
m: (Distributions.LogNormal{Float64}(μ=0.015077846825769316, σ=0.9744304178882592), Distributions.Pareto{Float64}(α=0.9907166680945756, θ=1.0004776976095284), Distributions.Gamma{Float64}(α=1.052760147593781, θ=0.9441034835961868), Distributions.Normal{Float64}(μ=0.021221643678072166, σ=0.9706193248521299))
)

Another possibility is to use an empirical copula and only fit the marginals:

other_fitted_model = fit(SklarDist{EmpiricalCopula,MyMargs},data)

SklarDist{EmpiricalCopula{4, Matrix{Float64}}, Tuple{Distributions.LogNormal{Float64}, Distributions.Pareto{Float64}, Distributions.Gamma{Float64}, Distributions.Normal{Float64}}}(
C: EmpiricalCopula{d}(4, 1000)
m: (Distributions.LogNormal{Float64}(μ=0.015077846825769316, σ=0.9744304178882592), Distributions.Pareto{Float64}(α=0.9907166680945756, θ=1.0004776976095284), Distributions.Gamma{Float64}(α=1.052760147593781, θ=0.9441034835961868), Distributions.Normal{Float64}(μ=0.021221643678072166, σ=0.9706193248521299))
)

This simple interface leverages the fit function from Distributions.jl. According to their documentation, this function is not supposed to use a particular method but to fit "quick and dirty" some distributions.

So you have to be careful: the fit method might not be the same for different copulas or different marginals. Copula families use maximum likelihood by default when they support it; here copula_method=:itau explicitly requests Kendall's-tau inversion, which is also more appropriate for this short, deterministic documentation example. Marginal estimators remain those provided by each distribution family in Distributions.jl.

Scatter of original data (first two dims)

using Plots
scatter(data[1,:], data[2,:]; ms=2, alpha=0.6, title="First two marginals (original scale)", legend=false)

Pseudo-observations vs simulated from fitted copula

U = Copulas.pseudos(data)              # pseudo-observations (uniforms)
Usim = rand(first(params(fitted_model)), size(data,2)) # simulate from fitted copula
P1 = scatter(U[1,:], U[2,:]; ms=2, alpha=0.6, xlim=(0,1), ylim=(0,1), title="Empirical uniforms", legend=false)
P2 = scatter(Usim[1,:], Usim[2,:]; ms=2, alpha=0.6, xlim=(0,1), ylim=(0,1), title="Fitted copula uniforms", legend=false)
plot(P1, P2; layout=(1,2), size=(850,350))