Vine copula models

Vine copula models are implemented in the class Vinecop. To use this class in your code, include the header vinecopulib/vinecop/class.hpp (or simply vinecopulib.hpp) at the top of your source file. This automatically enables all features for bivariate copula models.

Conditional simulation is covered in Conditional simulation, discrete variables in Discrete variables, and the score, gradient and Hessian of the log-likelihood in Scores, gradients and Hessians.

How to read the R-vine structure

The structure is held by RVineStructure, which stores it as a triangular array: rows are trees, and each column encodes one edge per tree. For

4 4 4 4
3 3 3
2 2
1

the pair copulas are

treeedgepair copula
00(1, 4)
1(2, 4)
2(3, 4)
10(1, 3; 4)
1(2, 3; 4)
20(1, 2; 3, 4)

Writing M[i, j] for the entry in row i and column j, the edge e of tree t in a d-dimensional vine is (M[d - 1 - e, e], M[t, e]; M[t - 1, e], ..., M[0, e]): start at the counter-diagonal element of the column for the first conditioned variable, jump up to row t for the second, and collect everything above it as the conditioning set.

RVineStructure documents the full list of validity conditions, the natural-order normalization, and how truncation is represented.

RVineStructure has two order-determined specializations, CVineStructure and DVineStructure, and can also generate random structures.

  // an R-vine structure from a matrix (columns are trees, the diagonal holds
  // the conditioned variable of each column)
  Eigen::Matrix<size_t, Eigen::Dynamic, Eigen::Dynamic> mat(4, 4);
  mat << 4, 4, 4, 4, 3, 3, 3, 0, 2, 2, 0, 0, 1, 0, 0, 0;
  RVineStructure structure(mat);

  // C- and D-vines are determined by an order alone
  CVineStructure cvine({ 1, 2, 3, 4 });
  DVineStructure dvine({ 1, 2, 3, 4 });

  // a random structure on 5 variables, truncated after 2 trees
  auto random = RVineStructure::simulate(5, false, { 7 });

  std::cout << "dimension: " << structure.get_dim()
            << ", truncation level: " << structure.get_trunc_lvl() << std::endl;
  std::cout << structure.str() << std::endl;

Set up a custom vine copula model

A model is a structure plus a std::vector<std::vector<Bicop>> of pair copulas, shaped like the structure's triangular array. Vinecop::make_pair_copula_store() builds one of the right shape, filled with independence copulas.

As for bivariate copulas, there are two ways of setting up a vine:

  • with known parameters,
  • from data (i.e., with estimated parameters).

The constructors without data take a dimension (giving a D-vine of independence copulas), or a structure (or an R-vine matrix) together with the pair copulas.

  // a 4-dimensional D-vine
  DVineStructure structure({ 1, 2, 3, 4 });

  // all pair copulas independence by default; fill in the ones you want
  auto pair_copulas = Vinecop::make_pair_copula_store(4);
  for (auto& tree : pair_copulas) {
    for (auto& pc : tree) {
      pc = Bicop(BicopFamily::clayton, 0, Eigen::VectorXd::Constant(1, 2.0));
    }
  }

  Vinecop model(structure, pair_copulas);

Fit and select a vine copula model

select() performs parameter estimation and model selection using the sequential procedure of Dissmann et al. (2013). It selects the structure as well as the families, unless the object already carries a structure, in which case that structure is kept. fit() re-estimates families and parameters on a fixed structure. The only mandatory argument is the data; the controls argument customizes the fit.

  size_t d = 5;
  auto data = dependent_data(500, d, 1);

  // select the structure along with the families
  Vinecop model(data);

  // keep a known structure and select only the families and parameters
  Vinecop model2(data, DVineStructure({ 1, 2, 3, 4, 5 }));

  // refit an existing model's families and parameters, structure unchanged
  model2.fit(data);

  std::cout << "loglik: " << model.get_loglik() << ", bic: " << model.get_bic()
            << std::endl;

