(Vine) Copulas for Dummies

Statistical Methods

QSAY PhD Student India introduces Copulas, a focus of her PhD

Author

India Richmond

Published

Mar 2026

(Vine) Copulas for Dummies

What are copulas?

A copula is a multivariate distribution that enables us to separate the marginal distributions of the response variable from the dependency structure between them (Haugh, 2016). Copulas are used to describe the relationships between random variables. They also allow for the computation of joint probabilities without caring about or even knowing what their marginal distributions are (Vandenberghe et al., 2010). Regular correlation metrics such as Pearson’s correlation can struggle to explain complex relationships, as it only measures linear correlations (Embrechts et al., 2001). Overall, copulas allow us to create multivariate distributions where each of the variables could, in theory, follow a different distribution.

Skylar’s Theorem

Skylar’s Theorem is the key theorem in copula modelling, and it states that a multivariate distribution \(F\) can be written as a copula \(C\) applied to uniform marginal univariate distributions. Mathematically this is (Czado, 2019): \[F(x_1,...,x_d)=C(F_1(x_1),...,F_d(x_d))\] with associated density or probability mass function \[f(x_1,...,x_d)=c(F_1(x_1),...,F_d(x_d))f_1(x_1)...f_d(x_d)\] where copula density or pmf is found by \[c(u_1,...,u_d) = \frac{f(F_1^{-1}(u_1),...,F_d^{-1}(u_d))}{f(F_1^{-1}(u_1))...f(F_d^{-1}(u_d))}.\] When your marginals are continuous this decomposition is unique everywhere, whereas when your marginals are discrete this is only unique in the set \[\text{Range}(F_1)\times...\times\text{Range}(F_d).\]

Bivariate Copula

Most often copulas are used between just two variables, due to the fact that most copula are determined by just one or two parameters, these are known as bivariate copula. Which bivariate copula to use comes down to a range of properties that each bivariate copula possess.

The most common of which is tail dependence. Upper tail sees clustering in the top right hand corner, while lower tail sees clustering in the bottom left. Tail dependence can be understood as correlation getting stronger as we move further into a tail. Examples of some copula and the kinds of tail dependence they can exhibit are given below:

Copula	Upper Tail	Lower Tail
Gaussian	No	No
Student t	Yes	Yes
Gumbel	Yes	No
Clayton	No	Yes
Frank	No	No
Joe	Yes	No
BB1	Yes	Yes
BB7	Yes	Yes

Understanding the behaviour of different copula is important, as choosing the wrong copula for your data will result in a poorly fit bivariate distribution.

Bivariate Copula Examples

Gaussian Bivariate Copula \[C_{\text{Gauss}}(u,v;\rho) = \phi_2(\phi^{-1}(u), \phi^{-1}(v);\rho)\] Frank’s Bivariate Copula

\[C_{\text{Frank}}(u,v) = -\kappa^{-1}\log\left(1-\frac{(1-e^{-\kappa u})(1-e^{-\kappa v})}{1-e^{-\kappa}}\right)\]

Clayton’s Bivariate Copula

\[C_{\text{Clayton}}(u,v) = \max\left((u^{-\kappa}+v^{-\kappa}-1)^{-1/\kappa},0\right)\]

\(\rho\) denotes the copula parameter for the Gaussian copula, and \(\kappa\) denote copula parameters for Frank’s and Clayton’s copula.

Vine Copulas

Often we don’t just want to couple together two variables, what if we had say six, or ten, or more? Then we could turn to pair copula constructions/vine copulas. Vine copulas were developed by Bedford and Cooke (2001, 2002), and give us a way of veiwing the multivariate copula decompostion in terms of trees, which they named vine copulas. These trees describe a sequence (vine) of bivariate copula. Joe (1996) gave the first pair copula construction (PCC) of a multivariate copula in terms of distribution functions.

Pair Copula Decompostion

The decompostion of a multivariate density in terms of bivariate comes from recursive factorisation (Czado, 2019): \[f(x_1,x_2,x_3) = f_{3|12}(x_3|x_1,x_2)f_{2|1}(x_2|x_1)f_1(x_1)\]

Consider each of \(f_{3\|12}(x_3|x_1,x_2)\), \(f_{2|1}(x_2|x_1)\), \(f_1(x_1)\) seperately. For the first term, consider instead \(f_{13|2}(x_1,x_3|x_2)\).

