Archimedeans

Project.toml

@@ -8,9 +8,9 @@ Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
 Cubature = "667455a9-e2ce-5579-9412-b964f529a492"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
-GSL = "92c85e6c-cbff-5e0c-80f7-495c94daaecd"
 LogExpFunctions = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
 MvNormalCDF = "37188c8d-bc69-4638-b057-733e744175ec"
+QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
@@ -24,12 +24,12 @@ Combinatorics = "1"
 Cubature = "1.5"
 Distributions = "0.25"
 ForwardDiff = "0.10"
-GSL = "1"
 HypothesisTests = "v0.11"
 InteractiveUtils = "1.6"
 LinearAlgebra = "1.6"
 LogExpFunctions = "0.3"
 MvNormalCDF = "0.2, 0.3"
+QuadGK = "2"
 Random = "1.6"
 Roots = "1, 2"
 SpecialFunctions = "2"

src/ArchimedeanCopula.jl

@@ -89,7 +89,6 @@ end
 
 williamson_dist(C::ArchimedeanCopula{d}) where d = WilliamsonTransforms.𝒲₋₁(t -> ϕ(C,t),d)
 function τ(C::ArchimedeanCopula)  
-    @show C
     return 4*Distributions.expectation(r -> ϕ(C,r), williamson_dist(C)) - 1
 end
 

src/ArchimedeanCopulas/FrankCopula.jl

@@ -44,10 +44,10 @@ end
 # Avoid type piracy by defiing it myself: 
 ϕ(  C::FrankCopula,       t::TaylorSeries.Taylor1) = C.θ > 0 ? -log(-expm1(LogExpFunctions.log1mexp(-C.θ)-t))/C.θ : -log1p(exp(-t) * expm1(-C.θ))/C.θ
 
-D₁ = GSL.sf_debye_1 # sadly, this is C code.
+# D₁ = GSL.sf_debye_1 # sadly, this is C code.
 # could be replaced by : 
 # using QuadGK
-# D₁(x) = quadgk(t -> t/(exp(t)-1), 0, x)[1]/x
+D₁(x) = QuadGK.quadgk(t -> t/(exp(t)-1), 0, x)[1]/x
 # to make it more general. but once gain, it requires changing the integrator at each evlauation, 
 # which is problematic. 
 # Better option is to try to include this function into SpecialFunctions.jl. 

src/ArchimedeanCopulas/GumbelBarnettCopula.jl

@@ -37,25 +37,25 @@ end
 ϕ(  C::GumbelBarnettCopula,       t) = exp((1-exp(t))/C.θ)
 ϕ⁻¹(C::GumbelBarnettCopula,       t) = log(1-C.θ*log(t))
 function τ(C::GumbelBarnettCopula)
-    # Define the function to integrate
-    f(x) = -x * (1 - C.θ * log(x)) * log(1 - C.θ * log(x)) / C.θ
-    result = Distributions.expectation(f,Distributions.Uniform(0,1))
-    # Calculate the integral using GSL
-    # result = 
-    # result, _ = GSL.integration_qags(f, 0.0, 1.0, [C.θ], 1e-7, 1000)
+    # Use a numerical integration method to obtain tau
+    result, _ = QuadGK.quadgk(x -> -((x-C.θ*x*log(x))*log(1-C.θ*log(x))/C.θ), 0, 1)
 
     return 1+4*result
 end
 function τ⁻¹(::Type{GumbelBarnettCopula}, tau)
     if tau == zero(tau)
         return tau
     end
+    if tau < 0
+        @warn "GumbelBarnettCopula cannot handle negative dependencies, returning independence..."
+        return zero(τ)
+    end
 
     # Define an anonymous function that takes a value x and computes τ 
     #for a GumbelBarnettCopula with θ = x
     τ_func(x) = τ(GumbelBarnettCopula(2,x))
 
     # Use the bisection method to find the root
-    x = Roots.find_zero(x -> τ_func(x) - tau, (0.0, 1.0))    
+    x = Roots.find_zero(x -> τ_func(x) - tau, (0.0, 1.0))
     return x
 end

src/ArchimedeanCopulas/GumbelCopula.jl

@@ -22,7 +22,6 @@ struct GumbelCopula{d,T} <: ArchimedeanCopula{d}
     θ::T
     function GumbelCopula(d,θ)
         if θ < 1
-            @show θ
             throw(ArgumentError("Theta must be greater than or equal to 1"))
         elseif θ == 1
             return IndependentCopula(d)

src/ArchimedeanCopulas/InvGaussianCopula.jl

@@ -39,26 +39,27 @@ end
 ϕ(  C::InvGaussianCopula,       t) = exp((1-sqrt(1+2*((C.θ)^(2))*t))/C.θ)
 ϕ⁻¹(C::InvGaussianCopula,       t) = ((1-C.θ*log(t))^(2)-1)/(2*(C.θ)^(2))
 function τ(C::InvGaussianCopula)
-    # Define the function to integrate
-    f(x) = -x * (((1 - C.θ * log(x))^2 - 1) / (2 * C.θ * (1 - C.θ * log(x))))
 
     # Calculate the integral using an appropriate numerical integration method
