Fitting and Selection

VineCopulas.jl uses sequential vine estimation. The package fits one tree at a time, computes conditional pseudo-observations through h-functions, and then uses those transformed observations to fit the next tree.

Sequential vine fitting

fit(CVineCopula, U), fit(DVineCopula, U), and fit(RVineCopula, U) return the fitted vine copula itself. method=:default and method=:sequential are equivalent; joint :mle fitting is not provided.

Quick fit

The simplest automatic R-vine fit is:

using VineCopulas
using Distributions: fit
using Random

truth = DVineCopula(
    [1, 2, 3],
    [[GaussianCopula(2, 0.55), ClaytonCopula(2, 1.4)],
     [FrankCopula(2, 2.0)]],
)

U = rand(MersenneTwister(7), truth, 250)
fitted = fit(RVineCopula, U)

(typeof(fitted), order(fitted), truncation(fitted))

(RVineCopula{3, 2, Tuple{Tuple{GaussianCopula{2, Float64}, SurvivalCopula{2, JoeCopula{2, Float64}}}, Tuple{FrankCopula{2, Float64}}}}, (1, 2, 3), 2)

The returned object is an ordinary RVineCopula: it can be evaluated, simulated, truncated, and passed to aic or bic.

Threaded fitting

Threading is opt-in. select_paircopula(U; threaded=true) fits candidate families and rotations concurrently, while fit(...; threaded=true) fits the independent edges of each vine tree concurrently. Both preserve the sequential statistical result and deterministic candidate ordering.

pair = select_paircopula(U[1:2, :]; threaded=true)
vine = fit(RVineCopula, U; threaded=true)

Benefits depend on the family set, sample size, vine dimension, truncation and available cores. Keep BLAS thread counts modest and avoid casually combining outer Threads.@threads loops with internally threaded fits. trace=true prioritizes readable deterministic output and may fit a tree's edges sequentially.

Statistical fitted models and inference

Use VineModel when the fitting recipe, observations, likelihood summaries, and coefficient vector are needed alongside the fitted distribution:

M = fit(VineModel, RVineCopula, U)

fitted_distribution(M)
fitting_method(M)       # :sequential
order(M); structure(M); truncation(M)
coef(M); coefnames(M)
StatsBase.aic(M); StatsBase.bic(M); mbicv(M)
edge_table(M)

edge_table(M) derives one row per active edge with its tree, edge index, conditioned and conditioning variables, family, rotation, and natural parameters. It does not store a second copy of the fitted structure.

StatsBase.aic(M) and StatsBase.bic(M) use the final vine likelihood and the number of active numerical pair-copula parameters. They are useful comparative scores for sequential fits to the same data, but the parameters are not a joint maximum-likelihood estimate; their usual joint-MLE asymptotic interpretation is therefore not asserted here.

VineCopulas.aicc(M) and VineCopulas.hqc(M) are currently internal convenience criteria computed from that same final likelihood and active numerical parameter count. They do not assert a special sequential-vine theory. mbicv(M) is the sparse-vine criterion (see Model comparison below); its psi0 defaults to the one the fit was called with, so a model fitted with trunc=:mbicv reports the score its truncation level was selected on.

Inference is a fixed-selection bootstrap:

I = infer(M; method=:bootstrap, nresamples=200, rng=MersenneTwister(1))
StatsBase.vcov(I)
StatsBase.stderror(I)
StatsBase.confint(I)

Conditional inference

Current vine covariance inference conditions on the selected truncation level, structure, pair families, and rotations. Every replicate resamples observations and sequentially refits all active pair parameters, propagating pseudo-observations through the trees. It does not include model-selection uncertainty, and it intentionally does not offer a joint-MLE Hessian.

Pair-family selection

Pair-copula selection chooses the best bivariate family for one edge. It can be used directly:

using VineCopulas
using Distributions: fit
using Random

C = ClaytonCopula(2, 1.8)
U2 = rand(MersenneTwister(9), C, 250)

pair = select_paircopula(
    U2;
    family_set=:default,
    selection_criterion=:bic,
    allow_rotations=true,
)

typeof(pair)

ClaytonCopula{2, Float64} (alias for ArchimedeanCopula{2, Copulas.ClaytonGenerator{Float64}})

