Goodness of fit and hypothesis tests

rcopula.gof.api

Goodness-of-fit testing for copulas.

Two ways to get a p-value out of a Cramer-von Mises statistic, because its null distribution depends on the unknown parameter and so cannot be tabulated.

Parametric bootstrap (simulation="pb"). Simulate from the fitted copula, refit, recompute the statistic, repeat. Correct and simple, and the reference procedure of Genest, Remillard & Beaudoin (2009) -- but it refits the copula on every replicate, so a thousand replicates cost a thousand fits.

Multiplier bootstrap (simulation="mult"). Replicate the limiting empirical process directly instead of resampling. Because :math:\sqrt n(\hat\theta - \theta) \approx H^{-1} n^{-1/2}\sum_i W_i for the same influence function :math:W_i the standard errors already use, a bootstrap replicate is obtained by reweighting with random multipliers :math:Z_i -- no refitting at all. Orders of magnitude faster, and the procedure that makes goodness-of-fit practical at large n.

The multiplier bootstrap has no implementation anywhere else in Python.

The p-value follows Pesarin's convention, :math:(0.5 + \#\{T_b \ge T\})/(N+1), which is strictly inside (0, 1) -- a p-value of exactly zero would be an artefact of finite N, not evidence.

References

Genest, C., Remillard, B. and Beaudoin, D. (2009). Goodness-of-fit tests for copulas: a review and a power study. Insurance: Mathematics and Economics 44(2), 199-213. Kojadinovic, I., Yan, J. and Holmes, M. (2011). Fast large-sample goodness-of-fit tests for copulas. Statistica Sinica 21(2), 841-871. The multiplier bootstrap implemented here. Remillard, B. and Scaillet, O. (2009). Testing for equality between two copulas. Journal of Multivariate Analysis 100(3), 377-386. Pesarin, F. (2001). Multivariate Permutation Tests. Wiley.

GofResult

Bases: NamedTuple

The result of a goodness-of-fit test: the score, its p-value, and how it was computed.

Returned by :func:gof_test and :func:gof_two_sample; you do not normally build one yourself. A small pvalue (say below 0.05) means the data are unlikely under the tested copula (or, for the two-sample test, that the two samples have different copulas).

It is a NamedTuple following the scipy.stats convention of exposing statistic and pvalue, so it unpacks like any other test result: stat, p, *rest = result. The constructor takes the attributes below, in order, as its arguments; copula defaults to None.

Attributes:

Name	Type	Description
statistic	float	The observed test statistic, zero or positive; larger means a worse fit.
pvalue	float	Bootstrap or permutation p-value, strictly between 0 and 1.
method	str	Name of the statistic used: "Sn", "Tn", "AnChisq", "AnGamma" or, for the two-sample test, "two-sample Sn".
simulation	str	How the null distribution was generated: "pb" (parametric bootstrap), "mult" (multiplier bootstrap) or "permutation".
n_rep	int	Number of bootstrap replicates or permutations used.
copula	Copula or None, default None	The fitted copula the data was tested against. None for a two-sample test, which compares two samples to each other with no model between them.

gof_test

gof_test(copula: Copula, data: ArrayLike, *, method: str = 'Sn', simulation: str = 'pb', estim_method: str = 'mpl', n_rep: int = 1000, random_state: Generator | int | None = None, ties_method: str = 'average') -> GofResult

Check whether a copula family is a plausible model for your data, and get a p-value.