-    result, _ = gsl_integration_qags(f, 0.0, 1.0, [C.θ], 1e-7, 1000)
+    result, _ = QuadGK.quadgk( x -> (x*((1-C.θ*log(x))^2-1))/(-2*C.θ*(1-C.θ*log(x))),0,1)
 
     return 1+4*result
 end
-function τ⁻¹(::Type{InvGaussianCopula}, τ)
-    if τ == zero(τ)
-        return τ
+function τ⁻¹(::Type{InvGaussianCopula}, tau)
+    if tau == zero(tau)
+        return tau
+    end
+    if tau < 0
+        @warn "GumbelBarnettCopula cannot handle negative dependencies, returning independence..."
+        return zero(τ)
     end
 
     # Define an anonymous function that takes a value x and computes τ for an InvGaussianCopula copula with θ = x
-    τ_func(x) = τ(InvGaussianCopula{2, Float64}(x))
+    τ_func(x) = τ(InvGaussianCopula(2,x))
 
     # Set an initial value for x₀ (adjustable)
-    x₀ = (1-τ)/4
-
-    return Roots.find_zero(x -> τ_func(x) - τ, x₀)
+    x = Roots.find_zero(x -> τ_func(x) - tau, (0.0, Inf))
+    return τ
 end
 
 williamson_dist(C::InvGaussianCopula{d,T}) where {d,T} = WilliamsonFromFrailty(Distributions.InverseGaussian(C.θ,1),d)

src/Copulas.jl

@@ -3,7 +3,6 @@ module Copulas
     import Base
     import Random
     import SpecialFunctions
-    import GSL
     import Roots
     import Distributions
     import StatsBase
@@ -14,6 +13,7 @@ module Copulas
     import WilliamsonTransforms
     import Combinatorics
     import LogExpFunctions
+    import QuadGK
 
     # Standard copulas and stuff. 
     include("utils.jl")

test/archimedean_tests.jl

@@ -10,8 +10,8 @@
         GumbelCopula,
         # IndependentCopula,
         # JoeCopula,
-        # GumbelBarnettCopula,
-        # InvGaussianCopula
+         GumbelBarnettCopula,
+         InvGaussianCopula
     )
         for τ in taus
             @test Copulas.τ(T(2,Copulas.τ⁻¹(T,τ))) ≈ τ
@@ -43,11 +43,11 @@ end
     using StableRNGs
     using Distributions
     rng = StableRNG(123)
-    C = ClaytonCopula(2,-1)
-    data = rand(rng,C,100)
-    @test all(pdf(C,data) .>= 0)
-    @test all(0 .<= cdf(C,data) .<= 1)
-    fit(ClaytonCopula,data)
+    C0 = ClaytonCopula(2,-1)
+    data0 = rand(rng,C0,100)
+    @test all(pdf(C0,data0) .>= 0)
+    @test all(0 .<= cdf(C0,data0) .<= 1)
+    fit(ClaytonCopula,data0)
     for d in 2:10
         for θ ∈ [-1/(d-1) * rand(rng),0.0,-log(rand(rng)), Inf]
             C = ClaytonCopula(d,θ)
@@ -66,17 +66,17 @@ end
     rand(rng,FrankCopula(2,-Inf),10)
     rand(rng,FrankCopula(2,log(rand(rng))),10)
 
-    C = FrankCopula(2,-Inf)
-    data = rand(rng,C,100)
-    @test all(pdf(C,data) .>= 0)
-    @test all(0 .<= cdf(C,data) .<= 1)
-    @test_broken fit(FrankCopula,data)
+    C0 = FrankCopula(2,-Inf)
+    data0 = rand(rng,C0,100)
+    @test all(pdf(C0,data0) .>= 0)
+    @test all(0 .<= cdf(C0,data0) .<= 1)
+    @test_broken fit(FrankCopula,data0)
 
-    C = FrankCopula(2,log(rand(rng)))
-    data = rand(rng,C,100)
-    @test all(pdf(C,data) .>= 0)
-    @test all(0 .<= cdf(C,data) .<= 1)
-    @test_broken fit(FrankCopula,data)
+    C1 = FrankCopula(2,log(rand(rng)))
+    data1 = rand(rng,C1,100)
+    @test all(pdf(C1,data1) .>= 0)
+    @test all(0 .<= cdf(C1,data1) .<= 1)
+    @test_broken fit(FrankCopula,data1)
 
 
     for d in 2:10
@@ -96,7 +96,6 @@ end
     rng = StableRNG(123)
     for d in 2:10
         for θ ∈ [1.0,1-log(rand(rng)), Inf]
-            @show d,θ
             C = GumbelCopula(d,θ)
             data = rand(rng,C,100)
             @test all(pdf(C,data) .>= 0)
@@ -144,7 +143,6 @@ end
     rng = StableRNG(123)
     for d in 2:10
         for θ ∈ [rand(rng),1.0, -log(rand(rng))]
-            @show d,θ
             C = InvGaussianCopula(d,θ)
             data = rand(rng,C,100)
             @test all(pdf(C,data) .>= 0)

test/margins_uniformity.jl

@@ -32,7 +32,6 @@
     for C in cops
         nfail = 0
         d = length(C)
-        @show C
         spl = rand(rng,C,n)
         @assert all(0 <= x <= 1 for x in spl)
         for i in 1:d