Bivariate copula models

Implemented bivariate copula families

Each family's parameter ranges, admissible rotations and tail behavior are documented on the BicopFamily enum itself.

Several convenience vectors of families live in the sub-namespace bicop_families, and are the intended way to constrain a search space:

  // print all available families
  std::cout << "Available families: ";
  for (auto family : bicop_families::all) {
    std::cout << get_family_name(family) << " ";
  }
  std::cout << std::endl;

Set up a custom bivariate copula model

There are essentially two ways of setting up a bivariate copula:

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

The constructor with known parameters takes four arguments, all defaulted:

• the copula family (default indep),
• the rotation, one of 0, 90, 180, 270 (default 0),
• the parameters (default: those of an independence copula),
• the variable types (default {"c", "c"}, i.e. both continuous).

  // 90 degree rotated Clayton with default parameter (independence)
  Bicop clayton(BicopFamily::clayton, 90);

  // Gauss copula with parameter 0.5
  Bicop gauss(BicopFamily::gaussian, 0, Eigen::VectorXd::Constant(1, 0.5));

The constructor from data takes the same arguments as select() and is described in the next section.

Fit and select a bivariate copula

You can either fit the parameters of a given Bicop object with fit(), or select the best fitting model from a set of families with select().

  // create a Gauss copula with parameter 0.5 and simulate 1000 observations
  Bicop model(BicopFamily::gaussian, 0, Eigen::VectorXd::Constant(1, 0.5));
  auto data = model.simulate(1000);

  // instantiate a Gaussian copula with default parameters and fit to data
  Bicop fitted(BicopFamily::gaussian);
  fitted.fit(data);
  std::cout << "estimated parameter: " << fitted.get_parameters() << std::endl;

  // assign another family to the same variable and fit to data
  fitted = Bicop(BicopFamily::student);
  fitted.fit(data);
  std::cout << "estimated parameters: " << fitted.get_parameters() << std::endl;

  // alternatively, select the family automatically
  fitted.select(data);
  std::cout << "family: " << fitted.get_family_name()
            << ", rotation: " << fitted.get_rotation() << std::endl;

As it is arguably the most important function of the Bicop class, it is worth understanding the second argument of select(), an object of the class FitControlsBicop:

• std::vector<BicopFamily> family_set — the set of families to select from. Takes a user-specified vector or any of the vectors above (default bicop_families::all).
• std::string parametric_method — "mle" (default) for maximum likelihood, or "itau" for Kendall's inversion (only for families in bicop_families::itau).
• std::string nonparametric_method — the degree of the density approximation for the transformation kernel estimator: "constant" (default), "linear" or "quadratic".
• double nonparametric_mult — a factor multiplying the smoothing parameters.
• size_t nonparametric_grid_size — the size of the interpolation grid backing tll (default 30).
• std::string selection_criterion — "loglik", "aic" (default), "bic" or "mbic".
• Eigen::VectorXd weights — optional observation weights.
• double psi0 — the prior probability of non-independence used by "mbic" (default 0.9).
• bool preselect_families — a heuristic preselection (default true) based on symmetry properties of the data; the unrotated Clayton is not preselected when the data show upper-tail dependence, for instance.
• bool allow_rotations — whether rotations may be selected (default true).
• size_t num_threads — threads used when fitting several families.

As mentioned above, the arguments of select() can be given to a constructor, instantiating a fitted object directly:

  // select the "best" archimedean family by the default criterion (AIC), with
  // parameters from maximum likelihood
  Bicop best_archimedean(data, FitControlsBicop(bicop_families::archimedean));
  std::cout << "family: " << best_archimedean.get_family_name()
            << ", rotation: " << best_archimedean.get_rotation() << ", "
            << best_archimedean.get_parameters() << std::endl;

  // select by BIC, with parameters from Kendall's tau inversion
  FitControlsBicop controls(bicop_families::itau, "itau");
  controls.set_selection_criterion("bic");
  Bicop best_itau(data, controls);
  std::cout << "family: " << best_itau.get_family_name()
            << ", rotation: " << best_itau.get_rotation() << ", "
            << best_itau.get_parameters() << std::endl;

Work with a bivariate copula model