Using Skylars theorem (associated density) we have: \[f_{13 \mid 2}\left(x_{1}, x_{3} \mid x_{2}\right)=c_{13 ; 2}\left(F_{1 \mid 2}\left(x_{1} \mid x_{2}\right), F_{3 \mid 2}\left(x_{3} \mid x_{2}\right) ; x_{2}\right) f_{1 \mid 2}\left(x_{1} \mid x_{2}\right) f_{3 \mid 2}\left(x_{3} \mid x_{2}\right) \space (*)\]

Note that this is equivelent to the bivariate case, just with each element conditioned on \(x_2\).

Consider the following: \(f_{1|2}(x_1\mid x_2) = c_{12}(F_1(x_1),F_2(x_2))f_2(x_2)\)

this is since \(f(x_1,x_2)=c_{12}(F_1(x_1),F_2(x_2))f_1(x_1)f_2(x_2)\), and therefore has conditional (on \(X_2\)) density \(f_{1\mid 2}(x_1 \mid x_2) = \frac{f(x_1,x_2)}{f_2(x_2)}\).

Following this same logic, to go from \(X_1\) and \(X_3\) conditional on \(X_2\) to \(X_3\) conditional on \(X_2\) and \(X_1\) we can apply this lemma be dividing by \(f_{1\mid 2}(x_1 \mid x_2)\). Then we have

\[f_{3 \mid 12}\left(x_{3} \mid x_{1}, x_{2}\right)=c_{13 ; 2}\left(F_{1 \mid 2}\left(x_{1} \mid x_{2}\right), F_{3 \mid 2}\left(x_{3} \mid x_{2}\right) ; x_{2}\right) f_{3 \mid 2}\left(x_{3} \mid x_{2}\right)\space \textbf{(1)}.\]

Applying the lemma directly again yeilds:

\[\begin{array}{l} f_{2 \mid 1}\left(x_{2} \mid x_{1}\right)=c_{12}\left(F_{1}\left(x_{1}\right), F_{2}\left(x_{2}\right)\right) f_{2}\left(x_{2}\right) \space \textbf{(2)}\\ f_{3 \mid 2}\left(x_{3} \mid x_{2}\right)=c_{23}\left(F_{2}\left(x_{2}\right), F_{3}\left(x_{3}\right)\right) f_{3}\left(x_{3}\right) \space \textbf{(3)} \end{array}\]

Subsitiuting in \((1)\), \((2)\) and \((3)\) into \((*)\) gives us the three dimensional pair copula decomposition:

\[\begin{aligned} f\left(x_{1}, x_{2}, x_{3}\right) & =c_{13 ; 2}\left(F_{1 \mid 2}\left(x_{1} \mid x_{2}\right), F_{3 \mid 2}\left(x_{3} \mid x_{2}\right) ; x_{2}\right) \times c_{23}\left(F_{2}\left(x_{2}\right), F_{3}\left(x_{3}\right)\right) \\ & \times c_{12}\left(F_{1}\left(x_{1}\right), F_{2}\left(x_{2}\right)\right) f_{3}\left(x_{3}\right) f_{2}\left(x_{2}\right) f_{1}\left(x_{1}\right) \end{aligned}\]

However, this decomposition is not unique! A reording of variables will result in a different vine decomposition. More on this later.

When modelling with vine copulas, and in most vine copula programming something known as the simplifying assumption is made. This assumes that the conditional bivariate copulas do not depend on the values of the conditioning variables:

\(C_{ij;D}(u_i, u_j \mid \textbf{u}_{D}) ≈ C_{ij;D}(u_i,u_j).\)

What this says is that we treat the dependence as constant, and that dependences does not change as \(u_D\) changes. Instead just the pseudo-observations account for the conditioning, found via h-functions.

h-functions of bivariate copula

Defined:

\[\begin{aligned} h_{1 \mid 2}\left(u_{1} \mid u_{2}\right) & :=\frac{\partial}{\partial u_{2}} C_{12}\left(u_{1}, u_{2}\right) \\ h_{2 \mid 1}\left(u_{2} \mid u_{1}\right) & :=\frac{\partial}{\partial u_{1}} C_{12}\left(u_{1}, u_{2}\right) . \end{aligned}\]

The exact form of the h-function will depend on which bivariate copula you are using. For example Gaussian Copula h-function \[h\left(u_{1} \mid u_{2}\right)=\Phi\left(\frac{\Phi^{-1}\left(u_{1}\right)-\rho \Phi^{-1}\left(u_{2}\right)}{\sqrt{1-\rho^{2}}}\right)\] where \(\rho\) is the copula parameter

