add rho TCopula and Gaussian CDF

Project.toml

@@ -8,12 +8,13 @@ Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 HCubature = "19dc6840-f33b-545b-b366-655c7e3ffd49"
+HypergeometricFunctions = "34004b35-14d8-5ef3-9330-4cdb6864b03a"
 LambertW = "984bce1d-4616-540c-a9ee-88d1112d94c9"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LogExpFunctions = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
-MvNormalCDF = "37188c8d-bc69-4638-b057-733e744175ec"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 PolyLog = "85e3b03c-9856-11eb-0374-4dc1f8670e7f"
+Primes = "27ebfcd6-29c5-5fa9-bf4b-fb8fc14df3ae"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
@@ -37,6 +38,7 @@ Combinatorics = "1"
 Distributions = "0.25"
 ForwardDiff = "0.10, 1"
 HCubature = "1"
+HypergeometricFunctions = "0.3"
 HypothesisTests = "v0.11"
 InteractiveUtils = "1"
 LambertW = "1.0.0"
@@ -46,6 +48,7 @@ MvNormalCDF = "0.2, 0.3"
 Optim = "1.13.2"
 Plots = "1"
 PolyLog = "2.6.0"
+Primes = "0.5.7"
 Printf = "1"
 QuadGK = "2"
 Random = "1"
@@ -65,9 +68,10 @@ Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 HypothesisTests = "09f84164-cd44-5f33-b23f-e6b0d136a0d5"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+MvNormalCDF = "37188c8d-bc69-4638-b057-733e744175ec"
 StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test", "InteractiveUtils", "LinearAlgebra", "HypothesisTests", "Aqua", "StableRNGs", "StatsBase"]
+test = ["Test", "InteractiveUtils", "LinearAlgebra", "HypothesisTests", "Aqua", "StableRNGs", "StatsBase", "MvNormalCDF"]

src/Copulas.jl

@@ -1,24 +1,25 @@
 module Copulas
 
     import Base
-    import Random
-    import SpecialFunctions
-    import Roots
+    import Combinatorics
     import Distributions
-    import Statistics
-    import StatsBase
-    import StatsFuns
     import ForwardDiff
     import HCubature
-    import MvNormalCDF
-    import Combinatorics
-    import LogExpFunctions
-    import QuadGK
-    import LinearAlgebra
-    import PolyLog
+    import HypergeometricFunctions
     import LambertW
+    import LinearAlgebra
+    import LogExpFunctions
     import Optim
+    import PolyLog
+    import Primes
     import Printf
+    import QuadGK
+    import Random
+    import Roots
+    import SpecialFunctions
+    import Statistics
+    import StatsBase
+    import StatsFuns
     import TaylorSeries
 
     # Main code

src/EllipticalCopula.jl

