lrnv · lrnv · Oct 8, 2025 · Oct 2, 2025 · Oct 3, 2025 · Oct 3, 2025
diff --git a/Copulas.jl b/Copulas.jl
diff --git a/Project.toml b/Project.toml
@@ -9,7 +9,6 @@ Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 HCubature = "19dc6840-f33b-545b-b366-655c7e3ffd49"
-InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 LambertW = "984bce1d-4616-540c-a9ee-88d1112d94c9"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LogExpFunctions = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
@@ -24,7 +23,7 @@ SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 StatsFuns = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
-WilliamsonTransforms = "48feb556-9bdd-43a2-8e10-96100ec25e22"
+TaylorSeries = "6aa5eb33-94cf-58f4-a9d0-e4b2c4fc25ea"
 
 [weakdeps]
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
@@ -59,9 +58,9 @@ StableRNGs = "1"
 Statistics = "1"
 StatsBase = "0.33, 0.34"
 StatsFuns = "0.9, 1.3"
+TaylorSeries = "0.20"
 Test = "1"
 TestItemRunner = "1"
-WilliamsonTransforms = "0.2"
 julia = "1"
 
 [extras]

diff --git a/README.md b/README.md
@@ -11,7 +11,6 @@
      <a href="https://github.com/lrnv/Copulas.jl/actions/workflows/CI.yml?query=branch%3Amain"><img src="https://github.com/lrnv/Copulas.jl/actions/workflows/CI.yml/badge.svg?branch=main" alt="Build Status" /></a>
      <a href="https://codecov.io/gh/lrnv/Copulas.jl"><img src="https://codecov.io/gh/lrnv/Copulas.jl/branch/main/graph/badge.svg"/></a>
      <a href="https://github.com/JuliaTesting/Aqua.jl"><img src="https://raw.githubusercontent.com/JuliaTesting/Aqua.jl/master/badge.svg" alt="Aqua QA" /></a>
-    <!-- <a href="https://benchmark.tansongchen.com/TaylorDiff.jl"><img src="https://img.shields.io/buildkite/2c801728055463e7c8baeeb3cc187b964587235a49b3ed39ab/main.svg?label=benchmark" alt="Benchmark Status" /></a> -->
 <br />
     <a href="https://opensource.org/licenses/MIT"><img src="https://img.shields.io/badge/license-MIT-blue.svg" alt="License: MIT" /></a>
     <a href="https://github.com/SciML/ColPrac"><img src="https://img.shields.io/badge/contributor's%20guide-ColPrac-blueviolet" alt="ColPrac: Contributor's Guide on Collaborative Practices for Community Packages" /></a>

diff --git a/docs/src/bestiary/archimedean.md b/docs/src/bestiary/archimedean.md
@@ -96,7 +96,7 @@ In this package, we implemented it through the [`WilliamsonGenerator`](@ref) cla
 
 `WilliamsonGenerator(X::UnivariateRandomVariable, d)`.
 