Clayton Copula h-function \[h_{1 \mid 2}\left(u_{1} \mid u_{2}\right)=\left(u_{1}^{-\kappa}+u_{2}^{-\kappa}-1\right)^{-1-1 / \kappa} u_{2}^{-\kappa-1}\] where \(\kappa\) is the copula parameter.

h-functions are used to filter out condtional dependence, therefore they are used at each stage of the decompostion.

Viewing PCC’s as Vines

As previously stated, a pair copula decomposition is not unique, a result of this is that it becomes important (particularly in high dimensions) that we both allow for different decompositions and easily differentiate between them. This is where Bedford and Cookes Vine copula method comes in, giving us a way to view the decompositions as vines (or matrices).

tree1diag.jpg

The above diagram shows an undirected graph, with no cycles but fully connected, known as a tree. A vine copula is made up of \(d-1\) trees, where \(d\) is the dimensions of the data. A tree is defined by its nodes \(N\), which above are the set \(\{1,2,3,4,5\}\) and its edges \(E\), above \(\{(1,2), (1,3), (3,4), (4,5)\}\). A tree in a vine copula will always have one less edge than it has nodes, and each tree will reduce in size by 1 for each subsequence tree. Nodes in tree 1 simply index the variables. Edges are defined by the nodes they connect, in tree one if they connect nodes \(a\) and \(b\), then the edge is \(a,b\). It is these edges where a bivariate copula is “assigned”, for edge \((4,5)\) for example we assign a bivariate copula that will describe the dependency between variables 4 and 5, which will be used in the corresponding PCC. In tree 2 nodes are now labelled as the previous trees edges:

tree2diag.jpg

This time edges are determined by taking the union of the two nodes it connects, then removing the intersection and condtioning on this intersection instead. This continues until just two nodes and one edge remains. Czado (2019) defines this specfically:

Regular R- vine tree sequence The set of trees \(\mathbb{V}=\{T_1,...,T_{d-1}\}\) is a regular vine tree sequence on d elements if:

\[\begin{array}{l} \text{(1) Each tree } T_{j}=\left(N_{j}, E_{j}\right) \text{ is connected, i.e. for all nodes } a, b \in T_{j}, j=1, \ldots, d-1 \text{,}\\ \quad \text{ there exists a path } n_{1}, \ldots, n_{k} \subset N_{j} \text{ with } a=n_{1}, b=n_{k} \text{.}\\ \text{(2) } T_{1} \text{ is a tree with node set } N_{1}={1, \ldots, d} \text{ and edge set } E_{1} \text{.}\\ \text{(3) For } j \geq 2, T_{j} \text{ is a tree with node set } N_{j}=E_{j-1} \text{ and edge set } E_{j} \text{.}\\ \text{(4) For } j=2, \ldots, d-1 \text{ and } {a, b} \in E_{j} \text{ it must hold that } \|a \cap b\|=1 \text{.}\end{array}\]

Property \((4)\) is known as the proximity condition, and ensures that for an edge to connect two nodes, then the two nodes must share a common node in the previous tree. This ensures there is always an intersection between two nodes on which to condition on. The example given above is what is known as a D-vine, where the nodes are connected in a line. Once the first tree is defined this way, all deeper trees are defined - no more choices to be made. For other vine sequences, known as R-vines, this is not the case. There may be multiple possible edges. For six variables there exists over 20,000 possible PCC’s and therefore equivalently many tree decompostions - this is how we keep track of the possible vines.

A joint distrbution \(F\) is then defined by a regular vine distrbution \((\mathbb{F}, \mathbb{V}, \mathbb{B})\), where \(\mathbb{F}\) is the variables marginal distributions, \(\mathbb{V}\) an R-vine tree sequence and \(\mathbb{B}\) the set of bivariate copulas, where each edge in \(V\) is assocated with a bivariate copula.

Simple Example

Hopefully we have given enough of the theoretical background, now time for a simple example. First we will introduce the package that will be used.

pyvinecopulib is the python interface to vinecopulib, a header-only C++ library for vine copula models. We will make use of this library to fit some vine copulas.

!pip install pyvinecopulib
import pyvinecopulib as pv       # Install Vine Copula Package
from pyvinecopulib import BicopFamily
import numpy as np
import pandas as pd
from scipy.stats import norm
from scipy.stats import kendalltau
from scipy.stats import norm, poisson, multivariate_normal
import matplotlib.pyplot as plt # Import matplotlib.pyplot for plotting