The main controls are:

Keyword	Meaning
family_set	Candidate families, such as :default, :all, or an explicit tuple
selection_criterion	:loglik, :aic, or :bic
allow_rotations	Whether to include rotated candidates
include_independence	Whether independence can be selected
preselect	Whether dependence sign may prune rotation candidates
pair_method	Estimator per candidate: :default (maximum likelihood) or :itau
independence_test	:none, or :kendall to accept independence on an edge the asymptotic Kendall test does not reject
independence_level	Size of that test, 0.05 by default

Tip

Use selection_criterion=:bic for a conservative default. Use selection_criterion=:loglik mainly for diagnostics or controlled comparisons where model dimension is fixed.

pair_method=:itau fits each candidate by fit(FT, U; method=:itau) of Copulas.jl, which inverts the sample Kendall tau instead of maximizing the likelihood. It is closed form for the Gaussian, Clayton, Gumbel, Frank and Joe families, and a one-dimensional profile in the degrees of freedom for the Student family, so a sequential vine fit is several times faster than under :mle. Student fitting delegates to Copulas.jl for both methods: finite estimates use the full valid ν > 0 domain. If its profile returns ν = Inf, VineCopulas represents that Gaussian endpoint canonically as GaussianCopula for selection and scoring. The BB families have two parameters, which tau alone does not identify, so they are fitted by :mle under :itau, as in vinecopulib; trace=true prints the method used by each candidate. A family outside the shipped sets that advertises no :itau fails as a candidate with the Copulas.jl error, which strict=true raises. A candidate whose estimate gives the sample zero density (a Clayton with θ < 0, which :itau returns on a negative tau when preselect=false lets the unrotated candidate run) scores Inf and loses; it is not an error under strict=true, because the estimator did not fail.

independence_test=:kendall is the usual companion of :itau: an edge whose sample tau the test does not reject at independence_level receives the independence copula before any family is fitted, whatever include_independence says, like the threshold gate.

Structure selection

Automatic R-vine fitting follows a Dissmann-style sequential procedure:

1. compute dependence weights for candidate edges;

2. build a maximum spanning tree;

3. fit pair-copulas on selected edges;

4. propagate conditional pseudo-observations;

5. build the next candidate graph under the proximity condition;

6. repeat until the requested truncation depth.

using VineCopulas
using Distributions: fit
using Random

truth = DVineCopula(
    [1, 2, 3, 4],
    [[GaussianCopula(2, 0.5), ClaytonCopula(2, 1.3), FrankCopula(2, 2.0)],
     [GumbelCopula(2, 1.2), JoeCopula(2, 1.3)]];
    trunc=2,
)

U = rand(MersenneTwister(10), truth, 250)

model = fit(
    RVineCopula,
    U;
    trunc=2,
    tree_criterion=:tau,
    tree_algorithm=:kruskal,
    family_set=:default,
)

(order = order(model), truncation = truncation(model), edge_count = sum(length, edges(model)))

(order = (1, 3, 2, 4), truncation = 2, edge_count = 5)

At tree the engine selects

where is the set of proximity-admissible candidate edges, its spanning trees, and the weight tree_criterion attaches to the candidate edge    from the two conditional pseudo-observation vectors , . tree_criterion changes only; the candidate set is the same for every criterion. The deterministic tree_algorithm=:kruskal path is the one used in the external parity tests.

Built-in tree criteria

Each built-in criterion is the absolute value of a documented dependence statistic, listed here with its estimator, its range and its reference (full citations on the References page):

tree_criterion		range	reference
:tau	, Kendall's tau-b (default)		Kendall (1938); Dißmann et al. (2013)
:rho	, Spearman's rho on average ranks		Spearman (1904)
:hoeffd	, Hoeffding's in the Hollander & Wolfe form, scaled by 30		Hoeffding (1948); Hollander, Wolfe & Chicken (2014)
:mcor	, the maximum correlation, estimated by ACE		Gebelein (1941); Rényi (1959); Breiman & Friedman (1985)
:joe	, the Pearson correlation of the normal scores		Joe (1989)
:cxi	, the symmetrised Chatterjee coefficient		Chatterjee (2021)