@@ -121,4 +121,172 @@ end
     Σ = L * L'
     Σ = (Σ + Σ')/2
     return Σ
-end
+end
+
+# ===================== Φ y Φ^{-1} (erfc/erfcinv) =====================
+@inline Φinv(p::T) where T = @fastmath sqrt(T(2)) * SpecialFunctions.erfcinv(T(2) * (T(1) - p))
+@inline Φ(x::T) where T = @fastmath T(0.5) * SpecialFunctions.erfc(-x / sqrt(T(2)))
+@inline φ(x::T) where T = @fastmath inv(sqrt(T(2π))) * exp(-x*x / T(2))
+
+# Richtmyer Generators (shared)
+@inline _δ(::Type{T}) where {T<:Real} = T(sqrt(eps(T)))
+@inline richtmyer_roots(T, n) = sqrt.(T.(Float64.(Primes.primes(1, max(n-1, Int(floor(5n*log(n+1)/4)))))))[1:n-1]
+function _chlrdr_orthant!(R::AbstractMatrix{T}, b::AbstractVector{T}) where {T}
+    @boundscheck (size(R,1) == size(R,2) == length(b)) || throw(DimensionMismatch())
+    n = length(b)
+    c = R           # trabajamos in-place
+    bp = b
+    y = zeros(T, n)
+
+    ϵ  = eps(T)
+    ep = T(1e-10)
+
+    @inbounds for k in 1:n
+        im = k
+        ckk = zero(T)
+        dem = one(T)
+        bm_tilde = zero(T)
+
+        # --- elegir pivote: min Φ(b̃) ---
+        for i in k:n
+            cii = c[i,i]
+            if cii > ϵ
+                cii_sqrt = sqrt(cii)
+                s = zero(T)
+                if k > 1
+                    @simd for j in 1:(k-1)
+                        s += c[i,j]*y[j]
+                    end
+                end
+                b_tilde = (bp[i] - s) / cii_sqrt
+                de = Φ(b_tilde)
+                if de <= dem
+                    dem = de; ckk = cii_sqrt; im = i; bm_tilde = b_tilde
+                end
+            end
+        end
+
+        # --- permutar im ↔ k ---
+        if im > k
+            c[im,im], c[k,k] = c[k,k], c[im,im]
+            bp[im],  bp[k]  = bp[k],  bp[im]
+            if k > 1
+                @simd for j in 1:(k-1)
+                    c[im,j], c[k,j] = c[k,j], c[im,j]
+                end
+            end
+            if im < n
+                @simd for j in (im+1):n
+                    c[j,im], c[j,k] = c[j,k], c[j,im]
+                end
+            end
+            if k < im - 1
+                for j in (k+1):(im-1)
+                    c[j,k], c[im,j] = c[im,j], c[j,k]
+                end
+            end
+        end
+
+        # anular parte superior de la fila k
+        if k < n
+            @simd for j in (k+1):n
+                c[k,j] = zero(T)
+            end
+        end
+
+        # --- actualización y y[k] ---
+        if ckk > k*ep
+            c[k,k] = ckk
+            inv_ckk = one(T)/ckk
+            for i in (k+1):n
+                c[i,k] *= inv_ckk
+            end
+            # actualización de rango-1 por doble bucle
+            for i in (k+1):n
+                cik = c[i,k]
+                @simd for j in (k+1):i
+                    c[i,j] -= cik * c[j,k]
+                end
+            end
+
+            # --- razón de Mills estable ---
+            if dem > ep  # dem ≈ Φ(b̃)
+                if bm_tilde < -T(10)
+                    # límite inferior: Φ(b̃) ~ 0 → usar aproximación de Mills
+                    y[k] = bm_tilde  # E[Z | Z ≤ b̃] ≈ b̃
+                else
+                    # fórmula exacta y estable para todo b̃ ∈ [-10, +∞)
+                    @fastmath y[k] = -φ(bm_tilde) / dem
+                end
+            else
+                # caso degenerado (Φ(b̃) ≈ 0)
+                y[k] = bm_tilde
+            end
+        else
+            # columna singular
+            @simd for i in k:n
+                c[i,k] = zero(T)
+            end
+            y[k] = zero(T)
+        end
+    end
+
+    return (c, bp)
+end
+# Generic kernel: assumes that rfill!(rvec) writes per-sample scales
+function qmc_orthant_core!(ch::AbstractMatrix{T}, bs::AbstractVector{T}; m::Integer=10_000, r::Integer=12,
+    rng::Random.AbstractRNG=Random.default_rng(), fill_w!::Function = (w::AbstractVector{T}, _j, _nv, _δ::T, _rng)->fill!(w, one(T))) where {T}
+
+    n = length(bs)
+    (size(ch,1) == size(ch,2) == n) || throw(DimensionMismatch())
+
+    q  = richtmyer_roots(T, n)               # Richtmyer √primes (n-1)
+    nv = max(div(m, r), 1)
+
+    y  = zeros(T, nv, n-1)
+    pv = Vector{T}(undef, nv)
+    dc = Vector{T}(undef, nv)
+    tv = Vector{T}(undef, nv)
+    w  = Vector{T}(undef, nv)
+
+    δ = _δ(T)
+    p = zero(T); e = zero(T)
+
+    @inbounds for j in 1:r
+        # generate scales w[k] for this replicate (t: random; normal: w=1)
+        fill_w!(w, j, nv, δ, rng)
+
+        # first coordinate: P(Z1 ≤ w*b1)
+        @simd for k in 1:nv
+            d1k = Φ((bs[1]/w[k]) / ch[1,1])
+            pv[k] = d1k
+            dc[k] = d1k
+        end
+
+        # dimensions 2..n
+        for i in 2:n
+            qi = q[i-1]; xr = rand(rng, T)
+            @inbounds @simd for k in 1:nv
+                t = k*qi + xr; t -= floor(t)
+                u = abs(2*t - one(T)) * dc[k]         # u ∈ (0,1)
+                p_safe = clamp(u, δ, one(T)-δ)
+                y[k, i-1] = Φinv(p_safe)
+            end
+            LinearAlgebra.mul!(tv, @view(y[:, 1:i-1]), @view(ch[i, 1:i-1]))
+            ci_inv = inv(ch[i,i]); bi = bs[i]
+            @inbounds @simd for k in 1:nv
+                val = Φ( (bi/w[k] - tv[k]) * ci_inv )
+                pv[k] *= val
+                dc[k]  = val
+            end
+        end
+
+        pj = sum(pv) / T(nv)
+        dm = (pj - p) / T(j)
+        p  += dm
+        e   = (T(j) - 2)*e/T(j) + dm*dm
+    end
+
+    return (p, r==1 ? zero(T) : 3*sqrt(e))
+end
+

src/EllipticalCopulas/GaussianCopula.jl

@@ -79,11 +79,6 @@ GaussianCopula{D, MT}(d::Int, ρ::Real) where {D, MT} = GaussianCopula(d, ρ)
 
 U(::Type{T}) where T<: GaussianCopula = Distributions.Normal()
 N(::Type{T}) where T<: GaussianCopula = Distributions.MvNormal
-function _cdf(C::CT,u) where {CT<:GaussianCopula}
-    x = StatsBase.quantile.(Distributions.Normal(), u)
-    d = length(C)
-    return MvNormalCDF.mvnormcdf(C.Σ, fill(-Inf, d), x)[1]
-end
 
 function rosenblatt(C::GaussianCopula, u::AbstractMatrix{<:Real})
     return Distributions.cdf.(Distributions.Normal(), inv(LinearAlgebra.cholesky(C.Σ).L) * Distributions.quantile.(Distributions.Normal(), u))
@@ -141,4 +136,17 @@ function _fit(CT::Type{<:GaussianCopula}, u, ::Val{:mle})
     Σ = Matrix(dd.Σ)
     return GaussianCopula(Σ), (; θ̂ = (; Σ = Σ))
 end
-_available_fitting_methods(::Type{<:GaussianCopula}, d) = (:mle, :itau, :irho, :ibeta)
+_available_fitting_methods(::Type{<:GaussianCopula}, d) = (:mle, :itau, :irho, :ibeta)
+
+function qmc_orthant_normal!(Σ::AbstractMatrix{T}, b::AbstractVector{T}; m::Integer=10_000, r::Integer=12,
+    rng::Random.AbstractRNG=Random.default_rng()) where {T}
+    (ch, bs) = _chlrdr_orthant!(Σ, b)    # ¡muta Σ y b!
+    qmc_orthant_core!(ch, bs; m=m, r=r, rng=rng)
+end
+function Distributions.cdf(C::GaussianCopula{d, MT}, u::AbstractVector; m::Integer = 2000*(d+1), r::Int = 8, rng = Random.default_rng()) where {d, MT}
+    x = StatsBase.quantile.(Distributions.Normal(), u)
+    Tx = eltype(x)
+    Σ_promoted = Tx.(copy(C.Σ))
+    p, _ = qmc_orthant_normal!(Σ_promoted, x; m=m, r=r, rng=rng)
+    return p
+end

src/EllipticalCopulas/TCopula.jl

@@ -49,7 +49,39 @@ end
 # Kendall tau of bivariate student: 
 # Lindskog, F., McNeil, A., & Schmock, U. (2003). Kendall’s tau for elliptical distributions. In Credit risk: Measurement, evaluation and management (pp. 149-156). Heidelberg: Physica-Verlag HD.
 τ(C::TCopula{2,MT}) where MT = 2*asin(C.Σ[1,2])/π 
-
+function τ(C::TCopula{d,MT}) where {d, MT}
+    T = (2/π) .* asin.(C.Σ)
+    @inbounds for i in 1:d
+        T[i,i] = 1.0
+    end
+    return LinearAlgebra.Symmetric(T, :U)
+end
+##############################
+function ρ(C::TCopula{2,ν,MT}) where {ν,MT}
+    ρ_ = C.Σ[1,2]
+    rtol = 1e-10
+    #  Normalization constant off_{Ṽ}
+    Cν = 2 * SpecialFunctions.gamma(ν)^2 * SpecialFunctions.gamma(3ν/2) / (SpecialFunctions.gamma(ν/2)^3 * SpecialFunctions.gamma(2ν))
+    f(v) = begin
+        # if we use HypergeometricFunctions.jl we can make:
+        F = HypergeometricFunctions.pFq((ν, ν), (2ν,), 1 - v^2)
+        # and if not... The implemented functions work well and in particular are quite fast.
+        # F = Copulas._Gauss2F1_hybrid(ν, 1 - v^2)
+        return asin(ρ_ * v) * Cν * v^(ν - 1) * (1 - v^2)^(ν/2 - 1) * F
+    end
+    try
+        val, _ = QuadGK.quadgk(f, 0.0, 1.0; rtol=rtol)
+        return (6/π) * val
+    catch err
+        if ν > 20
+            # asymptotic fallback (equivalent to normal copula)
+            ρ_norm = (6/π) * asin(ρ_/2)
+            return ρ_norm
+        else
+            rethrow(err)
+        end
+    end
+end
 # Conditioning colocated
 function DistortionFromCop(C::TCopula{D,ν,MT}, js::NTuple{p,Int}, uⱼₛ::NTuple{p,Float64}, i::Int) where {p,D,ν,MT}
     Σ = C.Σ; jst = js; ist = Tuple(setdiff(1:D, jst)); @assert i in ist
@@ -97,4 +129,37 @@ function _rebound_params(::Type{<:TCopula}, d::Int, α::AbstractVector{T}) where
     return (; ν = ν, Σ = Σ)
 end
 
-_available_fitting_methods(::Type{<:TCopula}, d) = (:mle,)
+_available_fitting_methods(::Type{<:TCopula}, d) = (:mle,)
+
+# t-ortant (copulates t with ν g.l.)
+function qmc_orthant_t!(R::AbstractMatrix{T}, b::AbstractVector{T}, ν::Integer; m::Integer = 10_000, r::Integer = 12,
+    rng::Random.AbstractRNG = Random.default_rng()) where T
+    # ¡muta R y b!
+    (ch, bs) = _chlrdr_orthant!(R, b)
+
+    # extra Richtmyer root for the radial dimension (χ²)
+    qχ  = richtmyer_roots(T, length(b) + 1)[end]
+    chi = Distributions.Chisq(ν)
+
+    # scale generator w[k] = √(ν / S_k), S_k ~ χ²_ν (quasi-random)
+    fill_w! = function (w::AbstractVector{T}, _j::Int, nv::Int, δ::T, rng_local)
+        xrχ = rand(rng_local, T)
+        @inbounds @simd for k in 1:nv
+            t = k*qχ + xrχ; t -= floor(t)
+            u = clamp(t, δ, one(T)-δ)                    # u ∈ (δ, 1-δ)
+            s = T(Distributions.quantile(chi, Real(u)))            # quantile χ²_ν
+            w[k] = sqrt(T(ν) / s)                       # radial scale
+        end
+        nothing
+    end
+
+    return qmc_orthant_core!(ch, bs; m=m, r=r, rng=rng, fill_w! = fill_w!)
+end
+
+function Distributions.cdf(C::TCopula{d,df,MT}, u::AbstractVector; m::Integer = 2000*(d+1), r::Int = 12, rng = Random.default_rng()) where {d,df,MT}
+    b = Distributions.quantile.(Distributions.TDist(df), u)
+    Tb = eltype(b)
+    Σ_promoted = Tb.(copy(C.Σ))
+    p, _ = qmc_orthant_t!(Σ_promoted, b, df; m=m, r=r, rng=rng)
+    return p
+end