Copulas.jl

Appearance

102

Fitting compound distributions ​

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

julia
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.0784939  0.997003  0.604455  …   1.82822   1.2783     0.446217
 76.4781     1.39965   3.77571       1.5642    4.16744   31.9205
  0.122785   0.634351  0.84995       2.41706   1.46092    0.0589099
  0.694961   0.331649  0.998912     -0.209601  0.263081   1.21281

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

julia
MyCop = SurvivalCopula{4,ClaytonCopula,(2,4)}
MyMargs = Tuple{LogNormal,Pareto,Gamma,Normal}
MyD = SklarDist{MyCop, MyMargs}
fitted_model = fit(MyD,data)
SklarDist{SurvivalCopula{4, ClaytonCopula{4, Float64}, (2, 4)}, Tuple{Distributions.LogNormal{Float64}, Distributions.Pareto{Float64}, Distributions.Gamma{Float64}, Distributions.Normal{Float64}}}(
C: SurvivalCopula(ClaytonCopula{4, Float64}(θ = 2.6778664734677733,))
m: (Distributions.LogNormal{Float64}(μ=-0.011822157090623123, σ=1.003768975981645), Distributions.Pareto{Float64}(α=1.0066863250981448, θ=1.001945355908185), Distributions.Gamma{Float64}(α=1.0117051489414404, θ=0.9176748986655682), Distributions.Normal{Float64}(μ=0.04878310534814991, σ=1.019810914820893))
)

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

julia
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.011822157090623123, σ=1.003768975981645), Distributions.Pareto{Float64}(α=1.0066863250981448, θ=1.001945355908185), Distributions.Gamma{Float64}(α=1.0117051489414404, θ=0.9176748986655682), Distributions.Normal{Float64}(μ=0.04878310534814991, σ=1.019810914820893))
)

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. For example, Archimedean copulas are fitted through inversion of the Kendall tau function, while the Gaussian copula is fitted by maximum likelihood.

Scatter of original data (first two dims) ​

julia
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 ​

julia
U = Copulas.pseudos(data)              # pseudo-observations (uniforms)
Usim = rand(fitted_model.C, size(data,2))  # simulate same length 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))