FitControlsVinecop inherits every member of FitControlsBicop (see Fit and select a bivariate copula) and adds the structure knobs:

  • size_t trunc_lvl — the number of trees that are actually fit. Pair copulas in trees above this level are set to independence. Defaults to no truncation.
  • std::string tree_criterion — the edge-weighting criterion used to grow each tree: "tau" (default) for Kendall's tau, "rho" for Spearman's rho, "hoeffd" for Hoeffding's D (suited to non-monotonic relationships), "mcor" for the maximum correlation, "cxi" for Chatterjee's xi symmetrized as $ \max\{\xi(X, Y), \xi(Y, X)\} $ (also suited to non-monotonic relationships), "joe" for the Gaussian-copula mutual information, or "custom" for a user-supplied callable set through set_tree_criterion_function(), which is always evaluated on the calling thread.
  • double threshold — pair copulas whose weight falls below this are set to independence (default 0, i.e. no thresholding).
  • std::string tree_algorithm — how the spanning tree is chosen: "mst_prim" (default) or "mst_kruskal" for maximum spanning trees, or "random_weighted" / "random_unweighted" for random ones.
  • bool select_trunc_lvl — select the truncation level automatically by mBICv (default false).
  • bool select_threshold — select the threshold automatically (default false).
  • bool select_families — whether families are selected at all, as opposed to keeping those already in place (default true).
  • std::vector<size_t> conditioning_set — variables to place at the tail of the order, so they can be conditioned on later (see Conditional simulation).
  • bool show_trace — print progress during fitting.
  • std::vector<int> seeds — seeds for the random tree algorithms.
  auto data = dependent_data(500, 5, 2);

  FitControlsVinecop controls;
  controls.set_family_set(bicop_families::itau);
  controls.set_parametric_method("itau");
  controls.set_trunc_lvl(2);           // force independence above tree 2
  controls.set_tree_criterion("mcor"); // edge weights for the spanning tree
  controls.set_threshold(0.05);        // weaker dependencies -> independence
  controls.set_tree_algorithm("mst_kruskal");
  controls.set_num_threads(2);

  Vinecop model(data, RVineStructure(), {}, controls);

  // or let the fit choose the truncation level and threshold by mBICv
  FitControlsVinecop selected;
  selected.set_select_trunc_lvl(true);
  selected.set_select_threshold(true);
  Vinecop sparse(data, RVineStructure(), {}, selected);
  std::cout << "selected truncation level: "
            << sparse.get_rvine_structure().get_trunc_lvl() << std::endl;

Work with a vine copula model

You can simulate from a model and evaluate its density, distribution function (by Monte Carlo), log-likelihood and information criteria, apply the Rosenblatt transform and its inverse, and truncate an existing model.

logpdf() is the log-density, and the accurate way to obtain it: a vine density is a product of one factor per edge, so in high dimensions or under strong dependence pdf() underflows to 0 while the log-density is still an ordinary double. loglik() sums it.

pdf_full() returns, alongside the density, the per-edge pair-copula densities and h-functions computed on the way; the derivative routines reuse them rather than recomputing the cascade.

  auto data = dependent_data(200, 4, 3);
  Vinecop model(data);

  // density, log-density, and the Monte-Carlo distribution function
  auto pdf = model.pdf(data);
  auto logpdf = model.logpdf(data);
  auto cdf = model.cdf(data, 1000);

  // pdf_full also returns the pair-copula densities and h-functions computed
  // along the way, which the derivative routines reuse
  auto full = model.pdf_full(data);

  // simulate, and the Rosenblatt transform with its inverse
  auto sim = model.simulate(100, false, 1, { 5 });
  auto z = model.rosenblatt(data);
  auto back = model.inverse_rosenblatt(z);

  // fit statistics
  std::cout << "loglik: " << model.get_loglik() << ", aic: " << model.get_aic()
            << ", bic: " << model.get_bic()
            << ", mbicv: " << model.get_mbicv(0.9) << std::endl;

  // a truncated copy
  auto truncated = model;
  truncated.truncate(1);

Inspecting the trees

get_trees() returns an RVineTrees: the structure as an explicit edge list, tree by tree. Each edge carries its conditioned pair (a, b), its conditioning set C, and its pair copula. RVineStructure can be constructed back from an RVineTrees, which makes this a round trip and a convenient way to build a structure edge by edge.

  auto data = dependent_data(300, 5, 4);
  Vinecop model(data);

  // the structure as an edge list, tree by tree: each edge carries the
  // conditioned pair (a, b), the conditioning set C, and its pair copula
  RVineTrees trees = model.get_trees();
  for (const auto& tree : trees.get_trees()) {
    for (const auto& edge : tree) {
      std::cout << edge.a << "," << edge.b;
      if (!edge.C.empty()) {
        std::cout << " | ";
        for (auto c : edge.C)
          std::cout << c << " ";
      }
      std::cout << "(" << edge.pair_copula.get_family_name() << ")  ";
    }
    std::cout << std::endl;
  }

  // and back again
  RVineStructure round_trip(trees);

Serialization

As for bivariate copulas, models can be written to and constructed from JSON or CBOR files and nlohmann::json objects. A .cbor extension selects the binary encoding; any other extension is JSON.

  auto data = dependent_data(200, 3, 6);
  Vinecop model(data);

  // JSON, via a file or an nlohmann::json object
  model.to_file(std::string("myvine.json"));
  Vinecop model2(std::string("myvine.json"));

  nlohmann::json j = model.to_json();
  Vinecop model3(j);

  // the .cbor extension selects binary CBOR automatically
  model.to_file(std::string("myvine.cbor"));
  Vinecop model4(std::string("myvine.cbor"));