We will now generate some correlated normal random variables data in three dimensions, \(X_1, X_2, X_3\), with their uniform counterparts denoted \(u_1, u_2, u_3\). \(X_1\), \(X_2\) and \(X_3\) are transformed onto the uniform distribution using the probability integral transform (PIT): \(u := F(x)\), this is equivalent to applying the cdf of a random variable.

np.random.seed(42)
n = 1000

# Correlated normal data
mean = np.zeros(3)
cov = [[1, 0.6, 0.5],
       [0.6, 1, 0.4],
       [0.5, 0.4, 1]]
X = np.random.multivariate_normal(mean, cov, size=n)

# probability integral transform to get uniform margins
# data must be uniform, according to Skylars theorem
U = norm.cdf(X)  # U is n x 3: each column is U_i = F_i(X_i)

We will fit just Gaussian Copulas, for the parameters we can use the following relation: \(\tau = \frac{2}{\pi}\arcsin(\rho)\) and use Kendall’s \(\tau\) between variables to determine the Gaussian copula parameter \(\rho\).

This relation can also be used when using student t copula. If instead you are using Clayton copula the relation would be \(\tau = \frac{\kappa}{\kappa +2}\) where \(\kappa\) would be the copula parameter. Different copula families have their own relations.

# Estimate Kendall's tau
tau_12 = kendalltau(U[:,0], U[:,1])[0] # Extract the correlation coefficients
tau_13 = kendalltau(U[:,0], U[:,2])[0]

# Convert to Gaussian copula parameter
def tau_to_rho(tau):
    return np.sin(np.pi / 2 * tau)

rho_12 = tau_to_rho(tau_12)
rho_13 = tau_to_rho(tau_13)

print(f"Estimated rho_12: {rho_12:.3f}, rho_13: {rho_13:.3f}")

Estimated rho_12: 0.601, rho_13: 0.469

# Create a bivariate Gaussian copula
gaussian_copula_12 = pv.Bicop(family=pv.BicopFamily.gaussian, parameters=np.array([[rho_12]]))
print(f"Fitted Gaussian Copula between u1 and u2: {gaussian_copula_12}")
# Create a bivariate Gaussian copula
gaussian_copula_13 = pv.Bicop(family=pv.BicopFamily.gaussian, parameters=np.array([[rho_13]]))
print(f"Fitted Gaussian Copula for u1 and u3: {gaussian_copula_13}")

Fitted Gaussian Copula between u1 and u2: <pyvinecopulib.Bicop> Bivariate copula: 
  family = Gaussian
  rotation = 0
  var_types = c,c
  parameters = 0.6

Fitted Gaussian Copula for u1 and u3: <pyvinecopulib.Bicop> Bivariate copula: 
  family = Gaussian
  rotation = 0
  var_types = c,c
  parameters = 0.47

In pyvinecopulib we can fit different copula families using their unique identifiers:

type	full name	string identifier
-	Independence	“indep”
Elliptical	Gaussian	“gaussian”
	Student t	“student”
Archimedean	Clayton	“clayton”
	Gumbel	“gumbel”
	Frank	“frank”
	Joe	“joe”
	Clayton-Gumbel (BB1)	“bb1”
	Joe-Gumbel (BB6)	“bb6”
	Joe-Clayton (BB7)	“bb7”
	Joe-Frank (BB8)	“bb8”
Extreme-Value	Tawn	“tawn”
Nonparametric	Transformation kernel	“tll”

We can also allow the package to determine the copula parameters for us:

U12 = np.column_stack((U[:,0], U[:,1]))
# Directly create a Bicop object with the Gaussian family
C12 = pv.Bicop(var_types=["c"] * 2,family=pv.BicopFamily.gaussian)
# Fit the Gaussian copula to the (uniform) data
C12.fit(U12)
print("Gaussian copula between u1 and u2:")
print(C12)

U13 = np.column_stack((U[:,0], U[:,2]))
# Directly create a Bicop object with the Gaussian family
C13 = pv.Bicop(var_types=["c"] * 2,family=pv.BicopFamily.gaussian)
# Fit the Gaussian copula to the (uniform) data
print("Gaussian copula between u1 and u3:")
C13.fit(U13)
print(C13)

Gaussian copula between u1 and u2:
<pyvinecopulib.Bicop> Bivariate copula: 
  family = Gaussian
  rotation = 0
  var_types = c,c
  parameters = 0.61