-This function computes the Williamson d-transform of the provided random variable $X$ using the [`WilliamsonTransforms.jl`](https://github.com/lrnv/WilliamsonTransforms.jl) package. See [williamson1955multiply, mcneil2009](@cite) for the literature. 
+This function computes the Williamson d-transform of the provided random variable $X$. See [williamson1955multiply, mcneil2009](@cite) for the literature. 
 
 !!! info "`max_monotony` of Williamson generators"
     The $d$-transform of a positive random variable is $d$-monotone but not $k$-monotone for any $k > d$. Its max monotony is therefore $d$. This has a few implications, one of the biggest is that the $d$-variate Archimedean copula that corresponds has no density.
@@ -115,15 +115,15 @@ The Williamson d-transform is a bijective transformation[^1] from the set of pos
 
     This bijection is to be taken carefuly: the bijection is between random variables *with unit scales* and generators *with common value at 1*, sicne on both rescaling does not change the underlying copula. 
 
-This transformation is implemented through one method in the Generator interface that is worth talking a bit about : `williamson_dist(G::Generator, d)`. This function computes the inverse Williamson d-transform of the d-monotone archimedean generator ϕ, still using the [`WilliamsonTransforms.jl`](https://github.com/lrnv/WilliamsonTransforms.jl) package. See [williamson1955multiply, mcneil2009](@cite).
+This transformation is implemented through one method in the Generator interface that is worth talking a bit about : `williamson_dist(G::Generator, d)`. This function computes the inverse Williamson d-transform of the d-monotone archimedean generator ϕ. See [williamson1955multiply, mcneil2009](@cite).
 
 To put it in a nutshell, for ``\phi`` a ``d``-monotone archimedean generator, the inverse Williamson-d-transform of ``\\phi`` is the cumulative distribution function ``F`` of a non-negative random variable ``R``, defined by : 
 
 ```math
 F(x) = 𝒲_{d}^{-1}(\phi)(x) = 1 - \frac{(-x)^{d-1} \phi_+^{(d-1)}(x)}{k!} - \sum_{k=0}^{d-2} \frac{(-x)^k \phi^{(k)}(x)}{k!}
 ```
 
-The [`WilliamsonTransforms.jl`](https://github.com/lrnv/WilliamsonTransforms.jl) package implements this transformation (and its inverse, the Williamson d-transfrom) in all generality. It returns this cumulative distribution function in the form of the corresponding random variable `<:Distributions.ContinuousUnivariateDistribution` from `Distributions.jl`. You may then compute : 
+It returns this cumulative distribution function in the form of the corresponding random variable `<:Distributions.ContinuousUnivariateDistribution` from `Distributions.jl`. You may then compute : 
 * The cdf via `Distributions.cdf`
 * The pdf via `Distributions.pdf` and the logpdf via `Distributions.logpdf`
 * Samples from the distribution via `rand(X,n)`.
@@ -184,14 +184,14 @@ This is why `williamson_dist(G::Generator,d)` is such an important function in t
 
 ```@example
 using Copulas: williamson_dist, FrankGenerator
-williamson_dist(FrankGenerator(7), Val{3}())
+williamson_dist(FrankGenerator(7), 3)
 ```
 
 For the Frank Copula, as for many classic copulas, the distribution used is known. We pull some of them from `Distributions.jl` but implement a few more, as this Logarithmic one. Another useful example are negatively-dependent Clayton copulas: 
 
 ```@example
 using Copulas: williamson_dist, ClaytonGenerator
-williamson_dist(ClaytonGenerator(-0.2), Val{3}())
+williamson_dist(ClaytonGenerator(-0.2), 3)
 ```
 
 for which the corresponding distribution is known but has no particular name, thus we implemented it under the `ClaytonWilliamsonDistribution` name.
@@ -209,7 +209,7 @@ for which the corresponding distribution is known but has no particular name, th
     We use this fraily approach for several generators, since sometimes it is faster, including e.g. the Clayton one with positive dependence:
     ```@example
     using Copulas: williamson_dist, ClaytonGenerator
-    williamson_dist(ClaytonGenerator(10), Val{3}())
+    williamson_dist(ClaytonGenerator(10), 3)
     ```
 
 

diff --git a/docs/src/bestiary/empirical.md b/docs/src/bestiary/empirical.md
@@ -216,7 +216,7 @@ Usage:
 
 - Build from data `u::d×n` (raw or pseudos): `Ĝ = EmpiricalGenerator(u; pseudo_values=true)`
 - Use directly in an Archimedean copula: `Ĉ = ArchimedeanCopula(d, Ĝ)`
-- Access the fitted radial law: `R̂ = williamson_dist(Ĝ, Val{d}())`
+- Access the fitted radial law: `R̂ = williamson_dist(Ĝ, d)`
 
 ```@docs; canonical=false
 EmpiricalGenerator

diff --git a/docs/src/examples/archimedean_radial_estimation.md b/docs/src/examples/archimedean_radial_estimation.md
@@ -205,7 +205,7 @@ R = Dirac(1.0)
 d, n = 3, 1000
 u = spl_cop(R, d, n)
 Ghat = EmpiricalGenerator(u)
-Rhat = Copulas.williamson_dist(Ghat, Val{d}())
+Rhat = Copulas.williamson_dist(Ghat, d)
 diagnose_plots(u, Rhat; R=R)
 ```
 
@@ -216,7 +216,7 @@ R = DiscreteNonParametric([1.0, 4.0, 8.0], fill(1/3, 3))
 d, n = 2, 1000
 u = spl_cop(R, d, n)
 Ghat = EmpiricalGenerator(u)
-Rhat = Copulas.williamson_dist(Ghat, Val{d}())
+Rhat = Copulas.williamson_dist(Ghat, d)
 diagnose_plots(u, Rhat; R=R)
 ```
 
@@ -227,7 +227,7 @@ R = DiscreteNonParametric([1.0, 4.0, 8.0], fill(1/3, 3))
 d, n = 3, 1000
 u = spl_cop(R, d, n)
 Ghat = EmpiricalGenerator(u)
-Rhat = Copulas.williamson_dist(Ghat, Val{d}())
+Rhat = Copulas.williamson_dist(Ghat, d)
 diagnose_plots(u, Rhat; R=R)
 ```
 
@@ -238,7 +238,7 @@ R = LogNormal(1, 3)
 d, n = 10, 1000
 u = spl_cop(R, d, n)
 Ghat = EmpiricalGenerator(u)
-Rhat = Copulas.williamson_dist(Ghat, Val{d}())
+Rhat = Copulas.williamson_dist(Ghat, d)
 diagnose_plots(u, Rhat; R=R, logged=true)
 ```
 
@@ -249,7 +249,7 @@ R = Pareto(1.0, 1/2)
 d, n = 10, 1000
 u = spl_cop(R, d, n)
 Ghat = EmpiricalGenerator(u)
-Rhat = Copulas.williamson_dist(Ghat, Val{d}())
+Rhat = Copulas.williamson_dist(Ghat, d)
 diagnose_plots(u, Rhat; R=R, logged=true)
 ```
 

diff --git a/docs/src/manual/fitting_interface.md b/docs/src/manual/fitting_interface.md
@@ -39,7 +39,7 @@ plot(Ĉ)
 ### Full Model (with metadata)
 
 ```@example fitting_interface
-M = fit(CopulaModel, GumbelCopula, U; method=:default, summaries=true)
+M = fit(CopulaModel, GumbelCopula, U; method=:default)
 M
 ```
 
@@ -48,15 +48,15 @@ Returns a `CopulaModel` with:
 * `result` (the fitted copula), `n`, `ll` (log-likelihood),
 * `method`, `converged`, `iterations`, `elapsed_sec`,
 * `vcov` (if available),
-* `method_details` (a named tuple with method metadata and, if `summaries=true`, **pairwise summaries**: means, deviations, minima, and maxima of empirical τ/ρ/β/γ).
+* `method_details` (a named tuple with method metadata).
 
 ---
 
 ## Behavior & conventions (important)
 
 - fit operates on types, not on pre-constructed parameterized instances. Always pass a Copula or SklarDist *type* to `fit`, e.g. `fit(GumbelCopula, U)` or `fit(CopulaModel, SklarDist{ClaytonCopula,Tuple{Normal,LogNormal}}, X)`. If you already have a constructed instance `C0`, re-estimate its parameters by calling `fit(typeof(C0), U)`.
 
-- Default method selection: each family exposes the list of available fitting strategies via `_available_fitting_methods(CT)`. When `method = :default` the first element of that tuple is used. Example: `Copulas._available_fitting_methods(MyCopula)`.
+- Default method selection: each family exposes the list of available fitting strategies via `_available_fitting_methods(CT, d)`. When `method = :default` the first element of that tuple is used. Example: `Copulas._available_fitting_methods(MyCopula, d)`.
 
 - `CopulaModel` is the full result object returned by the fits performed via `Distributions.fit(::Type{CopulaModel}, ...)`. The light-weight shortcut `fit(MyCopula, U)` returns only a copula instance; use `fit(CopulaModel, ...)` to get diagnostics and metadata.
 
@@ -80,26 +80,7 @@ The `CopulaModel{CT} <: StatsBase.StatisticalModel` type stores the result and s
 | `hqc(M)`           | Hannan–Quinn criterion          |
 
 
-Quick access to the contained copula: `_copula_of(M)` (returns the copula even if `result` is a `SklarDist`).
-
-
-### Pairwise summaries and `method_details`
-
-When you request `summaries=true` (default) the returned `CopulaModel` contains extra pre-computed pairwise statistics inside `M.method_details`. Typical keys are:
-
-- `:tau_mean`, `:tau_sd`, `:tau_min`, `:tau_max`
-- `:rho_mean`, `:rho_sd`, `:rho_min`, `:rho_max`
-- `:beta_mean`, `:beta_sd`, `:beta_min`, `:beta_max`
-- `:gamma_mean`, `:gamma_sd`, `:gamma_min`, `:gamma_max`
-
-Access example:
-
-```@example fitting_interface
-M = fit(CopulaModel, GumbelCopula, U; summaries=true)
-M.method_details.tau_mean  # average pairwise Kendall's tau
-```
-
-If `summaries=false` these keys will be absent and `method_details` will be smaller.
+By default, the returned `CopulaModel` contains a lot of extra statistics, that you can see by printing the model in the REPL. 
 
 ### `vcov` and inference notes
 
@@ -143,7 +124,7 @@ plot(Ŝ.result)
 The names and availiability of fitting methods depends on the model. You can check what is available with the following internal call : 
 
 ```@example fitting_interface
-Copulas._available_fitting_methods(ClaytonCopula)
+Copulas._available_fitting_methods(ClaytonCopula, 3)
 ```
 
 The first method in the list is the one used by default. 
@@ -168,7 +149,7 @@ When you add a new copula family, implement the following so the generic `fit` f
 1. `_example(CT, d)` — return a representative instance (used to obtain default params and initial values).
 2. `_unbound_params(CT, d, params)` — transform the family `NamedTuple` parameters to an unconstrained `Vector{Float64}` used by optimizers.
 3. `_rebound_params(CT, d, α)` — invert `_unbound_params`, returning a `NamedTuple` suitable for `CT(d, ...)` construction.
-4. `_available_fitting_methods(::Type{<:YourCopula})` — declare supported methods (examples:  `:mle, :itau, :irho, :ibeta, ...`).
+4. `_available_fitting_methods(::Type{<:YourCopula}, d::Int)` — declare supported methods (examples:  `:mle, :itau, :irho, :ibeta, ...`).
 5. `_fit(::Type{<:YourCopula}, U, ::Val{:mle})` (and other `Val{}` methods) — implement the method and return `(fitted_copula, meta::NamedTuple)`; include keys such as `:θ̂`, `:optimizer`, `:converged`, `:iterations` and optionally `:vcov`.
 
 Place this checklist and a minimal `_fit` skeleton in `docs/src/manual/developer_fitting.md` where contributors can copy/paste and adapt.

diff --git a/docs/src/manual/intro.md b/docs/src/manual/intro.md
@@ -61,9 +61,9 @@ You may define a copula object in Julia by simply calling its constructor:
 
 ```@example 1
 using Copulas
-d = 4 # The dimension of the model
+d = 3 # The dimension of the model
 θ = 7 # Parameter
-C = ClaytonCopula(4,7) # A 4-dimensional clayton copula with parameter θ = 7.
+C = ClaytonCopula(d,7) # A 3-dimensional clayton copula with parameter θ = 7.
 ```
 
 This object is a random vector, and behaves exactly as you would expect a random vector from `Distributions.jl` to behave: you may sample it with `rand(C,100)`, compute its pdf or cdf with `pdf(C,x)` and `cdf(C,x)`, etc:
@@ -84,6 +84,18 @@ plot(C, :logpdf)
 
 See [the visualizations page](@ref viz_page) for details on the visualisations tools. It’s often useful to get an intuition by looking at scatter plots.
 
+!!! example "Independence"
+    To give another example, the function 
+
+    $\Pi : \boldsymbol x \mapsto \prod_{i=1}^d x_i = \boldsymbol x^{\boldsymbol 1}$ is a copula, corresponding to independent random vectors.
+
+    This copula can be constructed using the [`IndependentCopula(d)`](@ref IndependentCopula) syntax as follows: 
+
+    ```@example 1
+    Π = IndependentCopula(d) # A 4-variate independence structure.
+    nothing # hide
+    ```
+
 One of the reasons that makes copulas so useful is the bijective map from the Sklar Theorem [sklar1959](@cite):
 
 !!! theorem "Sklar (1959)"
@@ -95,25 +107,14 @@ One of the reasons that makes copulas so useful is the bijective map from the Sk
 
 This result allows to decompose the distribution of $\boldsymbol X$ into several components: the marginal distributions on one side, and the copula on the other side, which governs the dependence structure between the marginals. This object is central in our work, and therefore deserves a moment of attention. 
 
-!!! example "Independence"
-    The function 
-
-    $\Pi : \boldsymbol x \mapsto \prod_{i=1}^d x_i = \boldsymbol x^{\boldsymbol 1}$ is a copula, corresponding to independent random vectors.
 
-The independence copula can be constructed using the [`IndependentCopula(d)`](@ref IndependentCopula) syntax as follows: 
-
-```@example 1
-Π = IndependentCopula(d) # A 4-variate independence structure.
-nothing # hide
-```
 
 We can then leverage the Sklar theorem to construct multivariate random vectors from a copula-marginals specification. The implementation we have of this theorem allows building multivariate distributions by specifying separately their marginals and dependence structures as follows:
 
 
 ```@example 1
 X₁, X₂, X₃ = Gamma(2,3), Pareto(), LogNormal(0,1) # Marginals
-C = ClaytonCopula(3,0.7) # A 3-variate Clayton Copula with θ = 0.7
-D = SklarDist(C, (X₁,X₂,X₃)) # The final distribution
+D = SklarDist(C, (X₁,X₂,X₃)) # The final distribution, using the previous copula C. 
 plot(D, scale=:sklar)
 nothing # hide
 ```
@@ -132,7 +133,7 @@ Sklar's theorem can be used the other way around (from the marginal space to the
 
 !!! info "Independent random vectors"
 
-    `Distributions.jl` provides the [`product_distribution`](https://juliastats.org/Distributions.jl/stable/multivariate/#Product-distributions) function to create independent random vectors with given marginals. `product_distribution(args...)` is essentially equivalent to `SklarDist(Π, args)`, but our approach generalizes to other dependence structures.
+    `Distributions.jl` provides the [`product_distribution`](https://juliastats.org/Distributions.jl/stable/multivariate/#Product-distributions) function to create independent random vectors with given marginals. `product_distribution(args...)` is essentially equivalent to `SklarDist(IndependentCopula(d), args)`, but our approach generalizes to other dependence structures.
 
 Copulas are bounded functions with values in [0,1] since they correspond to probabilities. But their range can be bounded more precisely, and [lux2017](@cite) gives us:
 
@@ -309,3 +310,5 @@ The documentation of this package aims to combine theoretical information and re
 Pages = [@__FILE__]
 Canonical = false
 ```
+
+
diff --git a/src/ArchimaxCopula.jl b/src/ArchimaxCopula.jl
@@ -155,7 +155,7 @@ function _rebound_params(CT::Type{<:ArchimaxCopula{2, <:Generator, <:Tail}}, d,
     NamedTuple{all_names}(all_vals)
 end
 
-_available_fitting_methods(::Type{<:ArchimaxCopula}) = (:mle,)
+_available_fitting_methods(::Type{<:ArchimaxCopula}, d) = (:mle,)
 
 # Fast conditional distortion binding (bivariate)
 DistortionFromCop(C::ArchimaxCopula{2}, js::NTuple{1,Int}, uⱼₛ::NTuple{1,Float64}, ::Int) = BivArchimaxDistortion(C.gen, C.tail, Int8(js[1]), float(uⱼₛ[1]))
@@ -177,7 +177,7 @@ end
 
 # --- log-PDF stable ---
 function Distributions._logpdf(C::ArchimaxCopula{2, TG, TT}, u) where {TG, TT}
-    T = promote_type(Float64, eltype(u))
+    T = typeof(A(C.tail, one(ϕ(C.gen, one(eltype(u))))/2))
     @assert length(u) == 2
     u1, u2 = u
     (0.0 < u1 ≤ 1.0 && 0.0 < u2 ≤ 1.0) || return T(-Inf)
@@ -201,7 +201,7 @@ function Distributions._logpdf(C::ArchimaxCopula{2, TG, TT}, u) where {TG, TT}
 
     s    = S * A0
     φp   = ϕ⁽¹⁾(C.gen, s)            # < 0
-    φpp  = ϕ⁽ᵏ⁾(C.gen, Val(2), s)            # > 0
+    φpp  = ϕ⁽ᵏ⁾(C.gen, 2, s)            # > 0
 
     base = su*sv + (φp/φpp)*suv
     base > 0 || return T(-Inf)