Three notes on the less common ones.

• :hoeffd is Hoeffding's statistic; its population value is   , which is 0 under independence and 1 under any monotone functional relation. A non-monotone functional relation scores below 1 (a symmetric parabola gives about ). Ties are counted strictly. The estimator needs at least five observations.

• :mcor is the maximum correlation coefficient of Gebelein and Rényi, estimated by the alternating conditional expectations algorithm of Breiman & Friedman with a running-mean smoother of half-width , at most 10 inner and 100 outer iterations, and stopping tolerances and   on the change of the mean squared difference between the two transforms. These are the constants of vinecopulib's mcor, so the two agree on the same data. For a Gaussian pair  , and a functional relation, monotone or not, scores near 1.

• :joe is the mutual information of a Gaussian copula with correlation , Joe's relative-entropy dependence measure under the Gaussian assumption. It is not bounded, and it is strictly increasing in : a maximum spanning tree depends only on the order of its weights, so :joe selects the same trees as on normal scores, and its scale is visible only through threshold. It is Inf on a pair whose normal scores are exactly linearly dependent.

• :cxi is Chatterjee's , which measures how well one variable is a measurable function of the other and is therefore asymmetric; the criterion takes the larger of the two directions. The population value is in ; the estimate can be slightly negative, and the absolute value is taken.

The rank correlations and :joe do not see a symmetric non-monotone relation (a parabola in has  ), while :hoeffd, :mcor and :cxi do; the last two rank a functional relation above a noisy monotone one.

A custom tree criterion

tree_criterion also accepts a function. It is called once per candidate edge   as f(u_a, u_b, a, b, D), where u_a and u_b are the two conditional pseudo-observation vectors and a, b, D are variable labels, and its value is . The same function drives C-vine root choice and child edges (there D is the roots chosen so far), D-vine order selection, and the R-vine spanning trees. The built-in :tau is (ua, ub, a, b, D) -> abs(VineCopulas._kendall_tau_b(ua, ub)).

The contract:

• the value must be a finite Real; NaN, ±Inf, missing, or anything that is not a Real throws an ArgumentError naming the edge and the value;

• the value is used as is: no absolute value is taken, and the tree maximises it, so the sign is the caller's choice;

• the function should be deterministic in its arguments. Each edge is scored once per tree, so a random criterion is not refused, but the fit is then not reproducible without control of its randomness.

# A precomputed weight matrix for the first tree, |tau| above it.
W = my_tree1_weights(U)
model = fit(RVineCopula, U;
    tree_criterion=(ua, ub, a, b, D) -> isempty(D) ? W[a, b] : abs(VineCopulas._kendall_tau_b(ua, ub)))

# The Pearson correlation of the normal scores, the monotone image of :joe.
normal_r = (ua, ub, a, b, D) -> abs(cor(quantile.(Normal(), ua), quantile.(Normal(), ub)))
model = fit(RVineCopula, U; tree_criterion=normal_r)

threshold follows the criterion's scale

threshold forces an independence pair-copula on every candidate edge with   threshold, on the criterion's own scale, and it does not remove the edge from the candidate set. Its admissible range therefore follows the criterion: for :tau, :rho, :hoeffd, :mcor and :cxi; for :joe; any finite value for a function.

With a custom criterion threshold is a cut on that function's scale and no longer a dependence threshold in the usual sense. For example (ua, ub, a, b, D) -> abs(VineCopulas._kendall_tau_b(ua, ub)) + 2 leaves deliberately; with it, threshold=2.3 is exactly :tau with threshold=0.3, and threshold=0.3 forces nothing.

Group-constrained first tree

groups is an optional structural constraint for automatic R-vine selection. It takes one integer group id per variable:

model = fit(RVineCopula, U; groups=[1, 1, 1, 2, 2, 2])

Let the variables be partitioned into groups . Among the admissible first-tree spanning trees, groups restricts selection to those for which every induces a connected subtree. Under the edge weights supplied by tree_criterion, tree 1 solves

This is a clustered spanning-tree connectivity constraint. For the complete candidate graph available at R-vine tree 1, the constrained maximum spanning tree decomposes into maximum spanning trees within the groups and a maximum spanning tree on the contracted graph of groups.

