Nested Archimedean copulas

Project.toml

@@ -36,6 +36,7 @@ CopulasPlotsExt = ["Plots", "RecipesBase"]
 ADTypes = "1.21.0"
 Aqua = "v0.8"
 Combinatorics = "1"
+DelimitedFiles = "1"
 Distributions = "0.25"
 ForwardDiff = "0.10,1"
 HCubature = "1"
@@ -64,6 +65,7 @@ julia = "1"
 
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
+DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
 HypothesisTests = "09f84164-cd44-5f33-b23f-e6b0d136a0d5"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -72,4 +74,4 @@ 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", "DelimitedFiles"]

docs/make.jl

@@ -32,6 +32,7 @@ makedocs(;
         "Bestiary" => [
             "Elliptical copulas"=>"bestiary/elliptical.md",
             "Archimedean copulas"=>"bestiary/archimedean.md",
+            "Nested Archimedean copulas"=>"bestiary/nested.md",
             "Extreme Value copulas"=>"bestiary/extremevalues.md",
             "Archimax copulas"=>"bestiary/archimax.md",
             "Empirical copulas"=>"bestiary/empirical.md",

docs/src/bestiary/nested.md

@@ -0,0 +1,260 @@
+```@meta
+CurrentModule = Copulas
+```
+
+# Nested Archimedean copulas
+
+A *nested* (hierarchical) Archimedean copula glues several Archimedean copulas
+together under an outer Archimedean generator, letting different blocks of
+variables share a stronger within-block dependence while still being coupled
+across blocks. This is the natural model for grouped or hierarchical
+dependence — for example several organ systems within a patient, or several
+assets within a sector.
+
+[`NestedArchimedeanCopula`](@ref) provides the density of such trees (and, via
+the standard [`condition`](@ref) / [`subsetdims`](@ref) framework, a
+per-variable censored / survival likelihood), following the algorithm of Yang &
+Li ([arXiv:2605.23134](https://arxiv.org/abs/2605.23134)).
+
+## Definition
+
+With an outer generator ``\phi_0`` over inner copulas ``C_1, \dots, C_m`` on
+disjoint coordinate blocks ``I_1, \dots, I_m`` (and possibly some bare
+coordinates attached directly to the root), the CDF is
+
+```math
+C(\mathbf u) = \phi_0\!\left(
+  \sum_{i \in \text{root leaves}} \phi_0^{-1}(u_i)
+  \;+\; \sum_{k=1}^m \phi_0^{-1}\bigl(C_k(\mathbf u_{I_k})\bigr)
+\right),
+```
+
+and each child ``C_k`` is itself Archimedean or, recursively, nested
+Archimedean — so trees nest to arbitrary depth.
+
+The density is the mixed partial of this CDF over the differentiated
+coordinates. Differentiating the composition of generators is exactly Faà di
+Bruno's formula; the implementation carries the partial Bell polynomials through
+truncated Taylor series over the generator tree, building only on the package's
+generator interface (`ϕ`, `ϕ⁻¹`, `ϕ⁽ᵏ⁾`).
+
+## Building a tree
+
+```@example nested
+using Copulas, Distributions
+using Copulas: ClaytonGenerator, JoeGenerator
+
+# Outer Clayton(2) over two inner Clayton panels on dims 1:2 and 3:4.
+C = NestedArchimedeanCopula(ClaytonGenerator(2.0);
+        children = [ClaytonCopula(2, 5.0), ClaytonCopula(2, 6.0)])
+logpdf(C, [0.3, 0.5, 0.4, 0.6])
+```
+
+Children are auto-placed on consecutive free dimension blocks in declaration
+order. You can instead pin a child to explicit dimensions with a `Pair`, and you
+can attach bare coordinates to the root with `leaves`:
+
+```@example nested
+# Root Clayton with a bare leaf on dim 1 and a Gumbel panel on dims 2:4.
+C2 = NestedArchimedeanCopula(ClaytonGenerator(2.0);
+         leaves = [1], children = [GumbelCopula(3, 2.0) => [2, 3, 4]])
+logpdf(C2, [0.25, 0.4, 0.55, 0.7])
+```
+
+Mixed families and arbitrary nesting depth are supported:
+
+```@example nested
+inner = NestedArchimedeanCopula(JoeGenerator(3.0); children = [JoeCopula(2, 4.0)])
+C3 = NestedArchimedeanCopula(ClaytonGenerator(1.5);
+         children = [GumbelCopula(2, 2.0), FrankCopula(2, 3.0), inner])
+logpdf(C3, [0.2, 0.3, 0.4, 0.5, 0.6, 0.7])
+```
+
+A purely flat declaration (only `leaves`, no `children`) returns the package's
+native [`ArchimedeanCopula`](@ref) so its fast specialised density is used.
+
+!!! note "Validity is the caller's responsibility"
+    The constructor does not check the nesting condition. For same-family
+    nestings the standard sufficient condition is that the inner generator be at
+    least as dependent as the outer one (e.g. for Clayton/Gumbel/Joe, the inner
+    parameter ``\ge`` the outer parameter). Mixed-family nestings are accepted
+    but must be validated by the user.
+
+## Fitting
+
+`fit` performs maximum-likelihood estimation of the generator parameters on a
+**fixed tree**: the leaf layout and the generator family at each node come from a
+template instance, and only the scalar ``\theta`` of each node is optimised. Pass
+the template and a `d×n` matrix of pseudo-observations (columns are observations).
+
+```@example nested
+using Random
+Ctrue = NestedArchimedeanCopula(ClaytonGenerator(2.0);
+            children = [ClaytonCopula(2, 6.0), ClaytonCopula(2, 8.0)])
+U = rand(Random.MersenneTwister(1), Ctrue, 2000)
+
+# Fit from a deliberately wrong same-shape template:
+Cstart = NestedArchimedeanCopula(ClaytonGenerator(1.0);
+             children = [ClaytonCopula(2, 3.0), ClaytonCopula(2, 3.0)])
+M = fit(CopulaModel, Cstart, U)
+M.result
+```
+
+The optimiser runs in an unconstrained space through a *parametrisation* — a map
+`α -> NestedArchimedeanCopula` decoupled from the generator objects. Above we fit
+a **template** tree. For full control, pass your own map and its initial point,
+`fit(CopulaModel, reparam, init, U)` — no template needed, since `reparam` builds
+the whole tree. This lets you share parameters across nodes, fit on a different
+scale, or encode a constraint.
+
+For instance, enforce the nesting condition by building each child's ``\theta`` as
+a non-negative increment over its parent's, so every step is a valid nesting:
+
+```@example nested
+softplus(x) = log1p(exp(-abs(x))) + max(x, zero(x))
+nest = α -> NestedArchimedeanCopula(ClaytonGenerator(exp(α[1]));
+    children = [ClaytonCopula(2, exp(α[1]) + softplus(α[2])),
+                ClaytonCopula(2, exp(α[1]) + softplus(α[3]))])
+Mn = fit(CopulaModel, nest, [0.0, 0.0, 0.0], U)
+(Mn.result.G.θ, Mn.result.children[1][1].G.θ)   # inner θ ≥ outer θ, by construction
+```
+
+Or share one ``\theta`` across the root and both panels — a single free parameter:
+
+```@example nested
+recon = α -> (θ = exp(α[1]);
+    NestedArchimedeanCopula(ClaytonGenerator(θ);
+        children = [ClaytonCopula(2, θ), ClaytonCopula(2, θ)]))
+Ms = fit(CopulaModel, recon, [0.0], U)
+(Ms.result.G.θ, Ms.result.children[1][1].G.θ)   # equal — the shared parameter
+```
+
+`fit(C0, U)` is a shorthand returning just the fitted copula; for the custom form
+use `fit(CopulaModel, reparam, init, U).result`.
+
+## Precision
+
+The recursion is generic in the value type. `logpdf` works on `Float64`
+out of the box; passing `BigFloat` (or `Double64`) coordinates carries that
+precision through the whole recursion, which is recommended for adversarial
+high-dimensional or deep-tail inputs where the alternating-sign Faà di Bruno sum
+can lose `Float64` precision:
+
+```@example nested
+logpdf(C, big.([0.3, 0.5, 0.4, 0.6]))
+```
+
+## Edge-composition method
+
+Each parent→child edge in the Faà di Bruno recursion needs the truncated Taylor
+expansion of the inner-to-outer link ``h = \phi^{-1}_{\text{outer}} \circ
+\phi_{\text{inner}}`` at the child's argument. It goes through the overloadable
+hook `composition_taylor(outer, inner, t₀, d)`, selected by dispatch exactly as
+you override `ϕ⁽ᵏ⁾` — most-specific method wins, no keyword or flag. Three methods
+are available:
+
+**1. Direct (default).** `composition_taylor_direct` puts a single jet through
+the explicit composition. Fast and accurate for ordinary inputs. It requires both
+``\phi`` and ``\phi^{-1}`` to accept a `Taylor1` argument.
+
+**2. Implicit.** `composition_taylor_implicit` solves
+``\phi_{\text{outer}}(h(t)) = \phi_{\text{inner}}(t)`` order-by-order, using only
+the scalar derivatives ``\phi^{(k)}`` of both generators and a single scalar
+``\phi^{-1}_{\text{outer}}`` — it never puts a `Taylor1` through ``\phi^{-1}``.
+Use it when a generator's ``\phi^{-1}`` has no `Taylor1` method (for instance, an
+inverse defined only through root-finding), where the direct jet cannot run.
+Select it globally by redefining the generic method:
+
+```julia
+Copulas.composition_taylor(o::Copulas.Generator, i::Copulas.Generator, t₀, d) =
+    Copulas.composition_taylor_implicit(o, i, t₀, d)
+```
+
+(Redefining the shipped default prints a benign "method overwritten" warning.)
+
+**3. Closed form (per generator pair).** Register a more-specific method when you
+know ``h`` analytically — the fastest and most robust option. The package ships
+Clayton/Clayton (in `Generator/ClaytonGenerator.jl`),
+``h(t) = ((1+\theta_{\text{in}} t)^{\theta_{\text{out}}/\theta_{\text{in}}} - 1)/
+\theta_{\text{out}}``; add your own the same way:
+
+```julia
+function Copulas.composition_taylor(outer::MyGenerator, inner::MyGenerator, t₀, d::Int)
+    # return [h'(t₀)/1!, …, h⁽ᵈ⁾(t₀)/d!]
+end
+```
+
+The method choice is about *availability and speed*, not accuracy: all three are
+exact in exact arithmetic. For deep-tail or very high-``d`` inputs where a
+`Float64` jet can overflow or lose precision, pass `BigFloat` coordinates (see
+[Precision](#Precision)) — that is the precision fix, independent of which
+composition method is in use.
+
+## Censored / survival likelihood
+
+Per-variable (right-)censoring is an **emergent capability** of the standard
+[`condition`](@ref) + [`subsetdims`](@ref) framework — there is no bespoke
+censored-likelihood function. Split the coordinates into an observed set ``O``
+(events) and a censored set ``C`` (the survival times we only know exceed their
+observation). The per-variable censored likelihood factorises as the *gist
+recipe*
+
+```math
+\underbrace{\log f_{O}(x_O)}_{\text{observed-marginal density}}
+\;+\;
+\underbrace{\log P(X_C \le x_C \mid X_O = x_O)}_{\text{conditional survival}},
+```
+
+In code these two terms are `logpdf(subsetdims(X, O), x_O)` and
+`logcdf(condition(X, O, x_O), x_C)`.
+
+which equals the observed-marginal density times the copula's mixed partial over
+the observed coordinates,
+
+```math
+\sum_{i\in O} \log f_i(x_i)
+\;+\; \log \frac{\partial^{|O|} C(\mathbf u)}{\prod_{i\in O}\partial u_i},
+\qquad \mathbf u = (F_1(x_1),\dots,F_d(x_d)),
+```
+
+because the denominator ``c_O`` in `condition` cancels against the
+`subsetdims` marginal density. Both factors route through the Faà di Bruno tree
+walk via the `subsetdims` / `condition` specialisations for this type — no
+ForwardDiff for the observed-marginal density nor for the conditional CDF (for
+any number of censored coordinates).
+
+```@example nested
+using Distributions
+Csurv = NestedArchimedeanCopula(ClaytonGenerator(2.0);
+            children = [ClaytonCopula(3, 5.0), ClaytonCopula(3, 6.0)])
+S = SklarDist(Csurv, ntuple(_ -> Exponential(1.0), 6))
+x = [0.7, 0.3, 0.9, 0.5, 0.4, 1.1]
+O = (1, 3, 4, 5)        # observed (events)
+C = (2, 6)              # right-censored
+logpdf(subsetdims(S, O), x[collect(O)]) +
+    log(cdf(condition(S, O, x[collect(O)]), x[collect(C)]))
+```
+
+When a *single* coordinate is censored, `condition(S, O, x_O)` returns a
+univariate conditional distribution and you use `logcdf(condition(...), x_C)`
+(a scalar `x_C`); when several are censored it returns a conditional joint
+distribution and you use `log(cdf(condition(...), x_C))` as above. With ``C``
+empty (all observed) the recipe reduces to the ordinary joint density
+`logpdf(S, x)`; with ``O`` empty (all censored) it reduces to
+``\log F(\mathbf x)``. The recipe works for any copula with `condition` /
+`subsetdims` support — flat [`ArchimedeanCopula`](@ref) as well as nested trees.
+
+!!! note "Multi-censored conditional CDF"
+    With two or more censored coordinates the conditional CDF is the mixed
+    partial of the nested CDF over the *observed* coordinates. The generic path
+    takes this by nesting one `ForwardDiff.derivative` per observed coordinate —
+    cost exponential in the number of observed dims, infeasible in high
+    dimension. A `_partial_cdf` specialisation routes it instead through the same
+    polynomial Faà di Bruno tree walk as the single-censored case (selected on the
+    conditional copula's concrete nested inner type), for any number of censored
+    coordinates.
+
+    At high differentiation order for fast-tail generators the `Float64` sum can
+    lose precision; pass `BigFloat` coordinates to recover the exact value (as for
+    the density). End-to-end `BigFloat` through `condition()` is not yet enabled —
+    upstream stores the conditioning values as `Float64`.

src/Copulas.jl

@@ -101,6 +101,14 @@ module Copulas
     include("Generator/InvGaussianGenerator.jl")
     include("Generator/JoeGenerator.jl")
 
+    # Nested (hierarchical) Archimedean copulas
+    include("nested/NestedArchimedeanDensity.jl")
+    include("NestedArchimedeanCopula.jl")
+    # Conditioning/subsetting fast paths (per-variable censoring is emergent here)
+    include("nested/NestedConditioning.jl")
+    # Maximum-likelihood fitting of a nested copula's generator params (fixed tree)
+    include("nested/NestedArchimedeanFitting.jl")
+
     #Extreme value copulas
     include("Tail.jl")
     include("ExtremeValueCopula.jl")
@@ -142,6 +150,7 @@ module Copulas
            SklarDist, # SklarDist to make multivariate models
            AMHCopula, # And a bunch of copulas.
            ArchimedeanCopula,
+           NestedArchimedeanCopula,
            AsymGalambosCopula,
            AsymLogCopula,
            AsymMixedCopula,

src/Generator/ClaytonGenerator.jl

@@ -57,7 +57,31 @@ max_monotony(G::ClaytonGenerator) = G.θ >= 0 ? Inf : (1 - 1/G.θ)
 ϕ⁽¹⁾(G::ClaytonGenerator, t) = (1+G.θ*t) ≤ 0 ? 0 : - (1+G.θ*t)^(-1/G.θ -1)
 ϕ⁻¹⁽¹⁾(G::ClaytonGenerator, t) = -t^(-G.θ-1)
 ϕ⁽ᵏ⁾(G::ClaytonGenerator, k::Int, t) = (1+G.θ*t) ≤ 0 ? 0 : (1 + G.θ * t)^(-1/G.θ - k) * prod(-1-ℓ*G.θ for ℓ in 0:k-1; init=1)
-ϕ⁽ᵏ⁾⁻¹(G::ClaytonGenerator, k::Int, t; start_at=t) = ((t / prod(-1-ℓ*G.θ for ℓ in 0:k-1; init=1))^(1/(-1/G.θ - k)) -1)/G.θ    
+ϕ⁽ᵏ⁾⁻¹(G::ClaytonGenerator, k::Int, t; start_at=t) = ((t / prod(-1-ℓ*G.θ for ℓ in 0:k-1; init=1))^(1/(-1/G.θ - k)) -1)/G.θ
+
+# Closed-form edge-composition override for a Clayton-over-Clayton nesting. Overrides the
+# default `composition_taylor` hook (nested/NestedArchimedeanDensity.jl) by dispatch, and
+# returns the Taylor coefficients [h'(t₀)/1!, …, h⁽ᵈ⁾(t₀)/d!] of the inner→outer change of
+# variables h(t) = ϕ⁻¹_outer(ϕ_inner(t)). With ϕ_θ(t)=(1+θt)^(-1/θ) and ϕ⁻¹_θ(u)=(u^(-θ)-1)/θ,
+# the link is h(t) = ((1+θ_in·t)^r − 1)/θ_out, r = θ_out/θ_in — a reparametrised power map
+# whose coefficients are a generalized binomial series, so it never touches the (ill-
+# conditioned) high-order derivatives of the inverse. NOTE: the θ live INSIDE ϕ here, so the
+# expansion base is 1+θ_in·t₀ (NOT 1+t₀); θ promotes into T so Float64/BigFloat stay exact.
+function composition_taylor(outer::ClaytonGenerator, inner::ClaytonGenerator, t₀::T, d::Int) where {T}
+    θ_out = T(outer.θ)
+    θ_in  = T(inner.θ)
+    r     = θ_out / θ_in
+    base  = 1 + θ_in * t₀
+    h = Vector{T}(undef, d)
+    binom = one(T)                                     # generalized binomial C(r,k), incremental
+    θ_in_pow = one(T)                                  # θ_in^k
+    for k in 1:d
+        binom    *= (r - (k - 1)) / k                  # C(r,k) = C(r,k-1)·(r-k+1)/k
+        θ_in_pow *= θ_in
+        h[k] = (θ_in_pow / θ_out) * binom * base^(r - k)
+    end
+    return h
+end
 
 τ(G::ClaytonGenerator) = ifelse(isfinite(G.θ), G.θ/(G.θ+2), 1)
 τ⁻¹(::Type{<:ClaytonGenerator},τ) = ifelse(τ == 1,Inf,2τ/(1-τ))

src/Generator/FrankGenerator.jl

@@ -55,6 +55,17 @@ function ϕ⁽ᵏ⁾(G::FrankGenerator, k::Int, t)
     return (-1)^k * (1 / G.θ) * PolyLog.reli(-(k - 1), -expm1(-G.θ) * exp(-t))
 end
 ϕ⁻¹(G::FrankGenerator, t) = G.θ > 0 ? LogExpFunctions.log1mexp(-G.θ) - LogExpFunctions.log1mexp(-t*G.θ) : -log(expm1(-t*G.θ)/expm1(-G.θ))
+
+# Taylor1-compatible ϕ / ϕ⁻¹: the θ>0 scalar forms above use `log1mexp`, which has
+# no `TaylorSeries.Taylor1` method, so the nested direct edge composition (a jet
+# over ϕ⁻¹_outer ∘ ϕ_inner) would fail when Frank appears in the tree. These add
+# the equivalent θ≠0 closed forms built only from `expm1`/`log1p` (all Taylor1-
+# compatible) — exactly the θ<0 branch forms, on a Taylor1 argument; the scalar
+# path is untouched. (A future TaylorSeries⊕LogExpFunctions extension giving
+# `log1mexp(::Taylor1)` would make these unnecessary; the implicit edge-composition
+# method already avoids them, using only the closed-form ϕ⁽ᵏ⁾.)
+ϕ(G::FrankGenerator, t::TaylorSeries.Taylor1{T}) where {T} = -log1p(exp(-t) * expm1(-T(G.θ))) / T(G.θ)
+ϕ⁻¹(G::FrankGenerator, t::TaylorSeries.Taylor1{T}) where {T} = -log(expm1(-T(G.θ) * t) / expm1(-T(G.θ)))
 𝒲₋₁(G::FrankGenerator, d::Int) = G.θ > 0 ? WilliamsonFromFrailty(Logarithmic(-G.θ), d) : @invoke 𝒲₋₁(G::Generator, d)
 frailty(G::FrankGenerator) = G.θ > 0 ? Logarithmic(-G.θ) : throw("The frank copula has no frailty when θ < 0")