First draft of william:son implementation

Project.toml

@@ -14,19 +14,21 @@ Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 TaylorSeries = "6aa5eb33-94cf-58f4-a9d0-e4b2c4fc25ea"
+WilliamsonTransforms = "48feb556-9bdd-43a2-8e10-96100ec25e22"
 
 [compat]
 Cubature = "1.5"
 Distributions = "0.25"
 ForwardDiff = "0.10"
 GSL = "1"
 MvNormalCDF = "0.2, 0.3"
+Random = "1.6"
 Roots = "1, 2"
 SpecialFunctions = "2"
 StatsBase = "0.33, 0.34"
 TaylorSeries = "0.12, 0.13, 0.14, 0.15"
+WilliamsonTransforms = "0.0"
 julia = "1.6"
-Random = "1.6"
 
 [extras]
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

src/ArchimedeanCopulas/WilliamsonCopula.jl

@@ -0,0 +1,28 @@
+struct WilliamsonCopula{d,Tϕ,TX} <: ArchimedeanCopula{d}
+    ϕ::Tϕ
+    X::TX
+end
+function WilliamsonCopula(X::Distributions.UnivariateDistribution, d)
+    ϕ = WilliamsonTransforms.𝒲(X,d)
+    return WilliamsonCopula{d,typeof(ϕ),typeof(X)}(ϕ,X)
+end
+function WilliamsonCopula(ϕ::Function, d)
+    X = WilliamsonTransforms.𝒲₋₁(ϕ,d)
+    return WilliamsonCopula{d,typeof(ϕ),typeof(X)}(ϕ,X)
+end
+function WilliamsonCopula(ϕ::Function, X::Distributions.UnivariateDistribution, d)
+    return WilliamsonCopula{d,typeof(ϕ),typeof(X)}(ϕ,X)
+end
+function Distributions._rand!(rng::Distributions.AbstractRNG, C::CT, x::AbstractVector{T}) where {T<:Real, CT<:WilliamsonCopula}
+    r = rand(rng,C.X)
+    Random.rand!(rng,x)
+    for i in 1:length(C)
+        x[i] = -log(x[i])
+    end
+    sx = sum(x)
+    for i in 1:length(C)
+        x[i] = ϕ(C,r * x[i]/sx)
+    end
+    return x
+end
+ϕ(  C::WilliamsonCopula,      t) = C.ϕ(t)

src/Copulas.jl

@@ -18,6 +18,7 @@ end Copulas
     import ForwardDiff
     import Cubature
     import MvNormalCDF
+    import WilliamsonTransforms
 
     # Standard copulas and stuff. 
     include("utils.jl")
@@ -59,12 +60,13 @@ end Copulas
     include("ArchimedeanCopulas/GumbelCopula.jl")
     include("ArchimedeanCopulas/FrankCopula.jl")
     include("ArchimedeanCopulas/AMHCopula.jl")
+    include("ArchimedeanCopulas/WilliamsonCopula.jl")
     export IndependentCopula, 
            ClaytonCopula,
            JoeCopula,
            GumbelCopula,
            FrankCopula,
-           AMHCopula
-
+           AMHCopula,
+           WilliamsonCopula
 
 end

test/williamson_test.jl

@@ -0,0 +1,18 @@
+
+@testitem "williamson test" begin
+    using Distributions, Random
+    Random.seed!(123)
+    taus = [0.0, 0.1, 0.5, 0.9, 1.0]
+
+    ϕ_clayton(x, θ) = max((1 + θ * x),zero(x))^(-1/θ)
+
+    Cops = (
+        WilliamsonCopula(Dirac(1),10),
+        WilliamsonCopula(x -> exp(-x),10),
+        WilliamsonCopula(x -> ϕ_clayton(x,2),2),
+        WilliamsonCopula(x -> ϕ_clayton(x,-0.3),2)
+    )
+    for C in Cops
+        x = rand(C,1000)
+    end
+end