Gaussian copula between u1 and u3:
<pyvinecopulib.Bicop> Bivariate copula: 
  family = Gaussian
  rotation = 0
  var_types = c,c
  parameters = 0.49

Notice the estimated parameters are approximately the same.

Once these copulas have been fit we can plot them.

# Plot the fitted bivariate copula
C13.plot()

Next we use h-functions for Gaussian copula to transform the data, removing \(u_3\)’s and \(u_2\)’s dependence on \(u_1\)

# By hand (see form of gaussian h-function above)
def h_gaussian(u_j, u_i, rho):
    q_j = norm.ppf(u_j)
    q_i = norm.ppf(u_i)
    return norm.cdf((q_j - rho * q_i) / np.sqrt(1 - rho**2))

u2_given_1 = h_gaussian(U[:,1], U[:,0], rho_12)
u3_given_1 = h_gaussian(U[:,2], U[:,0], rho_13)

# use hfunc2 as we want to remove dependence on the second element
# use hfunc1 to remove dependence on the first
u2_given_1_package = C12.hfunc2(np.column_stack((U[:,1], U[:,0])))
u3_given_1_package = C13.hfunc2(np.column_stack((U[:,2], U[:,0])))

# You can compare the results from the manual calculation and the package
print("u2 | u1 (manual):", u2_given_1[:5])
print("u2 | u1 (package):", u2_given_1_package[:5])
print("u3 | u1 (manual):", u3_given_1[:5])
print("u3 | u1 (package):", u3_given_1_package[:5])

u2 | u1 (manual): [0.65697832 0.25546435 0.09449687 0.35259542 0.25691763]
u2 | u1 (package): [0.6625114  0.25866312 0.09512119 0.35202629 0.25044428]
u3 | u1 (manual): [0.48577661 0.18648372 0.41339318 0.21803681 0.00974834]
u3 | u1 (package): [0.49143995 0.19028045 0.42238898 0.21661331 0.00866923]

Small differences likely down to rounding and the slight difference in parameters (C12/13 used parameters estimated by package rather than tau association directly)

Then all thats left to do is fit the copula between \(u_2\) and \(u_3\) conditional on \(u_1\)

tau_23_given_1 = kendalltau(u2_given_1, u3_given_1)[0]
rho_23_given_1 = tau_to_rho(tau_23_given_1)

print(f"Estimated rho_23|1: {rho_23_given_1:.3f}")

Estimated rho_23|1: 0.134

# fit the final copula using package directly
U23g1 = np.column_stack((u2_given_1_package, u3_given_1_package))
# Directly create a Bicop object with the Gaussian family
C23g1 = pv.Bicop(var_types=["c"] * 2,family=pv.BicopFamily.gaussian)
# Fit the Gaussian copula to the (uniform) data
C23g1.fit(U23g1)
print(C23g1)

<pyvinecopulib.Bicop> Bivariate copula: 
  family = Gaussian
  rotation = 0
  var_types = c,c
  parameters = 0.15

To get the final density we need: \(c(u_1,u_2) * c(u_1, u_3) * c(u_2 \mid u_1, u_3 \mid u_1)\)

# .pdf evaluates the copula density
c12_density = C12.pdf(U12)
c13_density = C13.pdf(U13)
c23g1_density = C23g1.pdf(U23g1)

vine_density = c12_density * c13_density * c23g1_density
# sum over these values and take the log for the final log-likelihood
print(np.sum(np.log(vine_density)))

357.5165404356702

We can then check our by hand compuation with using the package directly:

# limit to only the Gaussian Copula
gaussian_family = [
    pv.BicopFamily.gaussian
]
chosen_controls = pv.FitControlsVinecop(
    family_set=gaussian_family)

# using the package to choose structure and parameters
fitted_vine = pv.Vinecop.from_data(U, var_types=["c"] * 3, controls=chosen_controls)
print("log-likelihood of the given vine when using the package to fit: ", fitted_vine.loglik(U))

log-likelihood of the given vine when using the package to fit:  357.51654043567135

We can see that the package and the by hand compuations agree.

This illustrates the steps taken to fit a vine copula in Python. In practice we don’t need to do all the steps by hand (as seen above). Using the package to do this for us is by far the simpler method. We proceed with another example, this time using simulated correlated discrete data.

np.random.seed(42)
n = 1000
d = 5
R = np.full((d, d), 0.5)
np.fill_diagonal(R, 1.0)

