Transforms and the Kendall function

rcopula.transforms

Conditional copula distributions and the Rosenblatt transform.

The central object is the h-function, the conditional distribution of one coordinate given the others:

.. math:: h(u_1 \mid u_2) = P(U_1 \le u_1 \mid U_2 = u_2) = \frac{\partial C(u_1, u_2)}{\partial u_2}.

It does a surprising amount of work. It is the sampler for every copula without a closed-form generator (draw w, invert h); it is the building block of vine copulas; it is the Rosenblatt transform that turns dependent data into independent uniforms for goodness-of-fit; and in pairs trading it is the signal -- :math:h(u_1 \mid u_2) near zero says asset 1 is unusually cheap given where asset 2 sits, which is a statement no correlation can make.

R exposes this as cCopula and supports Archimedean and elliptical families only. Analytic forms are used here for those, and numerical differentiation elsewhere, so every family in the package is covered.

References

Rosenblatt, M. (1952). Remarks on a multivariate transformation. Annals of Mathematical Statistics 23(3), 470-472. Aas, K., Czado, C., Frigessi, A. and Bakken, H. (2009). Pair-copula constructions of multiple dependence. Insurance: Mathematics and Economics 44(2), 182-198. The h-function as the workhorse of vine copulas. Joe, H. (2014). Dependence Modeling with Copulas. Chapman & Hall/CRC.

conditional_cdf

conditional_cdf(copula: Copula, u: ArrayLike, given: int = 1) -> NDArray[np.float64]

Probability that one variable is at most its value, given the other variable's value.

This is the conditional distribution of a bivariate copula, often called the h-function: :math:P(U_j \le u_j \mid U_{given} = u_{given}). It answers "given where variable B is, how unusual is variable A?" -- a value near 0 means A is unusually low for that B, near 1 unusually high. It is also the building block of sampling, vine copulas and the Rosenblatt transform.

Parameters:

Name	Type	Description	Default
copula	Copula	A bivariate (dim=2) copula with its parameters set. Archimedean, Gaussian and Student-t families use exact formulas; every other family uses numerical differentiation of its CDF.	required
u	array_like of float, shape (n, 2)	Points in the unit square, one per row. A single point may be passed as a 1-D sequence of length 2. Values are clipped to [1e-12, 1 - 1e-12].	required
given	(0, 1)	Index of the column to condition on. given=1 returns :math:h(u_1 \mid u_2) = \partial C/\partial u_2; given=0 returns the other conditional, :math:P(U_2 \le u_2 \mid U_1 = u_1).	0

Returns:

Type	Description
numpy.ndarray of float, shape (n,)	Values in [0, 1]. Under the true copula these are uniform, which is what makes them usable as a standardised signal.

Raises:

Type	Description
ValueError	If u does not have exactly two columns, or given is not 0 or 1.

Examples:

Under independence, conditioning changes nothing:

>>> import numpy as np
>>> from rcopula import IndependenceCopula
>>> from rcopula.transforms import conditional_cdf
>>> u = np.array([[0.3, 0.8], [0.7, 0.2]])
>>> conditional_cdf(IndependenceCopula(2), u)
array([0.3, 0.7])

Under positive dependence, a low value of the conditioning asset makes the other look high by comparison:

>>> from rcopula import GaussianCopula
>>> float(round(conditional_cdf(GaussianCopula(0.8), [[0.5, 0.05]])[0], 6))
0.985851

The output is uniform when the copula is right -- the property the Rosenblatt transform and every h-function-based test rely on:

>>> from scipy import stats
>>> from rcopula import ClaytonCopula
>>> cop = ClaytonCopula(3.0)
>>> h = conditional_cdf(cop, cop.rvs(20_000, random_state=0))
>>> bool(stats.kstest(h, "uniform").pvalue > 0.01)
True

conditional_ppf

conditional_ppf(copula: Copula, w: ArrayLike, cond: ArrayLike, given: int = 1) -> NDArray[np.float64]

Find the value of one variable that has a given conditional probability, given the other.

The inverse of :func:conditional_cdf (the inverse h-function): returns x with :math:h(x \mid \text{cond}) = w. Its main use is simulation -- draw w uniform, invert, and the pair (x, cond) follows the copula -- including conditional scenarios where cond is fixed.

Parameters:

Name	Type	Description	Default
copula	Copula	A bivariate (dim=2) copula with its parameters set.	required
w	float or array_like of float, shape (n,)	Target conditional probabilities, each in :math:[0, 1].	required
cond	float or array_like of float, shape (n,)	Values of the conditioning variable, each in :math:[0, 1]. w and cond are broadcast against each other, so either may be a scalar.	required
given	(0, 1)	Which coordinate cond occupies. given=1 solves for the first coordinate given the second; given=0 solves for the second given the first.	0

Returns:

Type	Description
numpy.ndarray of float, shape (n,)	The solved values x, each in (0, 1).

Raises:

Type	Description
ValueError	If any w or cond lies outside :math:[0, 1], if given is not 0 or 1, or if copula is not bivariate (raised by :func:conditional_cdf).

Notes

Solved by 60 steps of vectorised bisection on [1e-12, 1 - 1e-12], which handles every family uniformly and takes the bracket below double precision. Out-of-range targets are rejected rather than quietly converging to the nearest end of the bracket, which used to return a plausible-looking but meaningless value. NaN in w or cond gives NaN in that position.

Examples:

Round-trips against :func:conditional_cdf:

>>> import numpy as np
>>> from rcopula import ClaytonCopula
>>> from rcopula.transforms import conditional_cdf, conditional_ppf
>>> cop = ClaytonCopula(2.5)
>>> w = np.array([0.1, 0.5, 0.9])
>>> cond = np.array([0.3, 0.6, 0.8])
>>> x = conditional_ppf(cop, w, cond)
>>> recovered = conditional_cdf(cop, np.column_stack([x, cond]))
>>> bool(np.allclose(recovered, w, atol=1e-8))
True

rosenblatt

rosenblatt(copula: Copula, u: ArrayLike) -> NDArray[np.float64]

Turn dependent copula data into independent uniform numbers, if the copula is right.

