Archimedean Copula Inference via Taylor-Mode AD

# Archimedean Copula Inference via Taylor-Mode AD
Source: [https://arxiv.org/html/2605.23134](https://arxiv.org/html/2605.23134)
Cambridge Yang Independent research@camyang\.com&Dongdong Li Harvard Medical School dongdongli@hsph\.harvard\.edu

###### Abstract

No existing nested Archimedean copula tool handles all three of \(a\) arbitrary per\-variable \(right\-\)censoring in survival analysis, \(b\) arbitrary nesting trees, and \(c\) exact parameter gradients\. Existing implementations handle only bivariate problems, low dimensional \(i\.e\.,d≤10d\\leq 10\) cases, two layers of nesting, or only hand\-derived copula nestings\. We presentacopula, a JAX\-native framework that, given any Archimedean generator—classical or neural—evaluates exact nested\-copula likelihoods and parameter gradients under arbitrary censoring masks in polynomial time\. The mechanism is polynomial powering of Taylor\-mode automatic differentiation output, which replaces per\-family hand\-derived partial Bell polynomial tables with a single differentiable computation that any user\-defined generator can drive\. We conduct extensive simulations to verify the correctness ofacopula\. We then demonstrate \(a\) per\-variable censoring on85,22985\{,\}229MIMIC\-IV ICU admissions in high dimensions withd=53d\{=\}53, fit by both classical Archimedean families and nested neural Archimedean copulas; \(b\) an 11\-sector hierarchical model on S&P 500 daily returns atd=98d\{=\}98; \(c\) family\-agnostic censored MLE across ten families, five of them with no prior implementation, on a retinopathy study; and \(d\) a∼650×\{\\sim\}650\\timesper\-density speedup over R’snacLLatd=35d\{=\}35, scaling quadratically tod=8,000d\{=\}8\{,\}000\.

## 1Introduction

Hierarchical multivariate data with partial observation is the default in modern survival, financial, and clinical risk settings\. Patients in a survival study with paired organs may never experience the event during follow\-up\(Huster et al\.,[1989](https://arxiv.org/html/2605.23134#bib.bib20)\)\. Stock returns cluster by industry sector and tail\-co\-move during crashes\. Intensive care unit \(ICU\) lab measurements group by clinical panel, and the analyst right\-censors each lab independently at hospital discharge\(Johnson et al\.,[2023](https://arxiv.org/html/2605.23134#bib.bib22)\)\. Each setting demands a model that captures the hierarchy with interpretable parameters—tail dependence and concordance—and accommodates an arbitrary censoring pattern per observation\.

##### Archimedean copulas\.

AcopulaC:\[0,1\]d→\[0,1\]C:\[0,1\]^\{d\}\\to\[0,1\]is a multivariate distribution on the unit cube withdduniform marginals; by Sklar’s theorem\(Sklar,[1959](https://arxiv.org/html/2605.23134#bib.bib39)\), every joint distribution on𝐮=\(u1,…,ud\)∈\[0,1\]d\\mathbf\{u\}=\(u\_\{1\},\\ldots,u\_\{d\}\)\\in\[0,1\]^\{d\}factors uniquely into its one\-dimensional marginals and a copula that captures their dependence\. This decomposition lets the analyst pair domain\-specific marginal models—e\.g\., Weibull\(Weibull,[1951](https://arxiv.org/html/2605.23134#bib.bib43)\)for survival times, GARCH\(Bollerslev,[1986](https://arxiv.org/html/2605.23134#bib.bib4)\)for returns—with an interpretable dependence structure whose parameters directly encode tail behaviour and concordance\.Nested Archimedean copulas\(McNeil,[2008](https://arxiv.org/html/2605.23134#bib.bib30); Hofert and Pham,[2013](https://arxiv.org/html/2605.23134#bib.bib17)\)extend copulas to hierarchical dependence by composing generator functions in a tree: variables under the same internal node share stronger dependence than variables across branches\.

##### The inference barrier\.

Fitting a nested copula requires its density—a mixed partial derivative of the copula cumulative distribution function \(CDF\), one differentiation per observed variable—through a tree of composed generators\(Hofert and Pham,[2013](https://arxiv.org/html/2605.23134#bib.bib17)\); three obstacles block existing tools\.

- •Per\-variable censoring\.We consider per\-variable right\-censoring under the survival setting\. Each subset of censored variables changes which mixed partial derivative of the copula CDF defines the likelihood, so ad=53d\{=\}53ICU dataset with thousands of distinct censoring patterns requires the framework to handle every per\-observation maskδ∈\{0,1\}d\\delta\\in\\\{0,1\\\}^\{d\}uniformly, whereδj=1\\delta\_\{j\}=1marks variablejjas observed\. Prior tools cover bivariate problems\(Sun and Ding,[2020](https://arxiv.org/html/2605.23134#bib.bib41); Emura and Chen,[2018](https://arxiv.org/html/2605.23134#bib.bib10)\)or shared censoring\(Liu et al\.,[2025b](https://arxiv.org/html/2605.23134#bib.bib29)\); none reachd\>30d\>30per\-variable\.
- •Arbitrary nesting and any generator\.R’scopulapackage\(Hofert et al\.,[2012](https://arxiv.org/html/2605.23134#bib.bib18)\)hard\-codes per\-family Stirling\-number tables for five generators and supports only two nesting layers; adding a sixth family—or any neural\-based generator such as the Deep Archimedean Copula family\(Ling et al\.,[2020](https://arxiv.org/html/2605.23134#bib.bib27); Liu et al\.,[2025b](https://arxiv.org/html/2605.23134#bib.bib29)\)—demands fresh hand derivation\.
- •Exact automatic parameter gradients\.Without automatic gradients, fitting ad=98d\{=\}98model relies on poorly\-scaling derivative\-free or finite\-difference methods, or Monte Carlo approximators\(Hofert et al\.,[2013](https://arxiv.org/html/2605.23134#bib.bib19); Ng et al\.,[2021](https://arxiv.org/html/2605.23134#bib.bib34)\)that are consistent but biased with finite samples\.

##### Our approach\.

We presentacopula, a JAX\-native framework that resolves all three obstacles\. The user specifies at minimum a copula generator—classical or neural—andacopulacomputes the nested copula likelihood for any tree shape and any censoring mask, with exact parameter gradients automatically\. The mechanism is polynomial powering of Taylor coefficients: Taylor\-mode AD\(Griewank and Walther,[2008](https://arxiv.org/html/2605.23134#bib.bib13); Radul et al\.,[2023](https://arxiv.org/html/2605.23134#bib.bib37)\)returns the higher\-order derivatives of each generator composition in one forward pass, and raising those coefficient vectors to successive powers via convolution recovers the partial Bell polynomials that the density formula ofHofert and Pham \([2013](https://arxiv.org/html/2605.23134#bib.bib17)\)demands without per\-family hand derivation, for any generator satisfying the standarddcd\_\{c\}\-monotonicity nesting condition\(McNeil,[2008](https://arxiv.org/html/2605.23134#bib.bib30)\)\.

##### Contributions\.

We make the following contributions:

1. 1\.Afamily\-agnostic, gradient\-based likelihood inference algorithmfor nested Archimedean copulas: given only a classical or neural generator \(e\.g\.,Deep Archimedean Copulas\(Ling et al\.,[2020](https://arxiv.org/html/2605.23134#bib.bib27); Liu et al\.,[2025b](https://arxiv.org/html/2605.23134#bib.bib29)\)\), the framework derives densities and gradients automatically\.
2. 2\.Per\-variable censoringas a first\-class operation: the likelihood algorithm handles arbitrary per\-observation masksδ∈\{0,1\}d\\delta\\in\\\{0,1\\\}^\{d\}directly through its core recurrence, lifting the bivariate\-only or single\-shared\-time restrictions of prior censoring\-aware tools\.
3. 3\.Acomplexity result:O​\(d3\)O\(d^\{3\}\)worst\-case for fixed\-depth trees andO​\(d2\)O\(d^\{2\}\)ford\\sqrt\{d\}\-decomposedfixed\-depth trees, where every non\-root subtree holds at mostO​\(d\)O\(\\sqrt\{d\}\)uncensored leaves; validated empirically up tod=8,000d\{=\}8\{,\}000\.
4. 4\.Anopen\-source JAX implementation,acopula, with a simple API for nested\-copula density and gradients, plus a per\-edgedcd\_\{c\}\-monotonicity diagnostic for the nesting\-validity condition\.
5. 5\.As a byproduct, conditional\-distribution sampling\(Hofert,[2008](https://arxiv.org/html/2605.23134#bib.bib15)\)of these copulas becomes practical, previously blocked by the high\-order derivatives; due to space, we defer to[Appendix˜F](https://arxiv.org/html/2605.23134#A6)\.

##### Related work\.

We positionacopulaagainst four neighbouring lines of work, each summarised below;[Table˜1](https://arxiv.org/html/2605.23134#S1.T1)aggregates the resulting capability gaps in tabular form\.

Table 1:Capability comparison across nested\-copula and copula\-survival tools \(Rcopula\(Hofert et al\.,[2012](https://arxiv.org/html/2605.23134#bib.bib18),[2013](https://arxiv.org/html/2605.23134#bib.bib19)\), HAC\(Okhrin and Ristig,[2014](https://arxiv.org/html/2605.23134#bib.bib35)\), Copulas\.jl\(Laverny and Jimenez,[2024](https://arxiv.org/html/2605.23134#bib.bib23)\), HACSurv†\(Liu et al\.,[2025b](https://arxiv.org/html/2605.23134#bib.bib29)\), DCSurvival\(Foomani et al\.,[2023](https://arxiv.org/html/2605.23134#bib.bib11)\), CopulaCenR\(Sun and Ding,[2020](https://arxiv.org/html/2605.23134#bib.bib41)\), and the Vine implementationrvinecopulib\(Nagler and Vatter,[2023](https://arxiv.org/html/2605.23134#bib.bib32)\)\)\. Testedddreports each tool’s largest successful attempt at a nested real\-data example, including in this work\.RcopulaHACCopulas\.jlHACSurv†DCSurvivalCopulaCenRVineOursUser\-defined classical generator×\\times×\\times✓\\checkmarkN/A×\\times×\\timesN/A✓\\checkmarkDeep Archimedean generators×\\times×\\times×\\times✓\\checkmark✓\\checkmark×\\timesN/A✓\\checkmarkSampling✓\\checkmark✓\\checkmark✓\\checkmarkN/A✓\\checkmark✓\\checkmark✓\\checkmark✓\\checkmarkNesting depth \(exact density\)2anyN/A2N/AN/AN/AanyAutomatic density gradients×\\times×\\times×\\times✓\\checkmark✓\\checkmark×\\times×\\times✓\\checkmarkFamily\-agnostic density×\\times×\\times×\\times×\\times✓\\checkmarkN/A×\\times✓\\checkmarkPer\-variable censoringnonenonenonecompetingbivar\.bivar\.noneper\-variableTesteddd\(nested real\-data\)3535∼10\{\\sim\}10N/A66222298989898- •Classical R\-package ecosystem\.The Rcopula\(Hofert et al\.,[2012](https://arxiv.org/html/2605.23134#bib.bib18),[2013](https://arxiv.org/html/2605.23134#bib.bib19)\)and HAC\(Okhrin and Ristig,[2014](https://arxiv.org/html/2605.23134#bib.bib35)\)packages compute generator derivatives via hand\-coded closed forms for five families, restrict every edge to same\-family pairs, and offer no analytic parameter gradients\. The simulated maximum likelihood estimator\(Hofert et al\.,[2013](https://arxiv.org/html/2605.23134#bib.bib19)\)relies on Laplace\-transform frailty sampling, with Monte Carlo variance that grows with dimension; our exact\-likelihood approach attains∼100×\{\\sim\}100\\timesspeedup at matched accuracy \([Appendix˜Q](https://arxiv.org/html/2605.23134#A17)\), and our deterministic\-density comparison against R’snacLLreaches∼650×\{\\sim\}650\\timesatd=35d\{=\}35before R aborts \([Section˜5\.2](https://arxiv.org/html/2605.23134#S5.SS2)\)\.
- •Censoring\-aware copula models\.CopulaCenR\(Sun and Ding,[2020](https://arxiv.org/html/2605.23134#bib.bib41)\)fits censored bivariate and trivariate copulas, including interval\-censoring via2d2^\{d\}\-corner inclusion–exclusion of the survival copula;Li et al\. \([2020](https://arxiv.org/html/2605.23134#bib.bib24),[2023](https://arxiv.org/html/2605.23134#bib.bib25)\)handle right\-censoring using brute\-force symbolic differentiation\. Both restrict to low ambient dimension\. HACSurv\(Liu et al\.,[2025b](https://arxiv.org/html/2605.23134#bib.bib29)\)fits hierarchical Archimedean copulas for competing risks where each observation reduces to one event time and the likelihood is a single first\-order partial of the nested CDF\. We handle the high\-order mixed partials that arise with per\-variable right\-censoring in the survival sense, and discuss interval\-censoring in[Section˜6](https://arxiv.org/html/2605.23134#S6)\.
- •Deep and probabilistic ML copulas\.Recent neural\-copula work\(Foomani et al\.,[2023](https://arxiv.org/html/2605.23134#bib.bib11); Liu et al\.,[2025a](https://arxiv.org/html/2605.23134#bib.bib28)\)learns flexible generators for bivariate dependent censoring and discrete\-diffusion joint sampling; we are complementary, targeting exact nested likelihoods at moderate\-to\-highddwith per\-variable censoring\. Vine copulas\(Aas et al\.,[2009](https://arxiv.org/html/2605.23134#bib.bib1); Nagler and Vatter,[2023](https://arxiv.org/html/2605.23134#bib.bib32); Stöber et al\.,[2013](https://arxiv.org/html/2605.23134#bib.bib40)\)and hierarchical Kendall copulas\(Brechmann,[2014](https://arxiv.org/html/2605.23134#bib.bib6)\)are density\-first: the censored likelihood requires a per\-observation multivariate integral over the censored block, exponential in its size\. METIC\(Chen et al\.,[2025](https://arxiv.org/html/2605.23134#bib.bib7)\)handles this at smallddvia a stochastic\-numerical scheme; no vine implementation reaches per\-variable censoring atd=53d\{=\}53\([Appendices˜U](https://arxiv.org/html/2605.23134#A21)and[V](https://arxiv.org/html/2605.23134#A22)\)\.
- •Higher\-order AD and differentiable probabilistic systems\.Our truncated\-Taylor\-series computation uses JAX’s Taylor\-mode AD primitive\(Bettencourt et al\.,[2019](https://arxiv.org/html/2605.23134#bib.bib3)\), also formalized in the Weil\-algebraic framework ofBetancourt \([2018](https://arxiv.org/html/2605.23134#bib.bib2)\): order\-ddderivatives in a single forward pass on a truncated polynomial algebra \([Section˜3\.2](https://arxiv.org/html/2605.23134#S3.SS2.SSS0.Px1)\)\.Lin \([2024](https://arxiv.org/html/2605.23134#bib.bib26)\)extends the JVP\-rule protocol to functional derivatives, complementary to the Taylor\-mode AD we use here\. NumPyro\(Phan et al\.,[2019](https://arxiv.org/html/2605.23134#bib.bib36)\)exposes first\-order but not higher\-order gradients required for nested copulas; Gen\-AC\(Ng et al\.,[2021](https://arxiv.org/html/2605.23134#bib.bib34)\)learns Laplace transforms with Monte Carlo variance that grows with derivative order\.

##### Roadmap\.

[Sections˜2](https://arxiv.org/html/2605.23134#S2),[3](https://arxiv.org/html/2605.23134#S3)and[4](https://arxiv.org/html/2605.23134#S4)formulate the problem, present the key algorithm, and describe the JAX implementation\.[Section˜5](https://arxiv.org/html/2605.23134#S5)validates correctness on simulation, characterises wall\-clock scaling against R’snacLLwith a∼650×\{\\sim\}650\\timesspeedup atd=35d\{=\}35before R aborts, and runs four real\-data experiments: MIMIC\-IV atd=53d\{=\}53with classical and neural Archimedean generators under per\-variable censoring; S&P 500 atd=98d\{=\}98—a nested likelihood beyond R’s reach; diabetic retinopathy across ten Archimedean families, five without any prior censored implementation; and end\-to\-end reproduction of HACSurv’s Framingham competing\-risks pipeline\.[Section˜6](https://arxiv.org/html/2605.23134#S6)addresses validity under mixed families and amortized JIT compile cost across observations\.

## 2Background and Notation

10ψ0\\psi\_\{0\}10ψ1\\psi\_\{1\}10ψ2\\psi\_\{2\}u1u\_\{1\}u2u\_\{2\}u4u\_\{4\}u5u\_\{5\}u6u\_\{6\}u3u\_\{3\}strongstrongFigure 1:Nested Archimedean copula on six leaves\. Nodes carry generators forming copula CDFs of their children; internal subtrees capture within\-group dependence\. Leaves are variables with marginal densities; the hatched leafu3u\_\{3\}is censored\.AcopulaC:\[0,1\]d→\[0,1\]C\\colon\[0,1\]^\{d\}\\to\[0,1\]captures the dependence amongdduniform random variables, independent of their marginal distributions\(Sklar,[1959](https://arxiv.org/html/2605.23134#bib.bib39); Embrechts et al\.,[2002](https://arxiv.org/html/2605.23134#bib.bib9); Nelsen,[2006](https://arxiv.org/html/2605.23134#bib.bib33)\)\. AnArchimedean copulaencodes that dependence through a singlegeneratorψ:\[0,∞\)→\[0,1\]\\psi\\colon\[0,\\infty\)\\to\[0,1\]—a continuous, strictly decreasing function withψ​\(0\)=1\\psi\(0\)\{=\}1that isdd\-monotone, equivalent to convex ford=2d\{=\}2and ensured by complete monotonicity for anydd\(McNeil and Nešlehová,[2009](https://arxiv.org/html/2605.23134#bib.bib31)\), producing the CDF:

C​\(𝐮\)=ψ​\(∑j=1dψ−1​\(uj\)\),𝐮∈\[0,1\]d\.C\(\\mathbf\{u\}\)=\\psi\\\!\\Bigl\(\\sum\_\{j=1\}^\{d\}\\psi^\{\-1\}\(u\_\{j\}\)\\Bigr\),\\qquad\\mathbf\{u\}\\in\[0,1\]^\{d\}\.\(2\)Nested Archimedean copulas\(McNeil,[2008](https://arxiv.org/html/2605.23134#bib.bib30)\)relax the flat exchangeability of[Equation˜2](https://arxiv.org/html/2605.23134#S2.E2)by arranging generators in a tree \([Figure˜1](https://arxiv.org/html/2605.23134#S2.F1)\)\. Each nodevvcarries its own generatorψv\\psi\_\{v\}and the CDF follows a recursive composition over its childrencc:

Cv​\(𝐮v\)=ψv​\(∑c∈children​\(v\)ψv−1​\(Cc​\(𝐮c\)\)\)\.C\_\{v\}\(\\mathbf\{u\}\_\{v\}\)=\\psi\_\{v\}\\\!\\Bigl\(\\sum\_\{c\\in\\mathrm\{children\}\(v\)\}\\psi\_\{v\}^\{\-1\}\\\!\\bigl\(C\_\{c\}\(\\mathbf\{u\}\_\{c\}\)\\bigr\)\\Bigr\)\.\(3\)[Equation˜3](https://arxiv.org/html/2605.23134#S2.E3)defines a valid copula whenever each compositionψouter−1∘ψinner\\psi\_\{\\mathrm\{outer\}\}^\{\-1\}\\circ\\psi\_\{\\mathrm\{inner\}\}is aBernstein function—one with completely monotone derivative—a sufficient condition due toMcNeil \([2008](https://arxiv.org/html/2605.23134#bib.bib30)\);[Section˜6](https://arxiv.org/html/2605.23134#S6)covers same\-family closed\-form ordering and the mixed\-family sign\-alternation diagnostic\.

##### Likelihood with per\-variable censoring\.

Fitting requires the density, so we move from[Equation˜3](https://arxiv.org/html/2605.23134#S2.E3)to its mixed partials and the censoring patterns they must accommodate\. Each observationiicarries acensoring maskδ\(i\)∈\{0,1\}d\\delta^\{\(i\)\}\\in\\\{0,1\\\}^\{d\}withδj\(i\)=1\\delta^\{\(i\)\}\_\{j\}=1marking variablejjuncensored and\|δ\|:=∑jδj\|\\delta\|:=\\sum\_\{j\}\\delta\_\{j\}counting uncensored variables\. Without loss of generality, we consider right\-censoring in the survival sense where marginals areuj=Sj​\(tj\):=Pr⁡\(Tj\>tj\)u\_\{j\}=S\_\{j\}\(t\_\{j\}\):=\\Pr\(T\_\{j\}\>t\_\{j\}\)and the copulaCCencodesPr⁡\(T1\>t1,…,Td\>td\)\\Pr\(T\_\{1\}\>t\_\{1\},\\ldots,T\_\{d\}\>t\_\{d\}\), with right\-censored observationsTj\>cjT\_\{j\}\>c\_\{j\}entering asuj=Sj​\(cj\)u\_\{j\}=S\_\{j\}\(c\_\{j\}\)\. The likelihood is the order\-\|δ\|\|\\delta\|mixed partial in the uncensored variables:

cδ​\(𝐮\)=∂\|δ\|C∏j:δj=1∂uj​\(𝐮\),ℓ​\(𝜽\)=∑i=1nlog⁡cδ\(i\)​\(𝐮\(i\);𝜽\),c\_\{\\delta\}\(\\mathbf\{u\}\)=\\frac\{\\partial^\{\|\\delta\|\}C\}\{\\prod\_\{j:\\delta\_\{j\}=1\}\\partial u\_\{j\}\}\(\\mathbf\{u\}\),\\qquad\\ell\(\\boldsymbol\{\\theta\}\)=\\sum\_\{i=1\}^\{n\}\\log c\_\{\\delta^\{\(i\)\}\}\(\\mathbf\{u\}^\{\(i\)\};\\,\\boldsymbol\{\\theta\}\),\(4\)with sample sizenn, observation𝐮\(i\)∈\[0,1\]d\\mathbf\{u\}^\{\(i\)\}\\in\[0,1\]^\{d\}, and𝜽\\boldsymbol\{\\theta\}stacking all generator parameters\. Two prior routes evaluatecδc\_\{\\delta\}: hand\-deriving a closed form perδ\\delta—infeasible at2d2^\{d\}patterns—or nesting\|δ\|\|\\delta\|first\-order AD calls atO​\(2\|δ\|\)O\(2^\{\|\\delta\|\}\)from intermediate trace blow\-up\(Betancourt,[2018](https://arxiv.org/html/2605.23134#bib.bib2)\); neither scales pastd≈10d\\approx 10\. R’scopulapackage, the dominant alternative, uses hand\-derived Stirling\-number tables for five families\(Hofert et al\.,[2012](https://arxiv.org/html/2605.23134#bib.bib18); Hofert and Pham,[2013](https://arxiv.org/html/2605.23134#bib.bib17)\), avoiding the exponential cost but covering only single\-family two\-level nesting, without censored likelihood or analytic gradients\. We assume covariate\-conditional independent censoring throughout;[Section˜6](https://arxiv.org/html/2605.23134#S6)returns to identifiability\.

##### TheHofert and Pham \([2013](https://arxiv.org/html/2605.23134#bib.bib17)\)density\.

Hofert and Pham \([2013](https://arxiv.org/html/2605.23134#bib.bib17)\)expressed the nested densitycδc\_\{\\delta\}as a tree\-structured sum involvingpartial Bell polynomialsBn,kB\_\{n,k\}of higher\-order derivatives of each parent–child edge compositionhv​c=ψv−1∘ψch\_\{vc\}=\\psi\_\{v\}^\{\-1\}\\circ\\psi\_\{c\}\(writingψ\(k\)\\psi^\{\(k\)\}for thekk\-th derivative throughout\):

cδ​\(𝐮\)=\(∑k=1\|δ\|βk\(r\)⋅ψr\(k\)​\(tr\)\)⋅∏ℓ∈leavesδℓ=1\(ψpar​\(ℓ\)−1\)′​\(uℓ\),c\_\{\\delta\}\(\\mathbf\{u\}\)=\\biggl\(\\sum\_\{k=1\}^\{\|\\delta\|\}\\beta\_\{k\}^\{\(r\)\}\\cdot\\psi\_\{r\}^\{\(k\)\}\(t\_\{r\}\)\\biggr\)\\cdot\\\!\\\!\\prod\_\{\\begin\{subarray\}\{c\}\\ell\\in\\mathrm\{leaves\}\\\\ \\delta\_\{\\ell\}=1\\end\{subarray\}\}\\\!\\\!\\bigl\(\\psi\_\{\\mathrm\{par\}\(\\ell\)\}^\{\-1\}\\bigr\)^\{\\prime\}\(u\_\{\\ell\}\),\(5\)withrrthe root,par​\(ℓ\)\\mathrm\{par\}\(\\ell\)the parent of leafℓ\\ell,tv:=∑c∈children​\(v\)ψv−1​\(Cc​\(𝐮c\)\)t\_\{v\}:=\\sum\_\{c\\in\\mathrm\{children\}\(v\)\}\\psi\_\{v\}^\{\-1\}\(C\_\{c\}\(\\mathbf\{u\}\_\{c\}\)\)the inverse\-generator argument at nodevv, and weightsβk\(r\)\\beta\_\{k\}^\{\(r\)\}defined recursively in[Section˜3\.2](https://arxiv.org/html/2605.23134#S3.SS2)—[Section˜A\.1](https://arxiv.org/html/2605.23134#A1.SS1)gives theBn,kB\_\{n,k\}formula\. Evaluating these Bell polynomials for user\-defined generators without hand\-derived tables is the obstacle[Section˜3](https://arxiv.org/html/2605.23134#S3)resolves by powering Taylor coefficients\.

##### Taylor\-mode AD\.

Our algorithm needs univariate higher\-order derivatives of each edge compositionh=ψv−1∘ψch=\\psi\_\{v\}^\{\-1\}\\circ\\psi\_\{c\}\.Taylor\-modeAD\(Griewank and Walther,[2008](https://arxiv.org/html/2605.23134#bib.bib13)\), orjet propagation, returns the factorial\-scaled derivative vectorJet​\(f,x,n\)=\(f​\(x\),f′​\(x\)/1\!,…,f\(n\)​\(x\)/n\!\)\\textsc\{Jet\}\(f,x,n\)=\\bigl\(f\(x\),f^\{\\prime\}\(x\)/1\!,\\ldots,f^\{\(n\)\}\(x\)/n\!\\bigr\)in one forward pass atO​\(n2\)O\(n^\{2\}\), versusO​\(2n\)O\(2^\{n\}\)fornnnested first\-order AD calls\. Our implementation builds on a JIT\-friendly variant of JAX’sjet\(Radul et al\.,[2023](https://arxiv.org/html/2605.23134#bib.bib37)\);[Section˜3](https://arxiv.org/html/2605.23134#S3)treats it as a black box\.

## 3Method

We walk through a concrete example exhibiting every algorithmic ingredient \([Section˜3\.1](https://arxiv.org/html/2605.23134#S3.SS1)\), state the bottom\-up recursion \([Section˜3\.2](https://arxiv.org/html/2605.23134#S3.SS2); full Bell\-polynomial derivation in[Sections˜A\.1](https://arxiv.org/html/2605.23134#A1.SS1)and[A\.3](https://arxiv.org/html/2605.23134#A1.SS3)\), give theO​\(d3\)O\(d^\{3\}\)complexity, and treat the per\-variable censoring mechanism\.

ψ0\\psi\_\{0\}ψ1\\psi\_\{1\}u1u\_\{1\}u2u\_\{2\}u3u\_\{3\}### 3\.1Worked example:d=3d\{=\}3nested copula with one censored leaf

Consider the tree on the right: a rootψ0\\psi\_\{0\}with sectorψ1\\psi\_\{1\}over uncensored leavesu1,u2u\_\{1\},u\_\{2\}and a right\-censored siblingu3u\_\{3\}\. The censoring maskδ=\(1,1,0\)\\delta=\(1,1,0\)gives an order\-2 mixed partial∂2C/∂u1​∂u2\\partial^\{2\}C/\\partial u\_\{1\}\\partial u\_\{2\}at𝐮\\mathbf\{u\}\. Each leaf carries anα\\alpha\-vector:α\(u1\)=α\(u2\)=\[0,1\]\\alpha^\{\(u\_\{1\}\)\}=\\alpha^\{\(u\_\{2\}\)\}=\[0,1\]andα\(u3\)=\[1\]\\alpha^\{\(u\_\{3\}\)\}=\[1\], the Cauchy\-product identity\. Four operations propagate these vectors upward and yield the density\.

1. 1\.Cauchy product at the sector\.Combine the uncensored leaves:β\(ψ1\)=α\(u1\)∗α\(u2\)=\[0,0,1\]\\beta^\{\(\\psi\_\{1\}\)\}=\\alpha^\{\(u\_\{1\}\)\}\\ast\\alpha^\{\(u\_\{2\}\)\}=\[0,0,1\], whereβ2\(ψ1\)=1\\beta\_\{2\}^\{\(\\psi\_\{1\}\)\}=1encodes the two differentiationsψ1\\psi\_\{1\}contributes upstream\.
2. 2\.Bell transform across the edgeψ0←ψ1\\psi\_\{0\}\\leftarrow\\psi\_\{1\}\.Convertβ\(ψ1\)\\beta^\{\(\\psi\_\{1\}\)\}into a child\-wiseα\\alpha\-vector at the root using higher\-order derivatives of the edge compositionh=ψ0−1∘ψ1h=\\psi\_\{0\}^\{\-1\}\\circ\\psi\_\{1\}att1=ψ1−1​\(u1\)\+ψ1−1​\(u2\)t\_\{1\}=\\psi\_\{1\}^\{\-1\}\(u\_\{1\}\)\+\\psi\_\{1\}^\{\-1\}\(u\_\{2\}\):αk\(ψ0,ψ1\)=∑j=k2βj\(ψ1\)​Bj,k​\(h′​\(t1\),h′′​\(t1\),…\)\\alpha\_\{k\}^\{\(\\psi\_\{0\},\\psi\_\{1\}\)\}=\\sum\_\{j=k\}^\{2\}\\beta\_\{j\}^\{\(\\psi\_\{1\}\)\}B\_\{j,k\}\(h^\{\\prime\}\(t\_\{1\}\),h^\{\\prime\\prime\}\(t\_\{1\}\),\\ldots\)\. Onlyj=2j=2contributes, givingα1\(ψ0,ψ1\)=h′′​\(t1\)\\alpha\_\{1\}^\{\(\\psi\_\{0\},\\psi\_\{1\}\)\}=h^\{\\prime\\prime\}\(t\_\{1\}\)andα2\(ψ0,ψ1\)=\(h′​\(t1\)\)2\\alpha\_\{2\}^\{\(\\psi\_\{0\},\\psi\_\{1\}\)\}=\(h^\{\\prime\}\(t\_\{1\}\)\)^\{2\}\.
3. 3\.Censored sibling\.The leafu3u\_\{3\}contributesα\(u3\)=\[1\]\\alpha^\{\(u\_\{3\}\)\}=\[1\], which leavesα\(ψ0,ψ1\)\\alpha^\{\(\\psi\_\{0\},\\psi\_\{1\}\)\}unshifted under the Cauchy product at the root:β\(ψ0\)=α\(ψ0,ψ1\)∗\[1\]=\[0,h′′​\(t1\),\(h′​\(t1\)\)2\]\\beta^\{\(\\psi\_\{0\}\)\}=\\alpha^\{\(\\psi\_\{0\},\\psi\_\{1\}\)\}\\ast\[1\]=\[0,h^\{\\prime\\prime\}\(t\_\{1\}\),\(h^\{\\prime\}\(t\_\{1\}\)\)^\{2\}\]\.
4. 4\.Density assembly\.The Hofert formula[Equation˜5](https://arxiv.org/html/2605.23134#S2.E5)gives cδ​\(𝐮\)=\[ψ0′​\(tr\)​h′′​\(t1\)\+ψ0′′​\(tr\)​\(h′​\(t1\)\)2\]⋅\(ψ1−1\)′​\(u1\)⋅\(ψ1−1\)′​\(u2\),c\_\{\\delta\}\(\\mathbf\{u\}\)=\\bigl\[\\psi\_\{0\}^\{\\prime\}\(t\_\{r\}\)\\,h^\{\\prime\\prime\}\(t\_\{1\}\)\+\\psi\_\{0\}^\{\\prime\\prime\}\(t\_\{r\}\)\\,\(h^\{\\prime\}\(t\_\{1\}\)\)^\{2\}\\bigr\]\\cdot\(\\psi\_\{1\}^\{\-1\}\)^\{\\prime\}\(u\_\{1\}\)\\cdot\(\\psi\_\{1\}^\{\-1\}\)^\{\\prime\}\(u\_\{2\}\),\(6\)withtr=ψ0−1​\(C1​\(u1,u2\)\)\+ψ0−1​\(u3\)t\_\{r\}=\\psi\_\{0\}^\{\-1\}\\\!\\bigl\(C\_\{1\}\(u\_\{1\},u\_\{2\}\)\\bigr\)\+\\psi\_\{0\}^\{\-1\}\(u\_\{3\}\)andC1C\_\{1\}the inner copula generated byψ1\\psi\_\{1\}\. The censored leafu3u\_\{3\}enterstrt\_\{r\}via its inverse\-generator term but drops out of the leaf product sinceδ3=0\\delta\_\{3\}=0\.

Algorithm 1Main algorithm\.- Require:Copula tree of valid nesting at rootrr, observations𝐮\\mathbf\{u\}
- Ensure:Log\-densitylog⁡c​\(𝐮\)\\log c\(\\mathbf\{u\}\)
- 1:Computetv=∑cψv−1​\(Cc\)t\_\{v\}=\\sum\_\{c\}\\psi\_\{v\}^\{\-1\}\(C\_\{c\}\)at every nodevv\(bottom\-up\)
- 2:foreach nodevv, from leaves to rootdo
- 3:foreach childccofvvdo
- 4:ifccis an uncensored leafthen
- 5:α\(c\)←\[0,1\]\\alpha^\{\(c\)\}\\leftarrow\[0,\\;1\] ⊳\\trianglerightone derivative, coefficient 1
- 6:elseifccis a censored leafthen
- 7:α\(c\)←\[1\]\\alpha^\{\(c\)\}\\leftarrow\[1\] ⊳\\trianglerightidentity under convolution
- 8:else
- 9:p←Jet​\(ψv−1∘ψc,tc,dc\)p\\leftarrow\\textsc\{Jet\}\\bigl\(\\psi\_\{v\}^\{\-1\}\\circ\\psi\_\{c\},\\;t\_\{c\},\\;d\_\{c\}\\bigr\) ⊳\\trianglerightTaylor coeffs ofhv​ch\_\{vc\}, computed via implicit equation solve in[Section˜A\.4](https://arxiv.org/html/2605.23134#A1.SS4)
- 10:q←pq\\leftarrow p ⊳\\trianglerightqqaccumulatespkp^\{k\}; starts atp1p^\{1\}
- 11:fork=1,…,dck=1,\\ldots,d\_\{c\}do
- 12:αk\(c\)←\(1/k\!\)​⟨β~\(c\),q⟩\\alpha^\{\(c\)\}\_\{k\}\\leftarrow\(1/k\!\)\\,\\langle\\tilde\{\\beta\}^\{\(c\)\},\\;q\\rangle ⊳\\trianglerightreads offBj,kB\_\{j,k\}
- 13:q←q∗pq\\leftarrow q\*p ⊳\\trianglerightq=pk\+1q=p^\{k\+1\}
- 14:β\(v\)←α\(c1\)∗α\(c2\)∗⋯\\beta^\{\(v\)\}\\leftarrow\\alpha^\{\(c\_\{1\}\)\}\*\\alpha^\{\(c\_\{2\}\)\}\*\\cdots ⊳\\trianglerightconvolve children
- 15:q1,…,qd←Jet​\(ψr,tr,d\)q\_\{1\},\\ldots,q\_\{d\}\\leftarrow\\textsc\{Jet\}\(\\psi\_\{r\},t\_\{r\},d\) ⊳\\trianglerightroot Taylor coeffs
- 16:returnlog⁡\(\(∑kβk\(r\)⋅k\!​qk\)⋅∏ℓ:δℓ=1\(ψpar​\(ℓ\)−1\)′​\(uℓ\)\)\\log\\Bigl\(\\bigl\(\\sum\_\{k\}\\beta\_\{k\}^\{\(r\)\}\\cdot k\!\\,q\_\{k\}\\bigr\)\\cdot\\prod\_\{\\ell:\\,\\delta\_\{\\ell\}=1\}\(\\psi\_\{\\mathrm\{par\}\(\\ell\)\}^\{\-1\}\)^\{\\prime\}\(u\_\{\\ell\}\)\\Bigr\)

TheBell polynomial valuesin[Equation˜6](https://arxiv.org/html/2605.23134#S3.E6)—h′′​\(t1\)h^\{\\prime\\prime\}\(t\_\{1\}\)atk=1k\{=\}1and\(h′​\(t1\)\)2\(h^\{\\prime\}\(t\_\{1\}\)\)^\{2\}atk=2k\{=\}2—equaln\!/k\!n\!/k\!times the coefficients ofεn\\varepsilon^\{n\}inp​\(ε\)kp\(\\varepsilon\)^\{k\}\(heren=2n\{=\}2, with2\!/1\!=22\!/1\!\{=\}2fork=1k\{=\}1and2\!/2\!=12\!/2\!\{=\}1fork=2k\{=\}2\), wherep​\(ε\)=h′​\(t1\)​ε\+h′′​\(t1\)​ε2/2\+⋯p\(\\varepsilon\)=h^\{\\prime\}\(t\_\{1\}\)\\,\\varepsilon\+h^\{\\prime\\prime\}\(t\_\{1\}\)\\,\\varepsilon^\{2\}/2\+\\cdots\([Section˜A\.3](https://arxiv.org/html/2605.23134#A1.SS3)\)\. Taylor\-mode AD applied tohhatt1t\_\{1\}returns these coefficients ofppin one forward pass; raisingppto thekk\-th power by truncated convolution yields all required Bell polynomial values without combinatorial enumeration or per\-family tables\. The Cauchy product, the Bell transform, and the polynomial powering of Taylor coefficients are the entire algorithm\.

### 3\.2Algorithm

The general procedure repeats the worked example bottom\-up at every internal node\. For each childccof nodevv, wheremv=\|children​\(v\)\|m\_\{v\}=\|\\mathrm\{children\}\(v\)\|, theα\\alpha\-vectorα\(v,c\)\\alpha^\{\(v,c\)\}encodescc’s contribution to the parent density, and theβ\\beta\-vectorβ\(v\)\\beta^\{\(v\)\}collects the combined contribution of children ofvv; uncensored leaves initialiseα=\[0,1\]\\alpha=\[0,1\]and censored leaves useα=\[1\]\\alpha=\[1\], the Cauchy\-product identity\. Letδv\\delta\_\{v\}be the count of uncensored leaves in the subtree atvv\. Two operations propagate these vectors upward, theBell transform\([Equation˜7](https://arxiv.org/html/2605.23134#S3.E7)\) across each edgev←cv\\leftarrow cand theCauchy product\([Equation˜8](https://arxiv.org/html/2605.23134#S3.E8)\) combining siblings:

∀k=1,…,\|δc\|:αk\(v,c\)=∑j=k\|δc\|βj\(c\)​Bj,k​\(hv​c′​\(tc\),hv​c′′​\(tc\),…\),\\ignorespaces\\ignorespaces\\forall k=1,\\ldots,\|\\delta\_\{c\}\|\\colon\\\\ \\alpha\_\{k\}^\{\(v,c\)\}=\\sum\_\{j=k\}^\{\|\\delta\_\{c\}\|\}\\beta\_\{j\}^\{\(c\)\}\\;B\_\{j,k\}\\\!\\bigl\(h^\{\\prime\}\_\{vc\}\(t\_\{c\}\),h^\{\\prime\\prime\}\_\{vc\}\(t\_\{c\}\),\\ldots\\bigr\),\(7\)
∀k=0,…,\|δv\|:βk\(v\)=∑k1\+⋯\+kmv=k∀ki≥0αk1\(v,c1\)​⋯​αkmv\(v,cmv\)\.\\ignorespaces\\ignorespaces\\forall k=0,\\ldots,\|\\delta\_\{v\}\|\\colon\\\\ \\beta\_\{k\}^\{\(v\)\}=\\\!\\\!\\sum^\{\\forall k\_\{i\}\\geq 0\}\_\{k\_\{1\}\+\\cdots\+k\_\{m\_\{v\}\}=k\}\\\!\\\!\\alpha\_\{k\_\{1\}\}^\{\(v,c\_\{1\}\)\}\\cdots\\alpha\_\{k\_\{m\_\{v\}\}\}^\{\(v,c\_\{m\_\{v\}\}\)\}\.\(8\)
The recursion followsHofert and Pham \([2013](https://arxiv.org/html/2605.23134#bib.bib17)\); our contribution is evaluating the Bell polynomialsBj,kB\_\{j,k\}in[Equation˜7](https://arxiv.org/html/2605.23134#S3.E7)via polynomial powering of Taylor coefficients\. Withε\\varepsilona formal indeterminate andpv​c​\(ε\)=∑m≥1hv​c\(m\)​\(tc\)​εm/m\!p\_\{vc\}\(\\varepsilon\)=\\sum\_\{m\\geq 1\}h\_\{vc\}^\{\(m\)\}\(t\_\{c\}\)\\,\\varepsilon^\{m\}/m\!the exponential generating function of the edge derivatives, a generating\-function identity givesBn,k=\(n\!/k\!\)⋅Coeff⁡\(pv​c​\(ε\)k;εn\)B\_\{n,k\}=\(n\!/k\!\)\\cdot\\operatorname\{Coeff\}\\\!\\bigl\(p\_\{vc\}\(\\varepsilon\)^\{k\};\\;\\varepsilon^\{n\}\\bigr\), whereCoeff⁡\(q;εn\)\\operatorname\{Coeff\}\(q;\\,\\varepsilon^\{n\}\)extracts the coefficient ofεn\\varepsilon^\{n\}in the polynomialqq\. Applying Taylor\-mode AD tohv​ch\_\{vc\}attct\_\{c\}returns the coefficients ofpv​cp\_\{vc\}directly; iterated truncated convolutionq←q⋅pv​cq\\leftarrow q\\cdot p\_\{vc\}producespv​ckp\_\{vc\}^\{k\}\. Substituting back,[Equation˜7](https://arxiv.org/html/2605.23134#S3.E7)reduces to a single inner productαk\(v,c\)=\(1/k\!\)​∑j=k\|δc\|β~j\(c\)​\(qk\)j\\alpha\_\{k\}^\{\(v,c\)\}=\(1/k\!\)\\,\\sum\_\{j=k\}^\{\|\\delta\_\{c\}\|\}\\tilde\{\\beta\}^\{\(c\)\}\_\{j\}\\,\(q\_\{k\}\)\_\{j\}withβ~j:=j\!​βj\\tilde\{\\beta\}\_\{j\}:=j\!\\,\\beta\_\{j\}the rescaledβ\\beta\-vector and\(qk\)j\(q\_\{k\}\)\_\{j\}theεj\\varepsilon^\{j\}coefficient ofpv​ckp\_\{vc\}^\{k\}\.

[Algorithm˜1](https://arxiv.org/html/2605.23134#alg1)summarises the resulting bottom\-up traversal\. Every step is differentiable in the parameters, so first\-order AD of the log\-likelihood enables gradient\-based optimisers\. For numerical stability, the algorithm performs operations in the log\-domain \([Appendix˜B](https://arxiv.org/html/2605.23134#A2)\)\. Moreover, at[˜9:](https://arxiv.org/html/2605.23134#S3.I2.i9), rather than feeding the explicit compositionψv−1∘ψc\\psi\_\{v\}^\{\-1\}\\circ\\psi\_\{c\}toJet, the algorithm computespv​cp\_\{vc\}by solving the implicit equationψv​\(hv​c​\(t\)\)=ψc​\(t\)\\psi\_\{v\}\(h\_\{vc\}\(t\)\)=\\psi\_\{c\}\(t\)for its Taylor coefficients to avoid the catastrophic cancellation that arises whenψc​\(t\)\\psi\_\{c\}\(t\)underflows beforeψv−1\\psi\_\{v\}^\{\-1\}can recover it \([Section˜A\.4](https://arxiv.org/html/2605.23134#A1.SS4)\)\. In principle, the algorithm can also use arbitrary\-precision arithmetic to avoid these tricks at the cost of efficiency\.

##### Complexity

Table 2:Asymptotic cost by tree shape, governed by the largest internal\-child subtree leaf count\|δv\|\|\\delta\_\{v\}\|\. Any shape withmax⁡\|δv\|≤O​\(d\)\\max\|\\delta\_\{v\}\|\\leq O\(\\sqrt\{d\}\)at every level reaches theO​\(d2\)O\(d^\{2\}\)root\-jet cost\.†Chains do not appear in real copula data to our knowledge\. Derivation in[Section˜A\.6](https://arxiv.org/html/2605.23134#A1.SS6)\.LetT​\(𝒯\)T\(\\mathcal\{T\}\)denote the total runtime of[Algorithm˜1](https://arxiv.org/html/2605.23134#alg1)on nesting tree𝒯\\mathcal\{T\}\. Each non\-root nodevvpaysO​\(\|δv\|3\)O\(\|\\delta\_\{v\}\|^\{3\}\)for the polynomial powering its parent applies to it, which dominates theO​\(\|δv\|2\)O\(\|\\delta\_\{v\}\|^\{2\}\)Cauchy product at any internalvv, giving

T​\(𝒯\)=O​\(d2\+∑v≠root\|δv\|3\),T\(\\mathcal\{T\}\)\\;=\\;O\\\!\\Bigl\(d^\{2\}\+\\textstyle\\sum\_\{v\\neq\\mathrm\{root\}\}\|\\delta\_\{v\}\|^\{3\}\\Bigr\),\(9\)whered2d^\{2\}comes from the order\-ddjet at the root\. The largest subtree therefore dominates: total cost reaches theO​\(d2\)O\(d^\{2\}\)root\-jet cost whenever∑v≠r\|δv\|3=O​\(d2\)\\sum\_\{v\\neq r\}\|\\delta\_\{v\}\|^\{3\}=O\(d^\{2\}\), a condition met ford\\sqrt\{d\}\-decomposedtrees—those where every non\-root subtree satisfies\|δv\|=O​\(d\)\|\\delta\_\{v\}\|=O\(\\sqrt\{d\}\)at constant tree depth—since per\-level disjointness gives∑v​at level\|δv\|3=O​\(d2\)\\sum\_\{v\\text\{ at level\}\}\|\\delta\_\{v\}\|^\{3\}=O\(d^\{2\}\)\. Any subtree exceedingd\\sqrt\{d\}raises the total by its\|δv\|3\|\\delta\_\{v\}\|^\{3\}contribution\.[Table˜2](https://arxiv.org/html/2605.23134#S3.T2)lists the common cases, and[Section˜A\.6](https://arxiv.org/html/2605.23134#A1.SS6)derives the bound rigorously\. Each edge’s implicit\-equation solve costsO​\(\|δc\|3\)O\(\|\\delta\_\{c\}\|^\{3\}\), matching the Bell powering it feeds, so it is absorbable into the total bound \([Section˜A\.4](https://arxiv.org/html/2605.23134#A1.SS4)\)\.

In contrast, prior tools make weaker trade\-offs: R’scopularelies on hard\-coded Stirling tables for five families at two nesting levels, and R’s symbolic differentiationD\(\)\(Li et al\.,[2020](https://arxiv.org/html/2605.23134#bib.bib24),[2023](https://arxiv.org/html/2605.23134#bib.bib25)\)incurs a combinatorial complexity tractable only ford≲6d\\lesssim 6\.[Appendix˜D](https://arxiv.org/html/2605.23134#A4)discusses FFT and scan\-based accelerations for the polynomial\-powering step\.

##### Handling per\-variable censoring

Censoring a leafℓ\\ellsetsα\(ℓ\)=\[1\]\\alpha^\{\(\\ell\)\}=\[1\]instead of\[0,1\]\[0,1\], which leaves the Cauchy product unchanged on its level and reduces the effective dimension\|δv\|\|\\delta\_\{v\}\|at every ancestor\. The reduction propagates automatically: shorter Cauchy product, lower Bell\-polynomial degree, and reducedJetorder; per\-node cost scales asO​\(\|δv\|3\)O\(\|\\delta\_\{v\}\|^\{3\}\), so heavy censoring runs faster \([Table˜8](https://arxiv.org/html/2605.23134#A1.T8)\)\.

## 4Implementation

@copula

classClayton:

theta:float

defgenerator\(self,t\):

return\(\(1\+t\)\*\*

\(\-1\.0/self\.theta\)\)

\\par@copula

classAMH:

theta:float

defgenerator\(self,t\):

e=jnp\.exp\(\-t\)

return\(\(1\-self\.theta\)\*e/

\(1\-self\.theta\*e\)\)

defnested\_amh\_clayton\(p,u\):

returnAMH\(p\.theta\_outer\)\(

Clayton\(p\.theta\_inner\)\(

marginal\(Weibull\(1\.5,1\.0\),

obs=u\[i,j\],

censored=\(\(i,j\)==\(1,3\)\)\)

forjinrange\(5\)\)

foriinrange\(4\)\)

\\parp0=Params\(theta\_outer=0\.5,theta\_inner=2\.0\)

obs=load\_data\(\)

m=compile\_model\(nested\_amh\_clayton,template=p0\)

m\.ll\_fn\(obs,m\.flatten\(p0\)\)\#logdensity

jax\.grad\(m\.ll\_fn,argnums=1\)\(obs,m\.flatten\(p0\)\)\#exactgrad

A nested AMH/Clayton copula with 4 sectors of 5 Weibull\-marginal leaves, per\-variable censoring, and end\-to\-end differentiable inference\.acopulais an open\-source JAX\(Bradbury et al\.,[2018](https://arxiv.org/html/2605.23134#bib.bib5)\)package\.111Available at[https://github\.com/thisiscam/acopula](https://github.com/thisiscam/acopula)\.The user supplies a generator class; the framework derives the inverse, all higher\-order derivatives, the Bell\-polynomial coefficients, the density values, and the parameter gradients automatically\.[Figure˜2](https://arxiv.org/html/2605.23134#S4.F2)shows two custom families and a 4\-sector nested model with per\-variable censoring fitted in fewer than twenty Python lines\.

The@copuladecorator registers parameters and derives the generator inverse via symbolic inversion\(Radul et al\.,[2023](https://arxiv.org/html/2605.23134#bib.bib37)\); when symbolic inversion fails,acopulafalls back to bisection with implicit\-function\-theorem \(IFT\) gradients \([Appendix˜E](https://arxiv.org/html/2605.23134#A5)\)\.acopulacompiles a single likelihood program that handles every per\-observation censoring maskδ\(i\)\\delta^\{\(i\)\}without recompilation \([Section˜A\.6](https://arxiv.org/html/2605.23134#A1.SS6.SSS0.Px7)\)\. As is standard in AD\-based frameworks, users may override any of the framework’s automatic derivations with a hand\-derived alternative — a custom generator inverse, a closed\-form edge composition, or a per\-family Taylor\-coefficient routine in place ofJet— for performance or numerical stability\.

## 5Experiments

Table 3:Replicated parameter recovery on nested Clayton copulas \(500500replications,500500samples per replication;LLdenotes nesting depth, EmpSE the across\-replication SE, HessSE the observed\-Fisher SE\)\. Bias stays under2\.7%2\.7\\%of the true value; the empirical SE matches the observed\-Fisher SE\. Results suggest standard Hessian\-based confidence intervals on our exact\-likelihood fits are reliable\.We verify correctness on simulated data, characterise wall\-clock scaling against R’snacLL, then evaluate three real\-data benchmarks that each exercise a capability prior tools cannot deliver: per\-variable censoring atd=53d\{=\}53on MIMIC\-IV, nested likelihood atd=98d\{=\}98on S&P 500, and family\-agnostic censored MLE on retinopathy\. All runs share a single JAX code base in float64\. The scaling sweep, the recovery study, R comparison, and S&P 500 and retinopathy fit run on a single AMD Threadripper 3960X CPU core in single\-process mode for the apples\-to\-apples R comparison; the MIMIC\-IV experiments run on an NVIDIA RTX 3090 Ti GPU \([Table˜27](https://arxiv.org/html/2605.23134#A18.T27)\)\.

### 5\.1Correctness validation

We verify the implementation in three ways before turning to real data\. Theparameter\-recoveryexperiment samples from a nested Clayton model with known true parameters, refits by maximum likelihood, and compares the estimate to the truth\.[Table˜3](https://arxiv.org/html/2605.23134#S5.T3)reports recovery withR=500R\{=\}500Monte Carlo replications,d∈\{10,20,50\}d\\in\\\{10,20,50\\\}, and tree depth up to three: every parameter recovers to within2\.7%2\.7\\%bias, and the empirical SE matches the observed\-Fisher SE within2\.9%2\.9\\%, confirming that Hessian\-based inference on the Bell\-polynomial likelihood is well\-calibrated\. We verify density agreement with R’snacLLto<10−14\{<\}10^\{\-14\}and with Mathematica’s symbolic derivatives up tod=60d\{=\}60\([Appendix˜J](https://arxiv.org/html/2605.23134#A10)\)\. AD gradients match central differences to<10−8\{<\}10^\{\-8\}in[Appendix˜K](https://arxiv.org/html/2605.23134#A11);[Section˜L\.1](https://arxiv.org/html/2605.23134#A12.SS1)extends parameter recovery to Frank and Gumbel\.

### 5\.2Scaling: head\-to\-head against RnacLL

10110^\{1\}10210^\{2\}10310^\{3\}10410^\{4\}10−110^\{\-1\}10110^\{1\}10310^\{3\}R stops hereO​\(d2\)O\(d^\{2\}\)dimensionddms per density evaluationFigure 3:Wall\-clock scaling\. R only runs at 2\-levelK=5K\{=\}5; aborts pastd=40d\{=\}40on Clayton/Gumbel\. Byd=35d\{=\}35R is∼650×\{\\sim\}650\\timesslower\.acopulaextends across topologies R cannot represent and reachesO​\(d2\)O\(d^\{2\}\)atd≥6,400d\\geq 6\{,\}400\.A wall\-clock comparison against R’snacLLon matched single\-thread CPU conditions atn=1n\{=\}1per evaluation showsacopulagrowing polynomially inddwhilenacLLgrows roughly exponentially in the number of sectors: byd=35d\{=\}35R is already∼650×\{\\sim\}650\\timesslower per density for both Clayton and Gumbel, aborts pastd=40d\{=\}40for both \([Figure˜3](https://arxiv.org/html/2605.23134#S5.F3)left\); on the same recovery task, our exact gradients also beat R’s frailty\-sampled SMLE by∼100×\{\\sim\}100\\timesat matched accuracy \([Appendix˜Q](https://arxiv.org/html/2605.23134#A17)\)\.acopulaextends across topologies that R cannot represent, sincenacLLrestricts to nesting depth two with same\-family pairs; both 3\-level fixed sector sizeK=5K\{=\}5leaves per inner node and balancedd\\sqrt\{d\}reach the theoreticalO​\(d2\)O\(d^\{2\}\)asymptote atd∈\[6,400,8,000\]d\\in\[6\{,\}400,\\,8\{,\}000\]\([Figure˜3](https://arxiv.org/html/2605.23134#S5.F3)right,[Appendix˜O](https://arxiv.org/html/2605.23134#A15)\)\.

### 5\.3MIMIC\-IV: classical nested copulas atd=53d\{=\}53

The MIMIC\-IV ICU cohort exercises per\-variable censoring at a dimension no prior nested\-copula tool reaches\. We fit two\-level nested Archimedean copulas on85,22985\{,\}229ICU admissions from MIMIC\-IV\(Johnson et al\.,[2023](https://arxiv.org/html/2605.23134#bib.bib22)\)atd=53d\{=\}53lab tests\. Each lab variable records the time for its value to become abnormal since ICU entry, right\-censored at discharge; per\-variable censoring rates span from11\.8%11\.8\\%for red blood cells to95\.8%95\.8\\%for specific gravity, and every admission censors a different subset, producing thousands of distinct masks \([Section˜R\.1](https://arxiv.org/html/2605.23134#A18.SS1)gives cohort details\)\. We group the labs into seven clinical panels—CBC, CMP, differential, liver, coagulation, ABG, and Other—and follow the semiparametric two\-stage recipe\(Genest et al\.,[1995](https://arxiv.org/html/2605.23134#bib.bib12); Joe and Xu,[1996](https://arxiv.org/html/2605.23134#bib.bib21); Li et al\.,[2020](https://arxiv.org/html/2605.23134#bib.bib24),[2023](https://arxiv.org/html/2605.23134#bib.bib25)\): Kaplan–Meier marginals at stage 1, joint copula by L\-BFGS\-B with exact gradients at stage 2\.

Table 4:MIMIC\-IV held\-out evaluation demonstratingacopula’s capability to fit nested Archimedean copulas\. Top: classical families \([Section˜5\.3](https://arxiv.org/html/2605.23134#S5.SS3)\)\. Bottom: nested neural Archimedean generators \([Section˜5\.4](https://arxiv.org/html/2605.23134#S5.SS4)\)\.Clayton/Clayton attains the best held\-out NLL among classical same\-family pairs with essentially zero train\-to\-test gap \([Table˜4](https://arxiv.org/html/2605.23134#S5.T4), top\)\. Frank/Frank’s symmetric tail beats Gumbel’s upper\-tail\-only dependence, consistent with ICU labs co\-degrading at both extremes during sepsis and recovery\. The fitted Kendall’sτ\\taurecovers panel\-level dependence that matches known mechanistic coupling: differential atτ=0\.72\\tau\{=\}0\.72, where WBC counts sum to100%100\\%; coagulation atτ=0\.61\\tau\{=\}0\.61, where PT and INR measure the same clotting factor; and ABG atτ=0\.60\\tau\{=\}0\.60, governed by Henderson–Hasselbalch equilibrium\. The outer level ties the seven panels atτouter=0\.15\\tau\_\{\\mathrm\{outer\}\}\{=\}0\.15\. The AMH/Clayton case exercises cross\-family nesting using Hofert\([2010](https://arxiv.org/html/2605.23134#bib.bib16)\)’s closed\-form nesting validity constraint\.

### 5\.4Capability demonstrations on nested neural Archimedean copulas

Our framework encompasses any Archimedean generator, including neural generators\(Ling et al\.,[2020](https://arxiv.org/html/2605.23134#bib.bib27); Ng et al\.,[2021](https://arxiv.org/html/2605.23134#bib.bib34); Liu et al\.,[2025b](https://arxiv.org/html/2605.23134#bib.bib29)\)\. We exercise this capability with two experiments: We fit ACNet and Gen\-AC fit on MIMIC\-IV and exactly reproduce HACSurv’s pipeline\.

##### Generators\.

ACNet\(Ling et al\.,[2020](https://arxiv.org/html/2605.23134#bib.bib27)\)parameterisesψ\\psias a softmax\-convex cascade of positive exponentials, enumerating exponentially many mixture atoms through polynomially many parameters \([Appendix˜I](https://arxiv.org/html/2605.23134#A9)\)\.Gen\-AC\(Ng et al\.,[2021](https://arxiv.org/html/2605.23134#bib.bib34)\)definesψ\\psias a Laplace transform of an MLP\-pushforward frailty,ψ​\(t\)=𝔼U​\[exp⁡\(−exp⁡\(MLP​\(U\)\)​t\)\]\\psi\(t\)=\\mathbb\{E\}\_\{U\}\[\\exp\(\-\\exp\(\\mathrm\{MLP\}\(U\)\)\\,t\)\]withU∼Uniform​\(0,1\)U\\sim\\mathrm\{Uniform\}\(0,1\); HACSurv\(Liu et al\.,[2025b](https://arxiv.org/html/2605.23134#bib.bib29)\)reuses this generator under a competing\-risks loss\. Two\-level trees compose each outer with an ACNet\-style Bernstein inner per panel, preserving complete monotonicity \([Section˜G\.1](https://arxiv.org/html/2605.23134#A7.SS1)\)\.

##### MIMIC\-IV results\.

Every neural row beats the best classical fit, with train\-to\-test gaps under0\.010\.01nats \([Table˜4](https://arxiv.org/html/2605.23134#S5.T4), bottom\)\. ACNet at\(w0,w1\)=\(4,4\)\(w\_\{0\},w\_\{1\}\)\{=\}\(4,4\)reaches4\.3104\.310with7070parameters, and wider or deeper ACNet variants give marginal further gains\. Gen\-AC at its HACSurv\-sized1,1421\{,\}142\-parameter configuration plateaus at4\.9144\.914, only0\.040\.04nats better than Clayton despite carrying16×16\\timesACNet’s parameters\. We isolate the cause as the smoothness prior that a small MLP imposes \([Appendix˜H](https://arxiv.org/html/2605.23134#A8)\)\. With Fourier\-feature\(Tancik et al\.,[2020](https://arxiv.org/html/2605.23134#bib.bib42)\)inputu↦\[cos⁡\(2​π​k​u\),sin⁡\(2​π​k​u\)\]k=1Ku\\mapsto\[\\cos\(2\\pi ku\),\\sin\(2\\pi ku\)\]\_\{k=1\}^\{K\}, Gen\-AC\-Fourier improves to4\.4934\.493at3\.6​k3\.6\\mathrm\{k\}parameters but still trails ACNet\. The Gen\-AC\-Fourier\-∫260\{\\int\}\_\{260\}row integrates over the frailty distribution by numerical integration instead of the Monte Carlo approximation used byNg et al\. \([2021](https://arxiv.org/html/2605.23134#bib.bib34)\), eliminating the estimation bias \([Sections˜6](https://arxiv.org/html/2605.23134#S6)and[Q](https://arxiv.org/html/2605.23134#A17)\)\.

##### Framingham reproduction\.

Table 5:Our port matchesLiu et al\. \([2025b](https://arxiv.org/html/2605.23134#bib.bib29), Table 7\); entries×102\\times 10^\{2\}\.HACSurv\(Liu et al\.,[2025b](https://arxiv.org/html/2605.23134#bib.bib29)\)restricts to competing risks with one observed event time, where the likelihood is a first\-order partial of the CDF, a strict subset of the mixed partials that our framework handles\.[Table˜5](https://arxiv.org/html/2605.23134#S5.T5)shows that we exactly reproduce their results on the Framingham dataset\.

### 5\.5S&P 500: nested likelihood atd=98d\{=\}98

Table 6:S&P 500 atd=98d\{=\}98: classical nested Archimedean families and an R\-vine\(Nagler and Vatter,[2023](https://arxiv.org/html/2605.23134#bib.bib32)\)for context\. Gumbel best of the three classical pairs \([Appendix˜U](https://arxiv.org/html/2605.23134#A21)\)\.The S&P 500 fit serves as a real\-data engineering check at a dimension prior tools cannot reach\. We model daily returns over1,2531\{,\}253trading days as an 11\-sector nested copulas atd∈\{50,75,98\}d\\in\\\{50,75,98\\\}with1212parameters; returns convert to pseudo\-observations via empirical\-CDF rank transform, so reported likelihoods isolate dependence from marginals\. Gumbel attains the best classical AIC atd=98d\{=\}98\([Table˜6](https://arxiv.org/html/2605.23134#S5.T6)\); per\-sectorθ\\thetatrack within\-sector dependence with Energy and Utilities strongest atθ=1\.88\\theta\{=\}1\.88and1\.751\.75\([Appendix˜T](https://arxiv.org/html/2605.23134#A20)\)\. We also ablate against R\-vine and show that it reaches lower NLL with per\-pair families and∼3,700\{\\sim\}3\{,\}700parameters\. However, it does not extend to per\-variable censoring \([Sections˜5\.3](https://arxiv.org/html/2605.23134#S5.SS3)and[U](https://arxiv.org/html/2605.23134#A21)\)\.

### 5\.6Retinopathy: family\-agnostic censored MLE

Table 7:Retinopathy held\-out NLL across all 10 generators \([Section˜R\.2](https://arxiv.org/html/2605.23134#A18.SS2)\);†\\daggermarks have no prior implementation\.The retinopathy paired\-eye study\(Huster et al\.,[1989](https://arxiv.org/html/2605.23134#bib.bib20)\)is a deliberately small benchmark withd=2d\{=\}2,n=197n\{=\}197, and60\.7%60\.7\\%censoring\. We use it to demonstrate the family\-agnostic API of our framework\. We use the same driver code to fit ten Archimedean generators, including five generators that no existing tool supports\. Nelsen12†achieves the best held\-out likelihood \([Table˜7](https://arxiv.org/html/2605.23134#S5.T7)\); on the two families CopulaCenR\(Sun and Ding,[2020](https://arxiv.org/html/2605.23134#bib.bib41)\)supports, our estimates agree to0\.6%0\.6\\%while CopulaCenR’s Frank optimiser fails on this dataset \([Appendix˜J](https://arxiv.org/html/2605.23134#A10)\)\.

## 6Discussions

We address concerns deferred from the main sections\.

Nesting validity\.[Equation˜3](https://arxiv.org/html/2605.23134#S2.E3)requires that each edge composition be a Bernstein function\. Every nested experiment in this paper uses a valid\-by\-construction parameterisation: the same\-family delta parameterisation ofMcNeil \([2008](https://arxiv.org/html/2605.23134#bib.bib30)\)for Clayton/Clayton, Frank/Frank, and Gumbel/Gumbel; Hofert\([2010](https://arxiv.org/html/2605.23134#bib.bib16)\)’s closed\-form shifted\-softplus constraint for AMH/Clayton; and the Bernstein\-composition construction of[Section˜G\.1](https://arxiv.org/html/2605.23134#A7.SS1)for the Gen\-AC and ACNet rows\. In each case, the nesting condition holds by construction without optimiser\-side enforcement\. The framework additionally returns a free sign\-alternation check on the edge\-composition Taylor coefficients per density evaluation, and[Appendix˜G](https://arxiv.org/html/2605.23134#A7)outlines a soft\-penalty method for cross\-family or custom generators that lack a closed\-form valid\-by\-construction parameterisation\.

Identifiability\.Identifiability is a structural limitation every copula\-based survival model inherits, including HACSurv\(Liu et al\.,[2025b](https://arxiv.org/html/2605.23134#bib.bib29)\)and DCSurvival\(Foomani et al\.,[2023](https://arxiv.org/html/2605.23134#bib.bib11)\)\. The practitioner response is to encode identifying assumptions through the choice of copula family\(Zheng and Klein,[1995](https://arxiv.org/html/2605.23134#bib.bib44); Foomani et al\.,[2023](https://arxiv.org/html/2605.23134#bib.bib11)\); the family becomes an explicit modelling choice rather than an implicit defect\. Our family\-agnostic inference makes that choice cheap: we enable the analyst to fit ten families with identical code and report sensitivity directly \([Section˜5\.6](https://arxiv.org/html/2605.23134#S5.SS6)\)\.

Censoring scope\.This work does not consider interval\-censored data\. In principle,acopulasupports this case by wrapping2m2^\{m\}calls to[Algorithm˜1](https://arxiv.org/html/2605.23134#alg1), wheremmis the number of interval\-censored variables, for the inclusion–exclusion over rectangle corners of the CDF\. The cost is exponential inmm, and we are not aware of any sub\-exponential algorithm to compute this CDF exactly\.

Stochastic density approximations\.Two prior approaches replace exact copula density evaluation with Monte Carlo: Simulated MLE \(SMLE\) samples frailties at the joint\-likelihood level\(Hofert et al\.,[2013](https://arxiv.org/html/2605.23134#bib.bib19)\); Gen\-AC samples latents at the generator\-derivative level\(Ng et al\.,[2021](https://arxiv.org/html/2605.23134#bib.bib34)\)\. Both yield consistent but biased estimators, with biasO​\(1/N\)O\(1/N\)from Jensen’s inequality onlog\\logcomposed with the nonlinear density\.[Appendix˜Q](https://arxiv.org/html/2605.23134#A17)empirically quantifies both:1\.SMLE atN=100N\{=\}100frailty samples per evaluation biases the inner Clayton parameter estimate by0\.360\.36;2\.Gen\-AC\-Fourier atN=32N\{=\}32frailty samples per evaluation biases the NLL by0\.090\.09–0\.120\.12nats per observation\.Our work enables exact density at high dimensions without finite\-sample approximation bias\.

Compile cost and memory\.JAX compilation takes 5–30 s for an S&P\-styled=98d\{=\}98nested likelihood; MIMIC\-IV hessianval\+gradatd=53d\{=\}53takes4040–5959s across families, broken down in[Table˜27](https://arxiv.org/html/2605.23134#A18.T27)\.acopulacaches compiled programs to disk, so reruns avoid compilation overhead; peak memory stays at1\.41\.4GB atd=98d\{=\}98and700700MB atd=53d\{=\}53\.

## 7Conclusion

acopulaunifies nested Archimedean copula inference under one differentiable pipeline\. Where prior tools offered only restricted slices—a few hand\-coded families, two\-level nesting, low\-dimensional mixed\-partial derivatives—acopulalifts all three to arbitrary families, arbitrary nesting depth, and high\-dimensional mixed\-partial derivatives, in one open\-source package, lowering the barrier to nested Archimedean copula inference for survival, finance, and reliability practitioners\. We hopeacopulacatalyzes expressive, interpretable censoring\-aware dependence models with principled likelihoods across survival, finance, and reliability practice\.

## References

- Aas et al\. \[2009\]Kjersti Aas, Claudia Czado, Arnoldo Frigessi, and Henrik Bakken\.Pair\-copula constructions of multiple dependence\.*Insurance: Mathematics and Economics*, 44\(2\), 2009\.
- Betancourt \[2018\]Michael Betancourt\.A geometric theory of higher\-order automatic differentiation\.*arXiv preprint arXiv:1812\.11592*, 2018\.
- Bettencourt et al\. \[2019\]Jesse Bettencourt, Matthew James Johnson, and David Duvenaud\.Taylor\-mode automatic differentiation for higher\-order derivatives in JAX\.In*Program Transformations for ML Workshop at NeurIPS*, 2019\.
- Bollerslev \[1986\]Tim Bollerslev\.Generalized autoregressive conditional heteroskedasticity\.*Journal of Econometrics*, 31\(3\), 1986\.
- Bradbury et al\. \[2018\]James Bradbury, Roy Frostig, Peter Hawkins, Matthew James Johnson, Chris Leary, Dougal Maclaurin, George Necula, Adam Paszke, Jake VanderPlas, Skye Wanderman\-Milne, and Qiao Zhang\.JAX: Composable transformations of Python\+NumPy programs, 2018\.URL[http://github\.com/jax\-ml/jax](http://github.com/jax-ml/jax)\.
- Brechmann \[2014\]Eike Christian Brechmann\.Hierarchical Kendall copulas: Properties and inference\.*Canadian Journal of Statistics*, 42\(1\), 2014\.
- Chen et al\. \[2025\]Xinyuan Chen, Yiwei Li, and Qian M\. Zhou\.Analysis of multivariate event times under informative censoring using vine copula\.*arXiv preprint arXiv:2502\.20608*, 2025\.
- Comtet \[1974\]Louis Comtet\.*Advanced Combinatorics: The Art of Finite and Infinite Expansions*\.D\. Reidel, 1974\.
- Embrechts et al\. \[2002\]Paul Embrechts, Alexander McNeil, and Daniel Straumann\.Correlation and dependence in risk management: properties and pitfalls\.In M\. A\. H\. Dempster, editor,*Risk Management: Value at Risk and Beyond*, pages 176–223\. Cambridge University Press, 2002\.
- Emura and Chen \[2018\]Takeshi Emura and Yi\-Hau Chen\.*Analysis of Survival Data with Dependent Censoring: Copula\-Based Approaches*\.Springer, 2018\.
- Foomani et al\. \[2023\]Ali Hossein Gharari Foomani, Michael Cooper, Russell Greiner, and Rahul G\. Krishnan\.Copula\-based deep survival models for dependent censoring\.In*Conference on Uncertainty in Artificial Intelligence*, 2023\.
- Genest et al\. \[1995\]Christian Genest, Kilani Ghoudi, and Louis\-Paul Rivest\.A semiparametric estimation procedure of dependence parameters in multivariate families of distributions\.*Biometrika*, 82\(3\), 1995\.
- Griewank and Walther \[2008\]Andreas Griewank and Andrea Walther\.*Evaluating Derivatives: Principles and Techniques of Algorithmic Differentiation*\.SIAM, 2008\.
- Higham \[2002\]Nicholas J\. Higham\.*Accuracy and Stability of Numerical Algorithms*\.SIAM, 2002\.
- Hofert \[2008\]Marius Hofert\.Sampling Archimedean copulas\.*Computational Statistics & Data Analysis*, 52\(12\), 2008\.
- Hofert \[2010\]Marius Hofert\.*Sampling Nested Archimedean Copulas: With Applications to CDO Pricing*\.PhD thesis, Ulm University, 2010\.
- Hofert and Pham \[2013\]Marius Hofert and David Pham\.Densities of nested Archimedean copulas\.*Journal of Multivariate Analysis*, 118, 2013\.
- Hofert et al\. \[2012\]Marius Hofert, Martin Mächler, and Alexander J\. McNeil\.Likelihood inference for Archimedean copulas in high dimensions under known margins\.*Journal of Multivariate Analysis*, 110, 2012\.
- Hofert et al\. \[2013\]Marius Hofert, Martin Mächler, and Alexander J\. McNeil\.Archimedean copulas in high dimensions: Estimators and numerical challenges motivated by financial applications\.*Journal de la Société Française de Statistique*, 154\(1\), 2013\.
- Huster et al\. \[1989\]William J\. Huster, Ron Brookmeyer, and Steven G\. Self\.Modelling paired survival data with covariates\.*Biometrics*, 45\(1\), 1989\.
- Joe and Xu \[1996\]Harry Joe and James J\. Xu\.The estimation method of inference functions for margins for multivariate models\.Technical Report 166, Department of Statistics, University of British Columbia, 1996\.
- Johnson et al\. \[2023\]Alistair E\. W\. Johnson, Lucas Bulgarelli, Lu Shen, Alvin Gayles, Ayad Shammout, Steven Horng, Tom J\. Pollard, Sicheng Hao, Benjamin Moody, Brian Gow, Li wei H\. Lehman, Leo A\. Celi, and Roger G\. Mark\.MIMIC\-IV, a freely accessible electronic health record dataset\.*Scientific Data*, 10\(1\), 2023\.
- Laverny and Jimenez \[2024\]Oskar Laverny and Santiago Jimenez\.Copulas\.jl: A fully distributions\.jl\-compliant copula package\.*Journal of Open Source Software*, 9\(94\):6189, 2024\.
- Li et al\. \[2020\]Dongdong Li, X\. Joan Hu, Mary L\. McBride, and John J\. Spinelli\.Multiple event times in the presence of informative censoring: modeling and analysis by copulas\.*Lifetime Data Analysis*, 26\(3\), 2020\.
- Li et al\. \[2023\]Dongdong Li, X\. Joan Hu, and Rui Wang\.Evaluating association between two event times with observations subject to informative censoring\.*Journal of the American Statistical Association*, 118\(542\), 2023\.
- Lin \[2024\]Min Lin\.Automatic functional differentiation in JAX\.In*International Conference on Learning Representations*, 2024\.
- Ling et al\. \[2020\]Chun Kai Ling, Fei Fang, and J\. Zico Kolter\.Deep Archimedean copulas\.In*Advances in Neural Information Processing Systems*, 2020\.
- Liu et al\. \[2025a\]Anji Liu, Oliver Broadrick, Mathias Niepert, and Guy Van den Broeck\.Discrete copula diffusion\.In*International Conference on Learning Representations*, 2025a\.
- Liu et al\. \[2025b\]Xin Liu, Weijia Zhang, and Min\-Ling Zhang\.HACSurv: A hierarchical copula\-based approach for survival analysis with dependent competing risks\.In*International Conference on Artificial Intelligence and Statistics*, 2025b\.
- McNeil \[2008\]Alexander J\. McNeil\.Sampling nested Archimedean copulas\.*Journal of Statistical Computation and Simulation*, 78\(6\), 2008\.
- McNeil and Nešlehová \[2009\]Alexander J\. McNeil and Johanna Nešlehová\.Multivariate Archimedean copulas,dd\-monotone functions andℓ1\\ell\_\{1\}\-norm symmetric distributions\.*The Annals of Statistics*, 37\(5B\), 2009\.
- Nagler and Vatter \[2023\]Thomas Nagler and Thibault Vatter\.*rvinecopulib: High Performance Algorithms for Vine Copula Modeling*, 2023\.URL[https://CRAN\.R\-project\.org/package=rvinecopulib](https://cran.r-project.org/package=rvinecopulib)\.R package, CRAN\.
- Nelsen \[2006\]Roger B\. Nelsen\.*An Introduction to Copulas*\.Springer, 2006\.
- Ng et al\. \[2021\]Yuting Ng, Ali Hasan, Khalil Elkhalil, and Vahid Tarokh\.Generative Archimedean copulas\.In*Uncertainty in Artificial Intelligence*, 2021\.
- Okhrin and Ristig \[2014\]Ostap Okhrin and Alexander Ristig\.Hierarchical Archimedean copulae: The HAC package\.*Journal of Statistical Software*, 58\(4\), 2014\.
- Phan et al\. \[2019\]Du Phan, Neeraj Pradhan, and Martin Jankowiak\.Composable effects for flexible and accelerated probabilistic programming in NumPyro\.*arXiv preprint arXiv:1912\.11554*, 2019\.
- Radul et al\. \[2023\]Alexey Radul, Adam Paszke, Roy Frostig, Matthew Johnson, and Dougal Maclaurin\.You only linearize once: Tangents transpose to gradients\.*Proceedings of the ACM on Programming Languages*, 7\(POPL\), 2023\.
- Rezapour \[2015\]Mohsen Rezapour\.On the construction of nested Archimedean copulas fordd\-monotone generators\.*Statistics & Probability Letters*, 101, 2015\.
- Sklar \[1959\]Abe Sklar\.Fonctions de répartition ànndimensions et leurs marges\.*Publications de l’Institut de Statistique de l’Université de Paris*, 8, 1959\.
- Stöber et al\. \[2013\]Jakob Stöber, Harry Joe, and Claudia Czado\.Simplified pair copula constructions—limitations and extensions\.*Journal of Multivariate Analysis*, 119, 2013\.
- Sun and Ding \[2020\]Tao Sun and Ying Ding\.CopulaCenR: Copula based regression models for bivariate censored data in R\.*The R Journal*, 12\(1\), 2020\.
- Tancik et al\. \[2020\]Matthew Tancik, Pratul P\. Srinivasan, Ben Mildenhall, Sara Fridovich\-Keil, Nithin Raghavan, Utkarsh Singhal, Ravi Ramamoorthi, Jonathan T\. Barron, and Ren Ng\.Fourier features let networks learn high frequency functions in low dimensional domains\.In*Advances in Neural Information Processing Systems*, 2020\.
- Weibull \[1951\]Waloddi Weibull\.A statistical distribution function of wide applicability\.*Journal of Applied Mechanics*, 18\(3\), 1951\.
- Zheng and Klein \[1995\]Ming Zheng and John P\. Klein\.Estimates of marginal survival for dependent competing risks based on an assumed copula\.*Biometrika*, 82\(1\), 1995\.

###### Appendix Contents

1. [1Introduction](https://arxiv.org/html/2605.23134#S1)
2. [2Background and Notation](https://arxiv.org/html/2605.23134#S2)
3. [3Method](https://arxiv.org/html/2605.23134#S3)1. [3\.1Worked example:d=3d\{=\}3nested copula with one censored leaf](https://arxiv.org/html/2605.23134#S3.SS1) 2. [3\.2Algorithm](https://arxiv.org/html/2605.23134#S3.SS2)
4. [4Implementation](https://arxiv.org/html/2605.23134#S4)
5. [5Experiments](https://arxiv.org/html/2605.23134#S5)1. [5\.1Correctness validation](https://arxiv.org/html/2605.23134#S5.SS1) 2. [5\.2Scaling: head\-to\-head against RnacLL](https://arxiv.org/html/2605.23134#S5.SS2) 3. [5\.3MIMIC\-IV: classical nested copulas atd=53d\{=\}53](https://arxiv.org/html/2605.23134#S5.SS3) 4. [5\.4Capability demonstrations on nested neural Archimedean copulas](https://arxiv.org/html/2605.23134#S5.SS4) 5. [5\.5S&P 500: nested likelihood atd=98d\{=\}98](https://arxiv.org/html/2605.23134#S5.SS5) 6. [5\.6Retinopathy: family\-agnostic censored MLE](https://arxiv.org/html/2605.23134#S5.SS6)
6. [6Discussions](https://arxiv.org/html/2605.23134#S6)
7. [7Conclusion](https://arxiv.org/html/2605.23134#S7)
8. [References](https://arxiv.org/html/2605.23134#bib)
9. [ADensity via Generator Derivatives: Full Derivation](https://arxiv.org/html/2605.23134#A1)1. [A\.1Partial Bell Polynomials](https://arxiv.org/html/2605.23134#A1.SS1) 2. [A\.2The Bottom\-Up Recursion](https://arxiv.org/html/2605.23134#A1.SS2) 3. [A\.3Computing Bell Polynomials via Polynomial Powering](https://arxiv.org/html/2605.23134#A1.SS3) 4. [A\.4Edge\-Taylor coefficients via implicit differentiation](https://arxiv.org/html/2605.23134#A1.SS4) 5. [A\.5Worked Example:2×22\\times 2Nested Copula](https://arxiv.org/html/2605.23134#A1.SS5) 6. [A\.6Complexity Analysis](https://arxiv.org/html/2605.23134#A1.SS6)
10. [BNumerical Stability of Taylor\-Mode Derivatives](https://arxiv.org/html/2605.23134#A2)
11. [CLog\-Space Taylor\-Mode Differentiation](https://arxiv.org/html/2605.23134#A3)1. [C\.1Why naive float64 fails for Frank\-type generators](https://arxiv.org/html/2605.23134#A3.SS1) 2. [C\.2Architecture:\(sign,log⁡\|c\|\)\(\\text\{sign\},\\log\|c\|\)throughout](https://arxiv.org/html/2605.23134#A3.SS2) 3. [C\.3Crossing the bell boundary without leaving log\-space](https://arxiv.org/html/2605.23134#A3.SS3) 4. [C\.4Validation](https://arxiv.org/html/2605.23134#A3.SS4) 5. [C\.5Cost and configurability](https://arxiv.org/html/2605.23134#A3.SS5)
12. [DPolynomial Multiplication Alternatives](https://arxiv.org/html/2605.23134#A4)
13. [ESymbolic Inversion and Numeric Fallback](https://arxiv.org/html/2605.23134#A5)
14. [FRosenblatt Sampling via Generator Derivatives](https://arxiv.org/html/2605.23134#A6)
15. [GValidity of Nested Constructions](https://arxiv.org/html/2605.23134#A7)1. [G\.1Validity of the nested Gen\-AC/ACNet construction](https://arxiv.org/html/2605.23134#A7.SS1)
16. [HGen\-AC MLP\-pushforward generator with Fourier features](https://arxiv.org/html/2605.23134#A8)
17. [IDeep ACNet Generator](https://arxiv.org/html/2605.23134#A9)
18. [JDetailed R Comparison](https://arxiv.org/html/2605.23134#A10)
19. [KAD Gradient Verification](https://arxiv.org/html/2605.23134#A11)
20. [LHigh\-Dimensional Parameter Recovery](https://arxiv.org/html/2605.23134#A12)1. [L\.1Frank and Gumbel Recovery](https://arxiv.org/html/2605.23134#A12.SS1)
21. [MCensoring Validation](https://arxiv.org/html/2605.23134#A13)
22. [NSupported Archimedean Families](https://arxiv.org/html/2605.23134#A14)
23. [OHigh\-Dimensional Nested Scaling](https://arxiv.org/html/2605.23134#A15)
24. [PMemory and Compilation Benchmarks](https://arxiv.org/html/2605.23134#A16)
25. [QExact MLE vs\. Simulated Maximum Likelihood](https://arxiv.org/html/2605.23134#A17)
26. [RPredictive Validation](https://arxiv.org/html/2605.23134#A18)1. [R\.1MIMIC\-IV](https://arxiv.org/html/2605.23134#A18.SS1)1. [R\.1\.1Ablation: independence and single\-layer baselines](https://arxiv.org/html/2605.23134#A18.SS1.SSS1) 2. [R\.1\.2Marginal misspecification](https://arxiv.org/html/2605.23134#A18.SS1.SSS2) 2. [R\.2Retinopathy](https://arxiv.org/html/2605.23134#A18.SS2)
27. [SHead\-to\-Head with HACSurv on Framingham](https://arxiv.org/html/2605.23134#A19)
28. [TS&P 500 Sector\-Level Tail Dependence](https://arxiv.org/html/2605.23134#A20)
29. [UVine Copula Comparison](https://arxiv.org/html/2605.23134#A21)
30. [VHierarchical Kendall Copulas](https://arxiv.org/html/2605.23134#A22)

## Appendix ADensity via Generator Derivatives: Full Derivation

This appendix establishes that the Bell polynomial recursion of[Section˜3\.2](https://arxiv.org/html/2605.23134#S3.SS2)computes the exact nested\-Archimedean density and admits a polynomial\-powering implementation whose cost matches the bounds we cite in the body\. We give the partial Bell polynomial background, the bottom\-up recursion at each internal node, the polynomial\-powering reformulation of the Bell\-transform step, a worked2×22\\times 2example, and a complexity analysis with empirical scaling\.

### A\.1Partial Bell Polynomials

For integersn≥k≥1n\\geq k\\geq 1, the partial Bell polynomialBn,k​\(x1,x2,…,xn−k\+1\)B\_\{n,k\}\(x\_\{1\},x\_\{2\},\\ldots,x\_\{n\-k\+1\}\)counts all partitions of\{1,…,n\}\\\{1,\\ldots,n\\\}into exactlykknon\-empty blocks, weighted by thexmx\_\{m\}values\. These polynomials arise in Faà di Bruno’s formula for the higher\-order chain rule:

dnd​tn​\[f​\(g​\(t\)\)\]=∑k=1nf\(k\)​\(g​\(t\)\)​Bn,k​\(g′​\(t\),g′′​\(t\),…,g\(n−k\+1\)​\(t\)\)\.\\frac\{d^\{n\}\}\{dt^\{n\}\}\\bigl\[f\(g\(t\)\)\\bigr\]=\\sum\_\{k=1\}^\{n\}f^\{\(k\)\}\(g\(t\)\)\\;B\_\{n,k\}\\bigl\(g^\{\\prime\}\(t\),\\,g^\{\\prime\\prime\}\(t\),\\,\\ldots,\\,g^\{\(n\-k\+1\)\}\(t\)\\bigr\)\.\(10\)They satisfy the exponential generating function identity\[Comtet,[1974](https://arxiv.org/html/2605.23134#bib.bib8), §3\.3\]

∑n≥kBn,k​\(x1,x2,…\)​tnn\!=1k\!​\(∑m≥1xmm\!​tm\)k,\\sum\_\{n\\geq k\}B\_\{n,k\}\(x\_\{1\},x\_\{2\},\\ldots\)\\,\\frac\{t^\{n\}\}\{n\!\}=\\frac\{1\}\{k\!\}\\left\(\\sum\_\{m\\geq 1\}\\frac\{x\_\{m\}\}\{m\!\}\\,t^\{m\}\\right\)^\{\\\!k\},\(11\)which provides the bridge to the polynomial\-powering algorithm\. For quick reference, the following table lists the values for dimensions up tod=4d=4:

### A\.2The Bottom\-Up Recursion

At each internal nodevvwith childrenc1,…,cmc\_\{1\},\\ldots,c\_\{m\}, the subtree coefficientsβk\(v\)\\beta\_\{k\}^\{\(v\)\}follow from three stages\.

##### Leaf base case\.

An uncensored leafℓ\\ellcontributesα1\(v,ℓ\)=1\\alpha\_\{1\}^\{\(v,\\ell\)\}=1\. A censored leaf contributesα0\(v,ℓ\)=1\\alpha\_\{0\}^\{\(v,\\ell\)\}=1, which acts as the identity in the subsequent Cauchy product\.

##### Bell polynomial step\.

For internal childccwith already\-computed subtree coefficientsβj\(c\)\\beta\_\{j\}^\{\(c\)\}, define the edge compositionhv​c​\(t\)=ψv−1​\(ψc​\(t\)\)h\_\{vc\}\(t\)=\\psi\_\{v\}^\{\-1\}\(\\psi\_\{c\}\(t\)\)\. The per\-child coefficients at the parent level are:

αk\(v,c\)=∑j=kdcβj\(c\)​Bj,k​\(hv​c′​\(tc​\(𝐮\)\),hv​c′′​\(tc​\(𝐮\)\),…\),k=1,…,dc\.\\alpha\_\{k\}^\{\(v,c\)\}=\\sum\_\{j=k\}^\{d\_\{c\}\}\\beta\_\{j\}^\{\(c\)\}\\;B\_\{j,k\}\\\!\\left\(h^\{\\prime\}\_\{vc\}\(t\_\{c\}\(\\mathbf\{u\}\)\),\\;h^\{\\prime\\prime\}\_\{vc\}\(t\_\{c\}\(\\mathbf\{u\}\)\),\\;\\ldots\\right\),\\qquad k=1,\\ldots,d\_\{c\}\.\(12\)

##### Cauchy product\.

The per\-childα\\alpha\-polynomials combine via polynomial multiplication:

βk\(v\)=∑k1\+⋯\+km=kki≥0∏i=1mαki\(v,ci\),k=0,…,dv\.\\beta\_\{k\}^\{\(v\)\}=\\sum\_\{\\begin\{subarray\}\{c\}k\_\{1\}\+\\cdots\+k\_\{m\}=k\\\\ k\_\{i\}\\geq 0\\end\{subarray\}\}\\prod\_\{i=1\}^\{m\}\\alpha\_\{k\_\{i\}\}^\{\(v,c\_\{i\}\)\},\\qquad k=0,\\ldots,d\_\{v\}\.\(13\)

##### Root assembly\.

At the rootrr, evaluate the root generator’s Taylor expansion attrt\_\{r\}up to orderddvia a single jet call, then form∑k=1dβk\(r\)​ψr\(k\)​\(tr\)\\sum\_\{k=1\}^\{d\}\\beta\_\{k\}^\{\(r\)\}\\psi\_\{r\}^\{\(k\)\}\(t\_\{r\}\)\.

### A\.3Computing Bell Polynomials via Polynomial Powering

We derive the polynomial\-powering algorithm that replaces explicit Bell polynomial evaluation with iterated polynomial convolution\.

##### Step 1\.

Letxm=hv​c\(m\)​\(tc​\(𝐮\)\)x\_\{m\}=h\_\{vc\}^\{\(m\)\}\(t\_\{c\}\(\\mathbf\{u\}\)\)and definep​\(t\)=∑m≥1\(xm/m\!\)​tmp\(t\)=\\sum\_\{m\\geq 1\}\(x\_\{m\}/m\!\)\\,t^\{m\}\.

##### Step 2\.

The coefficient oftnt^\{n\}inp​\(t\)kp\(t\)^\{k\}collects all ways to choosekkterms from copies ofppwith exponents summing tonn:

\[tn\]​p​\(t\)k=∑m1\+⋯\+mk=nmi≥1xm1m1\!​⋯​xmkmk\!\.\[t^\{n\}\]\\,p\(t\)^\{k\}=\\sum\_\{\\begin\{subarray\}\{c\}m\_\{1\}\+\\cdots\+m\_\{k\}=n\\\\ m\_\{i\}\\geq 1\\end\{subarray\}\}\\frac\{x\_\{m\_\{1\}\}\}\{m\_\{1\}\!\}\\cdots\\frac\{x\_\{m\_\{k\}\}\}\{m\_\{k\}\!\}\.\(14\)

##### Step 3\.

The generating function identity[Equation˜11](https://arxiv.org/html/2605.23134#A1.E11)givesBn,k=\(n\!/k\!\)​\[tn\]​p​\(t\)kB\_\{n,k\}=\(n\!/k\!\)\\,\[t^\{n\}\]\\,p\(t\)^\{k\}\.

##### Step 4\.

Substituting into[Equation˜12](https://arxiv.org/html/2605.23134#A1.E12):

αk=∑jβj⋅j\!k\!⋅\[tj\]​p​\(t\)k=1k\!​∑jj\!​βj⏟β~j⋅\[tj\]​p​\(t\)k⏟qk​\[j\]=1k\!​⟨β~,qk⟩\.\\alpha\_\{k\}=\\sum\_\{j\}\\beta\_\{j\}\\cdot\\frac\{j\!\}\{k\!\}\\cdot\[t^\{j\}\]\\,p\(t\)^\{k\}=\\frac\{1\}\{k\!\}\\sum\_\{j\}\\underbrace\{j\!\\,\\beta\_\{j\}\}\_\{\\tilde\{\\beta\}\_\{j\}\}\\cdot\\underbrace\{\[t^\{j\}\]\\,p\(t\)^\{k\}\}\_\{q\_\{k\}\[j\]\}=\\frac\{1\}\{k\!\}\\,\\langle\\tilde\{\\beta\},\\;q\_\{k\}\\rangle\.\(15\)

##### Step 5\.

Taylor\-mode AD returnsh\(m\)​\(tc​\(𝐮\)\)/m\!h^\{\(m\)\}\(t\_\{c\}\(\\mathbf\{u\}\)\)/m\!, which equals the coefficient oftmt^\{m\}inp​\(t\)p\(t\)\. The jet output*is*the coefficient vector ofpp\.

##### Step 6\.

Buildqk=pkq\_\{k\}=p^\{k\}iteratively via truncated polynomial convolution, implemented aslax\.scan\+lax\.conv\_general\_dilated:

q←p;for​k=1,…,dc:αk=⟨β~,q⟩k\!,q←conv​\(q,p\)\|≤dc\.q\\leftarrow p;\\qquad\\text\{for \}k=1,\\ldots,d\_\{c\}:\\;\\;\\alpha\_\{k\}=\\frac\{\\langle\\tilde\{\\beta\},\\,q\\rangle\}\{k\!\},\\;\\;q\\leftarrow\\mathrm\{conv\}\(q,p\)\\big\|\_\{\\leq d\_\{c\}\}\.\(16\)Each convolution of two length\-dcd\_\{c\}polynomials costsO​\(dc2\)O\(d\_\{c\}^\{2\}\), iterateddcd\_\{c\}times, givingO​\(dc3\)O\(d\_\{c\}^\{3\}\)per child\.

### A\.4Edge\-Taylor coefficients via implicit differentiation

The polynomial\-powering recipe of[Section˜A\.3](https://arxiv.org/html/2605.23134#A1.SS3)consumes the edge Taylor coefficientshv​c\(m\)​\(tc\)/m\!h\_\{vc\}^\{\(m\)\}\(t\_\{c\}\)/m\!\. The naive realisation appliesJetdirectly to the explicit compositionhv​c=ψv−1∘ψch\_\{vc\}=\\psi\_\{v\}^\{\-1\}\\circ\\psi\_\{c\}, propagating two univariate jets through the chain rule \(Faà di Bruno\) atO​\(dc2\)O\(d\_\{c\}^\{2\}\)\. The framework does not take this path by default: it is unstable for many generator pairs — Frank atd≥90d\\geq 90, Nelsen9 at moderateθ\\theta— becauseψc​\(t\)\\psi\_\{c\}\(t\)underflows to zero in float64 beforeψv−1\\psi\_\{v\}^\{\-1\}can recover it\.

The framework instead recoverspv​cp\_\{vc\}without ever evaluatingψv−1∘ψc\\psi\_\{v\}^\{\-1\}\\circ\\psi\_\{c\}as a single chain: the implicit solve reads only the separate jets ofψv\\psi\_\{v\}ath0h\_\{0\}and ofψc\\psi\_\{c\}attct\_\{c\}, both of which lie in the well\-conditioned interior of their respective generator domains\. This construction therefore avoids the composed\-chain underflow that defeats the naive path\. The same mechanism extends naturally to opaque generators \(e\.g\. a neural Archimedean family\[Ling et al\.,[2020](https://arxiv.org/html/2605.23134#bib.bib27), Liu et al\.,[2025b](https://arxiv.org/html/2605.23134#bib.bib29)\]\) without an analytic inverse: the algorithm usesψv−1\\psi\_\{v\}^\{\-1\}only at one point to computeh0h\_\{0\}, with all higher\-order Taylor information recovered algebraically fromψv\\psi\_\{v\}’s forward jet\. We start from the implicit equation

ψv​\(hv​c​\(t\)\)=ψc​\(t\)\.\\psi\_\{v\}\\bigl\(h\_\{vc\}\(t\)\\bigr\)\\;=\\;\\psi\_\{c\}\(t\)\.\(17\)Leth0:=ψv−1​\(ψc​\(tc\)\)h\_\{0\}:=\\psi\_\{v\}^\{\-1\}\(\\psi\_\{c\}\(t\_\{c\}\)\)and writehv​c​\(tc\+ε\)=h0\+Q​\(ε\)h\_\{vc\}\(t\_\{c\}\+\\varepsilon\)=h\_\{0\}\+Q\(\\varepsilon\)withQ​\(ε\)=∑m=1dcqm​εmQ\(\\varepsilon\)=\\sum\_\{m=1\}^\{d\_\{c\}\}q\_\{m\}\\,\\varepsilon^\{m\}\. Denote the Taylor coefficients of the two generators byam:=ψv\(m\)​\(h0\)/m\!a\_\{m\}:=\\psi\_\{v\}^\{\(m\)\}\(h\_\{0\}\)/m\!andbm:=ψc\(m\)​\(tc\)/m\!b\_\{m\}:=\\psi\_\{c\}^\{\(m\)\}\(t\_\{c\}\)/m\!\.[Equation˜17](https://arxiv.org/html/2605.23134#A1.E17)then expands as

∑m≥1am​Q​\(ε\)m=∑m≥1bm​εm,\\sum\_\{m\\geq 1\}a\_\{m\}\\,Q\(\\varepsilon\)^\{m\}\\;=\\;\\sum\_\{m\\geq 1\}b\_\{m\}\\,\\varepsilon^\{m\},\(18\)where the constant terms cancel becauseψv​\(h0\)=ψc​\(tc\)\\psi\_\{v\}\(h\_\{0\}\)=\\psi\_\{c\}\(t\_\{c\}\)by construction\. Matching the coefficient ofεk\\varepsilon^\{k\}yields a triangular system inq1,…,qdcq\_\{1\},\\ldots,q\_\{d\_\{c\}\}:

a1​qk\+∑m≥2am​\[εk\]​Q​\(ε\)m⏟depends only on​q1,…,qk−1=bk\.a\_\{1\}\\,q\_\{k\}\\;\+\\;\\underbrace\{\\sum\_\{m\\geq 2\}a\_\{m\}\\,\[\\varepsilon^\{k\}\]\\,Q\(\\varepsilon\)^\{m\}\}\_\{\\text\{depends only on \}q\_\{1\},\\ldots,q\_\{k\-1\}\}\\;=\\;b\_\{k\}\.\(19\)Each\[εk\]​Qm\[\\varepsilon^\{k\}\]\\,Q^\{m\}involvesqjq\_\{j\}only forj<kj<k, so we solve forqkq\_\{k\}in increasingkkonce we can read the correction term\.

##### Algorithm\.

The framework maintains the tableC​\[m,j\]:=\[εj\]​Q​\(ε\)m\+2C\[m,j\]:=\[\\varepsilon^\{j\}\]\\,Q\(\\varepsilon\)^\{m\+2\}asQQ’s coefficients fill in\. At outer stepkk, the algorithm has finalisedQQat slots1,…,k1,\\ldots,k, so columnsj≤kj\\leq kofCCare final and only columnj=k\+1j=k\+1needs updating\. The single\-column recurrence is

C​\[m,k\+1\]=∑i=1kqi⋅prevm​\[k\+1−i\],prevm=\{Qm=0,C​\[m−1,⋅\]m≥1\.C\[m,k\+1\]\\;=\\;\\sum\_\{i=1\}^\{k\}q\_\{i\}\\cdot\\mathrm\{prev\}\_\{m\}\[k\+1\-i\],\\qquad\\mathrm\{prev\}\_\{m\}\\;=\\;\\begin\{cases\}Q&m=0,\\\\ C\[m\-1,\\,\\cdot\]&m\\geq 1\.\\end\{cases\}\(20\)The readsC​\[m−1,k\+1−i\]C\[m\-1,k\+1\-i\]at stepkktouch columns≤k\\leq k— already final from prior steps — so all rowsmmin thej=k\+1j=k\+1column update in parallel\.[Algorithm˜2](https://arxiv.org/html/2605.23134#alg2)writes this out\.

Algorithm 2Implicit Taylor solve forpv​cp\_\{vc\}when the user has registered no closed\-form edge composition\. All steps are JIT\-compiled and differentiable in the generator parameters\.- Require:a=\(a1,…,adc\)a=\(a\_\{1\},\\ldots,a\_\{d\_\{c\}\}\)andb=\(b1,…,bdc\)b=\(b\_\{1\},\\ldots,b\_\{d\_\{c\}\}\), the rescaled Taylor coefficients ofψv\\psi\_\{v\}ath0h\_\{0\}and ofψc\\psi\_\{c\}attct\_\{c\}
- Ensure:q=\(q1,…,qdc\)q=\(q\_\{1\},\\ldots,q\_\{d\_\{c\}\}\)satisfying[Equation˜18](https://arxiv.org/html/2605.23134#A1.E18)
- q1←b1/a1q\_\{1\}\\leftarrow b\_\{1\}/a\_\{1\}
- C​\[m,j\]←0C\[m,j\]\\leftarrow 0for allm,jm,j;C​\[m,m\+2\]←q1m\+2C\[m,m\+2\]\\leftarrow q\_\{1\}^\{m\+2\}form=0,…,dc−2m=0,\\ldots,d\_\{c\}\-2 ⊳\\trianglerightInitialiseCCfromQ\(0\)=q1​εQ^\{\(0\)\}=q\_\{1\}\\,\\varepsilon
- fork=2,…,dck=2,\\ldots,d\_\{c\}do
- correction←∑m=0dc−2am\+2⋅C​\[m,k\]\\text\{correction\}\\leftarrow\\sum\_\{m=0\}^\{d\_\{c\}\-2\}a\_\{m\+2\}\\cdot C\[m,k\]
- qk←\(bk−correction\)/a1q\_\{k\}\\leftarrow\(b\_\{k\}\-\\text\{correction\}\)/a\_\{1\}
- Appendqkq\_\{k\}toQQat slotkk
- form=0,1,…,dc−2m=0,1,\\ldots,d\_\{c\}\-2do
- C​\[m,k\+1\]←∑i=1kqi⋅prevm​\[k\+1−i\]C\[m,k\+1\]\\leftarrow\\sum\_\{i=1\}^\{k\}q\_\{i\}\\cdot\\mathrm\{prev\}\_\{m\}\[k\+1\-i\], withprevm\\mathrm\{prev\}\_\{m\}as in[Equation˜20](https://arxiv.org/html/2605.23134#A1.E20)\.

##### Complexity\.

The correction read at stepkkis anO​\(dc\)O\(d\_\{c\}\)inner product\. The column update for slotk\+1k\+1doesdc−1d\_\{c\}\-1length\-kkdot products, costingO​\(dc⋅k\)O\(d\_\{c\}\\cdot k\)\. Summing overk=1,…,dc−1k=1,\\ldots,d\_\{c\}\-1givesO​\(dc3\)O\(d\_\{c\}^\{3\}\)per implicit\-solve call, matching the per\-child polynomial\-powering cost in[Section˜A\.3](https://arxiv.org/html/2605.23134#A1.SS3)\. Naively rebuilding the fullQ2,…,QdcQ^\{2\},\\ldots,Q^\{d\_\{c\}\}ladder at every outer step would inflate this toO​\(dc4\)O\(d\_\{c\}^\{4\}\);[Algorithm˜2](https://arxiv.org/html/2605.23134#alg2)keeps the ladder amortised acrosskkby writing one column per step\. At balancedd\\sqrt\{d\}topology the per\-edge cost is dominated by theO​\(dc3\)O\(d\_\{c\}^\{3\}\)Bell powering of[Section˜A\.3](https://arxiv.org/html/2605.23134#A1.SS3), so the headlineT​\(𝒯\)=O​\(d2\)T\(\\mathcal\{T\}\)=O\(d^\{2\}\)bound of[Equation˜24](https://arxiv.org/html/2605.23134#A1.E24)is unaffected by the implicit solve\.

##### Closed\-form override\.

Same\-family pairs \(Clayton, Gumbel, Joe, Frank, Nelsen9 withθv<θc\\theta\_\{v\}<\\theta\_\{c\}\) admit an algebraic simplification ofhv​c​\(t\)h\_\{vc\}\(t\)— for instance, Clayton/Clayton hash​\(t\)=\(1\+t\)θv/θc−1h\(t\)=\(1\+t\)^\{\\theta\_\{v\}/\\theta\_\{c\}\}\-1\. Users register such pairs viaregister\_compositionto bypass the implicit solve and apply a singleO​\(dc2\)O\(d\_\{c\}^\{2\}\)jet through the registered form, saving a factor ofdcd\_\{c\}per edge\. The framework ships built\-in registrations for the same\-family pairs above, so the classical single\-family experiments \(S&P 500, MIMIC\-IV with classical generators, retinopathy\) take this fast path\. The neural Archimedean rows in[Table˜4](https://arxiv.org/html/2605.23134#S5.T4)\(ACNet, Gen\-AC\) also bypass the implicit solve: their inner generator is constructed asψinner​\(t\)=ψouter​\(b​\(t\)\)\\psi\_\{\\text\{inner\}\}\(t\)=\\psi\_\{\\text\{outer\}\}\(b\(t\)\)withbba Bernstein function \([Section˜G\.1](https://arxiv.org/html/2605.23134#A7.SS1)\), so the edge composition ishv​c​\(t\)=b​\(t\)h\_\{vc\}\(t\)=b\(t\)in closed form\. Cross\-classical\-family pairs without a registered closed form \(e\.g\. AMH/Clayton\) fall through to the default implicit solve\.

### A\.5Worked Example:2×22\\times 2Nested Copula

We trace the polynomial\-powering recursion through the smallest non\-trivial nested copula—two sectors of two leaves—to expose how the per\-childα\\alpha\-vectors and the Cauchy product at the root combine into the closed\-form density\. Consider rootψ0\\psi\_\{0\}with two sectors, each governed byψ1\\psi\_\{1\}with two leaves\. At each sector,tv=ψ1−1​\(uv​1\)\+ψ1−1​\(uv​2\)t\_\{v\}=\\psi\_\{1\}^\{\-1\}\(u\_\{v1\}\)\+\\psi\_\{1\}^\{\-1\}\(u\_\{v2\}\),Cv=ψ1​\(tv\)C\_\{v\}=\\psi\_\{1\}\(t\_\{v\}\), andtr=ψ0−1​\(C1\)\+ψ0−1​\(C2\)t\_\{r\}=\\psi\_\{0\}^\{\-1\}\(C\_\{1\}\)\+\\psi\_\{0\}^\{\-1\}\(C\_\{2\}\)\.

Each sector has two leaf children, soβ2\(v\)=1\\beta\_\{2\}^\{\(v\)\}=1\. Defineh​\(t\)=ψ0−1​\(ψ1​\(t\)\)h\(t\)=\\psi\_\{0\}^\{\-1\}\(\\psi\_\{1\}\(t\)\)and run jet attvt\_\{v\}to getpvp\_\{v\}\. Polynomial powering givesα1\(r,v\)=h′′​\(tv\)\\alpha\_\{1\}^\{\(r,v\)\}=h^\{\\prime\\prime\}\(t\_\{v\}\)andα2\(r,v\)=\(h′​\(tv\)\)2\\alpha\_\{2\}^\{\(r,v\)\}=\(h^\{\\prime\}\(t\_\{v\}\)\)^\{2\}\. The Cauchy product at the root yields:

β2\(r\)\\displaystyle\\beta\_\{2\}^\{\(r\)\}=α1\(r,1\)​α1\(r,2\),\\displaystyle=\\alpha\_\{1\}^\{\(r,1\)\}\\,\\alpha\_\{1\}^\{\(r,2\)\},\(21\)β3\(r\)\\displaystyle\\beta\_\{3\}^\{\(r\)\}=α1\(r,1\)​α2\(r,2\)\+α2\(r,1\)​α1\(r,2\),\\displaystyle=\\alpha\_\{1\}^\{\(r,1\)\}\\,\\alpha\_\{2\}^\{\(r,2\)\}\+\\alpha\_\{2\}^\{\(r,1\)\}\\,\\alpha\_\{1\}^\{\(r,2\)\},\(22\)β4\(r\)\\displaystyle\\beta\_\{4\}^\{\(r\)\}=α2\(r,1\)​α2\(r,2\)\.\\displaystyle=\\alpha\_\{2\}^\{\(r,1\)\}\\,\\alpha\_\{2\}^\{\(r,2\)\}\.\(23\)The density isc​\(𝐮\)=\[β2​ψ0′′​\(tr\)\+β3​ψ0′′′​\(tr\)\+β4​ψ0\(4\)​\(tr\)\]⋅∏i=14\(ψ1−1\)′​\(ui\)c\(\\mathbf\{u\}\)=\[\\beta\_\{2\}\\psi\_\{0\}^\{\\prime\\prime\}\(t\_\{r\}\)\+\\beta\_\{3\}\\psi\_\{0\}^\{\\prime\\prime\\prime\}\(t\_\{r\}\)\+\\beta\_\{4\}\\psi\_\{0\}^\{\(4\)\}\(t\_\{r\}\)\]\\cdot\\prod\_\{i=1\}^\{4\}\(\\psi\_\{1\}^\{\-1\}\)^\{\\prime\}\(u\_\{i\}\)\.

### A\.6Complexity Analysis

We bound the runtime of the Bell\-polynomial recursion by summing the per\-node powering cost over the tree, then specialize the bound to common shapes and confirm it against measured wall times\.

##### General bound\.

At each non\-root nodevvthe dominant cost is polynomial powering:dvd\_\{v\}iterations ofO​\(dv2\)O\(d\_\{v\}^\{2\}\)truncated convolution applied byvv’s parent, givingO​\(dv3\)O\(d\_\{v\}^\{3\}\)\. The Cauchy product atvvcostsO​\(mv⋅dv⋅dmax​\(v\)\)O\(m\_\{v\}\\cdot d\_\{v\}\\cdot d\_\{\\max\}\(v\)\)withdmax​\(v\)=maxc⁡dcd\_\{\\max\}\(v\)=\\max\_\{c\}d\_\{c\}; becausemv≤dvm\_\{v\}\\leq d\_\{v\}, this Cauchy cost never exceedsO​\(dv3\)O\(d\_\{v\}^\{3\}\), so the powering term dominates\. Computing the edge Taylor vectorpv​cp\_\{vc\}that the powering step consumes addsO​\(dv3\)O\(d\_\{v\}^\{3\}\)via the default implicit solve of[Algorithm˜2](https://arxiv.org/html/2605.23134#alg2), orO​\(dv2\)O\(d\_\{v\}^\{2\}\)when a closed\-formhv​ch\_\{vc\}is registered \([Section˜A\.4](https://arxiv.org/html/2605.23134#A1.SS4)\); theO​\(dv3\)O\(d\_\{v\}^\{3\}\)powering bound absorbs both\. At the rootrr, the powering step is replaced by a single order\-ddjet call of costO​\(d2\)O\(d^\{2\}\)\. The total cost is therefore

T​\(𝒯\)=O​\(d2\+∑v≠rdv3\),T\(\\mathcal\{T\}\)=O\\\!\\left\(d^\{2\}\+\\sum\_\{v\\neq r\}d\_\{v\}^\{3\}\\right\),\(24\)wheredv≡\|δv\|d\_\{v\}\\equiv\|\\delta\_\{v\}\|counts uncensored leaves in the subtree ofvv, and the bound holds uniformly across the three edge\-Taylor implementations\.

##### Height bound\.

At each depth level the subtrees are disjoint, so thedvd\_\{v\}values at that level sum to at mostdd\. Since eachdv≤dd\_\{v\}\\leq d, we havedv3≤dv⋅d2d\_\{v\}^\{3\}\\leq d\_\{v\}\\cdot d^\{2\}, so the per\-level contribution satisfies∑v​at level​ℓdv3≤d2​∑vdv≤d3\\sum\_\{v\\,\\text\{at level\}\\,\\ell\}d\_\{v\}^\{3\}\\leq d^\{2\}\\sum\_\{v\}d\_\{v\}\\leq d^\{3\}\. Summing over thehhnon\-root internal levels gives∑v≠rdv3≤h⋅d3\\sum\_\{v\\neq r\}d\_\{v\}^\{3\}\\leq h\\cdot d^\{3\}, soT​\(𝒯\)=O​\(d2\+h⋅d3\)=O​\(h⋅d3\)T\(\\mathcal\{T\}\)=O\(d^\{2\}\+h\\cdot d^\{3\}\)=O\(h\\cdot d^\{3\}\)\. For the standard case of constant\-depth copula trees withh=1h=1or22, the algorithm runs inO​\(d3\)O\(d^\{3\}\)time\.

##### d\\sqrt\{d\}\-decomposed trees\.

When every non\-root subtree satisfiesdv≤c​dd\_\{v\}\\leq c\\sqrt\{d\}for a constantcc, the per\-level bound sharpens:dv2≤c2​dd\_\{v\}^\{2\}\\leq c^\{2\}dgivesdv3≤c2​d⋅dvd\_\{v\}^\{3\}\\leq c^\{2\}\\,d\\cdot d\_\{v\}, so∑v​at level​ℓdv3≤c2​d​∑vdv≤c2​d2\\sum\_\{v\\,\\text\{at level\}\\,\\ell\}d\_\{v\}^\{3\}\\leq c^\{2\}\\,d\\sum\_\{v\}d\_\{v\}\\leq c^\{2\}\\,d^\{2\}\. Summing over thehhnon\-root internal levels yields∑v≠rdv3≤h​c2​d2\\sum\_\{v\\neq r\}d\_\{v\}^\{3\}\\leq h\\,c^\{2\}\\,d^\{2\}, and the total cost collapses toT​\(𝒯\)=O​\(d2\)T\(\\mathcal\{T\}\)=O\(d^\{2\}\)for constant depth—matching theO​\(d2\)O\(d^\{2\}\)root\-jet cost\. More generally, theT​\(𝒯\)=O​\(d2\)T\(\\mathcal\{T\}\)=O\(d^\{2\}\)regime holds whenever∑v≠rdv3=O​\(d2\)\\sum\_\{v\\neq r\}d\_\{v\}^\{3\}=O\(d^\{2\}\), for whichdv=O​\(d\)d\_\{v\}=O\(\\sqrt\{d\}\)is a sufficient condition\.

*Two\-level symmetric illustration\.*As a concrete instance, take a two\-level tree withMMsectors ofK=d/MK=d/Mleaves each\. Substituting into[Equation˜24](https://arxiv.org/html/2605.23134#A1.E24)givesO​\(M​K3\+M2​K2\)=O​\(d3/M2\+d2\)O\(MK^\{3\}\+M^\{2\}K^\{2\}\)=O\(d^\{3\}/M^\{2\}\+d^\{2\}\); thed3/M2d^\{3\}/M^\{2\}term dominates whenMMgrows slower thand\\sqrt\{d\}, and the cost collapses toO​\(d2\)O\(d^\{2\}\)whenM=Ω​\(d\)M=\\Omega\(\\sqrt\{d\}\), equivalentlyK=O​\(d\)K=O\(\\sqrt\{d\}\)—the two\-level instance of the general\|δv\|=O​\(d\)\|\\delta\_\{v\}\|=O\(\\sqrt\{d\}\)condition above\. TakingM=2M=2sectors withK=d/2K=d/2leaves each givesO​\(d3\)O\(d^\{3\}\), notO​\(d2\)O\(d^\{2\}\)\.

##### Special cases\.

The following table summarises cost for common tree shapes\.

##### Comparison with alternative approaches\.

The following table contrasts our method with existing alternatives\.

##### Empirical scaling\.

TheO​\(d3\)O\(d^\{3\}\)worst case assumes the per\-sector leaf countKKgrows withdd\. In practice, the application domain fixes the tree branching factor—K=5K=5in our scaling benchmarks—so the number of sectorsM=d/KM=d/Kgrows instead\. The per\-sector polynomial powering then costsO​\(K3\)O\(K^\{3\}\), constant indd, and the total across all sectors isM⋅O​\(K3\)=O​\(d⋅K2\)=O​\(d\)M\\cdot O\(K^\{3\}\)=O\(d\\cdot K^\{2\}\)=O\(d\)\. The Cauchy product convolvesMMshort kernels of lengthKKwith add\-length accumulator atO​\(d⋅K⋅M\)=O​\(d2\)O\(d\\cdot K\\cdot M\)=O\(d^\{2\}\)\. Both terms carry small constants relative to the per\-scan\-iteration overhead inlax\.scan, so the wall\-time transition where the quadratic term begins to dominate sits atd∗≈103d^\{\*\}\\approx 10^\{3\}at fixedKK\.[Appendix˜O](https://arxiv.org/html/2605.23134#A15)reports the empirical slopes that follow this prediction: sub\-linear in the small\-ddregime, transient super\-quadratic during the L1/L2 cache transition, and the cleanO​\(d2\)O\(d^\{2\}\)asymptote atd≥4,000d\\geq 4\{,\}000for both the fixed\-K=5K\{=\}5and balancedd\\sqrt\{d\}topologies\.

##### Mask\-aware static\-shape JIT\.

acopulapads each node’sα\\alpha\-vector to the mask\-independent length\|Vv\|\+1\|V\_\{v\}\|\+1, whereVvV\_\{v\}is the full leaf set in the subtree ofvv\. Each child ofvvthen dispatches at runtime throughlax\.switchto one of three branches—censored leaf, uncensored leaf, or subtree—so a single JIT\-compiled program covers every per\-observation mask without recompilation\.

##### Censoring reduces runtime\.

The cost bound in[Equation˜24](https://arxiv.org/html/2605.23134#A1.E24)depends ondvd\_\{v\}, the number ofuncensoredleaves in each subtree—not on the total number of leaves\. Censoring a variable sets itsα\\alpha\-vector to\[1\]\[1\], which shortens the polynomial powering and Cauchy\-product steps at every ancestor\. To quantify this effect, we time thed=50d\{=\}50nested Clayton configuration with 10 sectors of 5 leaves onn=100n\{=\}100observations at four censoring rates on a single AMD Threadripper 3960X CPU core\.[Table˜8](https://arxiv.org/html/2605.23134#A1.T8)reports the results: 75% censoring reduces the effective dimension from 50 to 13 and delivers a13\.6×13\.6\\timesspeedup, in line with theO​\(dv3\)O\(d\_\{v\}^\{3\}\)scaling of the per\-vertex Bell recursion that dominates the cost at this dimension\. TheO​\(\|Vv\|2\)O\(\|V\_\{v\}\|^\{2\}\)Cauchy\-product floor \([Section˜3\.2](https://arxiv.org/html/2605.23134#S3.SS2.SSS0.Px2)\) remains a strict lower bound but does not bind atd=50d\{=\}50on CPU\.

Table 8:Wall\-clock time for a nested Clayton \(d=50d\{=\}50,n=100n\{=\}100\) log\-likelihood evaluation at varying censoring rates on a single AMD Threadripper 3960X CPU core\.deffd\_\{\\mathrm\{eff\}\}counts uncensored variables\. Median of 20 JIT\-compiled calls\.

## Appendix BNumerical Stability of Taylor\-Mode Derivatives

A naive float64 implementation of the Bell pipeline underflows long before reaching the dimensions our experiments target, because high\-order generator derivatives and the geometric powersp1jp\_\{1\}^\{j\}that the Bell transform multiplies in span hundreds of orders of magnitude\. This appendix specifies the log\-domain machinery that keeps every intermediate quantity representable and reports the stress tests that confirm its accuracy\.[Algorithm˜3](https://arxiv.org/html/2605.23134#alg3)states the log\-domain version of[Algorithm˜1](https://arxiv.org/html/2605.23134#alg1): per\-index rescaling to neutralize geometric decay in the polynomial powering, signed logsumexp for the Bell\-transform inner products and root assembly, and intermediate normalisation during the Cauchy product\.

Algorithm 3Log\-domain implementation of[Algorithm˜1](https://arxiv.org/html/2605.23134#alg1)\. All floating\-point computation uses float64\. We writeSignedLSEfor the signed logsumexp primitive that the subroutine below specifies\.- Require:Copula tree at rootrrwithdduncensored leaves, copula\-scale observations𝐮\\mathbf\{u\}withuℓ=Fℓ​\(xℓ\)u\_\{\\ell\}=F\_\{\\ell\}\(x\_\{\\ell\}\)for raw observationxℓx\_\{\\ell\}and marginal CDFFℓF\_\{\\ell\}
- Ensure:Joint log\-likelihoodlog⁡c​\(𝐮\)\+∑ℓ:δℓ=1log⁡fℓ​\(xℓ\)\\log c\(\\mathbf\{u\}\)\+\\sum\_\{\\ell:\\,\\delta\_\{\\ell\}=1\}\\log f\_\{\\ell\}\(x\_\{\\ell\}\)\(copula log\-density plus marginal log\-densities for uncensored leaves\)SubroutineNormalize​\(a\)\\textsc\{Normalize\}\(a\):⊳\\trianglerightPrevent underflow/overflow; keepa=eλ​a^a=e^\{\\lambda\}\\hat\{a\}
- s←maxj⁡\|aj\|s\\leftarrow\\max\_\{j\}\|a\_\{j\}\|;s←max⁡\(s,ϵmin\)s\\leftarrow\\max\(s,\\,\\sqrt\{\\epsilon\_\{\\min\}\}\) ⊳\\trianglerightϵmin\\epsilon\_\{\\min\}= smallest positive normal float64 \(≈2\.2×10−308\\approx 2\.2\{\\times\}10^\{\-308\}\); floor atϵmin≈1\.5×10−154\\sqrt\{\\epsilon\_\{\\min\}\}\\approx 1\.5\{\\times\}10^\{\-154\}protects∇\(a/s\)\\nabla\(a/s\)
- returna^←a/s\\hat\{a\}\\leftarrow a/s,λ←log⁡s\\lambda\\leftarrow\\log sSubroutineSignedLSE​\(ℓ1,…,ℓn;σ1,…,σn\)\\textsc\{SignedLSE\}\(\\ell\_\{1\},\\ldots,\\ell\_\{n\};\\;\\sigma\_\{1\},\\ldots,\\sigma\_\{n\}\):⊳\\trianglerightℓj=log⁡\|xj\|\\ell\_\{j\}=\\log\|x\_\{j\}\|,σj=sign⁡\(xj\)\\sigma\_\{j\}=\\operatorname\{sign\}\(x\_\{j\}\)
- m←maxj⁡ℓjm\\leftarrow\\max\_\{j\}\\ell\_\{j\}
- S←∑jσj​exp⁡\(ℓj−m\)S\\leftarrow\\sum\_\{j\}\\sigma\_\{j\}\\exp\(\\ell\_\{j\}\-m\)
- returnlog⁡\|S\|\+m\\log\|S\|\+m,sign⁡\(S\)\\operatorname\{sign\}\(S\) Step 1\. Bell transform\(α\(c\)\\alpha^\{\(c\)\}from childcc’sβ\(c\)\\beta^\{\(c\)\},dcd\_\{c\}uncensored descendants\):
- p0,p1,…,pdc←Jet​\(ψv−1∘ψc,tc,dc\)p\_\{0\},p\_\{1\},\\ldots,p\_\{d\_\{c\}\}\\leftarrow\\textsc\{Jet\}\\bigl\(\\psi\_\{v\}^\{\-1\}\\\!\\circ\\psi\_\{c\},\\;t\_\{c\},\\;d\_\{c\}\\bigr\) ⊳\\trianglerightpj=hv​c\(j\)​\(tc\)/j\!p\_\{j\}=h\_\{vc\}^\{\(j\)\}\(t\_\{c\}\)/j\!, sop0=hv​c​\(tc\)=uvp\_\{0\}=h\_\{vc\}\(t\_\{c\}\)=u\_\{v\}is the primal \(unused below\)Per\-index rescaling\(removes geometric decay ofp1jp\_\{1\}^\{j\}\):
- p^j←sign\(pj\)sign\(p1\)jexp\(log\|pj\|−jlog\|p1\|\)\\hat\{p\}\_\{j\}\\leftarrow\\operatorname\{sign\}\(p\_\{j\}\)\\,\\operatorname\{sign\}\(p\_\{1\}\)^\{j\}\\,\\exp\\\!\\bigl\(\\log\|p\_\{j\}\|\-j\\log\|p\_\{1\}\|\\bigr\)forj=1,…,dcj=1,\\ldots,d\_\{c\} ⊳\\trianglerightWriting𝒫​\(z\):=∑j=1dcpj​zj\\mathcal\{P\}\(z\):=\\sum\_\{j=1\}^\{d\_\{c\}\}p\_\{j\}\\,z^\{j\}, this rescaling produces𝒫^​\(z\)=𝒫​\(z/p1\)\\hat\{\\mathcal\{P\}\}\(z\)=\\mathcal\{P\}\(z/p\_\{1\}\)Polynomial powering\(all powers via scan\):
- Q^\(1\)←P^\\hat\{Q\}^\{\(1\)\}\\leftarrow\\hat\{P\}
- fork=2,…,dck=2,\\ldots,d\_\{c\}do
- Q^\(k\)←Q^\(k−1\)∗P^\\hat\{Q\}^\{\(k\)\}\\leftarrow\\hat\{Q\}^\{\(k\-1\)\}\*\\hat\{P\} ⊳\\trianglerightTruncated convolution to degreedcd\_\{c\}
- Signed logsumexp assembly\(incorporatesp1jp\_\{1\}^\{j\}factor in log domain\):
- wj←log⁡j\!\+log⁡\|βj\(c\)\|\+j​log⁡\|p1\|w\_\{j\}\\leftarrow\\log j\!\+\\log\|\\beta^\{\(c\)\}\_\{j\}\|\+j\\log\|p\_\{1\}\|,σj\(w\)←sign\(βj\(c\)\)sign\(p1\)j\\sigma\_\{j\}^\{\(w\)\}\\leftarrow\\operatorname\{sign\}\(\\beta^\{\(c\)\}\_\{j\}\)\\,\\operatorname\{sign\}\(p\_\{1\}\)^\{j\}forj=1,…,dcj=1,\\ldots,d\_\{c\} ⊳\\trianglerightPer\-jjlog weight; both arekk\-independent, computed once before the inner loop
- fork=1,…,dck=1,\\ldots,d\_\{c\}do
- ℓj←wj\+log⁡\|Q^j\(k\)\|\\ell\_\{j\}\\leftarrow w\_\{j\}\+\\log\|\\hat\{Q\}^\{\(k\)\}\_\{j\}\|;σj←σj\(w\)⋅sign⁡\(Q^j\(k\)\)\\sigma\_\{j\}\\leftarrow\\sigma\_\{j\}^\{\(w\)\}\\cdot\\operatorname\{sign\}\(\\hat\{Q\}^\{\(k\)\}\_\{j\}\)forj=1,…,dcj=1,\\ldots,d\_\{c\} ⊳\\trianglerightQ^j\(k\)\\hat\{Q\}^\{\(k\)\}\_\{j\}= coefficient ofzjz^\{j\}in𝒫^​\(z\)k\\hat\{\\mathcal\{P\}\}\(z\)^\{k\}
- Lk,Sk←SignedLSE​\(ℓ1,…,ℓdc;σ1,…,σdc\)L\_\{k\},S\_\{k\}\\leftarrow\\textsc\{SignedLSE\}\(\\ell\_\{1\},\\ldots,\\ell\_\{d\_\{c\}\};\\;\\sigma\_\{1\},\\ldots,\\sigma\_\{d\_\{c\}\}\)
- log⁡\|αk\(c\)\|←Lk−log⁡k\!\\log\|\\alpha^\{\(c\)\}\_\{k\}\|\\leftarrow L\_\{k\}\-\\log k\!;sign⁡\(αk\(c\)\)←Sk\\operatorname\{sign\}\(\\alpha^\{\(c\)\}\_\{k\}\)\\leftarrow S\_\{k\}
- 𝑟𝑒𝑓←maxk⁡log⁡\|αk\(c\)\|\\mathit\{ref\}\\leftarrow\\max\_\{k\}\\log\|\\alpha^\{\(c\)\}\_\{k\}\| ⊳\\trianglerightReference scale for this child
- αk\(c\)←Sk​exp⁡\(log⁡\|αk\(c\)\|−𝑟𝑒𝑓\)\\alpha^\{\(c\)\}\_\{k\}\\leftarrow S\_\{k\}\\exp\(\\log\|\\alpha^\{\(c\)\}\_\{k\}\|\-\\mathit\{ref\}\)fork=1,…,dck=1,\\ldots,d\_\{c\};α0\(c\)←0\\alpha^\{\(c\)\}\_\{0\}\\leftarrow 0 ⊳\\trianglerightInternal child: index\-0 entry unused in subsequent Cauchy product; differs from censored\-leaf base case whereα0=1\\alpha\_\{0\}=1acts as the product identity Step 2\. Cauchy product\(β\(v\)\\beta^\{\(v\)\}fromα\(c1\),…,α\(cM\)\\alpha^\{\(c\_\{1\}\)\},\\ldots,\\alpha^\{\(c\_\{M\}\)\}\):
- β\(v\)←α\(c1\)\\beta^\{\(v\)\}\\leftarrow\\alpha^\{\(c\_\{1\}\)\};Λ←0\\Lambda\\leftarrow 0 ⊳\\trianglerightCumulative log\-scale
- fori=2,…,Mi=2,\\ldots,Mdo
- β\(v\)←β\(v\)∗α\(ci\)\\beta^\{\(v\)\}\\leftarrow\\beta^\{\(v\)\}\*\\alpha^\{\(c\_\{i\}\)\} ⊳\\trianglerightTruncated polynomial product to degreedvd\_\{v\}
- β\(v\),λ←Normalize​\(β\(v\)\)\\beta^\{\(v\)\},\\,\\lambda\\leftarrow\\textsc\{Normalize\}\(\\beta^\{\(v\)\}\);Λ←Λ\+λ\\Lambda\\leftarrow\\Lambda\+\\lambda ⊳\\trianglerightRe\-normalize after each convolution
- Step 3\. Root assembly\(log\-density from root’sβ\(r\)\\beta^\{\(r\)\}\):
- q0,…,qd←Jet​\(ψr,tr,d\)q\_\{0\},\\ldots,q\_\{d\}\\leftarrow\\textsc\{Jet\}\(\\psi\_\{r\},\\,t\_\{r\},\\,d\) ⊳\\trianglerightqk=ψr\(k\)​\(tr\)/k\!q\_\{k\}=\\psi\_\{r\}^\{\(k\)\}\(t\_\{r\}\)/k\!
- ℓk←log⁡\|βk\(r\)\|\+log⁡\|qk\|\+log⁡k\!\\ell\_\{k\}\\leftarrow\\log\|\\beta^\{\(r\)\}\_\{k\}\|\+\\log\|q\_\{k\}\|\+\\log k\!fork=0,…,dk=0,\\ldots,d
- σk←sign⁡\(βk\(r\)\)​sign⁡\(qk\)\\sigma\_\{k\}\\leftarrow\\operatorname\{sign\}\(\\beta^\{\(r\)\}\_\{k\}\)\\,\\operatorname\{sign\}\(q\_\{k\}\)
- 𝑙𝑜𝑔𝑟𝑜𝑜𝑡,\_←SignedLSE​\(ℓ0,…,ℓd;σ0,…,σd\)\\mathit\{logroot\},\\,\\\_\\leftarrow\\textsc\{SignedLSE\}\(\\ell\_\{0\},\\ldots,\\ell\_\{d\};\\;\\sigma\_\{0\},\\ldots,\\sigma\_\{d\}\) Step 4\. Leaf product:
- 𝑙𝑜𝑔𝑙𝑒𝑎𝑓←∑ℓ:δℓ=1log⁡\|\(ψpar​\(ℓ\)−1\)′​\(uℓ\)\|\\mathit\{logleaf\}\\leftarrow\\sum\_\{\\ell:\\,\\delta\_\{\\ell\}=1\}\\log\\bigl\|\(\\psi\_\{\\mathrm\{par\}\(\\ell\)\}^\{\-1\}\)^\{\\prime\}\(u\_\{\\ell\}\)\\bigr\| ⊳\\trianglerightViajax\.gradofψ−1\\psi^\{\-1\}; vmap over uncensored leaves Step 5\. Log\-scale accumulation and final result\(Steps 1–2 apply at every internal node bottom\-up; each𝑟𝑒𝑓c\\mathit\{ref\}\_\{c\}andΛcauchy,v\\Lambda\_\{\\mathrm\{cauchy\},v\}produced below feeds a single global accumulatorΛtotal\\Lambda\_\{\\mathrm\{total\}\}\):
- Λtotal←∑v∈internal\(∑c∈ch​\(v\)𝑟𝑒𝑓c\+Λcauchy,v\)\\Lambda\_\{\\mathrm\{total\}\}\\leftarrow\\sum\_\{v\\in\\mathrm\{internal\}\}\\\!\\Bigl\(\\sum\_\{c\\in\\mathrm\{ch\}\(v\)\}\\mathit\{ref\}\_\{c\}\\;\+\\;\\Lambda\_\{\\mathrm\{cauchy\},v\}\\Bigr\) ⊳\\trianglerightSum every per\-child reference scale and every per\-node Cauchy normalisation across the whole tree
- return𝑙𝑜𝑔𝑟𝑜𝑜𝑡\+Λtotal\+𝑙𝑜𝑔𝑙𝑒𝑎𝑓\+∑ℓ:δℓ=1log⁡fℓ​\(xℓ\)\\mathit\{logroot\}\+\\Lambda\_\{\\mathrm\{total\}\}\+\\mathit\{logleaf\}\+\\sum\_\{\\ell:\\,\\delta\_\{\\ell\}=1\}\\log f\_\{\\ell\}\(x\_\{\\ell\}\) ⊳\\trianglerightAdd marginal log\-likelihoods for uncensored leaves

##### Generator derivative growth\.

Each family’s derivative growth rate sets the precision ceiling that the log\-domain machinery must accommodate\. For Clayton,ψ​\(t\)=\(1\+t\)−1/θ\\psi\(t\)=\(1\+t\)^\{\-1/\\theta\}, the Pochhammer identityψ\(k\)​\(t\)=\(−1\)k​\(1/θ\)k​\(1\+t\)−1/θ−k\\psi^\{\(k\)\}\(t\)=\(\-1\)^\{k\}\\,\(1/\\theta\)\_\{k\}\\,\(1\+t\)^\{\-1/\\theta\-k\}gives\|ψ\(k\)​\(t\)\|/k\!=O​\(k1/θ−1\)\|\\psi^\{\(k\)\}\(t\)\|/k\!=O\(k^\{1/\\theta\-1\}\)via Stirling, so the jet\-stored Taylor coefficients grow only polynomially inkkfor fixedθ\\thetaandtt\. Joe and other double\-exponential generators grow faster because the chain rule nests an exponential inside another\. Empirically, every family matches Mathematica to relative error below3×10−113\\times 10^\{\-11\}whereverD\[\]completes \([Table˜10](https://arxiv.org/html/2605.23134#A2.T10)\), and[Table˜9](https://arxiv.org/html/2605.23134#A2.T9)confirms finite log\-density evaluation up tod=60d=60for nine of the ten families and up tod=4d=4for Nelsen9, whosedd\-completely monotone parameter range tightens withdd\.

##### Edge composition safeguards\.

Theedge compositionh=ψouter−1∘ψinnerh=\\psi\_\{\\mathrm\{outer\}\}^\{\-1\}\\circ\\psi\_\{\\mathrm\{inner\}\}that the nested density requires admits three qualitatively different evaluation paths across our ten families \([Section˜A\.4](https://arxiv.org/html/2605.23134#A1.SS4)\), and Frank in particular needs a closed\-form rewrite to stay within float64\. Clayton, AMH, Joe, Gumbel, Nelsen9, and Nelsen12 provide closed\-form inverses, so the composition introduces no numerical root\-finding error\. InverseGaussian relies on Oryx’s symbolic inversion, which pattern\-matches individual primitives to their algebraic inverses\. Two families also expose a numerical hazard in thenestedcompositionψv−1​\(ψc​\(t\)\)\\psi\_\{v\}^\{\-1\}\(\\psi\_\{c\}\(t\)\)even with closed\-form inverses: Frank atd=98d=98overflows because intermediate log\-sum\-exp terms exceed float64, and Nelsen9 underflows becauseψc​\(t\)\\psi\_\{c\}\(t\)collapses to zero beforeψv−1\\psi\_\{v\}^\{\-1\}can recover it\. Both families ship a registered closed\-form edge composition \([Section˜A\.4](https://arxiv.org/html/2605.23134#A1.SS4)\) that algebraically cancels the offending exponentials, bypassing the default implicit\-differentiation solve of[Algorithm˜2](https://arxiv.org/html/2605.23134#alg2)\.

##### Precision at large sector size\.

Atd=60d\{=\}60with sector sizeK=30K\{=\}30,[Appendix˜J](https://arxiv.org/html/2605.23134#A10)reports a3\.5×10−53\.5\\times 10^\{\-5\}discrepancy against R’snacLL\. Mathematica’s symbolic differentiation of the copula CDF confirms that our framework matches the ground truth to machine epsilon atd=60d\{=\}60; the discrepancy originates in R’snacLL, not in our polynomial powering\.

##### Generator\-specific issues\.

Nelsen9’s double\-exponential structureψ​\(t\)=exp⁡\(\(1−et\)/θ\)\\psi\(t\)=\\exp\\bigl\(\(1\-e^\{t\}\)/\\theta\\bigr\)underflows to zero in float64 whenever the exponent\(1−et\)/θ\(1\-e^\{t\}\)/\\thetadrops below−745\-745, which occurs at moderatettfor smallθ\\theta\. We mitigate the underflow on three fronts\. A per\-family Taylor hook supplies the coefficients via the Touchard factorizationψ\(n\)​\(t\)=ψ​\(t\)⋅Tn​\(−et/θ\)\\psi^\{\(n\)\}\(t\)=\\psi\(t\)\\cdot T\_\{n\}\(\-e^\{t\}/\\theta\), where Faà di Bruno onψ=exp∘g\\psi=\\exp\\circ gwithg\(k\)​\(t\)=−et/θg^\{\(k\)\}\(t\)=\-e^\{t\}/\\thetafor everyk≥1k\\geq 1collapses the complete Bell polynomial at equal arguments to the Touchard polynomialTnT\_\{n\}\. The framework evaluateslog⁡\|ψ\(n\)​\(t\)\|=\(1−et\)/θ\+log⁡\|Tn​\(−et/θ\)\|\\log\|\\psi^\{\(n\)\}\(t\)\|=\(1\-e^\{t\}\)/\\theta\+\\log\|T\_\{n\}\(\-e^\{t\}/\\theta\)\|, which keeps the underflow\-prone exponential symbolic and confines the remaining float64 work to the polynomialTnT\_\{n\}, whose direct sum loses at mostnndecimal digits to alternating\-sign cancellation in thedd\-monotone\-valid range\. The hook bypasses jet on the closed\-formψ\\psi, whose JVP produces a NaN gradient whenψ\\psiunderflows\. The analyticalψ−1​\(u\)=log⁡\(1−θ​log⁡u\)\\psi^\{\-1\}\(u\)=\\log\(1\-\\theta\\log u\)replaces Oryx’s symbolic inverse, which we found to be numerically unstable at largett\. The closed\-form nested compositionh​\(t\)=t\+log⁡\(r\+\(1−r\)​e−t\)h\(t\)=t\+\\log\\\!\\bigl\(r\+\(1\-r\)\\,e^\{\-t\}\\bigr\)withr=θv/θcr=\\theta\_\{v\}/\\theta\_\{c\}replaces the naiveψv−1​\(ψc​\(t\)\)\\psi\_\{v\}^\{\-1\}\(\\psi\_\{c\}\(t\)\)chain that would underflow at moderatett, and the framework registers it as the closed\-form override for Nelsen9/Nelsen9 nestings in[Section˜A\.4](https://arxiv.org/html/2605.23134#A1.SS4)\.

##### Stress test across parameter regimes\.

[Table˜9](https://arxiv.org/html/2605.23134#A2.T9)reports a systematic stress test that sweeps the \(family, strength, dimension\) grid\. We map four target Kendall\-tau strengthsτ∈\{0\.10,0\.40,0\.70,0\.90\}\\tau\\in\\\{0\.10,\\,0\.40,\\,0\.70,\\,0\.90\\\}to each family’s parameter via numerical inversion when no closed form exists, build a two\-sector nested same\-family copula at sector sizeK∈\{2,5,10,15,20,25,30\}K\\in\\\{2,5,10,15,20,25,30\\\}for total dimensiond=2​Kd\{=\}2K, and evaluate the log\-density at100100randomuu\-vectors per cell\. A cellpasseswhen every evaluation returns a finite value\.

Table 9:Stability stress test: number ofKKvalues out of seven that pass at100%100\\%finiteness on100100randomuu\-vectors per cell\. Each row sweepsK∈\{2,5,10,15,20,25,30\}K\\in\\\{2,5,10,15,20,25,30\\\}for the listed Kendall\-tau strength\. AMH‡has bounded tauτ≤1/3\\tau\\leq 1/3, so itsτ=0\.40,0\.70,0\.90\\tau\{=\}0\.40,0\.70,0\.90columns collapse to the same parameter; we still report each column to keep the grid uniform\.Across all253253cells in the sweep that fall within each family’sdd\-completely\-monotone parameter range—ten generator families, four strengths, seven dimensions—every cell returns finite log\-densities on all100100random evaluation points, a100%100\\%pass rate\. Nelsen9 is a special case: thedd\-CM upper boundθmax​\(d\)\\theta\_\{\\max\}\(d\)contracts asddgrows \(θmax​\(4\)=0\.22\\theta\_\{\\max\}\(4\)\{=\}0\.22,θmax​\(10\)=0\.057\\theta\_\{\\max\}\(10\)\{=\}0\.057,θmax​\(60\)=0\.007\\theta\_\{\\max\}\(60\)\{=\}0\.007\), so only theτ=0\.10\\tau\{=\}0\.10,K=2K\{=\}2cell \(d=4d\{=\}4,θ≤0\.2\\theta\{\\leq\}0\.2\) lies inside the valid range; the other2727cells skip\.

We attribute the broad success to the log\-domain rescaling and signed\-logsumexp arithmetic of[Algorithm˜3](https://arxiv.org/html/2605.23134#alg3): per\-index rescaling absorbs the geometric blow\-up that the naive recursion would produce at high dimension, and signed logsumexp prevents catastrophic cancellation in the Bell\-transform inner products\.‡AMH’s bounded tauτ≤1/3\\tau\\leq 1/3collapses the three high\-strength columns to the same parameter, which is why all four columns repeat identically\.†Marks the five extra families that existing nested\-copula tools do not cover\.§Nelsen9 isdd\-completely monotone only forθ≤θmax​\(d\)=1/\|rd\|\\theta\\leq\\theta\_\{\\max\}\(d\)=1/\|r\_\{d\}\|, whererdr\_\{d\}is the most negative root of thedd\-th Touchard polynomialTd​\(y\)=∑jS​\(d,j\)​yjT\_\{d\}\(y\)=\\sum\_\{j\}S\(d,j\)\\,y^\{j\}; the other 27 cells in this grid violate that bound\.

##### Cross\-family accuracy verification\.

[Table˜10](https://arxiv.org/html/2605.23134#A2.T10)compares our framework against Mathematica’s symbolic differentiation of the copula CDF across all ten families\. Mathematica evaluates each cell at100100digit working precision; our framework runs in IEEE\-754float64\(≈16\\approx 16digits\) and matches the Mathematica ground truth to relative error below3×10−113\\times 10^\{\-11\}in every configuration throughd=60d\{=\}60whereD\[\]completes\.

The framework exposes a per\-family hook,Copula\.generator\_taylor\_coefficients\(t, k\), that returns the Taylor expansion\[ψ​\(t\),ψ′​\(t\)/1\!,…,ψ\(k\)​\(t\)/k\!\]\[\\psi\(t\),\\,\\psi^\{\\prime\}\(t\)/1\!,\\,\\ldots,\\,\\psi^\{\(k\)\}\(t\)/k\!\]directly, bypassing chain\-rule cancellations Taylor\-mode AD on the closed\-formψ\\psiwould produce at high derivative order\. The framework consults the hook everywhere it would otherwise call jet onψ\\psi: at the root density assembly and at every nested compositionh​\(t\)=ψouter−1∘ψinner​\(t\)h\(t\)=\\psi\_\{\\mathrm\{outer\}\}^\{\-1\}\\\!\\circ\\psi\_\{\\mathrm\{inner\}\}\(t\)\. For AMH \(ψ​\(t\)=\(1−θ\)/\(et−θ\)\\psi\(t\)=\(1\-\\theta\)/\(e^\{t\}\-\\theta\)\), Taylor\-mode AD on the closed form costs roughly two decimal digits per ten dimensions infloat64; the override sums the same\-sign Geometric\(1−θ\)\(1\-\\theta\)frailty seriesψ\(k\)​\(t\)=\(−1\)k​∑x≥1xk​\(1−θ\)​θx−1​e−t​x\\psi^\{\(k\)\}\(t\)=\(\-1\)^\{k\}\\sum\_\{x\\geq 1\}x^\{k\}\(1\-\\theta\)\\theta^\{x\-1\}e^\{\-tx\}in log\-domain via signed\-logsumexp, which has no cancellation and recovers machine precision\(≤2×10−15\)\(\\leq 2\\times 10^\{\-15\}\)at everyK≤30K\\leq 30\. Other completely monotone families \(Frank, Joe, Gumbel, InverseGaussian\) admit similar overrides via their respective frailty representations; we provide AMH because it was the only family whose chain\-rule path materially hit thefloat64cancellation wall in the verified range\.

Table 10:Numerical accuracy and symbolic\-differentiation scaling on two\-sector nested same\-family copulas\. “Max rel\. err\.” is the largest relative error between our framework and Mathematica’s symbolicD\[\]across the sector sizes where Mathematica completes; “maxKK” is the largest sector size Mathematica reaches within a 300 s timeout; “leaves” is theLeafCountof Mathematica’s symbolic derivative at thatKK; “D\[\]time” is the wall time to compute it\.The scaling column reveals the core problem with symbolic differentiation: the expression size depends on the generator’s algebraic complexity in ways that do not track with dimension alone\. Clayton’s power\-law generator and AMH’s rational generator remain tractable atK=30K\{=\}30with5\.1×1055\.1\{\\times\}10^\{5\}and2\.4×1072\.4\{\\times\}10^\{7\}leaves respectively, but every other family blows up exponentially—Joe reaches108108million leaves atK=15K\{=\}15, InverseGaussian164164million, and Nelsen127676million atK=5K\{=\}5\. Nelsen12 times out atK=10K\{=\}10; Nelsen13 and Nelsen17 reachK=15K\{=\}15before timing out at the same exponential rate\. Our polynomial\-powering algorithm avoids this family\-dependent cost: it processes all ten families inO​\(d3\)O\(d^\{3\}\)operations regardless of the generator’s algebraic form, because the cost depends only on the polynomial degreeKK, not on the symbolic complexity\.

## Appendix CLog\-Space Taylor\-Mode Differentiation

The log\-domain machinery of[Appendix˜B](https://arxiv.org/html/2605.23134#A2)keeps the Bell pipeline’s polynomial coefficients finite at highdd, but it operates on*outputs*ofjet\_array\.jetrather than inside the jet expansion itself\. For Frank and Joe generators on extreme\-tail leaves, that boundary is exactly where the failure now appears: the forward jet returns finite\-but\-denormal Taylor coefficients, and reverse\-mode AD then materialises1/c1/cfor those denormalcc, overflowing to∞\\inftyand producing∞⋅0\\infty\\cdot 0in the nextdot\_generalstep\. We resolve the failure by representing every Taylor coefficient as a\(sign,log⁡\|c\|\)\(\\text\{sign\},\\log\|c\|\)pair end\-to\-end; thelog\-space jet, opt\-in viacompile\_model\(…, log\_jet=True\)\.

### C\.1Why naive float64 fails for Frank\-type generators

For an Archimedean generator that decays asψ​\(t\)≈c​e−t\\psi\(t\)\\approx c\\,e^\{\-t\}at largett\(Frank, Joe\), thekk\-th Taylor coefficient at the composed argumenttrt\_\{r\}behaves likec​e−tr/k\!c\\,e^\{\-t\_\{r\}\}/k\!times a polynomial inkk\. Oncetr≳500t\_\{r\}\\gtrsim 500— which ford=75d\{=\}75flat copies ofu=10−3u\{=\}10^\{\-3\}givestr=−∑log⁡uℓ≈518t\_\{r\}=\-\\sum\\log u\_\{\\ell\}\\approx 518— the high\-order coefficients drop into float64 denormal range \(\|ck\|<2\.2×10−308\|c\_\{k\}\|<2\.2\\times 10^\{\-308\}\) byk≳60k\\gtrsim 60\.

The forward computation tolerates the denormals\. Backward does not: the upstream rules injet\_array\(notably\_log\_taylorand the body of\_div\_taylor\_rule\) recursively divide byccin their Faà di Bruno recurrences, and the Jacobian\-vector product of1/c1/catc∼10−310c\\sim 10^\{\-310\}is∼10310\\sim 10^\{310\}, beyondfloat64\.max\. The resulting∞\\inftythen meets a structural zero — for example, the constant divisor’s higher Taylor terms inu/θu/\\theta— anddot\_generalreturns NaN\. The minimal reproducer isd=75d\{=\}75flat\-Frank with all leaves atu=10−3u\{=\}10^\{\-3\}; raw forward log\-density is39\.9239\.92, raw gradient isNaN\\mathrm\{NaN\}\.

### C\.2Architecture:\(sign,log⁡\|c\|\)\(\\text\{sign\},\\log\|c\|\)throughout

We replace the raw float series injet\_array\.jetwith aLogSeriespair\(σ,λ\)\(\\sigma,\\lambda\)whose entries satisfyc=σ⋅eλc=\\sigma\\cdot e^\{\\lambda\}forσ∈\{−1,0,\+1\}\\sigma\\in\\\{\-1,0,\+1\\\}andλ∈ℝ\\lambda\\in\\mathbb\{R\}\. Structural zero encodes as\(σ,λ\)=\(0,−∞\)\(\\sigma,\\lambda\)=\(0,\-\\infty\), which acts as the additive identity𝟎\\mathbf\{0\}\. Multiplication and division become sign products and log\-magnitude addition or subtraction; addition uses signed logsumexp on the non\-zero support, with𝟎⊕x=x\\mathbf\{0\}\\oplus x=xby definition:

\(σa,λa\)⊕\(σb,λb\)=\(sign⁡\(S\),λM\+log⁡\|S\|\),λM=maxi∈A∗⁡λi,S=∑i∈A∗σi​eλi−λM,\(\\sigma\_\{a\},\\lambda\_\{a\}\)\\oplus\(\\sigma\_\{b\},\\lambda\_\{b\}\)\\;=\\;\\bigl\(\\operatorname\{sign\}\(S\),\\;\\lambda\_\{M\}\+\\log\|S\|\\bigr\),\\quad\\lambda\_\{M\}=\\max\_\{i\\in A^\{\*\}\}\\lambda\_\{i\},\\quad S=\\sum\_\{i\\in A^\{\*\}\}\\sigma\_\{i\}\\,e^\{\\lambda\_\{i\}\-\\lambda\_\{M\}\},\(25\)whereA∗=\{i∈\{a,b\}:σi≠0\}A^\{\*\}=\\\{i\\in\\\{a,b\\\}:\\sigma\_\{i\}\\neq 0\\\}is the active support; for emptyA∗A^\{\*\}the result is𝟎\\mathbf\{0\}, and a singleton reduces to its sole operand by inspection\. Inside each jet rule’s scan body — division, log, log1p, exp, expm1, integer power, and multiplication — the same recurrence runs, but every arithmetic operation dispatches to its log\-space sibling\. Coefficients with magnitude10±30010^\{\\pm 300\}now occupy normal float64 territory inλ\\lambda\(λ≈±690\\lambda\\approx\\pm 690\), so no JVP rule ever forms1/c1/cfor a denormalcc\.

### C\.3Crossing the bell boundary without leaving log\-space

A naiveraw on entry, log\-space inside, raw on exitwrapper re\-creates the failure:log\_to\_rawmaterialisesσ⋅eλ\\sigma\\cdot e^\{\\lambda\}at the boundary, and any genuinely\-tinyλ\\lambdabecomes a denormal float that the next downstreamjnp\.log\(\|c\|\)call differentiates as1/c→∞1/c\\to\\infty\. We extendjet\_array\.jetwithreturn\_log\_series=Trueso the output stays as aLogSeriespair, and rewritebell\.\_root\_assemblyto consume that pair directly: the existing signed\-logsumexp reduction readsλ\\lambdaaslog⁡\|c\|\\log\|c\|andσ\\sigmaas the sign, skipping thejnp\.log\(jnp\.abs\(⋅\)\)\\texttt\{jnp\.log\(jnp\.abs\}\(\\cdot\)\)step that was the residual NaN site\. The bell pipeline’s earlier stages already operate in log domain via\_normalize\_polyand the carriedtotal\_log\_scale, so no further bell changes are needed\.

### C\.4Validation

- •*Equivalence on benign inputs\.*Atd∈\{10,20,50\}d\\in\\\{10,20,50\\\}with moderateuu,log\_jet=Truematches the raw path on forward log\-density and gradient to relative error10−1010^\{\-10\}\.
- •*Stability where raw fails\.*Atd=75d\{=\}75flat\-Frank with allu=10−3u\{=\}10^\{\-3\}: raw gradientNaN\\mathrm\{NaN\}, log\-jet gradient29\.7029\.70\(finite\)\. At the nested 11\-by\-9 sp500\-style pattern with five extreme leaves per sector: raw gradientNaN\\mathrm\{NaN\}, log\-jet gradient finite; forward log\-density agrees to relative error10−710^\{\-7\}\.
- •*Production fit\.*S&P 500 d=98 Frank — the original failure that motivated this appendix — converges in1515L\-BFGS\-B function evaluations to NLL=−25749\.51=\-25749\.51withlog\_jet=True, improving on Clayton \(−25346\.78\-25346\.78\) at the same dimension while remaining behind Gumbel \(−26532\.01\-26532\.01\)\. Before the fix, the sp500 result table reportedconverged=Falsewith4444evaluations andNLL=−25897\.14\\mathrm\{NLL\}=\-25897\.14, a fictitious value the bell pipeline’s old NaN\-mask had silently swallowed \(see[Table˜6](https://arxiv.org/html/2605.23134#S5.T6)\)\.
- •*Regression\.*Theacopulatest suite passes withlog\_jet=Falseproducing identical bytes; the additionaltests/test\_log\_jet\.py\(six tests\) andtests/test\_log\_space\.py\(twenty\-two tests on theLogSeriesprimitives\) pass under JIT\.

### C\.5Cost and configurability

The log\-space rules add one scalar field to every Taylor coefficient and one sign\-tracking branch in each scan body, so memory doubles and per\-iteration arithmetic roughly doubles versus the raw path\. Compile time grows by3030–50%50\\%on the d=98 Frank fit \(3030seconds cold versus1818seconds for the raw d=75 case\)\. Steady\-state per\-evaluation runtime grows by a similar factor\. The kwarg defaults toFalse, so existing callers pay nothing; users of high\-ddFrank or Joe passlog\_jet=Trueonce and the rest of the pipeline adapts\.

## Appendix DPolynomial Multiplication Alternatives

We chose directO​\(K2\)O\(K^\{2\}\)convolution over the textbookO​\(K​log⁡K\)O\(K\\log K\)FFT path because FFT loses on both speed and precision in the regime copula applications occupy\. Our implementation multiplies truncated polynomials withlax\.conv\_general\_dilated, which executes a direct convolution\.[Table˜11](https://arxiv.org/html/2605.23134#A4.T11)confirms the choice empirically: direct convolution beats rfft\-based multiplication for every sector size up toK=192K\{=\}192and only loses atK=256K\{=\}256, far above the regime where copula applications operate\.

Table 11:Wall\-clock time per truncated polynomial\-powering scan \(p1,…,pKp^\{1\},\\ldots,p^\{K\}\) at sector sizeKK, JIT\-compiled in JAX on a single CPU thread\. Direct convolution wins by a factor of 2–5 for everyK≤192K\\leq 192; the FFT path’s relative error is unbounded pastK=64K\{=\}64, with the earlierK=16K\{=\}16spike \(∼0\.02\{\\sim\}0\.02\) andK=24K\{=\}24spike \(∼0\.56\{\\sim\}0\.56\) reflecting the alternating\-sign cancellation that destabilises the convolution before the magnitude blow\-up\.Three observations support this choice\. First, polynomial degrees in copula applications rarely exceedK=50K=50, where FFT’s constant overhead—complex\-arithmetic promotion, bit\-reversal permutations—wipes out its asymptotic savings\. Second, the short\-kernel optimisation in\_convolve\_short\_kernelalready reduces the Cauchy\-product step fromO​\(dv2\)O\(d\_\{v\}^\{2\}\)toO​\(dv⋅dmax​\(v\)\)O\(d\_\{v\}\\cdot d\_\{\\max\}\(v\)\), narrowing the gap further\. Third—decisively for our regime—FFT\-based polynomial multiplication satisfies only anorm\-wiseforward\-error bound\[Higham,[2002](https://arxiv.org/html/2605.23134#bib.bib14), Ch\. 24\],\|fl​\(a∗b\)k−\(a∗b\)k\|≤c​ϵ​log⁡K​‖a‖2​‖b‖2\|\\mathrm\{fl\}\(a\\ast b\)\_\{k\}\-\(a\\ast b\)\_\{k\}\|\\leq c\\,\\epsilon\\log K\\,\\\|a\\\|\_\{2\}\\\|b\\\|\_\{2\}, so every output entry inherits an absolute error proportional to thelargestcoefficient in the inputs\. Direct convolution, by contrast, preservescomponent\-wiseprecision: each output entry’s error scales with the coefficients that actually sum into it\. Bell polynomial coefficients within a singleβ\\beta\-vector span many orders of magnitude—factorial growth fromβ1\\beta\_\{1\}toβdv\\beta\_\{d\_\{v\}\}gives dynamic range above106010^\{60\}already atdv=50d\_\{v\}=50, far exceeding the∼1016\\sim\\\!10^\{16\}resolution of float64—so the small leading coefficients that the nested\-density formula multiplies againstψr\(k\)\\psi\_\{r\}^\{\(k\)\}at lowkklose all digits to FFT round\-off\. The repeated power scan compounds this loss at every step, so FFT\-based powering retains no usable precision well belowK=100K\{=\}100in float64 unless the implementation first rescales the polynomial into a narrow dynamic range—a rescaling that our log\-domain pipeline already performs for the direct path at zero asymptotic cost\. Should an application driveKKabove 200, swapping in an FFT path requires changing only two internal helpers, with no public\-API impact\.

##### Alternative Bell\-polynomial algorithms\.

Factorisation\-based algorithms compute partial Bell polynomials via recurrence relations or matrix factorisations, achieving the sameO​\(d3\)O\(d^\{3\}\)asymptotic complexity with different constants\. Our polynomial\-powering formulation sidesteps these alternatives because Taylor\-mode AD already produces jet coefficients in exactly the truncated\-power\-series representation that polynomial powering consumes, so no format conversion is needed\. Graph\-level AD optimisations such as cross\-country elimination\[Bettencourt et al\.,[2019](https://arxiv.org/html/2605.23134#bib.bib3)\]can reduce compilation overhead but remain orthogonal to the density\-level contribution of this work\.

## Appendix ESymbolic Inversion and Numeric Fallback

Two paths supply the inverse generatorψ−1\\psi^\{\-1\}that the density and sampler require: a symbolic derivation via Oryx for closed\-form generators, and a numeric bisection fallback for everything else\. Oryx\[Radul et al\.,[2023](https://arxiv.org/html/2605.23134#bib.bib37)\]derivesψ−1\\psi^\{\-1\}by tracing the generator through JAX’s intermediate representation and pattern\-matching each primitive against a table of known algebraic inverses—for example,exp↔log\\exp\\leftrightarrow\\log,x↦xa↔x↦x1/ax\\mapsto x^\{a\}\\leftrightarrow x\\mapsto x^\{1/a\}, andx↦x\+c↔x↦x−cx\\mapsto x\+c\\leftrightarrow x\\mapsto x\-c\. The procedure succeeds whenever the generator is a composition of individually invertible operations with no data\-dependent control flow\. All ten built\-in families satisfy this condition\.

Oryx’s pattern matcher cannot invert functions that contain data\-dependent branching—lax\.condorlax\.switch—because the branch the tracer follows at trace time may differ from the branch the runtime selects at evaluation time\. Functions that call numerical subroutines such as adaptive quadrature or iterative solvers also fail, because those routines expand into opaquewhile\_loopprimitives\. The matcher raises an error at trace time in either case, never silently returning an incorrect result\.

The framework probes each generator at initialisation by evaluating the Oryx\-derived inverse on a concrete test value\. When the probe succeeds, the symbolic inverse handles all subsequent calls\. When it fails, the framework switches transparently to a numeric path that solvesψ​\(s\)=t\\psi\(s\)=tforssviaoptimistix\.Bisectionwith tolerance10−1210^\{\-12\}\. Because generators are strictly decreasing on\[0,∞\)\[0,\\infty\), the bisection bracket always exists\. Gradients propagate through the root\-finding step via the implicit function theorem:optimistix\.ImplicitAdjointdifferentiates through the solution by solving∂ψ∂s\|s⋆​δ​s=δ​t\\frac\{\\partial\\psi\}\{\\partial s\}\\big\|\_\{s^\{\\star\}\}\\,\\delta s=\\delta t\. Our test suite verifies round\-trip accuracy—\|ψ​\(ψ−1​\(t\)\)−t\|<10−10\|\\psi\(\\psi^\{\-1\}\(t\)\)\-t\|<10^\{\-10\}—and gradient correctness against the analytic derivative\.

## Appendix FRosenblatt Sampling via Generator Derivatives

The same Bell\-polynomial machinery that delivers the density also samples from a nested Archimedean copula for free\. The Bell\-polynomial recursion of[Appendix˜A](https://arxiv.org/html/2605.23134#A1)supplies, as a byproduct, the high\-order generator\-composition derivatives that the textbookconditional\-distribution method—also called Rosenblatt sampling—requires\. The same sampler draws the synthetic data that validates the Frank and Gumbel parameter\-recovery results in[Section˜L\.1](https://arxiv.org/html/2605.23134#A12.SS1), and samples the InverseGaussian†, Nelsen9†, Nelsen12†, Nelsen13†, and Nelsen17†families used in the retinopathy experiment \([Section˜5\.6](https://arxiv.org/html/2605.23134#S5.SS6)\) with no per\-family implementation\.

##### The conditional\-distribution method\.

Sampling from an Archimedean copulaC​\(u1,…,ud\)=ψ​\(ψ−1​\(u1\)\+⋯\+ψ−1​\(ud\)\)C\(u\_\{1\},\\ldots,u\_\{d\}\)=\\psi\\bigl\(\\psi^\{\-1\}\(u\_\{1\}\)\+\\cdots\+\\psi^\{\-1\}\(u\_\{d\}\)\\bigr\)proceeds variable by variable\. Lettj=∑k=1jψ−1​\(uk\)t\_\{j\}=\\sum\_\{k=1\}^\{j\}\\psi^\{\-1\}\(u\_\{k\}\)denote the running sum of inverse\-generator transforms\. The conditional distribution ofuju\_\{j\}givenu1,…,uj−1u\_\{1\},\\ldots,u\_\{j\-1\}has the closed form

Cj∣1\.\.j−1​\(uj\)=ψ\(j−1\)​\(tj\)ψ\(j−1\)​\(tj−1\),C\_\{j\\mid 1\.\.j\-1\}\(u\_\{j\}\)\\;=\\;\\frac\{\\psi^\{\(j\-1\)\}\(t\_\{j\}\)\}\{\\psi^\{\(j\-1\)\}\(t\_\{j\-1\}\)\},\(26\)which depends only on the\(j−1\)\(j\{\-\}1\)\-th derivative of the generator evaluated at two arguments\. The sampler drawsVj∼Uniform​\(0,1\)V\_\{j\}\\sim\\mathrm\{Uniform\}\(0,1\)independently for each variable, then inverts[Equation˜26](https://arxiv.org/html/2605.23134#A6.E26)foruju\_\{j\}satisfyingCj∣1\.\.j−1​\(uj\)=VjC\_\{j\\mid 1\.\.j\-1\}\(u\_\{j\}\)=V\_\{j\}by bisection onuju\_\{j\}\. The nested case extends the same recursion through the tree by evaluating the appropriate edge\-composition derivative at each visited internal node\.[Algorithm˜4](https://arxiv.org/html/2605.23134#alg4)gives the flat\-copula sampler\. The nested case generalises by replacing the single\-generator derivativeψ\(j−1\)​\(⋅\)\\psi^\{\(j\-1\)\}\(\\cdot\)with the full Bell\-polynomial mixed partial of[Algorithm˜1](https://arxiv.org/html/2605.23134#alg1): each conditional CDFCj∣1\.\.j−1​\(uj\)C\_\{j\\mid 1\.\.j\-1\}\(u\_\{j\}\)is a ratio of two partial CDFs that differentiate variables1,…,j−11,\\ldots,j\{\-\}1and fix variablesj,…,dj,\\ldots,dto specific values, withuk=1u\_\{k\}=1marginalising the unsampled tail\. The censoring\-mask machinery that already drives likelihood evaluation supplies exactly this primitive, so the nested sampler reuses the density code path twice per variable \([Algorithm˜5](https://arxiv.org/html/2605.23134#alg5)\)\.

Algorithm 4Rosenblatt sampler for an Archimedean copula with generatorψ\\psiin dimensiondd\.- Require:Generatorψ\\psi, dimensiondd, RNG
- Ensure:Sampleu∈\(0,1\)du\\in\(0,1\)^\{d\}
- DrawV1,…,Vd∼Uniform​\(0,1\)V\_\{1\},\\ldots,V\_\{d\}\\sim\\mathrm\{Uniform\}\(0,1\)iid
- u1←V1u\_\{1\}\\leftarrow V\_\{1\};t1←ψ−1​\(u1\)t\_\{1\}\\leftarrow\\psi^\{\-1\}\(u\_\{1\}\)
- forj=2,…,dj=2,\\ldots,ddo
- DefineGj​\(u\):=ψ\(j−1\)​\(tj−1\+ψ−1​\(u\)\)/ψ\(j−1\)​\(tj−1\)G\_\{j\}\(u\)\\;:=\\;\\psi^\{\(j\-1\)\}\\\!\\bigl\(t\_\{j\-1\}\+\\psi^\{\-1\}\(u\)\\bigr\)\\,/\\,\\psi^\{\(j\-1\)\}\(t\_\{j\-1\}\) ⊳\\trianglerightψ\(j−1\)\\psi^\{\(j\-1\)\}via Taylor jet
- uj←u\_\{j\}\\leftarrowBisect\(Gj​\(u\)=VjG\_\{j\}\(u\)=V\_\{j\},u∈\(0,1\)u\\in\(0,1\)\)
- tj←tj−1\+ψ−1​\(uj\)t\_\{j\}\\leftarrow t\_\{j\-1\}\+\\psi^\{\-1\}\(u\_\{j\}\)
- return\(u1,…,ud\)\(u\_\{1\},\\ldots,u\_\{d\}\)

Algorithm 5Rosenblatt sampler for a nested Archimedean copula on treeTT\. Each inner call invokes the Bell\-polynomial mixed partial of[Algorithm˜1](https://arxiv.org/html/2605.23134#alg1)with a mask that differentiates only the already\-sampled variables\. An efficient implementation memoizes, at every internal nodevv, the Bell\-polynomial state of[Algorithm˜1](https://arxiv.org/html/2605.23134#alg1): the running sumtvt\_\{v\}, the countdvd\_\{v\}of leaves already sampled invv’s subtree, and the list of completed children’sα\(c\)\\alpha^\{\(c\)\}vectors; each per\-leaf step then traces only from the leaf’s parent up to the root\.- Require:Nested copula treeTT, leaf orderu1,…,udu\_\{1\},\\ldots,u\_\{d\}, RNG
- Ensure:Sampleu∈\(0,1\)du\\in\(0,1\)^\{d\}
- DrawV1,…,Vd∼Uniform​\(0,1\)V\_\{1\},\\ldots,V\_\{d\}\\sim\\mathrm\{Uniform\}\(0,1\)iid
- forj=1,…,dj=1,\\ldots,ddo
- δj←\(1,…,1⏟j−1,0,…,0\)\\delta\_\{j\}\\leftarrow\(\\underbrace\{1,\\ldots,1\}\_\{j\-1\},\\;0,\\ldots,0\) ⊳\\trianglerightdifferentiateu1,…,uj−1u\_\{1\},\\ldots,u\_\{j\-1\}only
- DefineDj​\(u\)←D\_\{j\}\(u\)\\leftarrowBellPartial\(T,\(u1,…,uj−1,u,1,…,1\),δjT,\\;\(u\_\{1\},\\ldots,u\_\{j\-1\},\\,u,\\,1,\\ldots,1\),\\;\\delta\_\{j\}\) ⊳\\trianglerightpartial CDF
- uj←u\_\{j\}\\leftarrowBisect\(Dj​\(u\)/Dj​\(1\)=VjD\_\{j\}\(u\)\\,/\\,D\_\{j\}\(1\)=V\_\{j\},u∈\(0,1\)u\\in\(0,1\)\)
- return\(u1,…,ud\)\(u\_\{1\},\\ldots,u\_\{d\}\)

##### Connection to prior work\.

Hofert \[[2008](https://arxiv.org/html/2605.23134#bib.bib15)\]identified the high\-order generator\-composition derivatives in[Equation˜26](https://arxiv.org/html/2605.23134#A6.E26)as the principal obstacle to applying this method beyond small dimensions and developed Marshall–Olkin frailty sampling as the practical alternative\. That alternative requires a known Laplace transform per family, which blocks cross\-family nesting and rules out any generator without a closed\-form mixing distribution—for example, the Gen\-AC generator ofNg et al\. \[[2021](https://arxiv.org/html/2605.23134#bib.bib34)\]adopted byLiu et al\. \[[2025b](https://arxiv.org/html/2605.23134#bib.bib29)\]\. The Bell\-polynomial recursion of[Appendix˜A](https://arxiv.org/html/2605.23134#A1)computes exactly these derivatives in closed form for any user\-supplied generator with the same complexity profile as the density evaluation, so sampling and likelihood share a single code path\.

## Appendix GValidity of Nested Constructions

A nested Archimedean copula defines a valid distribution only when every parent–child generator pair satisfies a sufficient nesting condition\.McNeil \[[2008](https://arxiv.org/html/2605.23134#bib.bib30)\]required the compositionh=ψouter−1∘ψinnerh=\\psi\_\{\\mathrm\{outer\}\}^\{\-1\}\\circ\\psi\_\{\\mathrm\{inner\}\}to be aBernstein function\(i\.e\., its derivative completely monotone\):\(−1\)k−1​h\(k\)​\(t\)≥0\(\-1\)^\{k\-1\}h^\{\(k\)\}\(t\)\\geq 0for allk≥1k\\geq 1andt≥0t\\geq 0\.McNeil and Nešlehová \[[2009](https://arxiv.org/html/2605.23134#bib.bib31)\]subsequently characteriseddd\-dimensional Archimedean copulas through the Williamsondd\-transform, weakening the complete\-monotonicity requirement on the generator itself todd\-monotonicity and setting the foundation for the results below\.Rezapour \[[2015](https://arxiv.org/html/2605.23134#bib.bib38)\]showed that a weaker condition suffices:hhneed only bedcd\_\{c\}\-monotone, that is,\(−1\)m−1​h\(m\)​\(t\)≥0\(\-1\)^\{m\-1\}h^\{\(m\)\}\(t\)\\geq 0form=1,…,dcm=1,\\ldots,d\_\{c\}and allt≥0t\\geq 0, wheredcd\_\{c\}counts uncensored leaves under the child node\. This finite condition matches the information our algorithm already computes: theJetcall at each edge returns exactlydcd\_\{c\}Taylor coefficients ofhh\.

For the classical families, closed\-form sufficient conditions exist: Clayton and Gumbel both requireθouter≤θinner\\theta\_\{\\mathrm\{outer\}\}\\leq\\theta\_\{\\mathrm\{inner\}\}, Frank requires the same ordering, and Joe requires1≤θouter≤θinner1\\leq\\theta\_\{\\mathrm\{outer\}\}\\leq\\theta\_\{\\mathrm\{inner\}\}\. These orderings guarantee complete monotonicity and thereforedcd\_\{c\}\-monotonicity\.

Our framework does*not*enforce these constraints automatically\. For same\-family pairs, a delta parameterisationθinner=θouter\+softplus⁡\(δ\)\\theta\_\{\\mathrm\{inner\}\}=\\theta\_\{\\mathrm\{outer\}\}\+\\operatorname\{softplus\}\(\\delta\)enforces the ordering by construction; all same\-family experiments use this approach \([Table˜12](https://arxiv.org/html/2605.23134#A7.T12)\)\.

Table 12:Nesting conditions for generator pairs in our experiments\.Same\-family rows enforceθouter≤θinner\\theta\_\{\\mathrm\{outer\}\}\\leq\\theta\_\{\\mathrm\{inner\}\}via the delta parameterisation\. The AMH/Clayton row uses Hofert\[[2010](https://arxiv.org/html/2605.23134#bib.bib16)\]’s closed\-form constraintθClayton≥1\\theta\_\{\\mathrm\{Clayton\}\}\\geq 1, enforced by shifted\-softplusθClayton=1\+softplus​\(δ\)\\theta\_\{\\mathrm\{Clayton\}\}=1\+\\mathrm\{softplus\}\(\\delta\); this is the single documented closed\-form cross\-family Archimedean pair with a McNeil\-2008 Bernstein\-function nesting proof\. InverseGaussian, Nelsen13, and Nelsen17 appear only in the retinopathy experiment, which uses single\-layer bivariate copulas; nesting conditions do not apply\.

##### Cross\-family and custom generators\.

The closed\-form constraints above cover the families used in our experiments, but a user nesting two different families—or supplying a custom generator—needs a more general route\. When no closed\-form parameter constraints exist, enforcing the nesting condition during maximum likelihood estimation becomes asemi\-infinite program\(SIP\)\. Theded\_\{e\}\-monotone condition must hold at everyt≥0t\\geq 0, producing infinitely many inequality constraints over the finite\-dimensional parameter vector𝜽∈ℝp\\boldsymbol\{\\theta\}\\in\\mathbb\{R\}^\{p\}stacking all per\-edge generator parameters\. Letℰ\\mathcal\{E\}denote the set of parent–child edges in the copula tree,ded\_\{e\}the number of uncensored leaves under the child of edgeee, andhe​\(t;𝜽\)=ψouter−1​\(ψinner​\(t\)\)h\_\{e\}\(t;\\,\\boldsymbol\{\\theta\}\)=\\psi\_\{\\mathrm\{outer\}\}^\{\-1\}\\\!\\bigl\(\\psi\_\{\\mathrm\{inner\}\}\(t\)\\bigr\)the edge composition, withhe≡hv​ch\_\{e\}\\equiv h\_\{vc\}whene=\(v,c\)e=\(v,c\)\. The constrained MLE is

max𝜽\\displaystyle\\max\_\{\\boldsymbol\{\\theta\}\}ℓ​\(𝜽;𝐗\)\\displaystyle\\ell\(\\boldsymbol\{\\theta\};\\,\\mathbf\{X\}\)\(27\)s\.t\.\(−1\)m−1​he\(m\)​\(t;𝜽\)≥0∀e∈ℰ,m=1,…,de,t≥0\.\\displaystyle\(\-1\)^\{m\-1\}\\,h\_\{e\}^\{\(m\)\}\(t;\\,\\boldsymbol\{\\theta\}\)\\;\\geq 0\\quad\\forall\\,e\\in\\mathcal\{E\},\\;\\;m=1,\\ldots,d\_\{e\},\\;\\;t\\geq 0\.

##### Grid discretisation\.

A standard approach to a SIP replaces the continuum of constraints with a finite grid\{t1,…,tG\}⊂\[0,tmax\]\\\{t\_\{1\},\\ldots,t\_\{G\}\\\}\\subset\[0,t\_\{\\max\}\]ofGGuser\-chosen points:

max𝜽\\displaystyle\\max\_\{\\boldsymbol\{\\theta\}\}ℓ​\(𝜽;𝐗\)\\displaystyle\\ell\(\\boldsymbol\{\\theta\};\\,\\mathbf\{X\}\)\(28\)s\.t\.\(−1\)m−1​he\(m\)​\(tg;𝜽\)≥0∀e∈ℰ,m=1,…,de,g=1,…,G\.\\displaystyle\(\-1\)^\{m\-1\}\\,h\_\{e\}^\{\(m\)\}\(t\_\{g\};\\,\\boldsymbol\{\\theta\}\)\\;\\geq 0\\quad\\forall\\,e\\in\\mathcal\{E\},\\;\\;m=1,\\ldots,d\_\{e\},\\;\\;g=1,\\ldots,G\.Checking each grid point at one edge requires a singleJetcall that returns theded\_\{e\}Taylor coefficients ofheh\_\{e\}; the signs of those coefficients encode theded\_\{e\}\-monotone condition\. Concretely, theJetcall returns the Taylor coefficient vectorpm=he\(m\+1\)​\(t\)/\(m\+1\)\!p\_\{m\}=h\_\{e\}^\{\(m\+1\)\}\(t\)/\(m\{\+\}1\)\!form=0,…,de−1m=0,\\ldots,d\_\{e\}\{\-\}1\. Because factorial scaling preserves sign, theded\_\{e\}\-monotone condition\(−1\)m−1​he\(m\)​\(t\)≥0\(\-1\)^\{m\-1\}h\_\{e\}^\{\(m\)\}\(t\)\\geq 0form=1,…,dem=1,\\ldots,d\_\{e\}reduces to checking that\(−1\)m​pm≥0\(\-1\)^\{m\}p\_\{m\}\\geq 0for eachmm—an alternating\-sign test on the Taylor coefficients\. Each density evaluation already performs one such call per edge att=tet=t\_\{e\}, soNNtraining observations contributeNNgrid points per edge at no extra cost; a denser grid requiresG−NG\-NadditionalJetcalls\.

##### Computing the penalty in principle\.

Theacopulaframework can compute the nesting penalty as a byproduct of the density evaluation\. At each edgeee, theJetcall that produces the Taylor coefficientsp0,p1,…,pdep\_\{0\},p\_\{1\},\\ldots,p\_\{d\_\{e\}\}for the Bell transform also provides the sign information for the diagnostic\. The penalty takes the form

P​\(𝜽;𝐱\)=∑e∈ℰ∑m=0de−1\[min⁡\(0,\(−1\)m​pm\)\]2,P\(\\boldsymbol\{\\theta\};\\,\\mathbf\{x\}\)=\\sum\_\{e\\in\\mathcal\{E\}\}\\;\\sum\_\{m=0\}^\{d\_\{e\}\-1\}\\bigl\[\\min\\\!\\bigl\(0,\\;\(\-1\)^\{m\}\\,p\_\{m\}\\bigr\)\\bigr\]^\{2\},\(29\)evaluated att=te​\(𝐱\)t=t\_\{e\}\(\\mathbf\{x\}\), the sector sum for the current observation𝐱\\mathbf\{x\}\. WhenP=0P=0, all Taylor coefficients satisfy the alternating\-sign condition at every evaluation point encountered during training\. A penalised optimiser would minimise−ℓ​\(𝜽\)\+ρ​P​\(𝜽\)\-\\ell\(\\boldsymbol\{\\theta\}\)\+\\rho\\,P\(\\boldsymbol\{\\theta\}\)with penalty weightρ≥0\\rho\\geq 0, but all experiments in this paper avoid this route: every nested fit uses one of the three valid\-by\-construction parameterisations summarised in[Table˜12](https://arxiv.org/html/2605.23134#A7.T12)—same\-family delta, AMH/Clayton shifted\-softplus, or the Gen\-AC/ACNet Bernstein composition \([Section˜G\.1](https://arxiv.org/html/2605.23134#A7.SS1)\)—which makes the penalty identically zero by construction\. We describe the penalty approach here as a conceptual extension for users who might want to experiment with cross\-family or custom generators that lack closed\-form constraints; we leave empirical validation of the penalised route to future work\.

##### Limitations\.

The grid relaxation[Equation˜28](https://arxiv.org/html/2605.23134#A7.E28)is necessary but not sufficient for the full SIP[Equation˜27](https://arxiv.org/html/2605.23134#A7.E27): violations can hide between grid points\. Certifying theded\_\{e\}\-monotone condition on all of\[0,∞\)\[0,\\infty\)for arbitrary generators remains an open problem\.

### G\.1Validity of the nested Gen\-AC/ACNet construction

The nested construction used by the Gen\-AC and ACNet rows of[Section˜5\.4](https://arxiv.org/html/2605.23134#S5.SS4)and by the HACSurv port of[Appendix˜S](https://arxiv.org/html/2605.23134#A19)pairs a shared outer generatorψout​\(s\)\\psi\_\{\\mathrm\{out\}\}\(s\)—a mixture of exponentialsLout−1​∑l=1Loutexp⁡\(−Ml​s\)L\_\{\\mathrm\{out\}\}^\{\-1\}\\sum\_\{l=1\}^\{L\_\{\\mathrm\{out\}\}\}\\exp\(\-M\_\{l\}s\)for Gen\-AC \(atLoutL\_\{\\mathrm\{out\}\}MC nodes\) and the softmax\-convex cascade of[Appendix˜I](https://arxiv.org/html/2605.23134#A9)for ACNet—with per\-panel innerψin,k​\(t\)=ψout​\(bk​\(t\)\)\\psi\_\{\\mathrm\{in\},k\}\(t\)=\\psi\_\{\\mathrm\{out\}\}\(b\_\{k\}\(t\)\)and inner Bernstein functionbk​\(t\)=eμk​t\+eβk​\(1−Lin−1​∑l=1Line−Nk,l​t\)b\_\{k\}\(t\)=e^\{\\mu\_\{k\}\}\\,t\+e^\{\\beta\_\{k\}\}\\bigl\(1\-L\_\{\\mathrm\{in\}\}^\{\-1\}\\sum\_\{l=1\}^\{L\_\{\\mathrm\{in\}\}\}e^\{\-N\_\{k,l\}t\}\\bigr\), whereLout,Lin∈ℕL\_\{\\mathrm\{out\}\},L\_\{\\mathrm\{in\}\}\\in\\mathbb\{N\}are the outer and inner mixture sizes,Ml,Nk,l≥0M\_\{l\},N\_\{k,l\}\\geq 0are exponential rates,kkindexes the per\-panel inner generator, andμk,βk∈ℝ\\mu\_\{k\},\\beta\_\{k\}\\in\\mathbb\{R\}are scalar log\-shift parameters\. The construction is valid for any\(Nk,μk,βk\)∈ℝ≥0Lin×ℝ2\(N\_\{k\},\\mu\_\{k\},\\beta\_\{k\}\)\\in\\mathbb\{R\}\_\{\\geq 0\}^\{L\_\{\\mathrm\{in\}\}\}\\times\\mathbb\{R\}^\{2\}, and the optimiser need not enforce any constraint\. The argument runs in three steps\. First,bkb\_\{k\}is a Bernstein function:bk​\(0\)=0b\_\{k\}\(0\)=0,bk\>0b\_\{k\}\>0, andbk′b\_\{k\}^\{\\prime\}is completely monotone \(CM\) as a sum of CM components\. Second, the composition of a CM function with a Bernstein function is itself CM\[McNeil,[2008](https://arxiv.org/html/2605.23134#bib.bib30), Rezapour,[2015](https://arxiv.org/html/2605.23134#bib.bib38)\], soψin,k\\psi\_\{\\mathrm\{in\},k\}is a valid generator\. Third,ψout−1∘ψin,k=bk\\psi\_\{\\mathrm\{out\}\}^\{\-1\}\\circ\\psi\_\{\\mathrm\{in\},k\}=b\_\{k\}is automatically Bernstein, so the parent–child pair satisfies the nesting condition\.

## Appendix HGen\-AC MLP\-pushforward generator with Fourier features

The Gen\-AC rows of[Table˜4](https://arxiv.org/html/2605.23134#S5.T4)implement the MLP\-pushforward generator ofNg et al\. \[[2021](https://arxiv.org/html/2605.23134#bib.bib34)\]with frailtyV=exp⁡\(MLP​\(U\)\)V=\\exp\(\\mathrm\{MLP\}\(U\)\),U∼Uniform​\(0,1\)U\\sim\\mathrm\{Uniform\}\(0,1\), so the Laplace transform isψ​\(t\)=𝔼U​\[exp⁡\(−V​t\)\]\\psi\(t\)=\\mathbb\{E\}\_\{U\}\[\\exp\(\-Vt\)\]\. HACSurv\[Liu et al\.,[2025b](https://arxiv.org/html/2605.23134#bib.bib29)\]uses this construction directly with a three\-layer MLP of hidden width1010and leaky\-ReLU activations; the row labelled Gen\-AC1\.1kreproduces that architecture with one outer MLP and a per\-panel inner MLP, composed with the Bernstein drift of[Section˜G\.1](https://arxiv.org/html/2605.23134#A7.SS1), for7⋅2\+8⋅141=1,1427\\cdot 2\+8\\cdot 141=1\{,\}142parameters\. We evaluate𝔼U\\mathbb\{E\}\_\{U\}atNz=32N\_\{z\}\{=\}32fixed equispaced nodes on\(0,1\)\(0,1\); equispaced nodes instead of per\-step Monte\-Carlo resamples remove the Jensen bias and shift held\-out NLL by under0\.020\.02nats\.

The raw\-uuinput of HACSurv’s MLP plateaus at held\-out NLL4\.9144\.914—barely better than Clayton’s4\.9494\.949—because a small MLP of Lipschitz constant≈1\\approx 1cannot place neighbouringlog⁡Mi\\log M\_\{i\}at widely separated rates: the 32 atoms collapse onto a narrow band\. We lift that prior by replacing the scalar input with Fourier features\[Tancik et al\.,[2020](https://arxiv.org/html/2605.23134#bib.bib42)\]atKKfrequencies:

u⟼\[cos⁡\(2​π​k​u\),sin⁡\(2​π​k​u\)\]k=1K∈ℝ2​K,u\\;\\longmapsto\\;\\bigl\[\\cos\(2\\pi ku\),\\,\\sin\(2\\pi ku\)\\bigr\]\_\{k=1\}^\{K\}\\;\\in\\;\\mathbb\{R\}^\{2K\},\(30\)passed into an otherwise identical three\-layer MLP of hidden widthhh\. The Fourier basis lets the network output sharply differentlog⁡Mi\\log M\_\{i\}at neighbouringuiu\_\{i\}without requiring a large Lipschitz constant, recovering the flexibility that direct trainablelog⁡M∈ℝ32\\log M\\in\\mathbb\{R\}^\{32\}would have\.[Table˜4](https://arxiv.org/html/2605.23134#S5.T4)reports two configurations: a parameter\-matched variant Gen\-AC\-Fourier260at\(K=2,h=3\)\(K\{=\}2,h\{=\}3\), and the full model Gen\-AC\-Fourier3\.6kat\(K=16,h=10\)\(K\{=\}16,h\{=\}10\)\. At matched capacity Gen\-AC remains less efficient than ACNet in this cohort—at260260parameters Gen\-AC reaches4\.8284\.828versus ACNet’s4\.3104\.310at7070parameters—and even the full Fourier model trails the best ACNet by∼0\.17\{\\sim\}0\.17nats once its parameter budget is released\.

## Appendix IDeep ACNet Generator

The ACNet rows of[Table˜4](https://arxiv.org/html/2605.23134#S5.T4)use the depth\-LLsoftmax\-convex\-combination generator ofLing et al\. \[[2020](https://arxiv.org/html/2605.23134#bib.bib27)\]as the outerψout\\psi\_\{\\mathrm\{out\}\}and the Bernstein\-composition innerψin,k​\(t\)=ψout​\(bk​\(t\)\)\\psi\_\{\\mathrm\{in\},k\}\(t\)=\\psi\_\{\\mathrm\{out\}\}\(b\_\{k\}\(t\)\)from[Section˜G\.1](https://arxiv.org/html/2605.23134#A7.SS1)\. Given widths\(w0,…,wL−1\)\(w\_\{0\},\\ldots,w\_\{L\-1\}\)and a queryt∈\[0,∞\)t\\in\[0,\\infty\), the generator is the terminal state of anLL\-layer cascade of positive exponential shifts interleaved with simplex\-constrained convex combinations, defined recursively as

Fi\(0\)\\displaystyle F^\{\(0\)\}\_\{i\}=e−s0,i​t,\\displaystyle\\;=\\;e^\{\-s\_\{0,i\}t\},\(31\)Fi\(ℓ\)\\displaystyle F^\{\(\\ell\)\}\_\{i\}=e−sℓ,i​t⋅\(Aℓ−1​F\(ℓ−1\)\)i​for​ℓ=1,…,L−1,\\displaystyle\\;=\\;e^\{\-s\_\{\\ell,i\}t\}\\cdot\\bigl\(A\_\{\\ell\-1\}F^\{\(\\ell\-1\)\}\\bigr\)\_\{i\}\\text\{ for \}\\ell=1,\\ldots,L\-1,ψoutACNet​\(t\)\\displaystyle\\psi\_\{\\mathrm\{out\}\}^\{\\mathrm\{ACNet\}\}\(t\)=AL−1​F\(L−1\),\\displaystyle\\;=\\;A\_\{L\-1\}F^\{\(L\-1\)\},wheresℓ=es~ℓ∈ℝwℓs\_\{\\ell\}=e^\{\\tilde\{s\}\_\{\\ell\}\}\\in\\mathbb\{R\}^\{w\_\{\\ell\}\}are per\-layer positive exponential rates ands~ℓ\\tilde\{s\}\_\{\\ell\}their unconstrained log parameters,Aℓ∈ℝwℓ\+1×wℓA\_\{\\ell\}\\in\\mathbb\{R\}^\{w\_\{\\ell\+1\}\\times w\_\{\\ell\}\}forℓ<L−1\\ell<L\-1andAL−1∈ℝ1×wL−1A\_\{L\-1\}\\in\\mathbb\{R\}^\{1\\times w\_\{L\-1\}\}are the row\-wise softmaxes of free matricesA~ℓ\\tilde\{A\}\_\{\\ell\}so that every row ofAℓA\_\{\\ell\}lies on the simplex\. The first layer’s convex combination acts on the scalar11, so the recursion absorbs it into the initialisationF\(0\)F^\{\(0\)\}\. Each per\-layer factore−sℓ,i​te^\{\-s\_\{\\ell,i\}t\}is a completely monotone function oftt, every row ofAℓA\_\{\\ell\}is a convex combination, and closure of complete monotonicity under convex combinations and under products of completely monotone functions \(applied at every step of the recursion\) yields a scalar\-valued completely monotone function, soψoutACNet\\psi\_\{\\mathrm\{out\}\}^\{\\mathrm\{ACNet\}\}is a valid Archimedean generator\. The parameterisation imposes no constraint, and the optimiser updatess~ℓ,A~ℓ\\tilde\{s\}\_\{\\ell\},\\tilde\{A\}\_\{\\ell\}freely\.

The inner generator follows the same Bernstein\-composition recipe as[Section˜G\.1](https://arxiv.org/html/2605.23134#A7.SS1),ψin,k​\(t\)=ψoutACNet​\(bk​\(t\)\)\\psi\_\{\\mathrm\{in\},k\}\(t\)=\\psi\_\{\\mathrm\{out\}\}^\{\\mathrm\{ACNet\}\}\(b\_\{k\}\(t\)\)withbk​\(t\)=eμk​t\+eβk​\(1−Lin−1​∑le−Nk,l​t\)b\_\{k\}\(t\)=e^\{\\mu\_\{k\}\}\\,t\+e^\{\\beta\_\{k\}\}\(1\-L\_\{\\mathrm\{in\}\}^\{\-1\}\\sum\_\{l\}e^\{\-N\_\{k,l\}t\}\); the same Schoenberg\-theorem argument showsψin,k\\psi\_\{\\mathrm\{in\},k\}is completely monotone and the parent–child compositionψout−1∘ψin,k=bk\\psi\_\{\\mathrm\{out\}\}^\{\-1\}\\circ\\psi\_\{\\mathrm\{in\},k\}=b\_\{k\}is Bernstein\. The compact\-parameter count in[Table˜4](https://arxiv.org/html/2605.23134#S5.T4)’s\|θ\|\|\\theta\|column sums the per\-layer unconstraineds~ℓ∈ℝwℓ\\tilde\{s\}\_\{\\ell\}\\in\\mathbb\{R\}^\{w\_\{\\ell\}\}andA~ℓ∈ℝwℓ\+1×wℓ\\tilde\{A\}\_\{\\ell\}\\in\\mathbb\{R\}^\{w\_\{\\ell\+1\}\\times w\_\{\\ell\}\}of the shared outer, the scalar log\-shift pair\(μk,βk\)\(\\mu\_\{k\},\\beta\_\{k\}\)per panel, and the per\-panel inner ratesNk∈ℝLinN\_\{k\}\\in\\mathbb\{R\}^\{L\_\{\\mathrm\{in\}\}\}\.

## Appendix JDetailed R Comparison

Our Bell polynomial log\-density agrees with R’snacLL\[Hofert and Pham,[2013](https://arxiv.org/html/2605.23134#bib.bib17)\]to machine precision on every configuration where R completes, establishing numerical correctness of the framework’s forward evaluation\. We compare the two implementations on 12 nested Clayton–Clayton configurations spanningd=4d=4to6060total dimensions with identical random input matrices,θouter=2\\theta\_\{\\mathrm\{outer\}\}=2,θinner=5\\theta\_\{\\mathrm\{inner\}\}=5, andn=5n=5observations per configuration\.[Table˜13](https://arxiv.org/html/2605.23134#A10.T13)reports the results\.

Table 13:Density correctness ofacopulaversus R’snacLLon nested Clayton–Clayton copulas\. Here,ssdenotes the number of sectors,KKthe leaves per sector, andd=s​Kd=sKthe total dimension\. The relative error is\|ℓours−ℓR\|/\|ℓR\|\|\\ell\_\{\\mathrm\{ours\}\}\-\\ell\_\{R\}\|/\|\\ell\_\{R\}\|\. Arbitrary\-precision verification attributes thed=60d\{=\}60discrepancy to R’snacLL, not our framework\.Eleven of 12 configurations agree to relative error below10−1410^\{\-14\}, matching float64 roundoff\. We attribute the singled=60d\{=\}60discrepancy to precision loss in R’snacLL, not our framework: Mathematica’s symbolic differentiation of the copula CDF viaD\[C, u1, …, ud\]provides an algorithm\-agnostic ground truth that agrees with our framework to machine epsilon at everyKKfrom 2 to 30, while R’snacLLdiverges atK≥31K\\geq 31\(see also[Appendix˜B](https://arxiv.org/html/2605.23134#A2)for the cross\-family Mathematica check\)\. The scriptsmathematica\_ground\_truth\.wlsandverify\_against\_mathematica\.pyreproduce this comparison\.

##### Controlled runtime comparison\.

Under matched single\-thread settings,acopularuns orders of magnitude faster than R’snacLLat moderate dimension and completes families thatnacLLcannot evaluate at all\.[Table˜14](https://arxiv.org/html/2605.23134#A10.T14)reports wall\-clock time for the log\-likelihood ofn=100n\{=\}100observations under an apples\-to\-apples single\-thread configuration that pins thread counts and the BLAS backend across both runtimes\. Both R and Python run withOMP\_NUM\_THREADS=OPENBLAS\_NUM\_THREADS=MKL\_NUM\_THREADS=1, read the same float64 input matrices from disk, and time identical nested same\-family configurations withθouter=2\\theta\_\{\\mathrm\{outer\}\}\{=\}2andθinner=5\\theta\_\{\\mathrm\{inner\}\}\{=\}5viamicrobenchmarkon the R side andtime\.perf\_counterafter a JIT warmup on the Python side\. We compare three families that R’snacLLpurports to support: Clayton, Frank, and Gumbel\.

Table 14:Controlled single\-thread runtime ofacopula, JIT\-compiled on CPU, versus RnacLLon nested same\-family copulas withn=100n\{=\}100observations,θouter=2\\theta\_\{\\mathrm\{outer\}\}\{=\}2, andθinner=5\\theta\_\{\\mathrm\{inner\}\}\{=\}5\. The “Compile” column reports the one\-time JIT cost; N/A indicates that R’snacLLaborts on the configuration\.Three findings emerge from[Table˜14](https://arxiv.org/html/2605.23134#A10.T14)\. First, our framework matches R’s log\-likelihood values to relative error below1\.7×10−151\.7\\times 10^\{\-15\}in every cell where R completes, ruling out a precision artefact in the timing gap\. Second, the speedup grows from∼4×\{\\sim\}4\\timesatd=4d\{=\}4to over700×700\\timesatd=30d\{=\}30for Gumbel: R’snacLLscales exponentially in the number of sectors via partition enumeration, whereas our polynomial powering scales asO​\(d3\)O\(d^\{3\}\), so the gap widens monotonically withdd\. Third, R’snacLLaborts on every Frank configuration with the errora\.coeff: Frank’s family is currently not supported; its documented support for nesting across the five classical families therefore fails for Frank, AMH, and Joe, becausenacLLhard\-codes the partial\-Bell coefficients only for Clayton and Gumbel\.[Table˜15](https://arxiv.org/html/2605.23134#A10.T15)aggregates the structural differences that drive these gaps\.

Table 15:Comparison ofacopulawith R’scopula/nacopulapackages\.Two practical caveats sharpen the comparison\. R’s closed\-form evaluations are faster per call than our jet\-basedkk\-step scan, but the Stirling and Eulerian numbers they rely on overflow 64\-bit floats at moderate orders\. Our JIT\-compiled pipeline pays a one\-time compilation cost of∼50\{\\sim\}50ms and thereafter evaluates in<0\.25<0\.25ms per observation\. TheCopulabase class also exposes an optionallog\_generator\_kth\_derivative\(t, k\)hook for users who have a closed\-formψ\(k\)\\psi^\{\(k\)\}, recovering R’s per\-call speed while retaining full framework compatibility\.

##### Polynomial scaling inss\.

Our Cauchy\-product algorithm replaces R’s exponential partition enumeration with polynomial multiplication, so cost in the number of sectorsssgrows polynomially rather than exponentially\. To quantify the difference, we time the JIT\-compiled log\-likelihood for nested Clayton copulas while varyingssfrom 2 to 50 at sector sizeK=2K=2, for whichd=2​sd=2s\. Execution time grows from 0\.4 ms ats=2s\{=\}2to 72 ms ats=50s\{=\}50, matching approximatelyO​\(s1\.7\)O\(s^\{1\.7\}\)scaling\. Repeating withK=5K=5up tos=30s\{=\}30, withd=5​sd=5s, yields a similar polynomial exponent\. By contrast, R’s partition count exceeds101510^\{15\}ats=50s=50, rendering exact enumeration infeasible\. Beyond raw scaling, the choice between exact and Monte Carloψ\(k\)\\psi^\{\(k\)\}evaluation involves further design tradeoffs that we discuss next\.

##### Extended Gen\-AC comparison\.

Our Taylor\-mode computation ofψ\(k\)\\psi^\{\(k\)\}and the Monte Carlo estimation in Gen\-AC\[Ng et al\.,[2021](https://arxiv.org/html/2605.23134#bib.bib34)\]occupy complementary points in the exactness–scalability tradeoff\. Taylor\-mode AD computes eachψ\(k\)\\psi^\{\(k\)\}exactly inO​\(k\)O\(k\)time per coefficient, yielding machine\-precision likelihoods that support gradient\-based parameter estimation, nested copula structures, and per\-variable censoring without approximation error\. Gen\-AC, in contrast, sidesteps generator differentiation by representingψ\(k\)\\psi^\{\(k\)\}as a Monte Carlo average over learned latent variables; this representation scales to very highkkwith bounded per\-sample cost but introduces sampling variance and restricts the generator to the learned Laplace class\. For the moderate dimensions typical of actuarial and survival\-analysis applications,d≤100d\\leq 100, the exact approach dominates: it produces unbiased gradients and handles arbitrary user\-defined generators\. Gen\-AC suits regimes where only samples—not likelihoods—matter, or whereddis so large that evenO​\(k\)O\(k\)exact derivative cost grows prohibitive\.

##### Bivariate censored comparison with CopulaCenR\.

acopulareproduces CopulaCenR’s parameter estimates on right\-censored bivariate data and succeeds where CopulaCenR’s finite\-difference optimiser fails\. We compare against CopulaCenR\[Sun and Ding,[2020](https://arxiv.org/html/2605.23134#bib.bib41)\], an R package that fits bivariate copulas to right\-censored data via joint MLE with parametric marginals\. On the retinopathy dataset, with197197patients at60\.7%60\.7\\%censoring, CopulaCenR’s Clayton fit yieldsη^=0\.895\\hat\{\\eta\}=0\.895versus our two\-stage Inference Functions for Margins \(IFM\)\[Joe and Xu,[1996](https://arxiv.org/html/2605.23134#bib.bib21), Genest et al\.,[1995](https://arxiv.org/html/2605.23134#bib.bib12)\]estimateθ^=0\.897\\hat\{\\theta\}=0\.897, a0\.3%0\.3\\%difference; its Gumbel fit yieldsη^=1\.256\\hat\{\\eta\}=1\.256versus ourθ^=1\.249\\hat\{\\theta\}=1\.249, a0\.6%0\.6\\%difference\. The CopulaCenR estimateη^\\hat\{\\eta\}uses the same Clayton/Gumbel parameterisation as ourθ^\\hat\{\\theta\}, so the two are directly comparable\. CopulaCenR’s Frank optimiser fails on this dataset with a finite\-difference overflow, while our AD\-based gradient computation converges toθ^=2\.255\\hat\{\\theta\}=2\.255\. The close parameter agreement on two families confirms that the copula likelihoods agree in the bivariate censored setting, despite two methodological differences: joint MLE versus IFM, and shared versus per\-margin Weibull marginals\.

## Appendix KAD Gradient Verification

JAX’s reverse\-mode AD produces parameter gradients that match central finite differences to relative error below10−810^\{\-8\}across families, dimensions, and censoring, validating the gradient path that drives parameter estimation\. We compare the AD gradientgAD=∇θℓg\_\{\\mathrm\{AD\}\}=\\nabla\_\{\\theta\}\\ell\(whereℓ\\ellis the total log\-likelihood across the configuration’s training sample\) thatjax\.gradreturns against the central finite\-difference approximationgFD,i=\[ℓ​\(θ\+ε​ei\)−ℓ​\(θ−ε​ei\)\]/\(2​ε\)g\_\{\\mathrm\{FD\},i\}=\\bigl\[\\ell\(\\theta\+\\varepsilon e\_\{i\}\)\-\\ell\(\\theta\-\\varepsilon e\_\{i\}\)\\bigr\]/\(2\\varepsilon\), withε=10−5\\varepsilon=10^\{\-5\}andei∈ℝpe\_\{i\}\\in\\mathbb\{R\}^\{p\}theii\-th standard basis vector for parameter dimensionpp\(theParamscolumn in[Table˜16](https://arxiv.org/html/2605.23134#A11.T16)\), across four nested copula configurations that exercise different families, dimensions, and censoring\.[Table˜16](https://arxiv.org/html/2605.23134#A11.T16)reports the maximum relative errormaxi⁡\|gAD,i−gFD,i\|/max⁡\(\|gAD,i\|,10−10\)\\max\_\{i\}\|g\_\{\\mathrm\{AD\},i\}\-g\_\{\\mathrm\{FD\},i\}\|/\\max\(\|g\_\{\\mathrm\{AD\},i\}\|,10^\{\-10\}\)\.

Table 16:AD gradient vs\. central finite differences \(ε=10−5\\varepsilon\{=\}10^\{\-5\}\)\. All configurations use 2\-level nesting with the stated sector structure\(M×K\)\(M\\times K\)denotingMMsectors ofKKleaves each\.All four configurations achieve relative error below10−810^\{\-8\}, confirming that JAX’s reverse\-mode AD produces correct gradients through the full Bell\-polynomial pipeline—including the log\-domain rescaling, signed logsumexp, and polynomial powering steps that[Appendix˜B](https://arxiv.org/html/2605.23134#A2)describes\. The censored configuration, with 55% of variables censored, shows no degradation relative to the fully observed cases, which validates the gradient path through the dynamic censoring mechanism\.

## Appendix LHigh\-Dimensional Parameter Recovery

The MLE recovers the true parameters with bias under2\.7%2\.7\\%and sampling variability that the observed\-Fisher SE calibrates to within2\.9%2\.9\\%on average, acrossN∈\{200,500,1000\}N\\in\\\{200,500,1000\\\}and 500 replications per cell\.[Table˜3](https://arxiv.org/html/2605.23134#S5.T3)in the main text reports replicated recovery atN=500N\{=\}500;[Table˜17](https://arxiv.org/html/2605.23134#A12.T17)extends that result to the fullNNgrid and reports bias, empirical SE, observed\-Fisher SE, RMSE, and the L\-BFGS\-B convergence rate per cell\. The empirical SE shrinks at the expected1/N1/\\sqrt\{N\}rate asNNgrows: the EmpSE ratio betweenN=200N\{=\}200andN=1000N\{=\}1000ranges from2\.032\.03to2\.332\.33across the seven parameters, all close to the theoretical5≈2\.24\\sqrt\{5\}\\approx 2\.24\. The observed\-Fisher SE matches the empirical SE within2\.9%2\.9\\%on average and within8\.3%8\.3\\%at worst, confirming that the Hessian\-based asymptotic approximation calibrates the actual sampling variability of the MLE on the Bell\-polynomial likelihood\. The bias remains uniformly small, under2\.7%2\.7\\%, and stable acrossNN, indicating that the residual bias arises from finite\-sample MLE bias rather than from convergence failures\.

Table 17:Replicated recovery on nested Clayton copulas withR=500R\{=\}500Monte Carlo replications per cell\. Each replication draws a fresh sample of sizeNNand refits the MLE\. EmpSE is the across\-replication standard deviation ofθ^\\hat\{\\theta\}; RMSE is its root\-mean\-squared error against the true value; Bias is the mean ofθ^−θtrue\\hat\{\\theta\}\-\\theta\_\{\\mathrm\{true\}\}\. HessSE isdiag​\(\(N​H^\)−1\)\\sqrt\{\\mathrm\{diag\}\(\(N\\hat\{H\}\)^\{\-1\}\)\}averaged across replications, whereH^=−N−1​∇2ℓ​\(θ^\)\\hat\{H\}=\-N^\{\-1\}\\nabla^\{2\}\\ell\(\\hat\{\\theta\}\)is the per\-observation observed information; every replication produces a positive\-definite Hessian\. Conv% is the L\-BFGS\-B convergence rate\.### L\.1Frank and Gumbel Recovery

The recovery result extends to Frank and Gumbel generators, confirming that the Bell\-polynomial likelihood and its AD gradients remain correct regardless of the generator choice\.[Table˜3](https://arxiv.org/html/2605.23134#S5.T3)uses Clayton exclusively because Clayton admits closed\-form Marshall–Olkin sampling;[Table˜18](https://arxiv.org/html/2605.23134#A12.T18)repeats the experiment with Frank and Gumbel generators and samples via the Rosenblatt transform\. L\-BFGS\-B recovers all parameters to within 5% atd=10d\{=\}10and within 2\.1% atd=20d\{=\}20for both families\. The Frankd=10d\{=\}10root parameter shows 13% error because only two sector\-level aggregates inform the root; atd=20d\{=\}20with four sectors this error drops below 1%\.

Table 18:Parameter recovery on nested Frank and Gumbel copulas \(N=500N\{=\}500, 2\-level nesting, L\-BFGS\-B optimiser\)\. The experiment samples data via the Rosenblatt transform, then recovers parameters using exact AD gradients\.

## Appendix MCensoring Validation

The per\-variable censoring mechanism reproduces the closed\-form mixed partial derivatives of the bivariate and trivariate Clayton and Gumbel CDFs to machine precision across every censoring pattern\. For the Clayton family with generatorψ​\(t\)=\(1\+t\)−1/θ\\psi\(t\)=\(1\+t\)^\{\-1/\\theta\}, thekk\-th order mixed partial of thedd\-variate CDF admits the closed form

∂kC∂u1​⋯​∂uk=∏j=0k−1\(1\+j​θ\)​∏i=1kui−\(θ\+1\)​\(∑i=1dui−θ−\(d−1\)\)−\(1/θ\+k\)\.\\frac\{\\partial^\{k\}C\}\{\\partial u\_\{1\}\\cdots\\partial u\_\{k\}\}=\\prod\_\{j=0\}^\{k\-1\}\(1\+j\\theta\)\\;\\prod\_\{i=1\}^\{k\}u\_\{i\}^\{\-\(\\theta\+1\)\}\\;\\Bigl\(\\textstyle\\sum\_\{i=1\}^\{d\}u\_\{i\}^\{\-\\theta\}\-\(d\-1\)\\Bigr\)^\{\-\(1/\\theta\+k\)\}\.\(32\)We compare the output of our Bell\-polynomial censoring implementation, which setsα\(0\)=\[1\]\\alpha^\{\(0\)\}=\[1\]for censored leaves, against this closed form ford∈\{2,3\}d\\in\\\{2,3\\\}with every possible censoring pattern running overk=0,…,dk=0,\\ldots,dobserved variables, 5–6 evaluation points per pattern, andθ∈\{0\.5,1,2,5\}\\theta\\in\\\{0\.5,1,2,5\\\}\. We perform analogous checks for the Gumbel family withψ​\(t\)=e−t1/θ\\psi\(t\)=e^\{\-t^\{1/\\theta\}\}andθ∈\{1\.5,2,3,5\}\\theta\\in\\\{1\.5,2,3,5\\\}, where nestedjax\.gradapplied to the known closed\-form CDF yields the reference values\.

Table 19:Censoring validation: framework output versus closed\-form expressions\. The test sweep covers every censoring pattern fromk=0k=0tok=dk=dobserved variables\.All 256 tests pass with a maximum absolute error of7\.1×10−157\.1\\times 10^\{\-15\}in log\-likelihood, confirming that the per\-variable censoring mechanism correctly reduces to the appropriate mixed partial derivative of the copula CDF for every censoring configuration in the sweep\.

## Appendix NSupported Archimedean Families

acopulaships with ten Archimedean families—the five classical generators plus five additional generators drawn from Nelsen’s catalogue\[Nelsen,[2006](https://arxiv.org/html/2605.23134#bib.bib33)\]that existing software cannot evaluate—and admits arbitrary user\-defined generators through a single decorator\.[Table˜20](https://arxiv.org/html/2605.23134#A14.T20)lists the included generators; the classical families above the rule are standard in the literature, and the five additional families below appear in our experiments\. Users can add further families by supplying only the generator function via the@copuladecorator, after whichacopuladerives all derivatives, inverses, and sampling procedures automatically\.

Table 20:One\-parameter Archimedean copula families included inacopula\. Classical families appear above the rule; additional families appear below\. Sampling uses numerical Laplace inversion for all families\.
## Appendix OHigh\-Dimensional Nested Scaling

This appendix establishes that theO​\(d2\)O\(d^\{2\}\)complexity bound of[Section˜3\.2](https://arxiv.org/html/2605.23134#S3.SS2.SSS0.Px1)holds in practice across three orders of magnitude indd, and that the prefactor stays small enough for single\-observation densities to evaluate in milliseconds atd=1,000d\{=\}1\{,\}000\. We measure single\-observation log\-likelihood wall\-clock for nested Clayton copulas atd∈\[10,8,000\]d\\in\[10,8\{,\}000\]on a single AMD Threadripper 3960X core under float64\. We scan three topologies:fixedK=5K\{=\}5two\-level \(root overM=d/5M\{=\}d/5Clayton sectors of five leaves\) and three\-level \(mid layer of five sectors per group\); andbalancedd\\sqrt\{d\}\(M=K=⌈d⌉M\{=\}K\{=\}\\lceil\\sqrt\{d\}\\,\\rceil\)\. All times are medians over five post\-warmup calls of a single JIT\-compiled XLA program;[Figure˜3](https://arxiv.org/html/2605.23134#S5.F3)plots the curves and[Table˜21](https://arxiv.org/html/2605.23134#A15.T21)lists selected dimensions\.

##### Two regimes\.

The empirical scaling exhibits two regimes separated byd≈1,000d\{\\approx\}1\{,\}000: a sub\-quadratic overhead\-dominated regime at the dimensions our experiments use, and the predictedO​\(d2\)O\(d^\{2\}\)asymptote pastd≈4,000d\{\\approx\}4\{,\}000\([Figure˜3](https://arxiv.org/html/2605.23134#S5.F3)\)\. In thereal\-data regimethat covers our experiments—retinopathy atd=2d\{=\}2, MIMIC\-IV atd=53d\{=\}53, and S&P 500 atd=98d\{=\}98—the empirical exponent is sub\-linear,α≈0\.55\\alpha\\approx 0\.55ond∈\[10,100\]d\\in\[10,100\]\. The Bell\-polynomial work scales asc⋅d2c\\cdot d^\{2\}withc≈10−5c\\approx 10^\{\-5\}ms, so per\-call XLA dispatch and the per\-sectorlax\.scanthat walks the tree dominate the small quadratic constant\. Across the measured regimed∈\[100,8,000\]d\\in\[100,8\{,\}000\]a global OLS log\-log fit givesα≈2\.03\\alpha\\approx 2\.03at 2\-level andα≈1\.92\\alpha\\approx 1\.92at 3\-level, both consistent with the theoreticalO​\(d2\)O\(d^\{2\}\)bound of[Section˜3\.2](https://arxiv.org/html/2605.23134#S3.SS2.SSS0.Px1)\. The local slopes in[Table˜21](https://arxiv.org/html/2605.23134#A15.T21)show a more nuanced picture: the 2\-level slope reachesαΔ=2\.06\\alpha\_\{\\Delta\}\{=\}2\.06atd=8,000d\{=\}8\{,\}000, settled at the asymptote, while the 3\-level slope atd=8,000d\{=\}8\{,\}000isαΔ=2\.89\\alpha\_\{\\Delta\}\{=\}2\.89, super\-quadratic\. This last value reflects L2 cache pressure as the per\-sector polynomial powering exceeds cache capacity, not an algorithmic shift; the pure\-Cauchy microbenchmark of[Section˜A\.6](https://arxiv.org/html/2605.23134#A1.SS6)isolates the cache transition from the algorithmic complexity\. The balancedd\\sqrt\{d\}topology hits the cache wall earlier because each sector containsK=⌈d⌉K\{=\}\\lceil\\sqrt\{d\}\\,\\rceilleaves—about8080atd=6,400d\{=\}6\{,\}400—inflating itsd≥1,000d\\geq 1\{,\}000slope toα≈3\.24\\alpha\\approx 3\.24\.

##### Head\-to\-head against RnacLL: balancedd\\sqrt\{d\}topology\.

The R wall atd≈40d\{\\approx\}40that[Figure˜3](https://arxiv.org/html/2605.23134#S5.F3)reports under the fixed\-K=5K\{=\}5topology persists under the harder balancedM=K=⌈d⌉M\{=\}K\{=\}\\lceil\\sqrt\{d\}\\,\\rceiltopology, where the per\-sector arity itself grows withdd\([Figure˜4](https://arxiv.org/html/2605.23134#A15.F4)\)\. The dominant per\-density cost in RnacLLis the partition enumeration inb\.coefffrom thecopulasource fileR/Rsource/dnac\.R, lines 95–145: for the Hofert–Pham formula, the inner loop isO​\(n⋅d0⋅∏vKv\)O\(n\\cdot d\_\{0\}\\cdot\\prod\_\{v\}K\_\{v\}\), whered0d\_\{0\}is the root arity andKvK\_\{v\}are the per\-sector sizes; with balanced sectors of sizeKKandM=d/KM\{=\}d/Ksectors this becomesO​\(KM\)O\(K^\{M\}\)per density, exponential in the number of sectors\. The classicalO​\(d2\)O\(d^\{2\}\)Stirling\-table closed form ofHofert et al\. \[[2013](https://arxiv.org/html/2605.23134#bib.bib19)\]applies only to flat single\-level Archimedean copulas; for nested copulas no Stirling shortcut exists onceM≥2M\\geq 2\. Timing R atn=1n\{=\}1requires a two\-line patch tonacLLthat addsdrop=FALSEatdnac\.Rline 189; upstreamnacLLonly handlesn≥2n\{\\geq\}2\.

10110^\{1\}10210^\{2\}10310^\{3\}10410^\{4\}10−110^\{\-1\}10010^\{0\}10110^\{1\}10210^\{2\}10310^\{3\}R walldimensionddwall\-clock per density evaluation \(ms\)Balancedd\\sqrt\{d\}head\-to\-head \(n=1n\{=\}1, CPU\)d2d^\{2\}referenceacopulaClaytonRnacLLClaytonRnacLLGumbelFigure 4:Head\-to\-head against RnacLLunder the balancedM=K=⌈d⌉M\{=\}K\{=\}\\lceil\\sqrt\{d\}\\,\\rceiltopology,n=1n\{=\}1single\-evaluation timings, single\-thread CPU, float64\. R aborts atd≥40d\{\\geq\}40for Clayton/Gumbel and on Frank for everydd\.acopulaextends tod=6,400d\{=\}6\{,\}400with the sameO​\(d2\)O\(d^\{2\}\)asymptote as the fixed\-K=5K\{=\}5topology \([Figure˜3](https://arxiv.org/html/2605.23134#S5.F3)\); the slope steepening ind∈\[1,600,6,400\]d\\in\[1\{,\}600,\\,6\{,\}400\]is a hardware cache transition as each sector sizeK=⌈d⌉K\{=\}\\lceil\\sqrt\{d\}\\,\\rceilapproaches8080leaves \([Section˜A\.6](https://arxiv.org/html/2605.23134#A1.SS6)\)\.Table 21:Single\-observation wall\-clock at selected dimensions on AMD Threadripper 3960X \(single thread, float64\)\.αΔ\\alpha\_\{\\Delta\}is the local log\-log slope from the previous row, illustrating the transition from overhead\-dominated \(αΔ≪2\\alpha\_\{\\Delta\}\\ll 2\) at smallddto theO​\(d2\)O\(d^\{2\}\)asymptote \(αΔ≈2\\alpha\_\{\\Delta\}\\approx 2\) at largedd; intermediate slopes also reflect L1/L2 cache transitions, not algorithmic shifts\.

## Appendix PMemory and Compilation Benchmarks

The compilation cost a user pays once at fit time grows linearly indd, the gradient program stays at the standard reverse\-mode∼3×\{\\sim\}3\\timesmultiple of the forward program, and the entire fused HLO fits in a few thousand instructions—small enough that the JIT cache amortises across the inner\-loop evaluations of any optimiser\.[Table˜22](https://arxiv.org/html/2605.23134#A16.T22)reports compilation overhead, execution time, and XLA program size for the Bell polynomial likelihood and its gradient at increasing dimensionsdd\. All measurements use a two\-level nested Clayton copula at sector sizeK=5K\{=\}5, on an AMD Threadripper 3960X CPU \(3\.8 GHz, 24 cores / 48 threads\) with JAX 0\.8\.3 in float64 mode\. We count HLO instructions in the single fused program XLA emits viajax\.jit\(f\)\.lower\(x\)\.compile\(\)\.as\_text\(\)\.

Table 22:Compilation and execution benchmarks for the Bell polynomial likelihood and gradient on CPU \(AMD Threadripper 3960X,n=1n\{=\}1observation\)\. “JIT” denotes one\-time compilation time; “Exec” denotes the median execution time over 10 post\-warmup runs\. “HLO ops” counts instructions in the compiled XLA program\.Three patterns stand out from[Table˜22](https://arxiv.org/html/2605.23134#A16.T22)\. First, JIT compilation time stays nearly flat inddbecause the stage\-group batched likelihood compiles each of the two nesting levels into one fusedlax\.scan; growingddextends the scan length but not the compiled program\. Second, the compiled HLO program grows minimally—from 1 700 to 1 800 forward ops and 5 000 to 5 300 gradient ops—confirming thatlax\.scanreuses the same loop body rather than unrolling\. Third, the gradient program contains roughly3×3\\timesthe forward operations, matching the standard reverse\-mode AD multiple for scan\-based programs\.

## Appendix QExact MLE vs\. Simulated Maximum Likelihood

Exact MLE through our Bell\-polynomial likelihood dominates thesimulated maximum likelihood estimation\(SMLE\) approach that Rnacopula\[Hofert et al\.,[2013](https://arxiv.org/html/2605.23134#bib.bib19)\]relies on, in both parameter accuracy and wall time, at every SMLE sample budget we test\. The SMLE estimator drawsNNfrailty samples from the Laplace\-transform hierarchy at each parameter evaluation, computes the conditional density at each observation given those frailties, and averages over theNNdraws to approximate the intractable marginal density\. Because the outer copula parameter enters only through the frailty distribution—not through the conditional density—the SMLE objective demands fresh frailty samples at every evaluation, producing a noisy objective that forces derivative\-free optimisation\. Our framework instead computes the exact density and exact gradients, opening the door to L\-BFGS\-B with analytic gradients\.

##### Setup\.

We sample 500 observations from a two\-level nested Clayton copula atd=10d\{=\}10with two sectors of sizeK=5K\{=\}5,θouter=1\.5\\theta\_\{\\mathrm\{outer\}\}\{=\}1\.5, andθinner=3\.0\\theta\_\{\\mathrm\{inner\}\}\{=\}3\.0\. For each of 30 replications, we fit four estimators on the same data: exact MLE via our Bell\-polynomial likelihood, and SMLE atN∈\{100,1,000,10,000\}N\\in\\\{100,\\;1\{,\}000,\\;10\{,\}000\\\}frailty samples per evaluation\. The exact MLE uses L\-BFGS\-B with analytic gradients; SMLE uses differential evolution—a stochastic global optimiser for noisy objectives—with frailty samples re\-drawn each evaluation\.

Table 23:Exact MLE vs\. SMLE on simulated nested Clayton \(d=10d\{=\}10, sample sizen=500n\{=\}500, 30 replications;NNdenotes the number of frailty samples per SMLE evaluation, distinct from the sample sizenn\)\. True parameters:θouter=1\.5\\theta\_\{\\mathrm\{outer\}\}\{=\}1\.5,θinner=3\.0\\theta\_\{\\mathrm\{inner\}\}\{=\}3\.0\. RMSEo, RMSEidenote the root\-mean\-squared error forθ^outer\\hat\{\\theta\}\_\{\\mathrm\{outer\}\}andθ^inner\\hat\{\\theta\}\_\{\\mathrm\{inner\}\}respectively\.[Table˜23](https://arxiv.org/html/2605.23134#A17.T23)reports the results\. The exact MLE recovers both parameters with near\-zero bias, RMSE≈0\.07\\approx 0\.07on each, in 1\.5 s per fit\. SMLE atN=100N\{=\}100underestimatesθinner\\theta\_\{\\mathrm\{inner\}\}by0\.360\.36on average, with an RMSE of0\.400\.40—roughly six times the exact MLE’s error at four times the wall time\. IncreasingNNto10,00010\{,\}000reduces the bias to below0\.020\.02and brings the RMSE to0\.070\.07onθinner\\theta\_\{\\mathrm\{inner\}\}, matching the exact MLE’s accuracy—but at153153\\,s per fit,one hundred times slowerthan exact inference on the same data\. For the moderate dimensions typical of actuarial and survival\-analysis applications,d≤100d\\leq 100, the exact approach dominates SMLE in both accuracy and wall time at every SMLE sample budget tested\. Differentiability compounds the advantage: L\-BFGS\-B with analytic gradients converges in far fewer iterations than the derivative\-free global search SMLE’s noisy objective requires\.

##### Generator\-derivative Monte Carlo: Gen\-AC ablation\.

The original Gen\-AC training procedure\[Ng et al\.,[2021](https://arxiv.org/html/2605.23134#bib.bib34)\]replacesψ\(k\)​\(t\)=𝔼U​\[\(−eMLP​\(U\)\)k​exp⁡\(−eMLP​\(U\)​t\)\]\\psi^\{\(k\)\}\(t\)=\\mathbb\{E\}\_\{U\}\[\(\-e^\{\\mathrm\{MLP\}\(U\)\}\)^\{k\}\\exp\(\-e^\{\\mathrm\{MLP\}\(U\)\}t\)\]with anNN\-sample Monte Carlo averageψ^N\(k\)\\hat\{\\psi\}^\{\(k\)\}\_\{N\}, mirroring SMLE’s role above but at the generator\-derivative level rather than the joint\-likelihood level\. The resulting log\-likelihood is biased by Jensen’s inequality onlog\\logcomposed with the nonlinear Hofert density formula inψ^\(k\)\\hat\{\\psi\}^\{\(k\)\}, with biasO​\(1/N\)O\(1/N\)—biased but consistent\. We measure this bias on Gen\-AC\-Fourier withK=2K\{=\}2Fourier features, hidden width33,262262parameters, MIMIC\-IVd=53d\{=\}53,N=32N\{=\}32frailty samples per evaluation, seed4242\.[Table˜24](https://arxiv.org/html/2605.23134#A17.T24)reports a four\-cell ablation: each of two trained models \(θstoch\\theta\_\{\\mathrm\{stoch\}\}trained with re\-sampled Monte Carlo gradients,θquad\\theta\_\{\\mathrm\{quad\}\}trained with deterministic quadrature gradients\) evaluated under both estimators \(the stochastic estimator averages5050independent draws to estimate𝔼​\[stoch NLL\]\\mathbb\{E\}\[\\text\{stoch NLL\}\]\)\. Eval\-time Jensen plus correlation bias accounts for0\.090\.09–0\.120\.12nats; train\-time optimisation differences contribute only0\.020\.02nats; quadrature evaluation eliminates both\.

Table 24:Gen\-AC\-Fourier260on MIMIC\-IV \(d=53d\{=\}53\): held\-out NLL under each combination of training and evaluation estimator, seed4242\. Lower is better\. Rows index the trained model, columns the evaluation method;θstoch\\theta\_\{\\mathrm\{stoch\}\}trains withN=32N\{=\}32Monte Carlo frailty draws per Adam step,θquad\\theta\_\{\\mathrm\{quad\}\}trains withN=32N\{=\}32Gauss–Legendre quadrature nodes per Adam step\.

## Appendix RPredictive Validation

We report held\-out predictive log\-likelihoods for the two real\-data experiments—MIMIC\-IV atd=53d\{=\}53and the diabetic retinopathy paired\-eye study atd=2d\{=\}2—to verify that the in\-sample family rankings of[Sections˜5\.3](https://arxiv.org/html/2605.23134#S5.SS3)and[5\.6](https://arxiv.org/html/2605.23134#S5.SS6)survive an 80/20 train/test split\.

### R\.1MIMIC\-IV

The MIMIC\-IV ablations break into three pieces: a held\-out evaluation that confirms the in\-sample family ranking, a compile\-cost audit that quantifies the one\-time XLA price of the Frank generator, and a Hessian\-based standard\-error check that confirms parameter identifiability atntrain=68,183n\_\{\\mathrm\{train\}\}\{=\}68\{,\}183\.

##### Cohort construction\.

The cohort draws from MIMIC\-IV 3\.1 on PhysioNet: one row perhadm\_idwith an ICU stay, joined oncharttimeto lab events that fall inside the ICU window\[intime, outtime\]\. For each admission\-lab pair, we recordTi​jT\_\{ij\}, the hours from ICU intime to the first abnormal reading of labjj, and right\-censor at the ICU length\-of\-stay when no abnormal reading occurs during the window\. The “abnormal” label uses MIMIC’s nativeflag=‘abnormal’column, which MIMIC curators set against each lab’s own reference range; we adopt this labelling verbatim and apply no additional reference range\. A lab enters the cohort only when at least30%30\\%of admissions have a measurement, at least5%5\\%of admissions have an abnormal reading, and at least90%90\\%of its rows carry a numericvaluenum\. These thresholds select thed=53d\{=\}53itemids in[Table˜25](https://arxiv.org/html/2605.23134#A18.T25), grouped into six named lab panels plus one catch\-all following the clinical groupings hard\-coded inexperiments/shared/mimic\.py\. Runningexperiments/data/build\_mimic\.pyon the MIMIC\-IV 3\.1 raw CSVs materialises the cohort as an\(n=85,229,d=53\)\(n\{=\}85\{,\}229,\\,d\{=\}53\)censored matrix inmimic\_full\_censored\.npz\.

Table 25:MIMIC\-IV 3\.1 cohort panels used in the nested copula fit\. Itemids come from MIMIC’sd\_labitemsdictionary; the “Other” panel collects the1414labs that pass the quality filters but fall outside the six named groupings, with the build\-time pipeline fixing its membership\.
##### Held\-out evaluation\.

We repeat the held\-out protocol on the MIMIC\-IV ICU lab dataset withd=53d\{=\}53labs in a seven\-panel nested hierarchy, using an 80/20 admission\-level split at seed 42\. Stage 1 estimates each lab’s marginal survivor function with Kaplan–Meier on training data only; Stage 2 fits copula parameters on the training pseudo\-observations via L\-BFGS\-B with analytic gradients\.[Table˜26](https://arxiv.org/html/2605.23134#A18.T26)reports per\-observation NLL on both splits \(ntrain=68,183n\_\{\\mathrm\{train\}\}\{=\}68\{,\}183,ntest=17,046n\_\{\\mathrm\{test\}\}\{=\}17\{,\}046\)\.

Table 26:MIMIC\-IV: held\-out predictive NLL per observation\. The symbol−ℓ¯\-\\bar\{\\ell\}denotes the per\-observation negative log\-likelihood; bold marks the best \(lowest\) test value\.ntrain=68,183n\_\{\\mathrm\{train\}\}\{=\}68\{,\}183,d=53d\{=\}53,88copula parameters\.The held\-out ranking matches the in\-sample AIC ranking: Clayton/Clayton attains the lowest test NLL, with essentially no train\-to\-test gap \(4\.9514\.951vs4\.9494\.949\) thanks to the fully non\-parametric Kaplan–Meier marginals that avoid the per\-lab tail\-misspecification bias\. Frank/Frank and Gumbel/Gumbel both fit worse than Clayton at every split\.[Table˜27](https://arxiv.org/html/2605.23134#A18.T27)reports one\-time fresh\-JIT costs on an RTX 3090 Ti for the MIMIC\-IV dynamic\-censoring pipeline atntrain=68,183n\_\{\\mathrm\{train\}\}\{=\}68\{,\}183\.acopula\.set\_compile\_cache\_dircaches every compile to disk, so subsequent runs read the compiled HLO program back in under a second\.

Table 27:MIMIC\-IV: one\-time XLA compile wall time on an RTX 3090 Ti\. Theval\+gradcolumn times JIT compilation ofjax\.value\_and\_gradon the chunked negative log\-likelihood; theHessian computecolumn times the chunked Hessian pass that supplies observed Fisher standard errors\.A warm cache reproduces the full[Table˜28](https://arxiv.org/html/2605.23134#A18.T28)in under five minutes \(L\-BFGS\-B fit plus chunked Hessian compute on already\-compiled kernels\)\.

##### Standard errors\.

jax\.hessianon the MIMIC\-IV fits yields positive\-definite Hessians and tight standard errors for all three nestings atntrain=68,183n\_\{\\mathrm\{train\}\}\{=\}68\{,\}183; see[Table˜28](https://arxiv.org/html/2605.23134#A18.T28)\. Every coefficient of variation \(CV\) stays below1\.5%1\.5\\%, so the L\-BFGS\-B optima are sharply identified\.

Table 28:MIMIC\-IV: Hessian\-based standard errors across three nestings, computed byjax\.hessianwith a delta\-method projection to the constrained space\. All three Hessians are positive definite, and the coefficient of variation stays below1\.5%1\.5\\%for every parameter\.ntrain=68,183n\_\{\\mathrm\{train\}\}\{=\}68\{,\}183\.
#### R\.1\.1Ablation: independence and single\-layer baselines

To quantify the contribution of the hierarchical nesting structure, we evaluate two simpler alternatives on the same held\-out split: an independence copula whose mixed partial ofC​\(𝐮\)=∏jujC\(\\mathbf\{u\}\)\{=\}\\prod\_\{j\}u\_\{j\}in the uncensored variables collapses to∏j:δj=0uj\\prod\_\{j:\\,\\delta\_\{j\}=0\}u\_\{j\}, contributing∑j:δj=0log⁡uj<0\\sum\_\{j:\\,\\delta\_\{j\}=0\}\\log u\_\{j\}<0per observation; and a single\-layer copula that places every lab under a single root node with one parameter, ignoring the panel grouping\.[Table˜29](https://arxiv.org/html/2605.23134#A18.T29)reports the results\. The hierarchical structure captures intra\-panel dependence—the near\-deterministic relationships within the coagulation panel and the differential count—that a single exchangeable copula cannot represent\. The eight\-parameter nested Clayton improves test NLL substantially over the best flat copula, because the panel structure absorbs heterogeneity across the5353labs into per\-panel parameters instead of forcing a single exchangeable parameter to average over the cohort\.

Table 29:MIMIC\-IV ablation: independence and single\-layer baselines versus the best nested model on the held\-out split \(ntrain=68,183n\_\{\\mathrm\{train\}\}\{=\}68\{,\}183,d=53d\{=\}53,88copula parameters, Kaplan–Meier marginals throughout\)\. Flat copulas use a single exchangeable parameter regardless of dimension\.
#### R\.1\.2Marginal misspecification

We refit each of the three classical pairs under two parametric marginals—per\-lab Weibull and per\-lab lognormal fitted by maximum likelihood with right\-censoring—in addition to the Kaplan–Meier marginals of the body, on the full cohort withn=85,229n\{=\}85\{,\}229\.[Table˜30](https://arxiv.org/html/2605.23134#A18.T30)reports per\-observation copula NLL and the fitted root parameterθ^outer\\hat\{\\theta\}\_\{\\mathrm\{outer\}\}, and two findings emerge\. First, theparametersof each family stay stable across the three marginals:θ^outer\\hat\{\\theta\}\_\{\\mathrm\{outer\}\}varies by at most30%30\\%for Clayton,9%9\\%for Frank, and stays essentially identical for Gumbel at1\.1701\.170–1\.1711\.171across all three marginals; the per\-panel inner parameters move by at most7%7\\%per panel\. Second, therankingof copula families changes with the marginal: Clayton wins under Kaplan–Meier and lognormal, while Weibull misspecification creates residual upper\-tail dependence that shifts the winner to Gumbel\. The framework’s parameter estimates are therefore robust to marginal choice, but downstream model selection should be cross\-validated against the marginal model rather than assumed transferable\.

Table 30:MIMIC\-IV marginal sensitivity \(n=85,229n\{=\}85\{,\}229,d=53d\{=\}53, full cohort,88copula parameters per fit\)\. Per\-observation copula NLL and root parameterθ^outer\\hat\{\\theta\}\_\{\\mathrm\{outer\}\}across three marginal families\. Per\-family parameter estimates are stable across marginals \(Gumbel especially\); the family ranking is not\.

### R\.2Retinopathy

Beyond the high\-dimensional MIMIC\-IV setting, we also evaluate the framework on the diabetic retinopathy paired\-eye study, which exercises the small\-sample heavy\-censoring regime withn=197n\{=\}197and censoring fraction60\.7%60\.7\\%, where parameter identifiability is the binding constraint\. To verify that the family rankings of[Section˜5\.6](https://arxiv.org/html/2605.23134#S5.SS6)are not artefacts of in\-sample overfitting, we perform held\-out predictive checks with an 80/20 train/test split at seed 42\.

Table 31:Retinopathy: held\-out copula log\-likelihood for all 10 generator families, withntrain=157n\_\{\\mathrm\{train\}\}\{=\}157,ntest=40n\_\{\\mathrm\{test\}\}\{=\}40, and seed 42\. We fit the Weibull marginals on training data alone; the copula negative log\-likelihood−ℓcop\-\\ell\_\{\\mathrm\{cop\}\}isolates the dependence contribution\. The dagger†marks generators with no prior censored\-copula implementation\.For retinopathy, the Nelsen12†generator—the Burr copula—achieves the best held\-out log\-likelihood among all1010families in[Table˜31](https://arxiv.org/html/2605.23134#A18.T31), although it converges to its lower parameter bound atθ^≈1\\hat\{\\theta\}\\approx 1\. Atθ=1\\theta\{=\}1the Nelsen12 generatorψ​\(t\)=\(1\+t1/θ\)−1\\psi\(t\)\{=\}\(1\{\+\}t^\{1/\\theta\}\)^\{\-1\}reduces to Clayton\(θ=1\)\(\\theta\{=\}1\), so the Nelsen12 fit is structurally Clayton constrained toθ≥1\\theta\{\\geq\}1; the unconstrained Clayton MLE in[Table˜31](https://arxiv.org/html/2605.23134#A18.T31)sits just below that bound atθ^=0\.61\\hat\{\\theta\}\{=\}0\.61, andLABEL:app:retino\-hessianreports a full\-data Clayton MLE ofθ^=0\.90±0\.31\\hat\{\\theta\}\{=\}0\.90\\pm 0\.31, both consistent with the Nelsen12 boundary fit being Clayton\-like\. Frank and InverseGaussian†form a tight cluster in second place, and the held\-out ranking closely mirrors the in\-sample AIC ordering\. Nelsen13†and Nelsen17†sit immediately behind that cluster with NLLs that differ in the third decimal \(21\.21621\.216vs\.21\.22321\.223\) at distinct fitted parametersθ^=1\.91\\hat\{\\theta\}\{=\}1\.91andθ^=2\.04\\hat\{\\theta\}\{=\}2\.04, confirming the two families are not algebraically equivalent\. The spread among the top nine families remains modest, withΔ−ℓtest≈1\.2\\Delta\{\-\}\\ell\_\{\\mathrm\{test\}\}\\approx 1\.2, consistent with the similar AIC values reported in the main table; Nelsen9†is the sole outlier, drifting toward the independence limitθ^→0\+\\hat\{\\theta\}\\to 0^\{\+\}and fitting poorly\.

##### Standard errors\.

We compute Hessian\-based standard errors from the observed Fisher information returned byjax\.hessian, with a delta\-method projection to the constrained parameter space; the top families show well\-identified parameters: Frankθ^=2\.25±0\.62\\hat\{\\theta\}\{=\}2\.25\\pm 0\.62at CV27%27\\%, Gumbelθ^=1\.25±0\.08\\hat\{\\theta\}\{=\}1\.25\\pm 0\.08at CV7%7\\%, and Claytonθ^=0\.90±0\.31\\hat\{\\theta\}\{=\}0\.90\\pm 0\.31at CV34%34\\%\. The higher CVs reflect the modest sample size ofn=197n\{=\}197and heavy censoring of60\.7%60\.7\\%; all Hessians remain positive definite, confirming MLEs at proper maxima\.

## Appendix SHead\-to\-Head with HACSurv on Framingham

HACSurv\[Liu et al\.,[2025b](https://arxiv.org/html/2605.23134#bib.bib29)\]is the closest neighbour to our framework in the competing\-risks regime \(up to five risks plus a censoring node\): a nested Gen\-AC Archimedean copula with neural density estimator \(NDE\) marginals, evaluated by an iterative inverse\-function\-theorem \(IFT\) solver for the generator inverse\. Our framework reproduces HACSurv’s per\-observation log copula partial to machine precision, isolating their density machinery from ours\. We verify the claim in two steps\. First, we replicate HACSurv’s published Framingham result\[Liu et al\.,[2025b](https://arxiv.org/html/2605.23134#bib.bib29)\]by running the authors’ code as released, with no modification beyond a CUDA\-device patch and a SciPy 1\.14simps\-API shim \([Table˜5](https://arxiv.org/html/2605.23134#S5.T5)\)\. Second, on the same trained checkpoints, ouracopulaframework—with the nested Gen\-AC generator of[Section˜G\.1](https://arxiv.org/html/2605.23134#A7.SS1)—computes the per\-observation copula log partial that enters HACSurv’s competing\-risks log\-likelihood to within machine precision across all five seeds and all three event classes \([Table˜32](https://arxiv.org/html/2605.23134#A19.T32)\); the comparison isolates the density machinery—HACSurv’s iterative IFT solve versus our exact compositional Bell polynomial powering—because both paths consume the same generator parametersMouter,Minner,μ,βM\_\{\\mathrm\{outer\}\},M\_\{\\mathrm\{inner\}\},\\mu,\\betaand the same NDE marginal evaluations\.

Table 32:Per\-observation log copula partiallog​∂jC​\(SE1,SE2,SC\)\\log\\partial\_\{j\}C\\\!\\left\(S\_\{\\mathrm\{E1\}\},S\_\{\\mathrm\{E2\}\},S\_\{\\mathrm\{C\}\}\\right\)from HACSurv’s iterative IFT solver compared with ouracopulaframework, on128128test observations per event class per seed\. Here,SE1,SE2S\_\{\\mathrm\{E1\}\},S\_\{\\mathrm\{E2\}\}are the NDE\-estimated marginal survivals at the two event times,SCS\_\{\\mathrm\{C\}\}the censoring survival, and∂jC:=∂C/∂uj\\partial\_\{j\}C:=\\partial C/\\partial u\_\{j\}the copula partial in variablejj\(the row’s event class\)\. Both paths share the same nested\-Archimedean generator and NDE marginals \(parametersMouter,Minner,μ,βM\_\{\\mathrm\{outer\}\},M\_\{\\mathrm\{inner\}\},\\mu,\\betaas in[Section˜G\.1](https://arxiv.org/html/2605.23134#A7.SS1)\), isolating the density\-evaluation step\.The largest seed\-wise discrepancy reaches2\.4×10−82\.4\\times 10^\{\-8\}on a single seed\-43 censored observation whose argument sits near the boundary of the unit cube, where HACSurv’s IFT iteration loses digits\[Liu et al\.,[2025b](https://arxiv.org/html/2605.23134#bib.bib29)\]; across all observations, the median agreement holds at machine precision near10−1510^\{\-15\}\. The replication establishes that ouracopulaframework reproduces HACSurv’s Framingham neural Archimedean copula likelihood\. The framework’s value\-add over HACSurv is the per\-variable censoring demonstration atd=53d\{=\}53in[Section˜5\.3](https://arxiv.org/html/2605.23134#S5.SS3), with the neural Archimedean generator extension in[Section˜5\.4](https://arxiv.org/html/2605.23134#S5.SS4)—a setting that HACSurv’s bivariate\-competing\-risks formulation cannot express\.

##### End\-to\-end JAX port with exact Bell polynomial density\.

Beyond evaluating the copula partial on HACSurv’s trained checkpoints, ouracopulaframework trains the full HACSurv competing\-risks model end\-to\-end without reusing any of their PyTorch code\. We port the HACSurv architecture to JAX/Equinox: the shared covariate embedding, three neural density estimator marginals with positive\-weight outcome networks for the two risks plus censoring, and the nested Gen\-AC Archimedean copula of[Section˜G\.1](https://arxiv.org/html/2605.23134#A7.SS1)\. The frailty isV=exp⁡\(MLP​\(U\)\)V\\\!=\\\!\\exp\(\\mathrm\{MLP\}\(U\)\)withU∼Uniform​\(0,1\)U\\sim\\mathrm\{Uniform\}\(0,1\), and the Laplace transformψ​\(t\)=𝔼U​\[exp⁡\(−V​t\)\]\\psi\(t\)\\\!=\\\!\\mathbb\{E\}\_\{U\}\[\\exp\(\-Vt\)\]is estimated at every optimiser step byNz=100N\_\{z\}\{=\}100Monte\-Carlo samples ofUU, matching the simulated\-maximum\-likelihood \(SMLE\) objective HACSurv trains\. The first\-order copula partial routes throughacopula’s Bell polynomial density path—the same machinery that powers thed=53d\{=\}53MIMIC\-IV experiment in[Section˜5\.4](https://arxiv.org/html/2605.23134#S5.SS4)—with per\-subject censoring masks selecting the single variable differentiated per competing\-risks observation\. Training reproduces HACSurv’s three\-stage pipeline: pairwise bivariate copula fits over the three event pairs, inner\-generator fitting on samples from the\(1,2\)\(1,2\)pair, and copula\-frozen marginal fine\-tuning with early stopping on the validation Conditional\-CIF C\-index\. We report Antolini’s time\-dependent C\-index and the IPCW\-weighted integrated Brier score on the same five HACSurv seeds in the Marginal\-SF prediction mode \([Table˜5](https://arxiv.org/html/2605.23134#S5.T5)\)\.

## Appendix TS&P 500 Sector\-Level Tail Dependence

The per\-sectorθinner\\theta\_\{\\mathrm\{inner\}\}values fitted atd=98d\{=\}98in[Section˜5\.5](https://arxiv.org/html/2605.23134#S5.SS5)reproduce the empirical tail behaviour of equity returns within each sector, confirming that the Gumbel/Gumbel hierarchy captures the upper\-tail piece of the joint distribution\. For every pair of stocks within a sector, withU,VU,Vthe copula\-scale pseudo\-observations of the two stocks, we compute the empirical upper\- and lower\-tail dependence coefficientsλ^U​\(q\)=Pr^​\(U\>1−q​∣V\>​1−q\)\\hat\{\\lambda\}\_\{U\}\(q\)=\\hat\{\\Pr\}\(U\>1\{\-\}q\\mid V\>1\{\-\}q\)andλ^L​\(q\)=Pr^​\(U​<q∣​V<q\)\\hat\{\\lambda\}\_\{L\}\(q\)=\\hat\{\\Pr\}\(U<q\\mid V<q\)atq∈\{0\.05,0\.10\}q\\in\\\{0\.05,0\.10\\\}, wherePr^\\hat\{\\Pr\}denotes the empirical relative\-frequency estimator over trading days; we then average over pairs within the sector and compare to the Gumbel\-implied upper\-tail dependenceλU​\(θ\)=2−21/θ\\lambda\_\{U\}\(\\theta\)=2\-2^\{1/\\theta\}\.[Table˜33](https://arxiv.org/html/2605.23134#A20.T33)reportsq=0\.05q\{=\}0\.05;q=0\.10q\{=\}0\.10shows the same pattern\.

Table 33:S&P 500d=98d\{=\}98Gumbel/Gumbel fit: per\-sectorθinner\\theta\_\{\\mathrm\{inner\}\}, theoretical upper\-tail dependence implied by Gumbel, and empirical upper/lower\-tail dependence at quantileq=0\.05q\{=\}0\.05, averaged over within\-sector stock pairs\.Two qualitative patterns emerge from the table\. First, the within\-sectorθinner\\theta\_\{\\mathrm\{inner\}\}values track empirical co\-movement: Energy and Utilities, the two sectors with the strongest physical co\-dependence, fit the highestθinner\\theta\_\{\\mathrm\{inner\}\}at1\.881\.88and1\.751\.75, while Communication Services and Consumer Discretionary—historically diverse mixes of independent business models—fit the lowest\. Second, every sector exhibits empirical lower\-tail dependence slightly above its empirical upper\-tail dependence \(average asymmetryλ^L−λ^U≈0\.03\\hat\{\\lambda\}\_\{L\}\-\\hat\{\\lambda\}\_\{U\}\\approx 0\.03atq=0\.05q\{=\}0\.05\), consistent with the well\-known crash\-driven correlation surge in equity returns\. Gumbel imposesλL=0\\lambda\_\{L\}=0by construction and so cannot represent the lower\-tail piece exactly; despite this limitation, Gumbel still attains the best AIC atd=98d\{=\}98in[Section˜5\.5](https://arxiv.org/html/2605.23134#S5.SS5), because the upper\-tail magnitudes it does capture explain more of the joint structure than Clayton or Frank reach with the same1212parameters\. A practitioner whose data exhibits lower\-tail asymmetry should swap Gumbel for Clayton \(λL\>0\\lambda\_\{L\}\>0,λU=0\\lambda\_\{U\}=0\) or for the180​°180°survival rotation of Gumbel/Joe; the framework treats all options with the same code path\.

## Appendix UVine Copula Comparison

Vine copulas trade parsimony for flexibility, so a head\-to\-head against our1212\-parameter nested Gumbel quantifies what we give up by collapsing pairwise dependence into a hierarchy\. Vine copulas\[Aas et al\.,[2009](https://arxiv.org/html/2605.23134#bib.bib1)\]model non\-exchangeable dependence via pair\-copula constructions;[Table˜34](https://arxiv.org/html/2605.23134#A21.T34)compares our nested Archimedean model—Gumbel, the best\-performing family in[Section˜5\.5](https://arxiv.org/html/2605.23134#S5.SS5)—against an R\-vine fitted byrvinecopulib\[Nagler and Vatter,[2023](https://arxiv.org/html/2605.23134#bib.bib32)\]on the S&P 500 pseudo\-observations\.

Table 34:Nested Archimedean vs\. R\-vine on S&P 500 data,n=1,253n\{=\}1\{,\}253trading days\.−ℓ¯\-\\bar\{\\ell\}is the per\-observation negative log\-likelihood, lower is better; AIC/n/nis the per\-observation Akaike Information Criterion \(\(2​\|θ\|−2​ℓ\)/n\\bigl\(2\|\\theta\|\-2\\ell\\bigr\)/n, lower is better\)\. All values in nats\.The R\-vine achieves lower negative log\-likelihood at every dimension, reflecting its greater flexibility withO​\(d2\)O\(d^\{2\}\)pair\-copula parameters versus our 12 hierarchical parameters\. This flexibility comes at a cost: the R\-vine atd=98d\{=\}98requires over4,7004\{,\}700pair\-copulas, each with its own dependence parameter, and those parameters lack the direct hierarchical interpretation that nested Archimedean parameters provide as within\-sector versus between\-sector dependence\. Standard vine implementations do not natively support censored observations at high dimension\.Stöber et al\. \[[2013](https://arxiv.org/html/2605.23134#bib.bib40)\]extends D\-vines to a single shared right\-censoring time per observation, sufficient for survival regression with one event but not for the per\-variable masks of[Section˜5\.3](https://arxiv.org/html/2605.23134#S5.SS3)\.Chen et al\. \[[2025](https://arxiv.org/html/2605.23134#bib.bib7)\]extends C\-vines to per\-variable right\-censoring under informative\-censoring assumptions, demonstrated on small dimension; the per\-observation cost is a multivariate numerical integration over the censored variables, infeasible atd=53d\{=\}53\. The computational pattern differs between the two model classes\. Vine is density\-first—c​\(𝐮\)=∏ece​\(Fie\|De,Fje\|De\)c\(\\mathbf\{u\}\)=\\prod\_\{e\}c\_\{e\}\(F\_\{i\_\{e\}\|D\_\{e\}\},F\_\{j\_\{e\}\|D\_\{e\}\}\)is closed\-form, but the censored likelihood requires the partial CDF∂\|δ\|C/∏j:δj=1∂uj\\partial^\{\|\\delta\|\}C/\\prod\_\{j:\\delta\_\{j\}=1\}\\partial u\_\{j\}, and recovering this from a density needs multivariate integration over the censored variables, repeated per observation since each mask defines a different integration domain\. Nested Archimedean is CDF\-first—C​\(𝐮\)=ψ​\(∑jψ−1​\(uj\)\)C\(\\mathbf\{u\}\)=\\psi\(\\sum\_\{j\}\\psi^\{\-1\}\(u\_\{j\}\)\)is closed\-form, the censored variables enter the inner sum as constants, and the partial is an algebraic Bell\-polynomial recurrence \([Section˜3\.2](https://arxiv.org/html/2605.23134#S3.SS2.SSS0.Px2)\)\. The two model classes therefore serve complementary roles: vines suit settings where flexible pairwise dependence matters more than parsimony and where censoring is absent or shared, while nested Archimedean copulas suit hierarchical data with arbitrary per\-variable censoring\.

## Appendix VHierarchical Kendall Copulas

Hierarchical Kendall copulas \(HKCs\) sit between vines and nested Archimedean copulas as a third hierarchical construction, but they aggregate clusters through a scalar Kendall transform rather than the full multivariate density—so they cannot supply the analytic gradients or per\-variable censoring that our framework provides on the same hierarchical structure\.Brechmann \[[2014](https://arxiv.org/html/2605.23134#bib.bib6)\]introducedhierarchical Kendall copulasas a construction intermediate between vines and nested Archimedean copulas\. HKCs aggregate each cluster of variables through theKendall distribution functionKC​\(t\)=Pr⁡\{C​\(𝐔\)≤t\}K\_\{C\}\(t\)=\\Pr\\\{C\(\\mathbf\{U\}\)\\leq t\\\}, building the joint distribution from cluster copulas linked by univariate Kendall transforms rather than from generator nesting\.Brechmann \[[2014](https://arxiv.org/html/2605.23134#bib.bib6)\]derives closed\-form HKC density expressions for the special case where every cluster copula is itself Archimedean, drawing on the analytical Kendall function of each generator; for non\-Archimedean cluster copulas the density requires multivariate integration, whichBrechmann \[[2014](https://arxiv.org/html/2605.23134#bib.bib6)\]approximates by sampling rather than evaluates in closed form\. Standard HKC implementations target sampling and pseudo\-likelihood estimation through L\-BFGS without analytic parameter gradients, and they do not support censored observations\. Ouracopulaframework therefore complements HKCs on the same hierarchical structure with two extra properties: exact analytic gradients through the full multivariate density, and per\-variable censoring\.