The family's parameters are first fitted to the data; the test then asks whether the gap between that fitted copula and the data's own empirical copula is bigger than chance alone would produce (R's gofCopula). Use it after choosing a family, to see whether the data clearly contradict it. A small p-value (say below 0.05) means "this family does not fit".

Parameters:

Name	Type	Description	Default
copula	Copula	The family to test, e.g. ClaytonCopula(). Its current parameter values are only a template; parameters are re-estimated from the data.	required
data	array_like of float, shape (n, d)	Observations, one row per observation and one column per variable, on any scale. They are always converted to ranks (pseudo-observations) first, so even data already in [0, 1] is re-ranked.	required
method	('Sn', 'Tn', 'AnChisq', 'AnGamma')	Which statistic to use (see :mod:rcopula.gof.statistics). "Sn" is the default and the recommended one.	"Sn"
simulation	('pb', 'mult')	How the p-value is computed. "pb" (parametric bootstrap) simulates and refits on every replicate; "mult" (multiplier bootstrap) reweights instead and is far faster. "mult" supports "Sn" only.	"pb"
estim_method	str	Estimation method passed to :func:~rcopula.fit, e.g. "mpl", "itau" or "irho".	"mpl"
n_rep	int	Number of bootstrap replicates, a positive integer. More replicates give a more precise p-value at proportionally more cost.	1000
random_state	(int, Generator or None)	Seed or generator, for reproducible p-values.	None
ties_method	str	How tied values are ranked; passed to :func:~rcopula.pseudo_obs. One of "average", "min", "max", "dense", "ordinal", "random".	"average"

Returns:

Type	Description
GofResult	With statistic and pvalue, as in scipy.stats, plus the fitted copula in copula.

Raises:

Type	Description
ValueError	If simulation is not "pb" or "mult", if simulation="mult" is combined with a method other than "Sn", if method is not a recognised statistic, or if estim_method is not a :func:~rcopula.fit method. All of these are checked before anything is fitted.

Notes

A large p-value is not evidence that the family is correct -- only that this test could not reject it at this sample size.

The "pb" cost is n_rep full model fits, which can be slow for large n or many parameters; "mult" needs one fit in total.

Examples:

The right family is not rejected; the wrong one is:

>>> import numpy as np
>>> from rcopula import ClaytonCopula, GumbelCopula, gof_test
>>> x = ClaytonCopula(4.0).rvs(300, random_state=0)
>>> right = gof_test(ClaytonCopula(), x, n_rep=100, random_state=1)
>>> wrong = gof_test(GumbelCopula(), x, n_rep=100, random_state=1)
>>> bool(right.pvalue > 0.05)
True
>>> bool(wrong.pvalue < 0.05)
True

The multiplier bootstrap reaches the same conclusion without refitting:

>>> fast = gof_test(GumbelCopula(), x, simulation="mult", n_rep=100, random_state=1)
>>> bool(fast.pvalue < 0.05)
True

gof_two_sample

gof_two_sample(x: ArrayLike, y: ArrayLike, n_rep: int = 1000, random_state: Generator | int | None = None, ties_method: str = 'average') -> GofResult

Check whether two datasets have the same dependence structure, ignoring their scales.

This is a two-sample copula test (R's gofT2stat). A small p-value (say below 0.05) means the way the variables move together differs between the two samples.

Every other test here compares one sample against a model. This compares two samples against each other, with no model in between -- which is the right question for "did the dependence change after the crisis?", "does this desk's book behave like that one's?", or "does my simulator reproduce the dependence it was fitted to?".

The statistic is a Cramer-von Mises distance between the two empirical copulas, evaluated at the pooled observations,

.. math::

T = \frac{n m}{n + m} \sum_{z} \bigl(C_n(z) - C_m(z)\bigr)^2,

and the null is generated by permuting the sample labels -- an exact randomisation test, valid at any sample size rather than asymptotically.

Parameters:

Name	Type	Description	Default
x	array_like of float, shape (n, d)	The first sample, on any scale. At least 2 rows.	required
y	array_like of float, shape (m, d)	The second sample, on any scale. At least 2 rows. It needs the same number of columns as x, not the same number of rows.	required
n_rep	int	Number of random relabellings (permutations) used to build the null distribution, a positive integer.	1000
random_state	(int, Generator or None)	Seed or generator, for reproducible p-values.	None
ties_method	str	How tied values are ranked; passed to :func:~rcopula.pseudo_obs for the observed statistic and for every permutation replicate alike, so the null is built the same way as the statistic it calibrates.	"average"

Returns:

Type	Description
GofResult	With method="two-sample Sn", simulation="permutation" and copula=None.

Raises:

Type	Description
ValueError	If x and y have different numbers of columns, or either has fewer than 2 rows.

Notes

Ranks are taken within each sample separately, which is the whole design and not an implementation detail: it makes the test blind to the margins, so two samples with wildly different scales and shapes but the same dependence are indistinguishable to it. If you want to compare the full distributions, this is the wrong test.

Examples:

Two samples from the same copula are not separated:

>>> import rcopula as rc
>>> from rcopula.gof import gof_two_sample
>>> a = rc.ClaytonCopula(2.0).rvs(600, random_state=0)
>>> b = rc.ClaytonCopula(2.0).rvs(600, random_state=1)
>>> bool(gof_two_sample(a, b, n_rep=200, random_state=0).pvalue > 0.05)
True

Two different copulas at the same Kendall's tau are:

>>> gumbel = rc.GumbelCopula.from_tau(0.5).rvs(600, random_state=1)
>>> clayton = rc.ClaytonCopula.from_tau(0.5).rvs(600, random_state=0)
>>> bool(gof_two_sample(clayton, gumbel, n_rep=200, random_state=0).pvalue < 0.05)
True

And changing only the margins changes nothing, which is the point:

>>> import numpy as np
>>> from scipy import stats
>>> same = rc.ClaytonCopula(2.0).rvs(600, random_state=1)
>>> rescaled = np.column_stack([stats.expon(scale=50).ppf(same[:, 0]),
...                             stats.norm(loc=-9, scale=0.01).ppf(same[:, 1])])
>>> bool(gof_two_sample(a, rescaled, n_rep=200, random_state=0).pvalue > 0.05)
True

rcopula.gof.statistics

Goodness-of-fit test statistics for copulas.

All of these compare a fitted parametric copula against the data's own empirical copula. They differ in what they compare and hence in what kind of departure they are sensitive to.

Sn Cramer-von Mises distance between the empirical and fitted copulas, :math:S_n = \sum_i \{C_n(\hat U_i) - C_{\hat\theta}(\hat U_i)\}^2. The blanket test: no tuning, good all-round power, and the one to reach for by default.

Tn The Kolmogorov-Smirnov analogue, :math:T_n = \sqrt n \max_i |C_n(\hat U_i) - C_{\hat\theta}(\hat U_i)|, the sup-norm of the same empirical process :math:\sqrt n (C_n - C_{\hat\theta}) whose squared :math:L^2 norm is Sn (Genest, Remillard and Beaudoin 2009). R's gofCopula has no Tn, so there is no R value to match; the scaling only multiplies the statistic by a constant at fixed n, so bootstrap p-values are unaffected by it. Less powerful than Sn in nearly every published comparison, but occasionally more sensitive to a single localised discrepancy.

AnChisq / AnGamma Anderson-Darling statistics applied after collapsing the Rosenblatt- transformed data to one dimension -- by a chi-squared or a gamma aggregation respectively. Cheap, and sensitive in the tails where Sn is weakest, but they only see the aggregated variable.

References

Genest, C., Remillard, B. and Beaudoin, D. (2009). Goodness-of-fit tests for copulas: a review and a power study. Insurance: Mathematics and Economics 44(2), 199-213. The definitive comparison; Sn is its recommendation. Genest, C. and Remillard, B. (2008). Validity of the parametric bootstrap for goodness-of-fit testing in semiparametric models. Annales de l'IHP Probabilites et Statistiques 44(6), 1096-1127.

empirical_copula_at

empirical_copula_at(u: ArrayLike, at: ArrayLike | None = None) -> NDArray[np.float64]

Compute, for each point, the share of observations that lie below it in every coordinate.

This is the empirical copula :math:C_n: the data's own, model-free estimate of the copula. Goodness-of-fit statistics compare it with a fitted parametric copula, and you can use it directly to see what the dependence in a sample looks like without assuming any family.

Parameters:

Name	Type	Description	Default
u	array_like of float, shape (n, d)	The observations defining the empirical copula, normally pseudo-observations in [0, 1] (see :func:~rcopula.pseudo_obs). A 1-D input is treated as a single row.	required
at	array_like of float, shape (m, d), or None	The points at which to evaluate it. None evaluates at the rows of u themselves, so m = n. Must have the same number of columns as u.	None

Returns:

Type	Description
numpy.ndarray of float, shape (m,)	For each point, the fraction of rows of u that are less than or equal to it in every coordinate; values lie in [0, 1].

Notes

Formally :math:C_n(\mathbf x) = \frac1n \sum_i \mathbf 1\{\mathbf U_i \le \mathbf x\}. Uses the 1/n scaling of Genest-Remillard-Beaudoin, matching R's C.n with offset = 0. Memory and time grow as m * n * d.

Examples:

>>> import numpy as np
>>> from rcopula.gof.statistics import empirical_copula_at
>>> u = np.array([[0.2, 0.3], [0.5, 0.6], [0.8, 0.9]])
>>> empirical_copula_at(u)
array([0.33333333, 0.66666667, 1.        ])

gof_statistic

gof_statistic(u: ArrayLike, copula: Copula | None = None, method: str = 'Sn') -> float

Score how far a fitted copula is from the data, where a bigger number means a worse fit.

This computes one goodness-of-fit statistic (R's gofTstat) and nothing else -- no p-value. Use it to compare candidate families on the same data, or as a building block; for a formal test with a p-value use :func:~rcopula.gof_test, which calibrates this number by bootstrap.

Parameters:

Name	Type	Description	Default
u	array_like of float, shape (n, d)	Pseudo-observations, values in [0, 1] (see :func:~rcopula.pseudo_obs). For "AnChisq"/"AnGamma" these should already be Rosenblatt-transformed with the fitted copula, as in R; this function does not do that transform for you.	required
copula	Copula or None	The fitted copula to compare against. Required for "Sn" and "Tn"; ignored by "AnChisq" and "AnGamma".	None
method	('Sn', 'Tn', 'AnChisq', 'AnGamma')	Which statistic to compute; see the module notes. "Sn" is the recommended all-round choice.	"Sn"

Returns:

Type	Description
float	The statistic, zero or positive. Larger values indicate worse fit. Values are only comparable between fits made with the same method and the same data size.

Raises:

Type	Description
ValueError	If method is not one of the four names above, or if method is "Sn" or "Tn" and copula is None.

Notes

Sn is :math:\sum_i \{C_n(\hat U_i) - C_{\hat\theta}(\hat U_i)\}^2. Tn is :math:\sqrt n \max_i |C_n(\hat U_i) - C_{\hat\theta}(\hat U_i)|, the Kolmogorov-Smirnov statistic of Genest, Remillard and Beaudoin (2009). (Versions up to 0.2 returned n rather than sqrt(n) times the maximum, contradicting the documentation; the p-value of :func:~rcopula.gof_test is the same under either scaling.) AnChisq and AnGamma collapse each row to one number (a sum of squared normal quantiles, or a sum of :math:-\log u), map it to (0, 1) with the chi-squared or gamma distribution function, and apply the Anderson-Darling uniformity statistic.

Examples:

A well-specified model gives a small Sn; a misspecified one does not:

>>> import numpy as np
>>> from rcopula import ClaytonCopula, GumbelCopula, fit, pseudo_obs
>>> from rcopula.gof import gof_statistic
>>> u = pseudo_obs(ClaytonCopula(4.0).rvs(500, random_state=0))
>>> right = gof_statistic(u, fit(ClaytonCopula(), u).copula)
>>> wrong = gof_statistic(u, fit(GumbelCopula(), u).copula)
>>> bool(wrong > 5 * right)
True

rcopula.htest.api

Specification tests for copulas.

These answer questions about the dependence structure that come before choosing a family, and each one narrows the search:

indep_test Is there any dependence at all? If not, nothing else matters. exch_test Is the copula exchangeable, :math:C(u,v) = C(v,u)? Every Archimedean and elliptical family is. Rejecting exchangeability rules all of them out at once and points at the Khoudraji device or a rotated family. rad_sym_test Is the copula radially symmetric, equal to its own survival copula? Gaussian, Student-t and Frank are; Clayton, Gumbel and Joe are not. Rejecting radial symmetry means the tails differ, which is exactly what a risk model needs to get right. ev_test Is the copula an extreme-value copula? The natural null when modelling componentwise maxima.

Null distributions are obtained by randomisation rather than by R's multiplier bootstrap. Each null hypothesis here states an invariance -- exchangeability says :math:(U,V) and :math:(V,U) have the same law, radial symmetry says :math:(U,V) and :math:(1-U,1-V) do -- so applying that transformation at random to each observation generates exact draws from the null. This is simpler than the multiplier process, needs no derivative estimates, and is valid in finite samples rather than asymptotically. The statistics match R's; the p-values are computed differently and agree only up to Monte-Carlo error.

References

Genest, C., Neslehova, J. and Quessy, J.-F. (2012). Tests of symmetry for bivariate copulas. Annals of the Institute of Statistical Mathematics 64(4), 811-834. Genest, C. and Neslehova, J. G. (2014). On tests of radial symmetry for bivariate copulas. Statistical Papers 55(4), 1107-1119. Genest, C., Remillard, B. (2004). Test of independence and randomness based on the empirical copula process. Test 13(2), 335-370. Kojadinovic, I. (2014). Some copula inference procedures adapted to the presence of ties. Overview of nonparametric tests of extreme-value dependence. International Statistical Review.

TestResult

Bases: NamedTuple

The result of a specification test: the score, its p-value, and what was tested.

Returned by :func:exch_test, :func:rad_sym_test, :func:indep_test and :func:ev_test; you do not normally build one yourself. A small pvalue (say below 0.05) is evidence against the property named in null.

It is a NamedTuple exposing statistic and pvalue, as scipy.stats tests do, so stat, p, *rest = result works. The constructor takes the attributes below, in order, as its arguments.

Attributes:

Name	Type	Description
statistic	float	The observed test statistic, zero or positive; larger means further from the null hypothesis.
pvalue	float	Randomisation or bootstrap p-value, strictly between 0 and 1.
null	str	Plain name of the null hypothesis tested, e.g. "independence", "exchangeability", "radial symmetry" or "extreme-value dependence".
n_rep	int	Number of randomisation or bootstrap replicates used.

DependogramResult dataclass

DependogramResult(subsets: list[tuple[int, ...]], statistics: NDArray[float64], pvalues: NDArray[float64], global_pvalue: float, n_rep: int)

Independence test results for every group of two or more variables, plus one overall p-value.

Returned by :func:dependogram and :func:serial_indep_test; you do not normally build one yourself. Read global_pvalue for "is anything dependent?", then :meth:significant or :meth:summary for "which groups?". The constructor takes the attributes below as keyword or positional arguments, in order. The object is immutable (a frozen dataclass).

Attributes:

Name	Type	Description
subsets	list of tuple of int	Every subset of column indices (0-based) of size two or more, in the order the statistics are given -- :math:2^d - d - 1 of them for d columns. For example (0, 2) is the pair of columns 0 and 2.
statistics	numpy.ndarray of float, shape (len(subsets),)	The Cramer-von Mises statistic for each subset; larger means more dependence.
pvalues	numpy.ndarray of float, shape (len(subsets),)	Permutation p-values, one per subset, strictly between 0 and 1. These are not corrected for testing many subsets at once.
global_pvalue	float	A single test combining all subsets, by the minimum p-value with a Bonferroni-style correction through the permutation distribution -- so it is exact rather than conservative.
n_rep	int	Number of permutations used.

significant

significant(level: float = 0.05) -> list[tuple[int, ...]]

List the groups of variables whose individual p-value is below level.

Parameters:

Name	Type	Description	Default
level	float	Significance threshold, between 0 and 1. A subset is reported when its p-value is strictly below it. No multiple-testing correction is applied; use global_pvalue for a single corrected answer.	0.05

Returns:

Type	Description
list of tuple of int	The significant subsets, in the same order as subsets. Empty if none is significant.

summary

summary(level: float = 0.05) -> str

Build a plain-text table of every subset's statistic and p-value, most significant first.

Parameters:

Name	Type	Description	Default
level	float	Subsets with a p-value strictly below this are marked with *.	0.05

Returns:

Type	Description
str	A multi-line report ending with the global p-value; pass it to print.

exch_test

exch_test(data: ArrayLike, n_rep: int = 1000, random_state: Generator | int | None = None, ties_method: str = 'average') -> TestResult

Check whether swapping the two variables leaves their dependence unchanged.

This is a test of exchangeability (R's exchTest): does the pair (X, Y) depend on each other in the same way as (Y, X)? A small p-value (say below 0.05) means the dependence is lopsided, which rules out most standard copula families at once.

Rejecting is informative out of proportion to its cost: every Archimedean and every elliptical copula is exchangeable, so a rejection eliminates all of them and directs you to asymmetric constructions such as the Khoudraji device.

Parameters:

Name	Type	Description	Default
data	array_like of float, shape (n, 2)	Bivariate observations on any scale; they are converted to ranks first. Exactly two columns.	required
n_rep	int	Number of random replicates used to build the null distribution, a positive integer.	1000
random_state	(int, Generator or None)	Seed or generator, for reproducible p-values.	None
ties_method	str	How tied values are ranked; passed to :func:~rcopula.pseudo_obs.	"average"

Returns:

Type	Description
TestResult	With null="exchangeability".

Raises:

Type	Description
ValueError	If data does not have exactly two columns.

Notes

:math:H_0: C(u,v) = C(v,u), tested with :math:S_n = \sum_i \{C_n(U_i, V_i) - C_n(V_i, U_i)\}^2. The null is generated by swapping the coordinates of a random half of the observations, which is an exact draw from the null when it holds.

Examples:

An exchangeable copula is not rejected:

>>> from rcopula import ClaytonCopula, exch_test
>>> x = ClaytonCopula(3.0).rvs(300, random_state=0)
>>> bool(exch_test(x, n_rep=200, random_state=1).pvalue > 0.05)
True

A deliberately asymmetric construction is:

>>> import numpy as np
>>> rng = np.random.default_rng(0)
>>> u = ClaytonCopula(6.0).rvs(400, random_state=2)
>>> u[:, 1] = np.where(u[:, 0] > 0.5, rng.uniform(size=400), u[:, 1])
>>> bool(exch_test(u, n_rep=200, random_state=1).pvalue < 0.05)
True

rad_sym_test

rad_sym_test(data: ArrayLike, n_rep: int = 1000, random_state: Generator | int | None = None, ties_method: str = 'average') -> TestResult

Check whether two variables are as strongly linked in joint crashes as in joint booms.

This is a test of radial symmetry (R's radSymTest): does the dependence look the same if you flip both variables upside down? A small p-value (say below 0.05) means the lower and upper tails behave differently.

This is the test that separates the families whose tails behave alike from those whose do not. Gaussian, Student-t and Frank are radially symmetric; Clayton, Gumbel and Joe are not. For risk work the asymmetry is the modelling question -- a copula symmetric in the tails cannot represent "crashes cluster but rallies do not".

Parameters:

Name	Type	Description	Default
data	array_like of float, shape (n, 2)	Bivariate observations on any scale; they are converted to ranks first. Exactly two columns.	required
n_rep	int	Number of random replicates used to build the null distribution, a positive integer.	1000
random_state	(int, Generator or None)	Seed or generator, for reproducible p-values.	None
ties_method	str	How tied values are ranked; passed to :func:~rcopula.pseudo_obs.	"average"

Returns:

Type	Description
TestResult	With null="radial symmetry".

Raises:

Type	Description
ValueError	If data does not have exactly two columns.

Notes

:math:H_0: :math:(U,V) and :math:(1-U,1-V) have the same distribution, equivalently :math:C equals its own survival copula. The null is generated by reflecting a random half of the observations.

Examples:

Frank is radially symmetric:

>>> from rcopula import FrankCopula, rad_sym_test
>>> x = FrankCopula(5.0).rvs(400, random_state=0)
>>> bool(rad_sym_test(x, n_rep=200, random_state=1).pvalue > 0.05)
True

Clayton, with its lower-tail dependence, is not:

>>> from rcopula import ClaytonCopula
>>> x = ClaytonCopula(6.0).rvs(400, random_state=0)
>>> bool(rad_sym_test(x, n_rep=200, random_state=1).pvalue < 0.05)
True

indep_test

indep_test(data: ArrayLike, n_rep: int = 1000, random_state: Generator | int | None = None, ties_method: str = 'average') -> TestResult

Check whether the variables are related to each other at all.

This is a rank-based test of independence (R's indepTest). It detects any kind of relationship, not only linear correlation. A small p-value (say below 0.05) means there is dependence worth modelling; a large one means the data give no evidence of it.

:math:H_0: C = \Pi, tested with a Cramer-von Mises distance between the empirical copula and the independence copula.

The null is generated by permuting each column independently, which destroys any dependence while leaving the margins untouched -- an exact randomisation test rather than an asymptotic approximation. R's indepTestSim instead simulates the limiting distribution, which costs :math:O(N n^2 p) time and :math:O(n^2 p) memory and becomes infeasible for large n; permutation has neither problem.

Works in any dimension, unlike the other tests in this module. To find out which variables are dependent, use :func:dependogram.

Parameters:

Name	Type	Description	Default
data	array_like of float, shape (n, d)	Observations on any scale, d >= 2 columns; they are converted to ranks first.	required
n_rep	int	Number of random permutations used to build the null distribution, a positive integer.	1000
random_state	(int, Generator or None)	Seed or generator, for reproducible p-values.	None
ties_method	str	How tied values are ranked; passed to :func:~rcopula.pseudo_obs.	"average"

Returns:

Type	Description
TestResult	With null="independence".

Examples:

>>> import numpy as np
>>> from rcopula import ClaytonCopula, IndependenceCopula, indep_test
>>> rng = np.random.default_rng(0)
>>> indep = IndependenceCopula(3).rvs(400, random_state=0)
>>> bool(indep_test(indep, n_rep=200, random_state=1).pvalue > 0.05)
True
>>> dep = ClaytonCopula(1.0, dim=3).rvs(400, random_state=0)
>>> bool(indep_test(dep, n_rep=200, random_state=1).pvalue < 0.05)
True

It detects dependence that Pearson correlation misses entirely -- the textbook demonstration that "uncorrelated" and "independent" are different claims:

>>> x = np.random.default_rng(7).normal(size=2000)
>>> quadratic = np.column_stack([x, x**2])
>>> bool(abs(np.corrcoef(x, x**2)[0, 1]) < 0.05)
True
>>> bool(indep_test(quadratic, n_rep=200, random_state=1).pvalue < 0.01)
True

ev_test

ev_test(data: ArrayLike, n_rep: int = 1000, random_state: Generator | int | None = None, ties_method: str = 'average') -> TestResult

Check whether the dependence has the shape typical of maxima, such as yearly peak floods.

This tests whether the copula is an extreme-value copula (R's evTestC). Use it before fitting an extreme-value family (Gumbel, Galambos, Husler-Reiss, ...) to block maxima or other componentwise extremes. A small p-value (say below 0.05) means an extreme-value copula is not a good description.

:math:H_0: :math:C is max-stable, :math:C(u^{1/m}, v^{1/m})^m = C(u,v).

The null distribution comes from a parametric bootstrap through the Gumbel family: Gumbel is extreme-value, so resampling from a Gumbel fitted to the data's Kendall tau gives draws satisfying the null with comparable dependence strength. This is cruder than R's multiplier construction and is best read as indicative.

Parameters:

Name	Type	Description	Default
data	array_like of float, shape (n, 2)	Bivariate observations on any scale; they are converted to ranks first. Exactly two columns. Extreme-value copulas cannot be negatively dependent; a sample Kendall's tau at or below zero triggers a warning (see Warns) and the null is built at near independence.	required
n_rep	int	Number of bootstrap replicates, a positive integer.	1000
random_state	(int, Generator or None)	Seed or generator, for reproducible p-values.	None
ties_method	str	How tied values are ranked; passed to :func:~rcopula.pseudo_obs.	"average"

Returns:

Type	Description
TestResult	With null="extreme-value dependence".

Raises:

Type	Description
ValueError	If data does not have exactly two columns.

Warns:

Type	Description
UserWarning	If the sample Kendall's tau is zero or negative. The Gumbel null is then built at tau = 1e-6 (near independence -- the only extreme-value copula without positive dependence), so the test can only tell whether the data look independent; clearly negative dependence already rules out an extreme-value copula.

Examples:

Gumbel is an extreme-value copula; Clayton is not:

>>> from rcopula import ClaytonCopula, GumbelCopula, ev_test
>>> x = GumbelCopula(3.0).rvs(400, random_state=0)
>>> bool(ev_test(x, n_rep=200, random_state=1).pvalue > 0.05)
True
>>> x = ClaytonCopula(5.0).rvs(400, random_state=0)
>>> bool(ev_test(x, n_rep=200, random_state=1).pvalue < 0.05)
True

dependogram

dependogram(data: ArrayLike, n_rep: int = 500, random_state: Generator | int | None = None, ties_method: str = 'average') -> DependogramResult

Find out which groups of variables are related, by testing every pair, triple and so on.

This runs an independence test on each subset of two or more columns (R's dependogram) and also reports one overall p-value.

:func:indep_test answers "are these independent?" with one number. This answers "and which of them are not?", by decomposing the empirical copula process into one component per subset (Genest and Remillard 2004). The components are asymptotically independent under the null, which is what makes reading them separately legitimate.

It routinely finds structure a global test cannot report: three variables can be pairwise independent and jointly dependent, and only the three-element subset shows it.

Parameters:

Name	Type	Description	Default
data	array_like of float, shape (n, d)	Observations on any scale, d >= 2 columns; they are converted to ranks first.	required
n_rep	int	Number of permutations, a positive integer. The null is generated by permuting each column independently, which is exact rather than asymptotic.	500
random_state	(int, Generator or None)	Seed or generator, for reproducible p-values.	None
ties_method	str	How tied values are ranked; passed to :func:~rcopula.pseudo_obs.	"average"

Returns:

Type	Description
DependogramResult	One statistic and p-value per subset, plus global_pvalue.

Raises:

Type	Description
ValueError	If data has fewer than 2 columns.

Notes

Cost is :math:O(n^2 2^d) per permutation, because each subset needs an :math:n \times n matrix product. That is fine to a few hundred observations in a handful of dimensions and prohibitive beyond; subsample rather than wait.

Ties cost power. The statistic is built from the empirical copula, and heavily tied margins make that coarse. The same pairwise-independent construction below is detected easily with continuous margins and not at all when the coordinates take two values each -- so for count or binary data, read a non-significant subset as uninformative rather than as evidence of independence.

Examples:

Pairwise independent and jointly dependent -- the case a global test can detect but not locate. With :math:Z = (X + Y) \bmod 1 every pair is exactly independent, while the triple is deterministic:

>>> import numpy as np
>>> from rcopula.htest import dependogram
>>> rng = np.random.default_rng(0)
>>> x, y = rng.uniform(size=(2, 400))
>>> data = np.column_stack([x, y, (x + y) % 1.0])
>>> result = dependogram(data, n_rep=200, random_state=1)
>>> result.significant()
[(0, 1, 2)]

Genuinely independent data leaves every subset alone:

>>> quiet = dependogram(rng.uniform(size=(300, 3)), n_rep=200, random_state=1)
>>> quiet.significant()
[]

serial_indep_test

serial_indep_test(x: ArrayLike, lags: int = 3, n_rep: int = 500, random_state: Generator | int | None = None, ties_method: str = 'average') -> DependogramResult

Check whether a time series' past values tell you anything about its next value.

This is a rank-based test of serial independence (R's serialIndepTest). A small global_pvalue (say below 0.05) means the observations are not independent over time -- for example, model residuals that still contain structure.

Embeds the series into its own lags and runs :func:dependogram on the result, so the answer names which lag structure carries the dependence rather than only whether some does.

The subsets are read as lag sets: (0, 1) is the dependence between consecutive observations, (0, 2) between observations two apart, and (0, 1, 2) a three-way structure that neither pair reveals -- which is exactly what a correlogram cannot show.

Being rank-based, it is invariant to any increasing transform of the series and needs no moments, so it applies to heavy-tailed returns where the autocorrelation function is not even well defined.

Parameters:

Name	Type	Description	Default
x	array_like of float, shape (n,)	A univariate series, in time order, oldest first. Any array is flattened. Needs more than lags values.	required
lags	int	Embedding dimension: how many consecutive values are looked at together, at least 2. 3 looks at lags 1 and 2. Cost grows as :math:2^{lags}.	3
n_rep	int	Number of permutations for the null, a positive integer.	500
random_state	(int, Generator or None)	Seed or generator, for reproducible p-values.	None
ties_method	str	How tied values are ranked; passed to :func:~rcopula.pseudo_obs.	"average"

Returns:

Type	Description
DependogramResult	Subsets are indices into the embedded vector, so 0 is the current observation and k the one k steps back.

Raises:

Type	Description
ValueError	If lags < 2, or if the series has lags or fewer values.

Notes

Permuting the embedded matrix column by column destroys serial dependence while preserving each column's marginal distribution, which is the exact randomisation this needs. It does not preserve the overlap between rows that embedding creates, so the test is mildly conservative -- the rows are not independent, and the permutation null treats them as if they were.

Examples:

White noise passes; an AR(1) does not:

>>> import numpy as np
>>> from rcopula.htest import serial_indep_test
>>> rng = np.random.default_rng(0)
>>> noise = rng.standard_normal(600)
>>> bool(serial_indep_test(noise, n_rep=200, random_state=1).global_pvalue > 0.05)
True
>>> series = np.zeros(600)
>>> for t in range(1, 600):
...     series[t] = 0.6 * series[t - 1] + rng.standard_normal()
>>> bool(serial_indep_test(series, n_rep=200, random_state=1).global_pvalue < 0.05)
True

A GARCH-like series is uncorrelated in level and dependent in magnitude, which is the case an autocorrelation function reports as clean:

>>> volatility = np.exp(0.5 * rng.standard_normal(2000).cumsum() / 30)
>>> returns = volatility * rng.standard_normal(2000)
>>> bool(abs(np.corrcoef(returns[:-1], returns[1:])[0, 1]) < 0.06)
True
>>> bool(serial_indep_test(np.abs(returns), n_rep=200, random_state=1).global_pvalue < 0.05)
True