X_latent = multivariate_normal.rvs(mean=np.zeros(d), cov=R, size=n)
U_copula = norm.cdf(X_latent)
# Let's pick different means for each marginal
means = [1.5, 2.0, 2.5, 3.0, 3.5]
X_discrete = np.stack([poisson.ppf(U_copula[:, i], mu=means[i]).astype(int) for i in range(d)], axis=1)

# when working with discrete data, we also need the values evaluated at the left hand side (loc=1)
# So our data needs to be stacked together with the original uniform variables (all 5) followed by location shifted
U_pseudo = np.stack([poisson.cdf(X_discrete[:, i], mu=means[i]) for i in range(d)], axis=1)
U_pseudoLHS = np.stack([poisson.cdf(X_discrete[:, i], mu=means[i], loc=1) for i in range(d)], axis=1)

# result is a 10x1000 dataframe
uniform_discrete = np.column_stack((U_pseudo, U_pseudoLHS))

Why do we need the location shifted data too?

Vine Copula Modelling with discrete variables

Probability mass functions (pmfs) for discrete rv’s with CDF’s \(F_1,...,F_d\) come from rectangle probailities. In the bivariate case this is (Joe, 2014):

\[\begin{aligned} f(a_1, a_2) &= P(a_1-1< Y_1 \leq a_1, a_2 -1 < Y_2 \leq a_2)\\ &= F(a_1, a_2) - F(a_1-1, a_2) - F(a_1, a_2 -1)+ F(a_1 -1, a_2-1)\\ &= C(F_1(a_1),F_2(a_2)) - C(F_1(a_1-1), F_2(a_2)) \\ &- C(F_1(a_1),F_2(a_2-1)) + C(F_1(a_1-1), F_2(a_2 -1)) \end{aligned}\]

The conditional distributions are defined from conditional probabilities applied to events.

\[\begin{aligned} f_{k j \mid S}\left(x_{k}, x_{j} \mid \boldsymbol{x}_{S}\right)= & P\left(X_{k}=x_{k}, X_{j}=x_{j} \mid \boldsymbol{X}_{S}=\boldsymbol{x}_{S}\right) \\ = & C_{k j ; S}\left(F_{k \mid S}^{+}, F_{j \mid S}^{+} ; \boldsymbol{x}_{S}\right)-C_{k j ; S}\left(F_{k \mid S}^{-}, F_{j \mid S}^{+} ; \boldsymbol{x}_{S}\right) \\ & -C_{k j ; S}\left(F_{k \mid S}^{+}, F_{j \mid S}^{-} ; \boldsymbol{x}_{S}\right)+C_{k j ; S}\left(F_{k \mid S}^{-}, F_{j \mid S}^{-} ; \boldsymbol{x}_{S}\right) \end{aligned}\]

# for this example we will set no controls, this will be explored later
# we allow the package to choose the structure, bivariate copula and their parameters for itself
# notice we change from "c" to "d" for discrete variable
discrete_vine =  pv.Vinecop.from_data(uniform_discrete, var_types=["d"] * 5)
print(discrete_vine)

<pyvinecopulib.Vinecop> Vinecop model with 5 variables
tree edge conditioned variables conditioning variables var_types   family rotation   parameters   df  tau 
   1    1                  2, 4                             d, d Gaussian        0         0.51  1.0 0.34 
   1    2                  3, 4                             d, d      TLL          [30x30 grid]  9.9 0.31 
   1    3                  1, 5                             d, d Gaussian        0         0.50  1.0 0.33 
   1    4                  4, 5                             d, d      TLL          [30x30 grid] 11.4 0.30 
   2    1                  2, 3                      4      d, d      TLL          [30x30 grid] 15.9 0.16 
   2    2                  3, 5                      4      d, d      BB1        0   0.24, 1.16  2.0 0.23 
   2    3                  1, 4                      5      d, d Gaussian        0         0.32  1.0 0.21 
   3    1                  2, 5                   3, 4      d, d Gaussian        0         0.23  1.0 0.15 
   3    2                  3, 1                   5, 4      d, d Gaussian        0         0.26  1.0 0.17 
   4    1                  2, 1                5, 3, 4      d, d Gaussian        0         0.21  1.0 0.13

Setting no controls will cause much slower compuation times, as it checks all possible bivariate copula, including non-parametric ones which are far slower to fit

# see log-likelihood and aic of the fitted vine, useful for model comparrison
print("Total log-likelihood: ", discrete_vine.loglik(uniform_discrete))
print("AIC: ", discrete_vine.aic(uniform_discrete))