This is the Rosenblatt transform (R's cCopula). Each column is replaced by its conditional probability given all the columns before it. If the data really came from copula, the output columns are independent and uniform on :math:[0, 1]; if not, they are visibly not. That turns a goodness-of-fit question into a simpler test of independence and uniformity.

Parameters:

Name	Type	Description	Default
copula	Copula	The copula to test against, with its parameters set. Archimedean, Gaussian and Student-t families work in any dimension; other families only for dim=2.	required
u	array_like of float, shape (n, d)	Observations on the copula scale (pseudo-observations), one row each. d must equal copula.dim. Values are clipped to [1e-12, 1 - 1e-12].	required

Returns:

Type	Description
numpy.ndarray of float, shape (n, d)	Transformed values in [0, 1]. The first column equals the (clipped) first column of u.

Raises:

Type	Description
ValueError	If the number of columns of u differs from copula.dim.
NotImplementedError	If copula is neither Archimedean nor elliptical and d > 2.

Notes

Maps dependent uniforms to independent uniforms by conditioning successively:

.. math:: Z_1 = U_1,\quad Z_2 = C(U_2 \mid U_1),\quad Z_3 = C(U_3 \mid U_1, U_2),\ \dots

If the copula is correct the result is a sample of independent standard uniforms, which turns any goodness-of-fit question into a test of independence and uniformity -- the basis of the SnB and SnC statistics.

Implemented analytically for Archimedean and elliptical families in any dimension, and for arbitrary bivariate copulas.

Examples:

The right copula gives independent uniforms; the wrong one does not:

>>> import numpy as np
>>> from scipy import stats
>>> from rcopula import ClaytonCopula
>>> from rcopula.transforms import rosenblatt
>>> cop = ClaytonCopula(3.0, dim=3)
>>> z = rosenblatt(cop, cop.rvs(5000, random_state=0))
>>> bool(all(stats.kstest(z[:, j], "uniform").pvalue > 0.01 for j in range(3)))
True
>>> bool(abs(np.corrcoef(z[:, 0], z[:, 1])[0, 1]) < 0.05)
True

Transformed with the wrong parameter, the output is visibly non-uniform:

>>> bad = rosenblatt(ClaytonCopula(0.2, dim=3), cop.rvs(5000, random_state=0))
>>> bool(stats.kstest(bad[:, 1], "uniform").pvalue < 1e-6)
True

inverse_rosenblatt

inverse_rosenblatt(copula: Copula, z: ArrayLike) -> NDArray[np.float64]

Turn independent uniform numbers into draws from a copula.

The inverse Rosenblatt transform, the reverse of :func:rosenblatt. Unlike copula.rvs, which makes its own randomness, this takes randomness you supply. That lets you push a quasi-random point set through a copula (see :mod:rcopula.sampling) or do conditional simulation: fix the first coordinates and vary the rest.

Parameters:

Name	Type	Description	Default
copula	Copula	The target copula, with its parameters set. Any family for which :func:rosenblatt works (Archimedean and elliptical in any dimension, others for dim=2).	required
z	array_like of float, shape (n, d)	Independent uniforms on :math:[0, 1], one row per draw. d must equal copula.dim. Values are clipped to [1e-12, 1 - 1e-12].	required

Returns:

Type	Description
numpy.ndarray of float, shape (n, d)	Draws from copula, each in (0, 1).

Raises:

Type	Description
ValueError	If the number of columns of z differs from copula.dim.
NotImplementedError	If :func:rosenblatt does not support copula (non-Archimedean, non-elliptical with d > 2).

Notes

Elliptical copulas invert analytically. Everything else is inverted by bisection against :func:rosenblatt itself, so the two are exact inverses by construction rather than by two derivations agreeing -- which is the failure mode this arrangement removes. The conditional CDF is monotone in its own argument, so bisection cannot fail; 60 steps takes the bracket below double precision.

Examples:

Round trip, in both directions:

>>> import numpy as np, rcopula as rc
>>> from rcopula.transforms import inverse_rosenblatt, rosenblatt
>>> cop = rc.ClaytonCopula(3.0, dim=3)
>>> u = cop.rvs(500, random_state=0)
>>> bool(np.max(np.abs(inverse_rosenblatt(cop, rosenblatt(cop, u)) - u)) < 1e-9)
True

Independent uniforms in, correctly dependent draws out:

>>> rng = np.random.default_rng(0)
>>> z = rng.uniform(size=(20000, 2))
>>> u = inverse_rosenblatt(rc.GumbelCopula(2.0), z)
>>> bool(abs(rc.cor_kendall(u)[0, 1] - 0.5) < 0.02)
True

radial_simplex

radial_simplex(copula: ArchimedeanCopula, u: ArrayLike) -> tuple[NDArray, NDArray]

Split Archimedean copula data into a "size" part and a "direction" part.

Every Archimedean copula sample can be written as a radial variable :math:R (one number per row, how far out the point is) times a direction :math:S that is uniform on the unit simplex. The direction part has the same distribution for every Archimedean copula, so if it does not look uniform, no Archimedean family fits the data. This is the decomposition of McNeil and Neslehova (2009).

Parameters:

Name	Type	Description	Default
copula	ArchimedeanCopula	Supplies the generator :math:\psi. Its parameter must be set.	required
u	array_like of float, shape (n, d)	Observations on the copula scale, strictly inside :math:(0, 1). d must equal copula.dim. A 1-D input is treated as a single row.	required

Returns:

Name	Type	Description
radial	numpy.ndarray of float, shape (n,)	:math:R = \sum_k \psi^{-1}(u_k), non-negative.
angular	numpy.ndarray of float, shape (n, d)	:math:S, whose rows sum to one. Where :math:R = 0 (every coordinate at 1, or so close that :math:\psi^{-1} underflows to 0) the direction is undefined -- 0/0 -- and is set to the simplex's barycentre, :math:1/d in every coordinate: the mean of the uniform angular law and the limit along the diagonal. Where :math:R = \infty (a coordinate at 0) the weight is shared equally among the infinite coordinates.

Raises:

Type	Description
TypeError	If copula is not an :class:ArchimedeanCopula.
ValueError	If the number of columns of u differs from copula.dim.

Notes

McNeil and Neslehova showed that :math:U \sim C for an Archimedean copula with generator :math:\psi if and only if

.. math:: \bigl(\psi^{-1}(U_1), \dots, \psi^{-1}(U_d)\bigr) = R\,S,

with :math:S uniform on the unit simplex and independent of the radial variable :math:R = \sum_k \psi^{-1}(U_k). All the family-specific information sits in the distribution of :math:R; the angular part is the same for every Archimedean copula in every dimension.

That is the whole basis for testing Archimedeanity rather than testing one particular family: if the angular part is not uniform on the simplex, no Archimedean copula fits, whatever generator is tried.

Examples:

The angular part is uniform on the simplex, so each of its coordinates has mean :math:1/d -- and this does not depend on the family or the parameter:

>>> import numpy as np, rcopula as rc
>>> from rcopula.transforms import radial_simplex
>>> for copula in (rc.ClaytonCopula(3.0, dim=4), rc.GumbelCopula(2.0, dim=4)):
...     u = copula.rvs(20000, random_state=0)
...     angular = radial_simplex(copula, u)[1]
...     print(bool(np.all(np.abs(angular.mean(axis=0) - 0.25) < 0.01)))
True
True

htrafo

htrafo(copula: ArchimedeanCopula, u: ArrayLike) -> NDArray[np.float64]

Turn Archimedean copula data into independent uniforms, even in high dimensions.

This is the Hering-Hofert transform (R's htrafo). It does the same job as :func:rosenblatt -- if the data came from copula, the output is independent and uniform -- but it avoids high-order derivatives of the generator, so it stays accurate at :math:d = 100 where the Rosenblatt transform breaks down. Use it for goodness-of-fit testing of Archimedean models with many variables.

Parameters:

Name	Type	Description	Default
copula	ArchimedeanCopula	The copula to test against, with its parameter set.	required
u	array_like of float, shape (n, d)	Observations on the copula scale, strictly inside :math:(0, 1). d must equal copula.dim.	required

Returns:

Type	Description
numpy.ndarray of float, shape (n, d)	Transformed values in [0, 1]. Columns 0 to d-2 come from the angular part, the last column from the radial part.

Raises:

Type	Description
TypeError	If copula is not an :class:ArchimedeanCopula.
ValueError	If the number of columns of u differs from copula.dim.

See Also

rosenblatt : the conditioning-based alternative, exact but derivative-hungry. radial_simplex : the decomposition this is built on.

Notes

The Rosenblatt transform conditions one coordinate at a time, which needs :math:d-1 derivatives of the generator and loses accuracy fast as :math:d grows. This transform goes through the simplex decomposition instead, so it needs no high-order derivatives at all -- which is why it exists.

Writing :math:S_j = \sum_{k \le j} \psi^{-1}(U_k),

.. math::

Y_j = \left(\frac{S_j}{S_{j+1}}\right)^{j}, \quad j = 1, \dots, d-1,
\qquad Y_d = K\bigl(\psi(S_d)\bigr),

with :math:K the Kendall distribution function. Under the null that u came from this copula the :math:Y_j are independent and uniform, so any test of multivariate uniformity becomes a goodness-of-fit test.

The first :math:d-1 components come from the angular part and the last from the radial part; since those are independent, so are the two groups.

Examples:

Data from the copula transforms to uniforms; data from a different copula does not:

>>> import numpy as np, rcopula as rc
>>> from scipy import stats
>>> from rcopula.transforms import htrafo
>>> copula = rc.ClaytonCopula(2.0, dim=5)
>>> y = htrafo(copula, copula.rvs(4000, random_state=0))
>>> bool(stats.kstest(y.ravel(), "uniform").pvalue > 0.01)
True
>>> wrong = htrafo(copula, rc.GumbelCopula(3.0, dim=5).rvs(4000, random_state=0))
>>> bool(stats.kstest(wrong.ravel(), "uniform").pvalue < 1e-6)
True

It still works where the Rosenblatt transform's derivatives would not:

>>> big = rc.GumbelCopula(2.0, dim=50)
>>> y = htrafo(big, big.rvs(500, random_state=0))
>>> bool(np.all(np.isfinite(y)) and y.shape == (500, 50))
True

radial_cdf

radial_cdf(copula: ArchimedeanCopula, x: ArrayLike) -> NDArray[np.float64]

Probability that the radial ("size") part of an Archimedean sample is at most x.

The distribution function of :math:R from :func:radial_simplex (R's pacR). Since the direction part is the same for every Archimedean copula, this curve carries all the family-specific information, which makes it the thing to compare between families or against data.

Parameters:

Name	Type	Description	Default
copula	ArchimedeanCopula	The copula, with its parameter set.	required
x	float or array_like of float, any shape	Radii at which to evaluate, each non-negative.	required

Returns:

Type	Description
numpy.ndarray of float, shape (m,)	:math:F_R(x) in :math:[0, 1], one value per element of x (flattened to 1-D).

Raises:

Type	Description
TypeError	If copula is not an :class:ArchimedeanCopula.
ValueError	If any element of x is negative.

Notes

:func:radial_simplex splits an Archimedean sample into a radial variable :math:R = \sum_j \psi^{-1}(U_j) and a point uniform on the simplex. The angular half is the same for every Archimedean copula in every dimension, so all the family-specific information lives in the distribution of :math:R -- which is this.

It needs no new machinery, because it is the Kendall distribution function in disguise. Since :math:C(\mathbf U) = \psi(R) and :math:\psi is decreasing,

.. math::

K(t) = P\{\psi(R) \le t\} = P\{R \ge \psi^{-1}(t)\},
\qquad\text{so}\qquad
F_R(x) = 1 - K(\psi(x)).

Examples:

>>> import numpy as np, rcopula as rc
>>> from rcopula.transforms import radial_cdf
>>> cop = rc.ClaytonCopula(2.0, dim=3)
>>> values = radial_cdf(cop, [1.5861, 11.6107, 346.2882])
>>> np.round(values, 3)
array([0.1  , 0.499, 0.9  ])

Which is what the sampled radial part actually does:

>>> radii = np.sum(cop.generator.ipsi(cop.rvs(50000, random_state=0), 2.0), axis=1)
>>> bool(np.all(np.abs(np.mean(radii[:, None] <= [1.5861, 11.6107], axis=0)
...                    - values[:2]) < 0.01))
True

radial_ppf

radial_ppf(copula: ArchimedeanCopula, q: ArrayLike) -> NDArray[np.float64]

Find the radius below which a given fraction of an Archimedean sample's radial part falls.

The quantile (inverse) function of :func:radial_cdf (R's qacR): returns x with :math:F_R(x) = q.

Parameters:

Name	Type	Description	Default
copula	ArchimedeanCopula	The copula, with its parameter set.	required
q	float or array_like of float, any shape	Probabilities, each in :math:[0, 1].	required

Returns:

Type	Description
numpy.ndarray of float, shape (m,)	Non-negative radii, one per element of q (flattened to 1-D).

Raises:

Type	Description
TypeError	If copula is not an :class:ArchimedeanCopula.
ValueError	If any element of q lies outside :math:[0, 1].

Notes

By the same identity as :func:radial_cdf: since :math:F_R(x) = 1 - K(\psi(x)), the quantile is :math:\psi^{-1}(K^{-1}(1 - q)), so it reduces to the Kendall function's own quantile and needs no root-finding of its own.

Examples:

>>> import numpy as np, rcopula as rc
>>> from rcopula.transforms import radial_cdf, radial_ppf
>>> cop = rc.GumbelCopula(2.0, dim=4)
>>> levels = np.array([0.1, 0.25, 0.5, 0.75, 0.9])
>>> bool(np.max(np.abs(radial_cdf(cop, radial_ppf(cop, levels)) - levels)) < 1e-6)
True

rcopula.kendall

The Kendall distribution function.

For a copula :math:C and :math:\mathbf U \sim C, the Kendall distribution function is the law of the copula evaluated at its own draw:

.. math:: K(t) = P\bigl(C(\mathbf U) \le t\bigr).

It is the univariate summary of a multivariate dependence structure: a single curve on :math:[0, 1] that determines Kendall's tau (:math:\tau = 3 - 4\int_0^1 K), identifies the generator of an Archimedean copula uniquely, and answers the question a joint quantile cannot.

That last point is what makes it useful rather than merely elegant. Ask how often is a flood this bad, and there is no single answer from the margins: the set :math:\{C(\mathbf u) > t\} is a region, not a point, and every point on its boundary is "equally extreme". :math:K measures that region. The Kendall return period :math:1/(1 - K(t)) is therefore the honest multivariate analogue of a return period, and it is systematically longer than the univariate one -- treating a compound event as if only its worst component mattered understates how rare it is.

============================= =============================================== :func:kendall_cdf :math:K(t); R's pK. :func:kendall_pdf Its density; R's dK. :func:kendall_ppf Its quantile function; R's qK. :func:kendall_rvs Draws of :math:C(\mathbf U); R's rK. :func:kendall_empirical The nonparametric estimator; R's Kn. :func:kendall_return_period :math:1/(1 - K(t)), and its inverse. ============================= ===============================================

For Archimedean copulas everything is closed form in any dimension. For other bivariate copulas :math:K is obtained by integrating the conditional distribution along the level curve of :math:C, which needs only the CDF.

References

Genest, C. and Rivest, L.-P. (1993). Statistical inference procedures for bivariate Archimedean copulas. JASA 88(423), 1034-1043. The bivariate form and the nonparametric estimator. Barbe, P., Genest, C., Ghoudi, K. and Remillard, B. (1996). On Kendall's process. Journal of Multivariate Analysis 58(2), 197-229. The multivariate generalisation used here. Salvadori, G. and De Michele, C. (2004). Frequency analysis via copulas. Water Resources Research 40, W12511. Kendall return periods, and why the univariate one is the wrong number. Nappo, G. and Spizzichino, F. (2009). Kendall distributions and level sets in bivariate exchangeable survival models. Information Sciences 179(17), 2878-2890.

kendall_cdf

kendall_cdf(copula: Copula, t: ArrayLike) -> NDArray[np.float64]

Probability that the copula, evaluated at a random draw from itself, is at most t.

This is the Kendall distribution function :math:K(t) = P(C(\mathbf U) \le t) (R's pK). It squeezes a whole multivariate dependence structure into one curve on :math:[0, 1]: the value :math:C(\mathbf u) says how "jointly low" a point is, and :math:K says how often points are at least that jointly low. Use it to compute joint (Kendall) return periods, to compare a fitted family with data via :func:kendall_empirical, or to recover Kendall's tau.

Parameters:

Name	Type	Description	Default
copula	Copula	The dependence model, with its parameters set. Archimedean families (Clayton, Gumbel, Frank, Joe, AMH, ...), :class:IndependenceCopula and :class:FrechetUpperCopula work in any dimension; every other family must be bivariate (dim=2).	required
t	float or array_like of float, any shape	Levels at which to evaluate :math:K. Values at or below 0 give 0 and values at or above 1 give 1.	required

Returns:

Type	Description
numpy.ndarray of float, shape (m,)	:math:K(t) for each level, in :math:[0, 1], where m is the number of elements in t (the input is flattened to 1-D; a scalar gives shape (1,)). :math:K(t) \ge t for every copula.

Raises:

Type	Description
NotImplementedError	If copula is neither Archimedean, independence nor comonotone and has dim other than 2.

Notes

For an Archimedean copula in :math:d dimensions (Barbe et al. 1996),

.. math:: K(t) = \sum_{k=0}^{d-1} \frac{(-1)^k}{k!}\, s^k\,\psi^{(k)}(s), \qquad s = \psi^{-1}(t),

which is a sum of non-negative terms once the alternating sign is absorbed into :math:|\psi^{(k)}|, so it never cancels. In two dimensions it reduces to the familiar :math:K(t) = t - s\,\psi'(s).

For any other bivariate copula, integrating along the level curve gives

.. math:: K(t) = t + \int_t^1 \frac{\partial C}{\partial u}\bigl(u, v_t(u)\bigr)\,du,

where :math:v_t(u) solves :math:C(u, v_t) = t. Below :math:u = t the event is certain, which is where the leading :math:t comes from.

Examples:

Independence has the closed form :math:K(t) = t(1 - \log t):

>>> import numpy as np
>>> from rcopula import IndependenceCopula
>>> from rcopula.kendall import kendall_cdf
>>> t = np.array([0.1, 0.3, 0.7])
>>> bool(np.allclose(kendall_cdf(IndependenceCopula(2), t), t * (1 - np.log(t))))
True

Comonotonicity leaves nothing to summarise -- :math:C(\mathbf U) = U, so :math:K is uniform:

>>> from rcopula import FrechetUpperCopula
>>> bool(np.allclose(kendall_cdf(FrechetUpperCopula(2), t), t))
True

K dominates the uniform for every copula, which is why the Kendall return period always exceeds the univariate one:

>>> from rcopula import ClaytonCopula
>>> bool(np.all(kendall_cdf(ClaytonCopula(3.0), t) >= t))
True

It recovers Kendall's tau through :math:\tau = 3 - 4\int_0^1 K:

>>> from scipy import integrate
>>> cop = ClaytonCopula(3.0)
>>> area = integrate.quad(lambda x: kendall_cdf(cop, x)[0], 0, 1)[0]
>>> bool(abs((3 - 4 * area) - cop.tau()) < 1e-6)
True

kendall_pdf

kendall_pdf(copula: Copula, t: ArrayLike) -> NDArray[np.float64]

Probability density of the Kendall distribution function (R's dK).

The slope of :func:kendall_cdf: where it is large, many draws of the copula land at that joint level. Useful for plotting the Kendall distribution or building likelihoods on it.

Parameters:

Name	Type	Description	Default
copula	Copula	The dependence model, with its parameters set. Same support as :func:kendall_cdf: Archimedean families in any dimension (closed form), other families bivariate only (numerical derivative).	required
t	float or array_like of float, any shape	Levels at which to evaluate the density. Points outside the open interval :math:(0, 1) give 0.	required

Returns:

Type	Description
numpy.ndarray of float, shape (m,)	Non-negative density values, one per element of t (flattened to 1-D; a scalar gives shape (1,)).

Raises:

Type	Description
NotImplementedError	If copula is not Archimedean and has dim other than 2 (raised by :func:kendall_cdf, which the numerical route calls). Independence and comonotone copulas use that numerical route.

Notes

Differentiating the Archimedean sum telescopes -- every term cancels against the next -- and only the last survives:

.. math:: K'(t) = \frac{s^{d-1}\,\bigl|\psi^{(d)}(s)\bigr|} {(d-1)!\,\bigl|\psi'(s)\bigr|}, \qquad s = \psi^{-1}(t),

manifestly non-negative, as a density must be. Elsewhere it is a central difference of :func:kendall_cdf with step 1e-6.

Examples:

Integrates to one, and matches the derivative of K:

>>> import numpy as np
>>> from scipy import integrate
>>> from rcopula import ClaytonCopula
>>> from rcopula.kendall import kendall_cdf, kendall_pdf
>>> cop = ClaytonCopula(3.0)
>>> total = integrate.quad(lambda x: kendall_pdf(cop, x)[0], 0, 1)[0]
>>> bool(abs(total - 1.0) < 1e-8)
True
>>> h = 1e-6
>>> slope = (kendall_cdf(cop, 0.4 + h)[0] - kendall_cdf(cop, 0.4 - h)[0]) / (2 * h)
>>> bool(abs(kendall_pdf(cop, 0.4)[0] - slope) < 1e-6)
True

kendall_ppf

kendall_ppf(copula: Copula, p: ArrayLike) -> NDArray[np.float64]

Find the joint level t below which a given fraction of copula draws fall.

The quantile (inverse) function of the Kendall distribution, R's qK: returns t with :math:K(t) = p. Solved by 60 steps of bisection, so it works for every family :func:kendall_cdf supports.

Parameters:

Name	Type	Description	Default
copula	Copula	The dependence model, with its parameters set (see :func:kendall_cdf for which families and dimensions are supported).	required
p	float or array_like of float, any shape	Probabilities, each in :math:[0, 1].	required

Returns:

Type	Description
numpy.ndarray of float, shape (m,)	Levels t in :math:[0, 1], one per element of p (flattened to 1-D). p = 0 maps to 0 and p = 1 maps to 1.

Raises:

Type	Description
ValueError	If any element of p lies outside :math:[0, 1].
NotImplementedError	If :func:kendall_cdf does not support copula.

Examples:

>>> import numpy as np
>>> from rcopula import GumbelCopula
>>> from rcopula.kendall import kendall_cdf, kendall_ppf
>>> cop = GumbelCopula(2.5)
>>> p = np.array([0.05, 0.5, 0.95])
>>> bool(np.allclose(kendall_cdf(cop, kendall_ppf(cop, p)), p, atol=1e-9))
True

kendall_rvs

kendall_rvs(copula: Copula, size: int = 1, random_state: Generator | int | None = None) -> NDArray[np.float64]

Draw random values from the Kendall distribution of a copula (R's rK).

It simulates points from the copula and evaluates the copula's CDF at each one, so the output is a sample of :math:C(\mathbf U).

Parameters:

Name	Type	Description	Default
copula	Copula	The dependence model, with its parameters set. Must support rvs and cdf (every family in the package does).	required
size	int	Number of draws; a positive integer.	1
random_state	(int, Generator or None)	Seed or generator for reproducibility. None uses fresh entropy.	None

Returns:

Type	Description
numpy.ndarray of float, shape (size,)	Independent draws, each in :math:[0, 1].

Notes

Exact by construction -- :math:K is the law of :math:C(\mathbf U) -- and cheaper than inverting :math:K.

Examples:

The draws reproduce the analytic K:

>>> import numpy as np
>>> from scipy import stats
>>> from rcopula import ClaytonCopula
>>> from rcopula.kendall import kendall_cdf, kendall_rvs
>>> cop = ClaytonCopula(3.0)
>>> w = kendall_rvs(cop, 20_000, random_state=0)
>>> p = stats.kstest(w, lambda x: kendall_cdf(cop, x)).pvalue
>>> bool(p > 0.01)
True

kendall_empirical

kendall_empirical(data: ArrayLike, t: ArrayLike | None = None) -> NDArray[np.float64]

Estimate the Kendall distribution function straight from data, with no model.

This is the nonparametric estimator of Genest and Rivest (1993), R's Kn. For each observation it counts the fraction of the other observations that are smaller in every coordinate, then takes the empirical distribution of those fractions. It uses only ranks, so the margins do not matter and no copula has to be fitted. Compare it with :func:kendall_cdf of a fitted family to see whether the family's joint behaviour matches the data.

Parameters:

Name	Type	Description	Default
data	array_like of float, shape (n, d)	Observations, one row per observation, on any scale (raw data or pseudo-observations). Needs n >= 2. A 1-D input is treated as a single row.	required
t	float, array_like of float or None	Levels at which to evaluate the estimate. None returns the sorted pseudo-values :math:W_i themselves instead.	None

Returns:

Type	Description
numpy.ndarray of float	If t is None: shape (n,), the sorted :math:W_i, each in :math:[0, 1]. Otherwise: shape (m,) with m the number of elements in t, the estimated :math:K at each level.

Raises:

Type	Description
ValueError	If data has fewer than two rows.

Notes

Genest & Rivest (1993) replace :math:C(\mathbf X_i) by the fraction of the sample it dominates,

.. math:: W_i = \frac{1}{n-1}\,#{j \ne i : \mathbf X_j < \mathbf X_i},

and take the empirical distribution of the :math:W_i. Being rank-based it needs no margins and no fitted copula, which makes it the natural check on whether a family's :math:K matches the data -- and the basis of the Sn^K goodness-of-fit statistics. Memory and time are :math:O(n^2).

Examples:

It tracks the true K without knowing the family:

>>> import numpy as np
>>> from rcopula import ClaytonCopula
>>> from rcopula.kendall import kendall_cdf, kendall_empirical
>>> cop = ClaytonCopula(3.0)
>>> x = cop.rvs(4000, random_state=0)
>>> grid = np.array([0.1, 0.3, 0.5, 0.8])
>>> gap = np.abs(kendall_empirical(x, grid) - kendall_cdf(cop, grid))
>>> bool(gap.max() < 0.03)
True

kendall_return_period

kendall_return_period(copula: Copula, t: ArrayLike, interval: float = 1.0) -> NDArray[np.float64]

How long, on average, until a joint event at least this severe happens again.

This is the Kendall return period :math:\text{interval}/(1 - K(t)): the multivariate counterpart of "the 100-year flood". Use it when an event is dangerous because several variables are extreme together (for example river flow and rainfall, or losses on several assets).

Parameters:

Name	Type	Description	Default
copula	Copula	The fitted dependence structure, with its parameters set (see :func:kendall_cdf for supported families and dimensions).	required
t	float or array_like of float, any shape	Critical level, in copula units: the value of :math:C(\mathbf u) at the event of interest, in :math:[0, 1].	required
interval	float	Mean time between observations, in whatever unit the answer should be in -- 1 for annual maxima. Must be positive.	1.0

Returns:

Type	Description
numpy.ndarray of float, shape (m,)	Return periods in the units of interval, one per element of t (flattened to 1-D). Exactly inf where :math:K(t) = 1 -- always at t = 1, an event that never happens, and possibly just below it once :math:K(t) rounds to 1 in floating point. This is the mathematically correct value, returned deliberately (no divide-by-zero warning is raised for it).

Raises:

Type	Description
ValueError	If interval is not positive, or any t lies outside :math:[0, 1] (or is NaN): no copula takes such a value, so there is no return period to report.
NotImplementedError	If :func:kendall_cdf does not support copula.

Notes

A univariate return period asks how often one variable is exceeded; the joint event :math:\{C(\mathbf U) > t\} is a region, and every point on its boundary is equally extreme, so no single margin can express it.

Because :math:K(t) \ge t for every copula, the Kendall return period is always longer than the univariate one at the same level. Reporting the univariate number for a compound event therefore overstates how often it happens.

Examples:

Comonotone variables move together exactly, so the joint event is the univariate one and the two return periods coincide. Everything else waits longer, and independence longest of all:

>>> from rcopula import FrechetUpperCopula, IndependenceCopula
>>> from rcopula.kendall import kendall_return_period
>>> float(round(kendall_return_period(FrechetUpperCopula(2), 0.99)[0], 1))
100.0
>>> float(round(kendall_return_period(IndependenceCopula(2), 0.99)[0], 1))
19933.2

The number that should give a practitioner pause: two copulas with the same Kendall's tau of 0.5 give 199 years and 6689 years. Gumbel has upper tail dependence and Clayton has none, and the critical layer is an upper-corner event -- so a rank correlation, however carefully estimated, does not determine the design life.

>>> from rcopula import ClaytonCopula, GumbelCopula
>>> float(round(kendall_return_period(GumbelCopula(2.0), 0.99)[0], 1))
199.0
>>> float(round(kendall_return_period(ClaytonCopula(2.0), 0.99)[0], 1))
6689.0

return_period_level

return_period_level(copula: Copula, period: ArrayLike, interval: float = 1.0) -> NDArray[np.float64]

Find the joint severity level that corresponds to a chosen return period.

The inverse of :func:kendall_return_period: given "the 100-year event", return the critical level :math:t (in copula units) whose Kendall return period is period. This is the quantity a design standard actually specifies -- the return period is fixed and the engineer needs the level that goes with it.

Parameters:

Name	Type	Description	Default
copula	Copula	The fitted dependence structure, with its parameters set (see :func:kendall_cdf for supported families and dimensions).	required
period	float or array_like of float, any shape	Target return periods, in the units of interval. Each must be at least interval: no event recurs, on average, more often than observations arrive. inf gives t = 1.	required
interval	float	Mean time between observations -- 1 for annual maxima. Must be positive.	1.0

Returns:

Type	Description
numpy.ndarray of float, shape (m,)	Critical levels t in :math:[0, 1], one per element of period (flattened to 1-D).

Raises:

Type	Description
ValueError	If interval is not positive, any period is not positive, or any period is shorter than interval (or is NaN). period == interval gives t = 0.

Examples:

>>> from rcopula import GumbelCopula
>>> from rcopula.kendall import kendall_return_period, return_period_level
>>> cop = GumbelCopula(2.0)
>>> t = return_period_level(cop, 100.0)
>>> bool(abs(kendall_return_period(cop, t)[0] - 100.0) < 1e-6)
True

rcopula.core.measures

Dependence measures by quadrature, for copulas with no closed form.

Structural copulas -- rotations, mixtures, Khoudraji devices -- generally have no analytic Kendall's tau. Both measures are integrals of the copula itself,

.. math:: \rho = 12\int_0^1!!\int_0^1 C(u,v)\,du\,dv - 3, \qquad \tau = 4\int_0^1!!\int_0^1 C(u,v)\,dC(u,v) - 1,

so a tensor-product rule evaluates them directly.

The rule is tanh-sinh (double exponential) rather than Gauss-Legendre, and the reason is specific: a copula with tail dependence has a density that blows up at a corner of the unit square. Gauss-Legendre assumes an analytic integrand and degrades to algebraic convergence when that fails -- for a Gumbel copula it reaches only 4e-4 with 128 nodes per axis and improves like :math:n^{-2}. The tanh-sinh substitution :math:u = \tfrac{1}{2}(1 + \tanh(\tfrac{\pi}{2}\sinh t)) sends the node spacing to zero doubly exponentially at each end, which absorbs endpoint singularities; the same Gumbel case reaches 5e-10 with 159 nodes.

References

Nelsen, R. B. (2006). An Introduction to Copulas, 2nd ed. Springer, Theorems 5.1.1 and 5.1.3 -- the integral representations used here. Takahasi, H. and Mori, M. (1974). Double exponential formulas for numerical integration. Publications of RIMS, Kyoto University 9(3), 721-741.

rho_by_quadrature

rho_by_quadrature(copula: Copula, level: int = _LEVEL) -> float

Compute Spearman's rho of a two-variable copula by numerical integration.

Spearman's rho is a rank correlation between -1 and 1: 0 means no monotone association, positive values mean the two variables tend to be large together. Use this for copulas that have no closed-form rho, such as rotations, mixtures and Khoudraji constructions.

Technically, it evaluates :math:\rho = 12\int\int C\,du\,dv - 3 with a tanh-sinh rule. It needs only the CDF, so it works for every bivariate copula -- including those with no density at all.

Parameters:

Name	Type	Description	Default
copula	Copula	A fully specified copula with dim == 2 (no NaN parameters).	required
level	int	Half-width of the tanh-sinh node index; larger means more nodes per axis and higher accuracy at more cost. The default gives roughly 150 usable nodes per axis and about 1e-9 accuracy.	140

Returns:

Type	Description
float	Spearman's rho, clipped to [-1, 1].

Raises:

Type	Description
ValueError	If copula.dim is not 2, or the copula has unspecified parameters.

Examples:

Reproduces the analytic value for families that have one:

>>> from rcopula import ClaytonCopula, GumbelCopula
>>> from rcopula.core.measures import rho_by_quadrature
>>> for cop in (ClaytonCopula(2.0), GumbelCopula(2.5)):
...     print(bool(abs(rho_by_quadrature(cop) - cop.rho()) < 1e-8))
True
True

tau_by_quadrature

tau_by_quadrature(copula: Copula, level: int = _LEVEL) -> float

Compute Kendall's tau of a two-variable copula by numerical integration.

Kendall's tau is a rank correlation between -1 and 1 based on how often pairs of observations move in the same direction. Use this for copulas with a density but no closed-form tau. If the copula has a singular part or very strong dependence, prefer :func:tau_by_partials.

Technically, it evaluates :math:\tau = 4\int\int C\,dC - 1 = 4\int\int C\,c\,du\,dv - 1 with a tanh-sinh rule. It needs the density as well as the CDF. Unlike Spearman's rho, tau is not linear in the copula -- which is exactly why mixtures and rotations cannot simply average their components' values and have to come here instead.

Parameters:

Name	Type	Description	Default
copula	Copula	A fully specified copula with dim == 2 that has a density (pdf).	required
level	int	Half-width of the tanh-sinh node index; larger means more nodes per axis and higher accuracy at more cost.	140

Returns:

Type	Description
float	Kendall's tau, clipped to [-1, 1].

Raises:

Type	Description
ValueError	If copula.dim is not 2, or the copula has unspecified parameters.

Examples:

>>> from rcopula import GumbelCopula, JoeCopula
>>> from rcopula.core.measures import tau_by_quadrature
>>> for cop in (GumbelCopula(2.5), JoeCopula(3.0)):
...     print(bool(abs(tau_by_quadrature(cop) - cop.tau()) < 1e-7))
True
True

tau_by_partials

tau_by_partials(copula: Copula, nodes: int = _PARTIAL_NODES) -> float

Compute Kendall's tau of a two-variable copula using only its CDF.

Kendall's tau is a rank correlation between -1 and 1. This is the most robust of the quadrature routines: use it when the copula has a singular part (no density everywhere) or very strong dependence.

Technically, it evaluates :math:\tau = 1 - 4\int\int \partial_1 C\,\partial_2 C\,du\,dv, with the partial derivatives taken by central differences and a Gauss-Legendre rule. The same quantity as :func:tau_by_quadrature, integrated a different way, and better in two situations that matter.

• It needs only the CDF. The :math:\int\int C\,dC form silently assumes the copula is absolutely continuous: a copula with a singular component -- Marshall-Olkin, or Clayton below zero -- carries mass that no density represents, and integrating the density alone misses it. Here the partial derivatives are conditional distribution functions, which exist regardless.
• The integrand is bounded in :math:[0, 1]. Under strong dependence the density concentrates into a near-delta on the diagonal and no quadrature can resolve it: at Plackett theta = 1e4 the density form returns 1.0 against a true 0.9756. The bounded form gets 0.97572.

The trade is a numerical derivative, so it carries about 4e-8 rather than 1e-8 -- immaterial next to what it buys.

Parameters:

Name	Type	Description	Default
copula	Copula	A fully specified copula with dim == 2.	required
nodes	int	Gauss-Legendre nodes per axis. The cost grows with nodes**2 CDF evaluations (times four, for the differences).	400

Returns:

Type	Description
float	Kendall's tau, clipped to [-1, 1].

Raises:

Type	Description
ValueError	If copula.dim is not 2, or the copula has unspecified parameters.

Examples:

Agrees with the closed forms:

>>> from rcopula import ClaytonCopula, GumbelCopula
>>> from rcopula.core.measures import tau_by_partials
>>> for cop in (ClaytonCopula(2.0), GumbelCopula(3.0)):
...     print(bool(abs(tau_by_partials(cop) - cop.tau()) < 1e-7))
True
True

And works where the density form cannot -- Marshall-Olkin puts mass on a curve, so its density integrates to less than one:

>>> from rcopula import MarshallOlkinCopula
>>> cop = MarshallOlkinCopula(0.7, 0.4)
>>> bool(abs(tau_by_partials(cop) - cop.tau()) < 1e-3)
True

rcopula.sampling

Quasi-random and variance-reduced copula sampling.

Every Monte Carlo answer in this package -- a CDO tranche spread, a basket option, a portfolio's expected shortfall -- has a standard error, and the default way to halve it is to quadruple the sample. These are the cheaper ways.

Quasi-random. A Sobol or Halton point set fills the unit cube more evenly than independent uniforms do, because it is designed to rather than merely tending to. Pushing such a set through the inverse Rosenblatt transform gives a sample from the copula whose integration error falls like :math:(\log n)^d / n rather than :math:n^{-1/2}. For a smooth payoff in a few dimensions that is the difference between four digits and two.

The catch, stated plainly because it is usually not: the improvement depends on the integrand being smooth and of low effective dimension. A payoff with a kink (any option), a high dimension, or a discontinuity (a default indicator) erodes it, and the rate can fall back to :math:n^{-1/2}. :func:variance_ratio measures what you actually got instead of assuming it.

Antithetic. Pair every draw :math:z with :math:1 - z. Free, and helps exactly when the integrand is monotone in the uniforms -- which a European call on a single asset is, and a straddle is not. It can increase variance for a symmetric payoff, so it is not a default.

Latin hypercube. Stratifies each margin into :math:n equal bins and permutes. Guarantees marginal coverage, does nothing about the joint structure, and is a good default when the answer depends mostly on the margins.

============================ ================================================ :func:quasi_rvs Sobol or Halton draws from a copula. :func:antithetic_rvs Draws paired with their reflections. :func:latin_hypercube_rvs Stratified draws. :func:variance_ratio What the variance reduction actually was. ============================ ================================================

Examples:

>>> import numpy as np, rcopula as rc
>>> from rcopula.sampling import quasi_rvs
>>> u = quasi_rvs(rc.ClaytonCopula(2.0), 1024, random_state=0)
>>> u.shape
(1024, 2)
>>> bool(abs(rc.cor_kendall(u)[0, 1] - 0.5) < 0.02)
True

References

Cambou, M., Hofert, M. and Lemieux, C. (2017). Quasi-random numbers for copula models. Statistics and Computing 27(5), 1307-1329. The paper this implements: low-discrepancy sequences through the inverse Rosenblatt transform. Sobol', I. M. (1967). On the distribution of points in a cube and the approximate evaluation of integrals. USSR Computational Mathematics and Mathematical Physics 7(4), 86-112. Owen, A. B. (1997). Scrambled net variance for integrals of smooth functions. Annals of Statistics 25(4), 1541-1562. Why scrambling is what makes an error estimate possible at all. McKay, M. D., Beckman, R. J. and Conover, W. J. (1979). A comparison of three methods for selecting values of input variables in the analysis of output from a computer code. Technometrics 21(2), 239-245. Latin hypercube sampling. Glasserman, P. (2003). Monte Carlo Methods in Financial Engineering. Springer. Chapter 4 on variance reduction, chapter 5 on quasi-Monte Carlo.

quasi_rvs

quasi_rvs(copula: Copula, size: int, *, sequence: Sequence = 'sobol', scramble: bool = True, random_state: Any = None) -> NDArray[np.float64]

Draw a sample from a copula whose points are spread more evenly than random ones.

This is quasi-random (quasi-Monte Carlo) sampling. Use it in place of copula.rvs when you are averaging something over the sample -- a price, a risk measure -- and want a smaller error for the same number of draws.

It generates a Sobol or Halton set in :math:[0,1]^d and pushes it through the copula's inverse Rosenblatt transform (Cambou, Hofert and Lemieux 2017). The result has the right distribution and covers the space more evenly than independent draws.

Parameters:

Name	Type	Description	Default
copula	Copula	Any fitted or constructed copula of dimension d; it must support the inverse Rosenblatt transform.	required
size	int	Number of points, at least 1. For Sobol, a power of two (512, 1024, ...) preserves the balance properties the construction is built on; scipy warns otherwise, and so does the note below.	required
sequence	('sobol', 'halton')	Which point set to use. Sobol is the usual choice and better in moderate dimensions; Halton is simpler and degrades faster past about ten.	"sobol"
scramble	bool	Owen scrambling. Keep it on: an unscrambled set is deterministic, so it gives one answer with no way to estimate its error, and its first point is the corner of the cube. Scrambling preserves the low discrepancy and restores the ability to replicate.	True
random_state	(int, Generator or None)	Seeds the scrambling. Pass an int for reproducible output.	None

Returns:

Type	Description
numpy.ndarray of float, shape (size, d)	Draws from the copula, every value strictly between 0 and 1.

Raises:

Type	Description
ValueError	If size is less than 1, or sequence is not "sobol" or "halton".

Notes

The gain depends on the integrand. Smooth and low-dimensional: large. Kinked (any option payoff), high-dimensional, or discontinuous (a default indicator): smaller, sometimes none. Measure it with :func:variance_ratio rather than assuming.

Examples:

The margins are uniform and the dependence is right, as for rvs:

>>> import numpy as np, rcopula as rc
>>> from rcopula.sampling import quasi_rvs
>>> u = quasi_rvs(rc.GumbelCopula(2.0, dim=3), 4096, random_state=0)
>>> bool(np.all(np.abs(u.mean(axis=0) - 0.5) < 0.01))
True
>>> bool(abs(rc.cor_kendall(u)[0, 1] - 0.5) < 0.02)
True

Coverage is visibly more even -- the largest gap between consecutive sorted values in a margin is smaller than for independent draws:

>>> plain = rc.GumbelCopula(2.0, dim=3).rvs(4096, random_state=0)
>>> gap = lambda a: np.max(np.diff(np.sort(a[:, 0])))
>>> bool(gap(u) < gap(plain))
True

antithetic_rvs

antithetic_rvs(copula: Copula, size: int, *, random_state: Any = None) -> NDArray[np.float64]

Draw a sample from a copula in mirror-image pairs, to reduce Monte Carlo error.

This is antithetic sampling. It is free and helps when the quantity you average goes up (or down) steadily as the draws go up, such as a call option payoff; it can hurt for symmetric payoffs (see the warning below).

Each independent uniform vector :math:z is used twice, as :math:z and :math:1-z, before the inverse Rosenblatt transform. The two resulting draws are negatively associated, so a Monte Carlo average over them has lower variance when the integrand is monotone in the uniforms.

Parameters:

Name	Type	Description	Default
copula	Copula	Any copula of dimension d that supports the inverse Rosenblatt transform.	required
size	int	Total number of draws, at least 1. Rounded up to an even number, since they come in pairs; the first half and second half are the partners of each other, so with m = len(u) // 2 row u[i] pairs with u[m + i].	required
random_state	(int, Generator or None)	Seed for the underlying uniforms. Pass an int for reproducible output.	None

Returns:

Type	Description
numpy.ndarray of float, shape (2 * ceil(size / 2), d)	Draws from the copula: size rows when size is even, and size + 1 rows when it is odd, so that every draw keeps its partner.

Raises:

Type	Description
ValueError	If size is less than 1.

Warnings

For a payoff symmetric about the middle of the distribution -- a straddle, an absolute deviation, a variance -- antithetic pairing can raise variance rather than lower it, because the two halves of a pair move the same way. Check with :func:variance_ratio before adopting it.

Examples:

>>> import numpy as np, rcopula as rc
>>> from rcopula.sampling import antithetic_rvs
>>> u = antithetic_rvs(rc.ClaytonCopula(2.0), 2000, random_state=0)
>>> u.shape
(2000, 2)

The pairing is exact on the first coordinate, which the transform leaves alone:

>>> half = 1000
>>> bool(np.allclose(u[:half, 0] + u[half:, 0], 1.0))
True

latin_hypercube_rvs

latin_hypercube_rvs(copula: Copula, size: int, *, random_state: Any = None) -> NDArray[np.float64]

Draw a sample from a copula whose margins are spread evenly by construction.

This is Latin hypercube sampling. Use it when the answer you are computing depends mostly on the individual margins rather than on how they move together.

Each of the :math:d uniform coordinates is stratified into size equal bins with exactly one point per bin, then the columns are permuted independently. The marginal coverage is therefore guaranteed rather than random, which is worth having when the answer depends mostly on the margins.

It says nothing about the joint structure -- that comes from the inverse Rosenblatt transform, as usual.

Parameters:

Name	Type	Description	Default
copula	Copula	Any copula of dimension d that supports the inverse Rosenblatt transform.	required
size	int	Number of draws, at least 1. Also the number of bins per margin.	required
random_state	(int, Generator or None)	Seed for the design. Pass an int for reproducible output.	None

Returns:

Type	Description
numpy.ndarray of float, shape (size, d)	Draws from the copula, every value strictly between 0 and 1.

Raises:

Type	Description
ValueError	If size is less than 1.

Examples:

Every margin covers its range by construction, so the sample mean of a uniform coordinate is almost exactly one half:

>>> import numpy as np, rcopula as rc
>>> from rcopula.sampling import latin_hypercube_rvs
>>> u = latin_hypercube_rvs(rc.FrankCopula(5.0, dim=3), 2000, random_state=0)
>>> bool(abs(u[:, 0].mean() - 0.5) < 0.005)
True
>>> bool(abs(rc.cor_kendall(u)[0, 1] - rc.FrankCopula(5.0).tau()) < 0.03)
True

variance_ratio

variance_ratio(copula: Copula, payoff: Callable[[NDArray[float64]], ArrayLike], size: int, *, method: str = 'sobol', replicates: int = 20, random_state: Any = None) -> dict[str, float]

Measure how much a variance-reduction sampler shrinks the Monte Carlo error for your payoff.

Use it before adopting :func:quasi_rvs, :func:antithetic_rvs or :func:latin_hypercube_rvs: it tells you, for your copula and your payoff, whether the method helps and by how much.

Runs the estimator replicates times under plain Monte Carlo and again under method, and compares the spread of the answers. This is the honest way to report a quasi-Monte Carlo gain: the theory gives a rate, not a number, and the number depends on the payoff.

Parameters:

Name	Type	Description	Default
copula	Copula	The copula to sample from, of dimension d.	required
payoff	callable	Function payoff(u) taking a numpy.ndarray of shape (n, d) of copula draws and returning an array_like of float, shape (n,): one value per row. The quantity estimated is its mean.	required
size	int	Draws per replicate, at least 1.	required
method	('sobol', 'halton', 'antithetic', 'lhs')	The sampler to compare against plain copula.rvs: :func:quasi_rvs with a Sobol or Halton set, :func:antithetic_rvs, or :func:latin_hypercube_rvs.	"sobol"
replicates	int	Independent repetitions, at least 2. Each uses a different scramble or seed, which is what makes a spread measurable at all -- an unscrambled quasi-random set would give the same answer every time and no error estimate.	20
random_state	(int, Generator or None)	Seed from which every replicate's seed is drawn.	None

Returns:

Type

Description

dict of str to float

"plain_mean" Average of the plain Monte Carlo estimates. "reduced_mean" Average of the estimates under method; should agree with "plain_mean" up to noise. "plain_se" Standard deviation of the plain estimates across replicates, i.e. the standard error of one plain estimate of size draws. "reduced_se" The same under method. "ratio" plain_se / reduced_se. Above 1 means the method helped; inf if reduced_se is zero. "equivalent_sample_factor" ratio ** 2: how many times more plain draws would be needed to match the method's error.

Raises:

Type	Description
ValueError	If replicates is less than 2, or method is not one of the four names above.

Examples:

A smooth payoff, where quasi-random sampling does well:

>>> import numpy as np, rcopula as rc
>>> from rcopula.sampling import variance_ratio
>>> smooth = lambda u: u[:, 0] * u[:, 1]
>>> out = variance_ratio(rc.ClaytonCopula(2.0), smooth, 1024, replicates=12,
...                      random_state=0)
>>> bool(out["ratio"] > 2.0)
True
>>> bool(abs(out["plain_mean"] - out["reduced_mean"]) < 0.01)
True