Sample from a clayton copula

clayton.jl

using Distributions
using SpecialFunctions
using Plots

"""
Generator
"""
function φ(x::Real, ϑ::Real)
    return (1+x)^(−1/ϑ)
end

"""
Inverse Generator
"""
function φ⁻(x::Real, ϑ::Real)
    return x^(-ϑ) - 1
end

"""
First generator derivative
"""
function φ¹(x::Real, ϑ::Real)
    α = 1/ϑ
    return gamma(1 + α)/gamma(α) * (1 + x)^(-(1 + α))
end

function sample(ϑ::Real, n::Int64)
    Eₖ = rand(Exponential(1), 2, n)
    M = rand(Gamma(1/ϑ, 1), 1, n)

    return φ.(Eₖ ./ M, ϑ)
end

function rosenblatt(M::AbstractMatrix, ϑ::Real)
    u1 = M[1, :]
    u2 = M[2, :]

    return [u1'; (φ¹.(φ⁻.(u1, ϑ) + φ⁻.(u2, ϑ), ϑ) ./ φ¹.(φ⁻.(u1, ϑ), ϑ))']
end

function inverse_rosenblatt(U::AbstractMatrix, ϑ::Real)
    u1 = U[1, :]
    u2 = U[2, :]

    return [u1'; ((u1.^(-ϑ) .* (u2.^(-ϑ/(ϑ + 1)) .- 1) .+ 1).^(-1/ϑ))']
end

ϑ = 5
M = sample(ϑ, 2500)
U = rosenblatt(M, ϑ)
M2 = inverse_rosenblatt(U, ϑ)

p1 = scatter(M[1, :], M[2, :], aspect_ratio = :equal, xlims=[0, 1], ylims=[0, 1], xlabel="Clayton", label=:none)
p2 = scatter(U[1, :], U[2, :], aspect_ratio = :equal, xlims=[0, 1], ylims=[0, 1], xlabel="Rosenblatt", label=:none)
p3 = scatter(M2[1, :], M2[2, :], aspect_ratio = :equal, xlims=[0, 1], ylims=[0, 1], xlabel="Inverse Rosenblatt", label=:none)

p = plot(p1, p2, p3)

savefig(p, "rosenblatt.png")