You can simulate from a model, and evaluate its density, distribution function, h-functions, inverse h-functions, log-likelihood and information criteria.

  // Gauss copula with parameter 0.5
  Bicop bicop(BicopFamily::gaussian, 0, Eigen::VectorXd::Constant(1, 0.5));

  // simulate 100 observations
  auto sim_data = bicop.simulate(100);

  // evaluate the pdf and cdf
  auto pdf = bicop.pdf(sim_data);
  auto cdf = bicop.cdf(sim_data);

  // evaluate the two h-functions and their inverses
  auto h1 = bicop.hfunc1(sim_data);
  auto h2 = bicop.hfunc2(sim_data);
  auto hi1 = bicop.hinv1(sim_data);
  auto hi2 = bicop.hinv2(sim_data);

  // evaluate the log-likelihood and the information criteria
  auto ll = bicop.loglik(sim_data);
  auto aic = bicop.aic(sim_data);
  auto bic = bicop.bic(sim_data);
  auto mbic = bicop.mbic(sim_data);
  std::cout << "loglik: " << ll << ", aic: " << aic << ", bic: " << bic
            << ", mbic: " << mbic << std::endl;

Dependence measures

get_tau() reports Kendall's of the current parameters, and parameters_to_tau() does the same for an arbitrary parametrization of the current family. The inverse, tau_to_parameters(), only exists where determines the parameters completely — see its documentation for the exact set of families.

get_beta() returns Blomqvist's beta, and get_taildep() a matrix whose and entries are the lower and upper tail-dependence coefficients (NaN for families with no closed form).

  Bicop bicop(BicopFamily::bb1, 0, (Eigen::VectorXd(2) << 1.0, 2.0).finished());

  // Kendall's tau of the current parameters
  double tau = bicop.get_tau();

  // Blomqvist's beta, and the 2 x 2 matrix of tail-dependence coefficients
  // (lower in (0, 0), upper in (1, 1))
  double beta = bicop.get_beta();
  Eigen::MatrixXd taildep = bicop.get_taildep();

  // The parameter <-> tau maps are only defined for the families in
  // bicop_families::itau, which BB1 is not.
  Bicop clayton(BicopFamily::clayton);
  Eigen::MatrixXd par = clayton.tau_to_parameters(0.5);
  double tau_again = clayton.parameters_to_tau(par);
  std::cout << "BB1 tau: " << tau << ", beta: " << beta
            << ", lower/upper tail dependence: " << taildep(0, 0) << "/"
            << taildep(1, 1) << "; Clayton tau round trip: " << tau_again
            << std::endl;

Per-observation parameters

For parametric families, pdf(), cdf(), hfunc1(), hfunc2(), hinv1(), hinv2(), loglik(), simulate() and the derivative methods also have overloads taking one parameter set per observation, as in conditional or covariate-dependent copula models.

  Bicop bicop(BicopFamily::gaussian, 0, Eigen::VectorXd::Constant(1, 0.5));

  // One parameter set per observation: simulate 100 draws, each from its own
  // Gauss copula, then evaluate at the same per-observation parameters.
  Eigen::MatrixXd parameters = Eigen::VectorXd::LinSpaced(100, 0.1, 0.8);
  auto sim = bicop.simulate(parameters);
  auto pdf = bicop.pdf(sim, parameters);
  auto loglik = bicop.loglik(sim, parameters);

  std::cout << "per-row draws: " << sim.rows() << " x " << sim.cols()
            << ", loglik: " << loglik << std::endl;

Serialization

Bivariate copula models can be written to and constructed from JSON or CBOR files and nlohmann::json objects. CBOR is a binary encoding of exactly the same logical representation: it is generally substantially smaller and faster to parse, while JSON remains useful for inspection and interchange. A .cbor extension selects CBOR; any other extension is JSON.

  // Gauss copula with parameter 0.5
  Bicop bicop(BicopFamily::gaussian, 0, Eigen::VectorXd::Constant(1, 0.5));

  // save as a nlohmann::json object
  nlohmann::json bicop_json = bicop.to_json();

  // write into a JSON file
  std::string filename = "myfile.json";
  std::ofstream file(filename);
  file << bicop_json << std::endl;
  file.close();

  // equivalently
  bicop.to_file(filename);

  // a new Bicop can be constructed from the nlohmann::json object
  Bicop bicop2(bicop_json);

  // or from the JSON file
  Bicop bicop3(std::string("myfile.json"));

  // the .cbor extension selects binary CBOR automatically
  bicop.to_file(std::string("myfile.cbor"));
  Bicop bicop4(std::string("myfile.cbor"));