The implementation obtains the same optimum in one deterministic Kruskal pass: within-group candidates are processed before cross-group candidates, while decreasing edge weight and the existing deterministic tie break are preserved inside each class.

Thus grouped tree-1 selection has the same asymptotic complexity as the unconstrained Kruskal selection. Its implementation is tested against brute force enumeration of all admissible spanning trees on small instances.

groups affects tree 1 only. From tree 2 onward, candidate edges are determined by the R-vine proximity condition and selection proceeds exactly as in the ordinary sequential fit.

Two limiting cases reproduce the unconstrained fit exactly:

groups = ones(Int, p)   # all variables in one group
groups = collect(1:p)   # one group per variable

groups composes with tree_criterion, threshold, trunc, and sampling_tail. It cannot be combined with structure, because a fixed structure already determines tree 1.

The partition should be interpreted as structural information supplied by the user rather than inferred from the data. Its usefulness therefore depends on whether the grouping is meaningful for the application.

The effect is particularly relevant for truncated vines, where the first trees determine which dependencies are represented explicitly. In the special full-depth Gaussian case, different R-vine structures can represent the same Gaussian copula; remaining differences are then due to finite-sample sequential estimation rather than to the group constraint itself.

Observation weights

Every fit method, select_paircopula and fit(VineModel, ...) take weights, a vector with one non-negative, finite entry per column of U (not all zero), or nothing for the unweighted fit:

using VineCopulas
using Distributions: fit
using Random

truth = DVineCopula(
    [1, 2, 3],
    [[GaussianCopula(2, 0.55), ClaytonCopula(2, 1.4)],
     [FrankCopula(2, 2.0)]],
)
U = rand(MersenneTwister(21), truth, 300)

# Exponentially decaying weights: recent columns count for more.
w = exp.(-(size(U, 2) .- (1:size(U, 2))) ./ 100)
weighted = fit(RVineCopula, U; weights=w)

(order(weighted), truncation(weighted))

((1, 2, 3), 2)

A weight is "how many observations this column counts for". The weights enter the fit in three places:

• the tree criterion is the weighted statistic: Kendall's tau (:tau) and Spearman's rho (:rho) with every pair (i, j) weighted by w[i] * w[j], and the Gaussian mutual information (:joe) of the weighted correlation of the normal scores;

• each pair-copula maximises the weighted pseudo-likelihood ∑ᵢ w[i] log c(U[:, i]);

• the selection criterion uses the weighted sample size ∑ᵢ w[i].

The weights are scaled internally so that they sum to the number of observations, so weights=fill(c, n) reproduces the unweighted fit for any c > 0, and the :bic penalty is unchanged. With integer weights the fit equals the unweighted fit of the sample in which column i is repeated w[i] times, up to optimiser tolerance. A zero weight removes the column. This is the contract of fit(::Type{<:Copula}, U; weights) in Copulas.jl, and the pair families that Copulas.jl fits (every family outside the default set, and every pair_method other than :mle) take the same vector through it.

The weights do not change the pseudo-observations: U is taken as given, and weighted ranks are the caller's responsibility.

fit(VineModel, VT, U; weights=w) stores the scaled weights in the recipe, and its loglikelihood, aic, bic and deviance are the weighted ones; nobs stays the number of columns. infer(M) resamples (observation, weight) pairs, so each bootstrap replicate is a weighted refit.

Criteria without a weighted form

:hoeffd, :mcor and :cxi have no weighted form here, so fit refuses them with weights instead of silently ignoring the weights. A custom tree_criterion function receives the full columns, aligned with the weights, and no weight argument: a weighted custom criterion closes over the weight vector.

Fixed-structure fitting

If the R-vine structure is part of the statistical design, pass it explicitly:

using VineCopulas
using Distributions: fit
using Random

ord = [3, 1, 4, 2]
S = ([1, 4, 2], [4, 2], [2])
st = RVineStructure(ord, S; trunc=3)

source = RVineCopula(
    st,
    [
        [GaussianCopula(2, 0.4), ClaytonCopula(2, 1.5), FrankCopula(2, 2.0)],
        [GaussianCopula(2, 0.3), ClaytonCopula(2, 1.2)],
        [FrankCopula(2, 1.2)],
    ],
)

