Merge pull request #40 from Santymax98/PlackettCopula

src/Copulas.jl

@@ -59,8 +59,9 @@ end Copulas
     include("MiscellaneousCopulas/MCopula.jl")
     include("MiscellaneousCopulas/WCopula.jl")
     include("MiscellaneousCopulas/SurvivalCopula.jl")
+    include("MiscellaneousCopulas/PlackettCopula.jl")
     export MCopula,
            WCopula,
-           SurvivalCopula
-
+           SurvivalCopula,
+           PlackettCopula
 end

src/MiscellaneousCopulas/PlackettCopula.jl

@@ -0,0 +1,69 @@
+#= Details about Plackett copulation are found in Joe, H. (2014). 
+   Dependence modeling with copulas. CRC press, Page.164
+==#
+
+#Create an instance of the Plackett copula
+struct PlackettCopula{P} <: Copula{2} # since it is only bivariate.
+    θ::P  # Copula parameter
+
+    function PlackettCopula(θ)
+        if θ == 1
+            return IndependentCopula(2)
+        elseif θ == 0
+            return MCopula(2)
+        elseif θ == Inf
+            return WCopula(2)
+        else
+            θ >= 0 || throw(ArgumentError("Theta must be non-negative"))
+            return new{typeof(θ)}(θ)
+        end
+    end
+end
+
+# Base.length(S::PlackettCopula{P}) where {P} = 2 # length should return the dimension of the copula, bnut i think it is already working without this definition.
+Base.eltype(S::PlackettCopula{P}) where {P} = P # this shuold be P. 
+
+# CDF calculation for bivariate Plackett Copula
+function Distributions.cdf(S::PlackettCopula{P}, uv) where {P}
+    u, v = uv
+    η = S.θ - 1
+    term1 = 1 + η * (u + v)
+    term2 = sqrt(term1^2 - 4 * S.θ * η * u * v)
+    return 0.5 * η^(-1) * (term1 - term2)
+end
+
+# PDF calculation for bivariate Plackett Copula
+function Distributions._logpdf(S::PlackettCopula{P}, uv) where {P}
+    u, v = uv
+    η = S.θ - 1
+    term1 = S.θ * (1 + η * (u + v - 2 * u * v))
+    term2 = (1+η*(u+v))^2-4*(S.θ)*η*u*v
+    return log(term1) - 3 * log(term2)/2 # since we are supposed to return the logpdf. 
+end
+import Random
+
+#=  Details about the algorithm to generate copula samples 
+    can be seen in the following references
+    Johnson, Mark E. Multivariate statistical simulation:
+    A guide to selecting and generating continuous multivariate distributions.
+    Vol. 192. John Wiley & Sons, 1987. Page 193.
+    Nelsen, Roger B. An introduction to copulas. Springer, 2006. Exercise 3.38.
+==#
+
+function Distributions._rand!(rng::Distributions.AbstractRNG, C::CT, x::AbstractVector{T}) where {T<:Real, CT<:PlackettCopula}
+    u = rand(rng)
+    t = rand(rng)
+    a = t * (1 - t)
+    b = C.θ + a * (C.θ - 1)^2
+    cc = 2a * (u * C.θ^2 + 1 - u) + C.θ * (1 - 2a)
+    d = sqrt(C.θ) * sqrt(C.θ + 4a * u * (1 - u) * (1 - C.θ)^2)
+    v = (cc - (1 - 2t) * d) / (2b)
+    x[1] = u
+    x[2] = v
+    return x
+end
+
+# Calculate Spearman's rho based on the PlackettCopula parameters
+function ρ(c::PlackettCopula{P}) where P
+    return (c.θ+1)/(c.θ-1)-(2*c.θ*log(c.θ)/(c.θ-1)^2)
+end

test/biv_Plackett.jl

@@ -0,0 +1,68 @@
+@testitem "PlackettCopula Constructor" begin
+    @test isa(PlackettCopula(1), IndependentCopula)
+    @test isa(PlackettCopula(Inf),WCopula) # should work in any dimenisons if theta is smaller than the bound.
+    @test isa(PlackettCopula(0),MCopula)
+    @test_throws ArgumentError PlackettCopula(-0.5)
+end
+
+@testitem "PlackettCopula CDF" begin
+    using Distributions
+    u = 0.1:0.18:1
+    v = 0.4:0.1:0.9
+    l1 = [
+        0.055377800527509735,
+        0.1743883734874062,
+        0.3166277269195278,
+        0.48232275012183223,
+        0.6743113969874872,
+        0.8999999999999999
+    ]
+    l2 = [
+        0.026208734813001233,  
+        0.10561162651259381,
+        0.23491134194308438,
+        0.4162573282722253,
+        0.6419254774317229,
+        0.9
+    ]
+
+    for i in 1:6
+        @test cdf(PlackettCopula(2.0), [u[i], v[i]]) ≈ l1[i]
+        @test cdf(PlackettCopula(0.5), [u[i], v[i]]) ≈ l2[i]
+    end
+end
+
+@testitem "PlackettCopula PDF" begin
+    using Distributions
+    u = 0.1:0.18:1
+    v = 0.4:0.1:0.9
+    l1 = [
+        1.0592107420343486,
+        1.023290881054283,
+        1.038466936984394,
+        1.1100773231007635,
+        1.2729591789643138,
+        1.652892561983471
+    ]
+    l2 = [
+        0.8446203068160272,
+        1.023290881054283,
+        1.0648914416282562,
+        0.9360170818943749,
+        0.7346611825055718,
+        0.5540166204986149
+    ]
+
+    for i in 1:6
+        @test pdf(PlackettCopula(2.0), [u[i], v[i]]) ≈ l1[i]
+        @test pdf(PlackettCopula(0.5), [u[i], v[i]]) ≈ l2[i]
+    end
+end
+
+@testitem "PlackettCopula Sampling" begin
+    using Random
+    n_samples = 100
+    C = PlackettCopula(0.8)
+    samples = rand(C, n_samples)
+    @test size(samples) == (2, n_samples)
+end