Total log-likelihood:  796.2077530438696
AIC:  -1501.8345076400144

Vine Structure

We can now explore these results

Starting with the structure, which we can be viewed as a series of trees:

discrete_vine.plot()  # optional, add labels by using argument vars_names =["A",...,"Z"]

The structure of the vine is typically chosen via a maximum spanning tree algorithm. Consider the following fully connected graph (note this is not the data used in the above example, this is some real football data collected from https://www.football-data.co.uk/:

Screenshot 2025-08-05 at 09.48.33.png

The weights on the edges between nodes denote the absolute values of Kendalls’ tau between those variables. The maximum spanning tree finds the set of edges which adhere to the vine copula rules (which we have already set out) that maximise the total possible (absolute value) of the correlations between variables. The edges marked in red represent the maximum spanning tree for the given connected graph, this then would be used as the first tree in the corresponding vine. This is repeated at deeper trees, after data has been transformed using h-functions. This is known as a top down approach, where the aim is to capture the most important (strongest) dependencies in the earlier trees.

Vine Copula Controls

If we wished to exhibit more control on the fit of the vine copula, we can use FitControlsVinecop. Using this allows us to limit which bivariate copula families can be used, the method for estimating copula parameters or the tree algorithm used, to name a few. We will set out an example of possible controls below and show their effect on the fitted vine.

# limit copula families
chosen_candidate_families = [
    pv.BicopFamily.gaussian,
    pv.BicopFamily.clayton,
    pv.BicopFamily.frank,
    pv.BicopFamily.indep
]
chosen_controls = pv.FitControlsVinecop(
    family_set=chosen_candidate_families,
    allow_rotations=False,         # do not allow rotations
    selection_criterion='loglik',  # change selection criteria to log-likelihood from aic (how bivariate copula are chosen)
    trunc_lvl = 3                  # truncates (sets all bivariate copula to independence) after 3rd tree
)

controls_vine =  pv.Vinecop.from_data(uniform_discrete, var_types=["d"] * 5, controls=chosen_controls)
print(controls_vine)

<pyvinecopulib.Vinecop> Vinecop model with 5 variables
tree edge conditioned variables conditioning variables var_types   family rotation parameters  df  tau 
   1    1                  2, 4                             d, d Gaussian        0       0.51 1.0 0.34 
   1    2                  3, 4                             d, d Gaussian        0       0.53 1.0 0.35 
   1    3                  1, 5                             d, d Gaussian        0       0.50 1.0 0.33 
   1    4                  4, 5                             d, d Gaussian        0       0.51 1.0 0.34 
   2    1                  2, 3                      4      d, d Gaussian        0       0.29 1.0 0.19 
   2    2                  3, 5                      4      d, d Gaussian        0       0.35 1.0 0.23 
   2    3                  1, 4                      5      d, d Gaussian        0       0.31 1.0 0.20 
   3    1                  2, 5                   3, 4      d, d Gaussian        0       0.23 1.0 0.15 
   3    2                  3, 1                   5, 4      d, d Gaussian        0       0.26 1.0 0.17

# see log-likelihood and aic of the fitted vine with controls
print("Total log-likelihood: ", controls_vine.loglik(uniform_discrete))
print("AIC: ", controls_vine.aic(uniform_discrete))

Total log-likelihood:  750.1378934915335
AIC:  -1482.275786983067

For full set of controls see https://vinecopulib.github.io/pyvinecopulib/_generate/pyvinecopulib.FitControlsVinecop.init.html

C/D- Vines

Above we have fit what is known as an R-vine. Within R-vines there exsists two subclasses of vines known as C-vines (canonical) and D-vines (drawable). C/D- vines can be fit by determining an order. For example, say we have the order \({1, 2, 3, 4, 5}\), for a C-vine the first tree in the vine connects variable 5 with all other variables, the second tree then connects variable 4 with all others, and so on. Whereas for a D-vine variables are connected in a sequence \(1-2-3-4-5\), all further trees are then decided due to the proximity condition. We can see examples of these:

# fit and plot a C-vine
cvineorder = pv.CVineStructure([1, 2, 3, 4, 5])  # determines C-vine order
CVine_fit = pv.Vinecop.from_data(uniform_discrete, var_types=["d"] * 5, structure = cvineorder, controls=chosen_controls)
CVine_fit.plot()

# fit and plot a D-vine
Dvineorder = pv.DVineStructure([1, 2, 3, 4, 5])  # determines D-vine order
DVine_fit = pv.Vinecop.from_data(uniform_discrete, var_types=["d"] * 5, structure = Dvineorder, controls=chosen_controls)
DVine_fit.plot()

We can also give the structure for any R-vine:

rand_Rvine_structure = pv.RVineStructure.simulate(d=5)   # this gives a random VALID R-vine structure
Rand_Rvine = pv.Vinecop.from_data(uniform_discrete, var_types=["d"] * 5, structure = rand_Rvine_structure, controls=chosen_controls)
Rand_Rvine.plot()

Note these vines have been truncated at the third level, resulting in just three trees (compared with five if we weren’t to truncate)

Vine copula simulations

We can also simulate from vine copulas. This works by doing the following:

Generate i.i.d. observations \(Z_1,...,Z_d ∼ \text{Uniform}(0,1)\)
Use the vine structure to apply the inverse h-function to each \(Z_i\) of the appropriate conditional copula. This transforms \(Z_i\) into the dependent variable \(U_i\), conditional on the previous ones.
After all transformations, the result is a matrix of \(n\times d\) dependent uniforms — consistent with the vine copula structure.

To put it simply, we work backwards, applying the inverse h-functions recursively to mimic the dependence captured by the vine copula. This could be a long process, but thankfully the package has us covered! Note the final simulations are still on the uniform scale, so in order to inspect them we need to apply the inverse CDF of their chosen distribution functions.

# can simulate from our fitted vine copula too, use seeds =[X] if need reproducable simulations
u_sims = controls_vine.simulate(10)   # number of simulations
print(u_sims)

[[0.16138634 0.43681965 0.15963652 0.31162637 0.88811142]
 [0.13300476 0.64368038 0.26510405 0.73835714 0.27482131]
 [0.20608251 0.50664383 0.19080473 0.61536064 0.58914311]
 [0.07239539 0.3286797  0.54564728 0.39091807 0.62388021]
 [0.99449737 0.75095777 0.57176362 0.62683536 0.66324885]
 [0.68761029 0.56641544 0.69888026 0.23104566 0.72399266]
 [0.97707278 0.97488406 0.94226797 0.99700456 0.96051241]
 [0.87749008 0.56329399 0.40452368 0.40300404 0.46622737]
 [0.51894709 0.28128173 0.4115985  0.42349014 0.23420049]
 [0.28688302 0.44279004 0.7390947  0.57051331 0.70053365]]

# they then just need transforming back into original scales
OS_sims1 = poisson.ppf(u_sims[:,0], means[0])
OS_sims2 = poisson.ppf(u_sims[:,1], means[1])
OS_sims3 = poisson.ppf(u_sims[:,2], means[2])
OS_sims4 = poisson.ppf(u_sims[:,3], means[3])
OS_sims5 = poisson.ppf(u_sims[:,4], means[4])
final_sims = np.column_stack((OS_sims1,OS_sims2,OS_sims3,OS_sims4,OS_sims5))
print(final_sims)

[[0. 2. 1. 2. 6.]
 [0. 2. 1. 4. 2.]
 [0. 2. 1. 3. 4.]
 [0. 1. 3. 2. 4.]
 [5. 3. 3. 3. 4.]
 [2. 2. 3. 2. 4.]
 [4. 5. 5. 9. 7.]
 [3. 2. 2. 2. 3.]
 [1. 1. 2. 3. 2.]
 [1. 2. 3. 3. 4.]]

Overall

We have learnt why vine copulas work, fit a simple vine copula from scratch, used pyvinecopulib to fit a vine for us and generated simulations!

References

Nelsen, R. B. (2006). An Introduction to Copulas. Springer. https://doi.org/10.1007/0-387-28678-0

Embrechts, P., Lindskog, F., & McNeil, A. (2001). Modelling dependence with copulas. Rapport technique, Département de mathématiques, Institut Fédéral de Technologie de Zurich, Zurich, 14:1–50.

Vandenberghe, S., Verhoest, N. E. C., & De Baets, B. (2010). Fitting bivariate copulas to the dependence structure between storm characteristics: A detailed analysis based on 105 year 10 min rainfall. Water Resources Research, 46(1). https://doi.org/10.1029/2009WR00785

Czado, C. (2019). Analyzing dependent data with vine copulas. Lecture Notes in Statistics, Springer, 222

Joe, H. (2014). Dependence modeling with copulas. CRC press.

Bedford, T., & Cooke, R. M. (2002). Vines–a new graphical model for dependent random variables. The Annals of statistics, 30(4), 1031-1068