U = rand(MersenneTwister(12), source, 250)

refit = fit(RVineCopula, U; structure=st, family_set=:default)
(order(refit), truncation(refit))

((3, 1, 4, 2), 3)

Fixed structure is not fixed family

Passing structure=st fixes the vine graph/order. It does not freeze the pair-copula families unless the API path you use also supplies explicit edges or a previously constructed vine.

Truncation

trunc=q limits fitting to the first q trees:

fit(RVineCopula, U; trunc=2)

For an existing model, use:

smaller = truncate(vine, 2)

Truncation returns a new model or structure and does not mutate the original. It cannot restore pair-copulas that are not present in the input.

Selecting the truncation level

trunc=:mbicv selects the level tree by tree with the modified BIC for vines of Nagler, Bumann and Czado (2019):

fit(RVineCopula, U; trunc=:mbicv)                 # search up to p - 1 trees
fit(RVineCopula, U; trunc=:mbicv, max_trunc=5)    # cap the search
fit(RVineCopula, U; trunc=:mbicv, psi0=0.8)       # a sparser prior

After each tree t is fitted, the model truncated at t is scored with mbicv. If the score does not improve on the score at t - 1, tree t is dropped and the fit stops; otherwise the search continues up to the ceiling, which is max_trunc or p - 1. The rule is greedy, so it costs at most one extra tree. It finds the first local improvement boundary along this sequential fit, not an exhaustive optimum over levels, structures, orders, or family assignments. psi0 ∈ (0, 1) is the prior probability that a tree-1 pair copula is not independence; the prior decays geometrically with the tree level, so deeper trees must earn their parameters. The published default is 0.9; a smaller value favours sparser models.

The selected level is truncation(fitted), and the returned vine is the one trunc=truncation(fitted) returns. trace=true prints the score of every level tried. The same keywords apply to CVineCopula, DVineCopula, and to an R-vine fitted on a fixed structure, where the ceiling is the structure's own truncation level. Through VineModel the selected level and its score are read off the model:

M = fit(VineModel, RVineCopula, U; trunc=:mbicv, psi0=0.8)
(truncation(M), mbicv(M))     # mbicv(M) reads psi0 = 0.8 from the recipe

Level zero

The mathematical case truncate(vine, 0) corresponds to multivariate independence. The current public truncation API starts at level 1 because the existing engines assume at least one active tree. Level-zero truncation is a roadmap item, not a hidden feature.

Model comparison

For explicit or fitted vine objects:

using VineCopulas
using Random

vine = DVineCopula(
    [1, 2, 3],
    [[GaussianCopula(2, 0.5), ClaytonCopula(2, 1.3)],
     [FrankCopula(2, 1.8)]],
)
U = rand(MersenneTwister(14), vine, 200)

(loglikelihood = loglikelihood(vine, U),
 npars = npars(vine),
 AIC = aic(vine, U),
 BIC = bic(vine, U),
 mBICV = mbicv(vine, U))

(loglikelihood = 95.27958396863679, npars = 3, AIC = -184.55916793727357, BIC = -174.66421583762946, mBICV = -173.82133171236686)

AIC, BIC, and mBICV are post-fit scores for the supplied model and data. They do not by themselves prove that the sequential fitting procedure found a global full-vine maximum likelihood optimum. mbicv adds to BIC a prior on the sparsity pattern: with p variables, qₜ non-independence pair copulas in tree t, and ψₜ = psi0^t,

Trees beyond the truncation level are independence, so they contribute (p - t) log(1 - ψₜ) each.

For weighted fits, mbicv uses VineCopulas' normalized weighted log-likelihood convention: weights sum to the number of observations, so the criterion retains that observation count in its BIC term. This is a package extension; the paper treats ordinary iid likelihoods.

Metadata and diagnostics

The fitted vine is deliberately a distribution, not a fitted-result container. Its order, truncation, selected edge families, and post-fit scores can be inspected directly, and VineModel carries the fitting recipe and the data beside it. The selected truncation level is truncation of the fitted vine or model, and mbicv scores it after the fact. Richer diagnostics such as selection traces and edge-level convergence details belong in VineModel rather than in the vine distribution object, and are planned separately.