diff --git a/Project.toml b/Project.toml index cc3e81fc6..ac15aa709 100644 --- a/Project.toml +++ b/Project.toml @@ -21,6 +21,7 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" +StatsAPI = "82ae8749-77ed-4fe6-ae5f-f523153014b0" StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" StatsFuns = "4c63d2b9-4356-54db-8cca-17b64c39e42c" TaylorSeries = "6aa5eb33-94cf-58f4-a9d0-e4b2c4fc25ea" @@ -60,6 +61,7 @@ Roots = "1, 2, 3" SpecialFunctions = "2" StableRNGs = "1" Statistics = "1" +StatsAPI = "1.8.0" StatsBase = "0.33, 0.34" StatsFuns = "0.9, 1.3, 2" TaylorSeries = "0.20, 0.21, 0.22" diff --git a/docs/make.jl b/docs/make.jl index 000a366f8..288c85f36 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -27,6 +27,7 @@ makedocs(; "Conditioning and subsetting"=>"manual/conditioning_and_subsetting.md", "Dependence metrics"=>"manual/dependence_measures.md", "Fitting"=>"manual/fitting_interface.md", + "Hypothesis testing" => "manual/hypothesis_testing.md", "Visualizations"=>"manual/visualizations.md", ], "Bestiary" => [ diff --git a/docs/src/assets/references.bib b/docs/src/assets/references.bib index bcf53f050..0618abdbd 100644 --- a/docs/src/assets/references.bib +++ b/docs/src/assets/references.bib @@ -1064,4 +1064,124 @@ @article{gudendorf2012multivariate year={2012}, publisher={Elsevier}, doi={10.1016/j.jspi.2012.05.007} -} +} + +@article{genest2004independence, + title = {Test of Independence and Randomness Based on the Empirical Copula Process}, + author = {Genest, Christian and Rémillard, Bruno}, + year = {2004}, + journal = {TEST}, + volume = {13}, + number = {2}, + pages = {335--369}, + doi = {10.1007/BF02595777}, + keywords = {copula, empirical copula process, independence, Cramer-von Mises} +} + +@article{fermanian2004empirical, + title = {Weak Convergence of Empirical Copula Processes}, + author = {Fermanian, Jean-David and Radulović, Dragan and Wegkamp, Marten}, + year = {2004}, + journal = {Bernoulli}, + volume = {10}, + number = {5}, + pages = {847--860}, + doi = {10.3150/bj/1099579158}, + keywords = {copula, empirical process, weak convergence} +} + +@article{genest2012symmetry, + title = {Tests of Symmetry for Bivariate Copulas}, + author = {Genest, Christian and Nešlehová, Johanna and Quessy, Jean-François}, + year = {2012}, + journal = {Annals of the Institute of Statistical Mathematics}, + volume = {64}, + number = {4}, + pages = {811--834}, + doi = {10.1007/s10463-011-0337-6}, + keywords = {copula, symmetry, exchangeability, empirical copula process} +} + +@article{harder2017exchangeability, + title = {Testing Exchangeability of Copulas in Arbitrary Dimension}, + author = {Harder, Michael and Stadtm{"u}ller, Ulrich}, + year = {2017}, + journal = {Journal of Nonparametric Statistics}, + volume = {29}, + number = {1}, + pages = {40--60}, + doi = {10.1080/10485252.2016.1253841}, + keywords = {copula, exchangeability, multiplier bootstrap, empirical copula} +} + +@article{beare2020symmetry, + title = {Randomization Tests of Copula Symmetry}, + author = {Beare, Brendan K. and Seo, Juwon}, + year = {2020}, + journal = {Econometric Theory}, + volume = {36}, + number = {6}, + pages = {1025--1063}, + doi = {10.1017/S0266466619000410}, + keywords = {copula, radial symmetry, exchangeability, randomization test} +} + +@article{kojadinovic2011extremevalue, + title = {Large-Sample Tests of Extreme-Value Dependence for Multivariate Copulas}, + author = {Kojadinovic, Ivan and Segers, Johan and Yan, Jun}, + year = {2011}, + journal = {Canadian Journal of Statistics}, + volume = {39}, + number = {4}, + pages = {703--720}, + doi = {10.1002/cjs.10110}, + keywords = {copula, extreme value, max-stability, multiplier bootstrap} +} + +@article{remillard2009equality, + title = {Testing for Equality Between Two Copulas}, + author = {Rémillard, Bruno and Scaillet, Olivier}, + year = {2009}, + journal = {Journal of Multivariate Analysis}, + volume = {100}, + number = {3}, + pages = {377--386}, + doi = {10.1016/j.jmva.2008.05.004}, + keywords = {copula, empirical process, multiplier bootstrap, Cramer-von Mises} +} + +@article{bucher2010bootstrap, + title = {A Note on Bootstrap Approximations for the Empirical Copula Process}, + author = {Bücher, Axel and Dette, Holger}, + year = {2010}, + journal = {Statistics \& Probability Letters}, + volume = {80}, + number = {23--24}, + pages = {1925--1932}, + doi = {10.1016/j.spl.2010.08.021}, + keywords = {copula, empirical process, multiplier bootstrap} +} + +@article{genest2008bootstrap, + title = {Validity of the Parametric Bootstrap for Goodness-of-Fit Testing in Semiparametric Models}, + author = {Genest, Christian and Rémillard, Bruno}, + year = {2008}, + journal = {Annales de l'Institut Henri Poincaré, Probabilités et Statistiques}, + volume = {44}, + number = {6}, + pages = {1096--1127}, + doi = {10.1214/07-AIHP148}, + keywords = {bootstrap, goodness of fit, copula, semiparametric model} +} + +@article{genest2009gof, + title = {Goodness-of-Fit Tests for Copulas: A Review and a Power Study}, + author = {Genest, Christian and Rémillard, Bruno and Beaudoin, David}, + year = {2009}, + journal = {Insurance: Mathematics and Economics}, + volume = {44}, + number = {2}, + pages = {199--213}, + doi = {10.1016/j.insmatheco.2007.10.005}, + keywords = {copula, goodness of fit, parametric bootstrap, Cramer-von Mises} +} diff --git a/docs/src/manual/developer_guide.md b/docs/src/manual/developer_guide.md index bb43a0e91..4192edd73 100644 --- a/docs/src/manual/developer_guide.md +++ b/docs/src/manual/developer_guide.md @@ -246,7 +246,435 @@ Once the above methods are implemented, your family becomes automatically compat the requires _fit function for :mle for example might already be provided by the package, test your case. +## 1.6 Hypothesis-testing interface + +The hypothesis-testing machinery is dispatch-oriented, similarly to the fitting +interface. The user-facing testing API is deliberately small; the underscore- +prefixed hooks described below are **internal contributor interfaces** used to +implement tests inside `Copulas.jl`. They are not public extension points and +are not covered by SemVer, so their signatures may evolve as the framework +develops. + +```text +Hypothesis × Statistic × Calibration + ↓ + CopulaTest +``` + +The role of each component is deliberately separated: + +* a `CopulaHypothesis` describes the mathematical null hypothesis; +* `_teststatistic` defines how a statistic is computed; +* `_calibrate` or one of the generic calibration hooks defines how its null distribution is approximated; +* `CopulaTest` stores the common result and implements the `StatsAPI.HypothesisTest` interface. + +The generic constructor knows nothing about independence, exchangeability, max-stability, goodness of fit, or any other particular statistical problem. Those properties enter exclusively through Julia dispatch. + +### Relationship with the fitting interface + +The design intentionally mirrors the existing fitting API. Fitting uses + +```text +_available_fitting_methods(CT, d) + ↓ + first = default + ↓ +_find_method(...) + ↓ +_fit(CT, U, Val(method)) +``` + +Hypothesis testing uses + +```text +_available_statistics(h) + ↓ + first = default + ↓ +_find_statistic(...) + ↓ +_teststatistic(h, Val(statistic), U) +``` + +and subsequently + +```text +_available_calibrations(h, Val(statistic)) + ↓ + first = default + ↓ +_find_calibration(...) + ↓ +_calibrate(h, Val(calibration), Val(statistic),...) +``` + +Within the internal testing machinery, **the order returned by `_available_statistics` and `_available_calibrations` is significant**. The first entry is the default used when the corresponding public keyword is `:default`. + +For example, + +```julia +_available_statistics(::MyHypothesis) = (:cvm, :ks) +``` + +means that both statistics are supported, while `:cvm` is the default. Likewise, + +```julia +_available_calibrations(::MyHypothesis, ::Val{:cvm},) = (:multiplier, :simulation) +``` + +means that both calibrations are valid for `:cvm`, with `:multiplier` selected by default. + +### Defining a new hypothesis + +A new null hypothesis starts by subtyping [`CopulaHypothesis`](@ref): + +```julia +struct MyHypothesis <: CopulaHypothesis end +``` + +The hypothesis should then define a human-readable test name and null hypothesis: + +```julia +Copulas.testname(::MyHypothesis) = "My copula hypothesis test" + +Copulas.nullhypothesis(::MyHypothesis) = "The null hypothesis holds." +``` + +`testname` and `nullhypothesis` are extension hooks used by the generic display machinery. `testname` is intentionally not exported, so contributors should extend it through the qualified name `Copulas.testname`. + +Next declare the available statistics: + +```julia +Copulas._available_statistics(::MyHypothesis) = (:Sn, :ks) +``` + +and the available calibrations for each statistic: + +```julia +Copulas._available_calibrations(::MyHypothesis, ::Val{:Sn},) = (:simulation,) + +Copulas._available_calibrations(::MyHypothesis, ::Val{:ks},) = (:simulation,) +``` + +Finally implement the mathematical statistics: + +```julia +function Copulas._teststatistic(::MyHypothesis, ::Val{:Sn}, U::AbstractMatrix; kwargs...,) + # Compute Sₙ. +end + +function Copulas._teststatistic(::MyHypothesis, ::Val{:ks}, U::AbstractMatrix; kwargs...,) + # Compute the KS statistic. +end +``` + +At this point no constructor-specific branch has been introduced. Selection is performed entirely by dispatch on `Val{:Sn}` or `Val{:ks}`. + +### Adding another statistic to an existing hypothesis + +Suppose an existing hypothesis currently declares + +```julia +_available_statistics(::SomeHypothesis) = (:cvm,) +``` + +and a contributor implements a new Kolmogorov--Smirnov statistic. The capability declaration becomes + +```julia +_available_statistics(::SomeHypothesis) = (:cvm, :ks) +``` + +and the new statistic is added independently: + +```julia +function _teststatistic(::SomeHypothesis, ::Val{:ks}, U::AbstractMatrix; kwargs...,) + # ... +end +``` + +Its valid calibrations are then declared: + +```julia +_available_calibrations(::SomeHypothesis, ::Val{:ks},) = (:simulation,) +``` + +The existing public constructor automatically accepts + +```julia +SomeCopulaTest(U; statistic=:ks,) +``` + +without any modification to `CopulaTest`, `show`, or a central routing table. If the new statistic should become the default, simply change the order: + +```julia +_available_statistics(::SomeHypothesis) = (:ks, :cvm) +``` + +The same convention is already used by `_available_fitting_methods`. + +### Generic calibration engines + +A contributor should normally reuse one of the generic calibration engines rather than reimplementing resampling loops. + +The currently available engines are: + +| Calibration | Required hypothesis-specific hook | Purpose | +| ----------------------- | ---------------------------------------- | ------------------------------------------------ | +| `:simulation` | `_simulation_sample(h, U, rng)` | Generate a sample directly under `H_0` | +| `:randomization` | `_randomization_sample(h, U, rng)` | Apply a null-invariant random transformation | +| `:multiplier` | `_multiplier_representation(h, stat, U)` | Supply the empirical-process representation | +| `:parametric_bootstrap` | `_bootstrap_copula(h)` | Supply the copula used for parametric simulation | + +For parametric-bootstrap tests, a hypothesis may additionally implement + +```julia +_bootstrap_hypothesis(h, Ustar) +``` + +when something must be recomputed from every bootstrap sample. For example, composite goodness-of-fit testing uses this hook to refit the model in every bootstrap replicate. + +--- + +### Simulation calibration + +If the null distribution can be simulated directly, declare + +```julia +_available_calibrations(::MyHypothesis, ::Val{:Sn},) = (:simulation,) +``` + +and implement only + +```julia +function _simulation_sample(::MyHypothesis, U::AbstractMatrix, rng::Distributions.AbstractRNG,) + d, n = size(U) + + # Generate a d × n sample under H₀. + sample = ... + + return sample +end +``` + +The generic engine handles: + +* repetition over `N` resamples; +* conversion to pseudo-observations; +* recomputation of the selected statistic; +* exceedance counting; +* p-value construction; +* storage of the actual number of resamples. + +--- + +### Randomization calibration + +For a hypothesis characterized by a transformation group that leaves the null distribution invariant, implement + +```julia +function _randomization_sample(::MyHypothesis, U::AbstractMatrix, rng::Distributions.AbstractRNG,) + # Randomly transform U under H₀. +end +``` + +Optional information printed or stored with the result can be provided through + +```julia +_randomization_details(::MyHypothesis) = (; some_property=value) +``` + +The radial-symmetry test is an example: every observation is independently kept or radially reflected with probability `1/2`. + +--- + +### Multiplier calibration + +Statistics based on an empirical-copula process can use the generic multiplier engine by implementing + +```julia +function _multiplier_representation(h::MyHypothesis, ::Val{:Sn}, U::AbstractMatrix,) + # Construct the empirical-process representation. + return (matrices=..., scale=..., details=(; ...),) +end +``` + +The returned `NamedTuple` may contain + +| Field | Meaning | +| ------------ | ----------------------------------------------------- | +| `matrices` | Linear representations used by the multiplier process | +| `scale` | Final scale applied to the bootstrap statistic | +| `weights` | Optional observation weights | +| `strict` | Whether exceedance means `>` instead of `≥` | +| `correction` | Optional Monte Carlo p-value correction | +| `details` | Metadata stored in the resulting `CopulaTest` | + +The generic multiplier engine generates centered exponential multipliers and performs the repeated linear-algebra operations. + +The hypothesis-specific method should therefore describe the **mathematical empirical-process representation**, not the mechanics of bootstrap iteration. + +--- + +### Parametric bootstrap calibration + +For a parametric null hypothesis, implement + +```julia +_bootstrap_copula(h::MyHypothesis) = ... +``` + +and declare + +```julia +_available_calibrations(::MyHypothesis, ::Val{:Sn},) = (:parametric_bootstrap,) +``` + +The generic calibration engine repeatedly generates + +```math +\boldsymbol U_1^\star,\ldots,\boldsymbol U_n^\star +\sim C_0, +``` + +converts them to pseudo-observations, constructs the bootstrap hypothesis, and recomputes the statistic. + +The default bootstrap hypothesis is unchanged: + +```julia +_bootstrap_hypothesis(h::CopulaHypothesis, Ustar::AbstractMatrix,) = h +``` + +Override this only when the null model must be re-estimated. For example, a composite goodness-of-fit hypothesis uses + +```julia +function _bootstrap_hypothesis(h::GoodnessOfFitHypothesis{<:CopulaModel}, Ustar::AbstractMatrix,) + # Refit the same copula family on Ustar. +end +``` + +This keeps parameter-estimation uncertainty inside the bootstrap rather than treating fitted parameters as fixed. + +--- + +### Implementing a new calibration mechanism + +If none of the reusable engines is appropriate, a new calibration can be added through dispatch: + +```julia +function _calibrate(h::MyHypothesis, ::Val{:my_calibration}, stat::Val, U::AbstractMatrix, observed::Real; kwargs...,) + # ... + return p, n_resamples, details +end +``` + +The return contract is always + +```julia +(p, n_resamples, details) +``` + +where + +* `p` is the resulting p-value; +* `n_resamples` is the actual number of resampling replicates used; +* `details` is a `NamedTuple` containing calibration-specific metadata. + +A non-resampling calibration, such as a future analytical or asymptotic calibration, can therefore return + +```julia +return p, 0, (;) +``` + +without pretending that the constructor's `N` value was used. Validation of `N` belongs to resampling calibrations, not to the generic `CopulaTest` constructor. + +--- + +### Result interface + +Every procedure ultimately returns the same type: + +```julia +CopulaTest{H<:CopulaHypothesis, S<:Real, P<:Real, D<:NamedTuple} +``` + +The common API is + +```julia +teststatistic(test) +pvalue(test) +StatsBase.nobs(test) +Copulas.testname(test) +Copulas.nullhypothesis(test) +``` + +and the result contains + +```julia +test.hypothesis +test.dimension +test.statistic +test.calibration +test.n_resamples +test.details +``` + +This common result type is the reason a contributor should generally add behavior through `CopulaHypothesis`, statistic, and calibration dispatch rather than introducing a new result struct. + +--- + +### Generic display + +`Base.show(::MIME"text/plain", test::CopulaTest)` is generic. + +A new hypothesis therefore receives the common output automatically from + +```julia +testname(h) +nullhypothesis(h) +teststatistic(test) +pvalue(test) +``` + +without modifying `show.jl`. + +When a hypothesis needs additional output, extend the corresponding display hook rather than adding hypothesis-specific branching to the generic printer. + +The guiding rule is: + +> Mathematical differences should be expressed by dispatch; common mechanics +> should remain generic. + +--- + +### Minimal extension checklist + +To add a new hypothesis using an existing calibration engine, the usual minimum +is: + +```julia +struct MyHypothesis <: CopulaHypothesis end + +Copulas.testname(::MyHypothesis) = "..." +Copulas.nullhypothesis(::MyHypothesis) = "..." + +Copulas._available_statistics(::MyHypothesis) = (:my_statistic,) + +Copulas._available_calibrations(::MyHypothesis, ::Val{:my_statistic},) = (:simulation,) + +Copulas._teststatistic(::MyHypothesis, ::Val{:my_statistic}, U; kwargs...,) = ... + +Copulas._simulation_sample(::MyHypothesis, U, rng,) = ... +``` + +After those methods are defined, + +```julia +CopulaTest(MyHypothesis(), U; N=1000,) +``` + +uses the complete common infrastructure automatically. +No modification of the generic constructor is required. # 2. Specific sub-APIs Some families of copulas in `Copulas.jl` have additional internal structures or specific mathematical representations. diff --git a/docs/src/manual/hypothesis_testing.md b/docs/src/manual/hypothesis_testing.md new file mode 100644 index 000000000..44ad69df7 --- /dev/null +++ b/docs/src/manual/hypothesis_testing.md @@ -0,0 +1,1026 @@ +```@meta +CurrentModule = Copulas +``` + +# [Hypothesis testing](@id hypothesis_testing) + +`Copulas.jl` provides a common interface for rank-based hypothesis tests on copulas. The framework separates three distinct ingredients: + +1. the **null hypothesis** being tested; +2. the **test statistic** used to measure departures from the null; +3. the **calibration method** used to obtain a p-value. + +Conceptually, + +```text +CopulaHypothesis + × + Statistic + × + Calibration + ↓ + CopulaTest +``` + +This separation makes it possible to reuse the same calibration machinery across different hypotheses and to introduce new statistics or hypotheses without modifying the generic test constructor. + +The current implementation includes tests of: + +* mutual independence; +* exchangeability; +* radial symmetry; +* extreme-value dependence (max-stability); +* goodness of fit for a specified copula; +* goodness of fit for a fitted copula family. + +The procedures are based on empirical-copula processes, resampling, multiplier methods, and parametric bootstrap ideas developed throughout the copula-testing literature; see, among others, + +[genest2004independence](@cite), +[fermanian2004empirical](@cite), +[remillard2009equality](@cite), +[bucher2010bootstrap](@cite), and +[genest2009gof](@cite). + +--- + +## Data convention + +As elsewhere in `Copulas.jl`, observations are represented by a `d\times n` matrix + +```math +U = +\begin{pmatrix} +U_{11} & \cdots & U_{1n}\\ +\vdots & & \vdots\\ +U_{d1} & \cdots & U_{dn} +\end{pmatrix}, +``` + +where each column + +```math +\boldsymbol U_i = (U_{1i},\ldots,U_{di})^\top +``` + +is one `d`-dimensional observation. + +Hypothesis tests are rank based. When + +```julia +pseudo_values=false +``` + +the input matrix is transformed internally with [`pseudos`](@ref). If the data already consist of pseudo-observations in `[0,1]^d`, use + +```julia +pseudo_values=true +``` + +to avoid ranking them again. + +::: warning Continuous margins and ties + +The currently implemented copula hypothesis tests assume continuous margins and therefore require tie-free observations in every margin. Tied or discrete data are rejected with an `ArgumentError`. + +This is intentional: ordinal ranking would otherwise assign distinct ranks to tied observations and could produce apparently valid p-values without the tie-aware empirical-process or bootstrap theory required for such data. Tie-aware procedures are outside the scope of the current implementation. + +::: + +Given pseudo-observations $\boldsymbol U_1,\ldots,\boldsymbol U_n$, the empirical copula is + +```math +C_n(\boldsymbol u) += +\frac{1}{n} +\sum_{i=1}^{n} +\mathbf 1 +\left( +\boldsymbol U_i\le\boldsymbol u +\right), +``` + +where the inequality is understood componentwise. + +Empirical-copula processes and their weak convergence form the theoretical basis for many of the statistics and multiplier approximations used below [fermanian2004empirical](@cite). + +--- + +## Common interface + +All tests return a [`CopulaTest`](@ref), which implements `StatsAPI.HypothesisTest`. + +For example, + +```@example hypothesis_testing +using Copulas, Distributions, Random, StatsBase + +U = rand(Xoshiro(123), ClaytonCopula(2, 3.0), 80) + +test = IndependenceCopulaTest(U; N=49, rng=Xoshiro(123),) + +nothing # hide +``` + +The common result interface is + +```@example hypothesis_testing +teststatistic(test) +``` + +```@example hypothesis_testing +pvalue(test) +``` + +```@example hypothesis_testing +nobs(test) +``` + +A test also records: + +* `test.statistic`: the statistic used; +* `test.calibration`: the calibration method; +* `test.n_resamples`: number of resampling replicates; +* `test.dimension`: dimension of the copula; +* `test.details`: test-specific metadata. + +Printing the object gives a summary of the hypothesis, statistic, calibration, p-value, and relevant test-specific information. + +```@example hypothesis_testing +test +``` + +::: note Number of resamples + +The small values of `N` used in the documentation keep the examples fast. +For statistical work, substantially larger values should generally be used, +depending on the desired Monte Carlo precision. + +::: + +--- + +# Mutual independence + +## Null hypothesis + +Let `C` denote the copula of the random vector. Mutual independence is equivalent to the product copula + +```math +\Pi(\boldsymbol u) += +\prod_{j=1}^{d}u_j. +``` + +Thus + +```math +H_0: +C(\boldsymbol u) += +\Pi(\boldsymbol u) +\qquad +\text{for every } +\boldsymbol u\in[0,1]^d. +``` + +Rank-based independence tests constructed from the empirical copula process are studied in [genest2004independence](@cite). + +## Cramér--von Mises statistic + +The statistic currently available in `Copulas.jl` is `:cvm`. The implementation evaluates the squared discrepancy between the empirical copula and the product copula at the observed pseudo-observations: + +```math +S_n^{\mathrm{ind}} += +\sum_{i=1}^{n} +\left[ +C_n(\boldsymbol U_i) +- +\prod_{j=1}^{d}U_{ji} +\right]^2. +``` + +Large values indicate departure from mutual independence. + +## Calibration + +Under `H_0`, the coordinates are independent uniforms. The default calibration is therefore `:simulation`: + +1. generate `n` observations from the `d`-dimensional product copula; +2. transform the generated sample to pseudo-observations; +3. recompute `S_n^{\mathrm{ind}}`; +4. repeat the procedure `N` times; +5. compare the observed statistic with its simulated null distribution. + +## Usage + +```@example hypothesis_testing +Uind = rand(Xoshiro(1), IndependentCopula(3), 100) + +tind = IndependenceCopulaTest(Uind; N=49, rng=Xoshiro(2),) + +(tind.statistic, tind.calibration, pvalue(tind)) +``` + +The current defaults are + +```text +statistic = :cvm +calibration = :simulation +``` + +--- + +# Exchangeability + +## Null hypothesis + +A copula `C` is exchangeable when it is invariant under permutations of its coordinates. + +For a permutation + +```math +\pi: +\{1,\ldots,d\} +\longrightarrow +\{1,\ldots,d\}, +``` + +write + +```math +\boldsymbol u_\pi += +(u_{\pi(1)},\ldots,u_{\pi(d)}). +``` + +Full exchangeability means + +```math +H_0: +C(\boldsymbol u) += +C(\boldsymbol u_\pi) +``` + +for every `\boldsymbol u\in[0,1]^d` and every coordinate permutation `\pi`. + +Empirical-copula tests for bivariate symmetry were developed in [genest2012symmetry](@cite) and extended to arbitrary dimension by [harder2017exchangeability](@cite). + +## Statistic + +For a collection `\mathcal G` of non-identity permutations, the implemented statistic is + +```math +S_n^{\mathrm{ex}} += +\sum_{\pi\in\mathcal G} +\sum_{i=1}^{n} +\left[ +C_n(\boldsymbol U_i) +- +C_n(\boldsymbol U_{i,\pi}) +\right]^2 +w_\pi(\boldsymbol U_i). +``` + +The default weight is `weight=:wm2`. + +Let + +```math +m(\boldsymbol u) += +\min_{1\le j\le d}u_j, +``` + +and + +```math +b(\boldsymbol u) += +d-1+m(\boldsymbol u)-\sum_{j=1}^{d}u_j. +``` + +For a transposition exchanging coordinates `a` and `b`, define + +```math +\omega_\pi(\boldsymbol u) += +|u_a-u_b|. +``` + +For a general permutation, let + +```math +u_{(1)}\le\cdots\le u_{(d)} +``` + +denote the ordered coordinates and define the implementation's permutation separation term by + +```math +\omega_\pi(\boldsymbol u) += +\sum_{k=\lceil d/2\rceil+1}^{d} +\left( +u_{(k)}-m(\boldsymbol u) +\right). +``` + +The `:wm2` weight is then + +```math +w_\pi(\boldsymbol u) += +\left[ +\max +\left\{ +0, +\min +\left( +m(\boldsymbol u), +\omega_\pi(\boldsymbol u), +b(\boldsymbol u) +\right) +\right\} +\right]^2. +``` + +Alternatively, + +```julia +weight=:none +``` + +sets `w_\pi(\boldsymbol u)=1`. + +## Permutation generators + +The keyword `permutations` controls the set `\mathcal G`. + +### `permutations=:G2` + +This is the default. + +For `d=2`, it contains the only nontrivial transposition, + +```math +(12). +``` + +For `d>2`, it uses the transposition + +```math +(12) +``` + +together with the cyclic left shift + +```math +(12\cdots d). +``` + +### `permutations=:G1` + +Uses the transpositions + +```math +(12),(13),\ldots,(1d). +``` + +### `permutations=:all` + +Uses all non-identity permutations. + +Because the current multiplier implementation materializes one dense `n × n` matrix for every selected permutation, the resolved permutation collection is protected by a memory-cost guard. For `permutations=:all`, the factorial-sized collection is checked before it is materialized; explicit custom collections are checked after validation using their actual resolved size. Problems whose estimated matrix storage exceeds the safety limit raise an `ArgumentError`. For larger dimensions or samples, use `:G1`, `:G2`, or a smaller custom collection of permutations. + +A custom permutation or collection of permutations can also be supplied directly. + +## Multiplier calibration + +The default calibration is `:multiplier`. + +The empirical-copula process has a nontrivial correction caused by replacing the unknown margins with ranks. The implementation therefore constructs the corresponding multiplier representation, including finite-difference estimates of the partial derivatives of $C_n$. + +The derivative bandwidth is + +```math +h_n=n^{-1/2}. +``` + +For coordinate $j$, the derivative is approximated by a boundary-corrected finite difference of the form + +```math +\dot C_{n,j}(\boldsymbol u) +\approx +\frac{ +C_n(\boldsymbol u+h_n\boldsymbol e_j) +- +C_n(\boldsymbol u-h_n\boldsymbol e_j) +}{ +\text{effective width} +}. +``` + +Independent exponential multipliers are generated and centered before applying the empirical-process representation. This type of multiplier approximation is closely related to the methods discussed in [remillard2009equality](@cite), [bucher2010bootstrap](@cite), and [harder2017exchangeability](@cite). + +## Usage + +```@example hypothesis_testing +Uex = rand(Xoshiro(4), GumbelCopula(3, 2.0), 80) + +tex = ExchangeabilityCopulaTest(Uex; permutations=:G2, weight=:wm2, N=49, rng=Xoshiro(5),) + +(tex.statistic, tex.calibration) +``` + +The current defaults are + +```text +statistic = :Sn +calibration = :multiplier +``` + +--- + +# Radial symmetry + +## Null hypothesis + +A copula is radially symmetric when + +```math +\boldsymbol U +\overset{d}{=} +\boldsymbol 1-\boldsymbol U. +``` + +Equivalently, if $C^{\mathrm{rad}}$ denotes the copula of $\boldsymbol 1-\boldsymbol U$, then + +```math +H_0: +C += +C^{\mathrm{rad}}. +``` + +Nonparametric tests of copula symmetry and randomization procedures based on the corresponding invariance group are studied in [beare2020symmetry](@cite). + +## Statistic + +Let $C_n$ denote the empirical copula of the original pseudo-observations and let $\bar{C_n}$ denote the empirical copula constructed from + +```math +\boldsymbol 1-\boldsymbol U_1, +\ldots, +\boldsymbol 1-\boldsymbol U_n. +``` + +The implemented statistic is + +```math +S_n^{\mathrm{rad}} += +\sum_{i=1}^{n} +\left[ +C_n(\boldsymbol U_i) +- +\bar C_n(\boldsymbol U_i) +\right]^2. +``` + +Large values indicate radial asymmetry. + +## Randomization calibration + +Under radial symmetry, an observation and its radial reflection are distributionally equivalent. For every observation $i$, independently generate + +```math +B_i\sim\operatorname{Bernoulli}(1/2), +``` + +and construct + +```math +\boldsymbol U_i^\star += +\begin{cases} +\boldsymbol U_i, +& +B_i=0,\\[2mm] +\boldsymbol 1-\boldsymbol U_i, +& +B_i=1. +\end{cases} +``` + +The randomized sample is converted back to pseudo-observations before the statistic is evaluated. Because radial reflection of rank-grid values can create exact ties even when the original sample is tie-free, ties induced by the randomization itself are reranked using average ranks. This does not relax the requirement that the original input margins be tie-free. + +Thus the default reflection probability is exactly + +```math +\Pr(B_i=1)=\frac12. +``` + +The procedure exploits the group invariance associated with radial symmetry, following the randomization-testing principle developed in [beare2020symmetry](@cite). + +## Usage + +```@example hypothesis_testing +Urad = rand(Xoshiro(6), GaussianCopula([1.0 0.6; 0.6 1.0]), 80) + +trad = RadialSymmetryCopulaTest(Urad; N=49, rng=Xoshiro(7),) + +(trad.statistic, trad.calibration, trad.details.reflection_probability) +``` + +The current defaults are + +```text +statistic = :Sn +calibration = :randomization +``` + +--- + +# Extreme-value dependence + +## Max-stability + +Extreme-value copulas are characterized by max-stability. For any $r>0$, + +```math +C(u_1^r,\ldots,u_d^r) += +C(u_1,\ldots,u_d)^r. +``` + +Equivalently, for $r>1$, + +```math +C(\boldsymbol u) += +C(\boldsymbol u^{1/r})^r, +``` + +where + +```math +\boldsymbol u^{1/r} += +(u_1^{1/r},\ldots,u_d^{1/r}). +``` + +This characterization provides a direct way to test + +```math +H_0: +C\text{ belongs to the extreme-value class}. +``` + +Large-sample tests based on this max-stability identity, the empirical copula, and multiplier approximations are developed by [kojadinovic2011extremevalue](@cite). + +## Statistic + +For a finite collection of powers + +```math +\mathcal R += +\{r_1,\ldots,r_K\}, +\qquad +r_k>1, +``` + +the implemented statistic is + +```math +S_n^{\mathrm{EV}} += +\sum_{r\in\mathcal R} +\sum_{i=1}^{n} +\left[ +C_n(\boldsymbol U_i^{1/r})^r +- +C_n(\boldsymbol U_i) +\right]^2. +``` + +The default powers are + +```math +\mathcal R=\{3,4,5\}. +``` + +They can be changed through the `powers` keyword. + +## Multiplier calibration + +Approximate p-values are obtained from a multiplier representation of the empirical-copula process, following the max-stability testing strategy in +[kojadinovic2011extremevalue](@cite). + +As in the exchangeability test, the finite-difference bandwidth used for the empirical partial derivatives is + +```math +h_n=n^{-1/2}. +``` + +The multiplier variables are exponential and centered before the bootstrap process is evaluated. + +## Usage + +```@example hypothesis_testing +Uev = rand(Xoshiro(8), GumbelCopula(2, 2.5), 80) + +tev = ExtremeValueCopulaTest(Uev; powers=3:5, N=49, rng=Xoshiro(9),) + +(tev.statistic, tev.calibration, tev.details.powers) +``` + +A single power is also allowed: + +```julia +ExtremeValueCopulaTest(U; powers=2) +``` + +All supplied powers must be finite and strictly larger than one. + +The current defaults are + +```text +statistic = :Sn +calibration = :multiplier +``` + +--- + +# Goodness of fit + +Copula goodness-of-fit procedures compare the empirical dependence structure with a proposed parametric copula model. Empirical-process and Cramér--von Mises procedures of this form are reviewed extensively in [genest2009gof](@cite). + +`Copulas.jl` distinguishes a **simple** null hypothesis from a **composite** null hypothesis. + +--- + +## Simple goodness of fit + +Suppose that a fully specified copula `C_0` is given, including all its parameters. + +The null hypothesis is + +```math +H_0: +C=C_0. +``` + +The implemented Cramér--von Mises-type statistic is + +```math +S_n^{\mathrm{GOF}} += +\sum_{i=1}^{n} +\left[ +C_n(\boldsymbol U_i) +- +C_0(\boldsymbol U_i) +\right]^2. +``` + +Use + +```@example hypothesis_testing +C0 = ClaytonCopula(2, 3.0) +Ugof = rand(Xoshiro(10), C0, 80) + +tsimple = GOFCopulaTest(C0, Ugof; N=49, rng=Xoshiro(11),) + +tsimple.hypothesis.kind +``` + +which produces a `:simple` goodness-of-fit hypothesis. + +### Parametric bootstrap + +For every bootstrap replicate: + +1. generate `n` observations from `C_0`; +2. transform the sample to pseudo-observations; +3. compute the same goodness-of-fit statistic; +4. compare the bootstrap statistic with the observed value. + +No parameters are re-estimated because `C_0` is fully specified. + +--- + +## Composite goodness of fit + +Suppose instead that + +```math +\mathcal C += +\{ +C_\theta:\theta\in\Theta +\} +``` + +is a parametric copula family and that $\widehat\theta$ is estimated from the data. + +The null hypothesis becomes + +```math +H_0: +C\in\mathcal C, +``` + +and the observed statistic is + +```math +S_n^{\mathrm{GOF}} += +\sum_{i=1}^{n} +\left[ +C_n(\boldsymbol U_i) +- +C_{\widehat\theta}(\boldsymbol U_i) +\right]^2. +``` + +Because $\widehat\theta$ is estimated, the uncertainty introduced by fitting must also be reproduced in the bootstrap. Parametric-bootstrap validity for this type of semiparametric goodness-of-fit problem is studied in [genest2008bootstrap](@cite); practical copula GOF procedures and their finite sample behavior are discussed in [genest2009gof](@cite). + +In `Copulas.jl`, every composite bootstrap replicate performs the following steps: + +```math +\boldsymbol U_1^\star,\ldots,\boldsymbol U_n^\star +\sim +C_{\widehat\theta}, +``` + +then refits using the **same estimator specification**, + +```math +\widehat\theta^\star += +\operatorname{fit} +\left( +\boldsymbol U_1^\star,\ldots,\boldsymbol U_n^\star +\right), +``` + +and computes + +```math +S_n^\star += +\sum_{i=1}^{n} +\left[ +C_n^\star(\boldsymbol U_i^\star) +- +C_{\widehat\theta^\star}(\boldsymbol U_i^\star) +\right]^2. +``` + +Thus parameter estimation is repeated inside every bootstrap replicate rather than treating the fitted parameters as fixed. + +The fitting procedure itself is also reproduced. In particular, estimator-defining runtime information such as the fitting method, method-specific keywords, and copula structure is retained whenever the fitted model records a reproducible fitting specification. +This matters for models whose fitting procedure cannot be reconstructed from the fitted copula type alone. + +## Usage + +First fit a model: + +```@example hypothesis_testing +M = fit(CopulaModel, ClaytonCopula, Ugof; vcov=false,) + +nothing # hide +``` + +Then run the test directly from the fitted model: + +```@example hypothesis_testing +tcomposite = GOFCopulaTest(M; N=49, rng=Xoshiro(12),) + +(tcomposite.hypothesis.kind, pvalue(tcomposite)) +``` + +`GOFCopulaTest(M)` tests the data used to fit `M`. The stored fitting input is preprocessed consistently with the original fit before the observed statistic is computed. The fitted model `M` supplies both the estimated null model and the estimator specification that is replayed in every bootstrap replicate. + +A separate sample can be tested with + +```julia +GOFCopulaTest(M, U) +``` + +In this form, `M` is interpreted as an **estimator specification**, not as a fixed set of parameter estimates. The same fitting procedure is first reapplied to `U`, so the observed statistic uses parameters estimated from the sample being tested. Each parametric-bootstrap replicate then repeats that same fitting procedure. + +Consequently, the observed statistic and every bootstrap statistic are based on the same estimation rule. If the original fitting procedure cannot be reproduced safely, composite goodness-of-fit testing raises an `ArgumentError` rather than silently replacing it by a different estimator. + +The current defaults are + +```text +statistic = :Sn +calibration = :parametric_bootstrap +``` + +--- + +# Statistics and calibrations + +The public keywords + +```julia +statistic=:default +calibration=:default +``` + +are resolved through capability declarations. + +::: warning Internal dispatch hooks + +The underscore-prefixed functions shown below are implementation interfaces +for contributors to `Copulas.jl`. They are not public extension points and are +not covered by SemVer. + +::: + +For a hypothesis `h`, the available statistics are declared by + +```julia +Copulas._available_statistics(h) +``` + +and the available calibrations for a statistic `s` by + +```julia +Copulas._available_calibrations(h, Val(s)) +``` + +The **first element** of each returned tuple is the default. + +For example, the independence hypothesis declares conceptually + +```julia +_available_statistics(::IndependenceHypothesis) = (:cvm,) + +_available_calibrations(::IndependenceHypothesis, ::Val{:cvm},) = (:simulation,) +``` + +while the extreme-value hypothesis declares + +```julia +_available_statistics(::ExtremeValueHypothesis) = (:Sn,) + +_available_calibrations(::ExtremeValueHypothesis, ::Val{:Sn},) = (:multiplier,) +``` + +This convention deliberately mirrors the fitting interface: + +```text +_available_fitting_methods + ↓ + first = default + ↓ +_fit(..., Val(method)) +``` + +and, for hypothesis tests, + +```text +_available_statistics + ↓ + first = default + ↓ +_teststatistic(..., Val(statistic)) +``` + +followed by + +```text +_available_calibrations + ↓ + first = default + ↓ +_calibrate(..., Val(calibration), Val(statistic)) +``` + +The generic `CopulaTest` constructor therefore does not need to know which statistics are implemented by any particular hypothesis. + +::: info Why use `Val` internally? + +Users interact with ordinary symbols such as `:Sn`, `:cvm`, `:simulation`, +and `:multiplier`. Internally those symbols are converted to `Val` objects, +allowing Julia's multiple dispatch to select the appropriate mathematical +implementation without central `if`/`elseif` tables. + +::: + +--- + +# Calibration engines + +The framework currently provides four reusable calibration mechanisms. + +| Calibration | Principle | Typical use | +| ----------------------- | ------------------------------------------------- | ----------------------------------------- | +| `:simulation` | Generate directly under `H_0` | Independence | +| `:randomization` | Exploit invariance under `H_0` | Radial symmetry | +| `:multiplier` | Approximate an empirical-copula process | Exchangeability, extreme-value dependence | +| `:parametric_bootstrap` | Simulate from a parametric fitted/specified copula | Goodness of fit | + +The empirical-copula multiplier methodology is related to [remillard2009equality](@cite) and [bucher2010bootstrap](@cite), while the parametric-bootstrap framework for composite goodness-of-fit hypotheses is studied in [genest2008bootstrap](@cite). + +The concrete hypothesis only provides the mathematical ingredients required by the selected engine. The mechanics of repeated simulation, randomization, multiplier generation, or parametric bootstrap remain centralized. + +--- + +# Monte Carlo p-values + +Let $T_n$ be the observed statistic and let + +```math +T_n^{(1)},\ldots,T_n^{(N)} +``` + +denote resampled statistics. + +For calibrations using the finite-sample correction in the generic engine, `Copulas.jl` computes + +```math +\widehat p += +\frac{ +1/2+ +\sum_{b=1}^{N} +\mathbf 1 +\left\{ +T_n^{(b)}\ge T_n +\right\} +}{ +N+1 +}. +``` + +Specific calibration methods may override the comparison convention when their theoretical construction requires it. In particular, the exchangeability multiplier implementation uses strict exceedances and its corresponding uncorrected empirical proportion. + +Accordingly, $N$ controls Monte Carlo precision rather than the definition of the test statistic itself. + +--- + +# Internal contributor extension mechanism + +The internal hypothesis-testing machinery is designed so that contributors to +`Copulas.jl` can add new procedures while reusing the common constructor and +existing calibration engines. + +A new hypothesis starts with + +```julia +struct MyHypothesis <: CopulaHypothesis end +``` + +and then declares its name, null hypothesis, supported statistics, and calibrations: + +```julia +Copulas.testname(::MyHypothesis) = "My copula hypothesis test" + +Copulas.nullhypothesis(::MyHypothesis) = "The null hypothesis holds." + +Copulas._available_statistics(::MyHypothesis) = (:Sn, :ks) + +Copulas._available_calibrations(::MyHypothesis, ::Val{:Sn},) = (:simulation,) +``` + +The statistic is added through dispatch: + +```julia +function Copulas._teststatistic(::MyHypothesis, ::Val{:Sn}, U; kwargs...,) + # Compute and return the observed statistic. +end +``` + +If the generic simulation engine is appropriate, the hypothesis only needs to specify how to generate data under its null: + +```julia +function Copulas._simulation_sample(::MyHypothesis, U, rng,) + # Return a d × n sample generated under H₀. +end +``` + +The generic constructor then works automatically: + +```julia +test = CopulaTest(MyHypothesis(), U; N=999,) +``` + +No change to `CopulaTest`, the generic result type, or the display machinery is required. + +For a more complete contributor-facing description of these internal hooks, see the [Developer Guide](@ref developer_fitting). + +--- + +## References + +```@bibliography +Pages = [@__FILE__] +Canonical = false +``` diff --git a/src/CopulaTest.jl b/src/CopulaTest.jl new file mode 100644 index 000000000..a1c52a65c --- /dev/null +++ b/src/CopulaTest.jl @@ -0,0 +1,786 @@ +""" + CopulaHypothesis + +Abstract supertype for hypotheses about copulas. + +Subtypes describe what is being tested. The generic [`CopulaTest`](@ref) +constructor combines a hypothesis, a statistic, and a calibration method to +produce a standard `StatsAPI.HypothesisTest` result. +""" +abstract type CopulaHypothesis end + +""" + CopulaTest{H<:CopulaHypothesis} <: HypothesisTest + +Result of a copula hypothesis test. + +The hypothesis stores the mathematical null being tested; `CopulaTest` stores +the common result fields: sample size, dimension, observed statistic, p-value, +resampling method, and details useful for display or reproducibility. +""" +struct CopulaTest{H<:CopulaHypothesis,S<:Real,P<:Real,D<:NamedTuple} <: HypothesisTest + hypothesis::H + n::Int + dimension::Int + statistic_value::S + p::P + n_resamples::Int + statistic::Symbol + calibration::Symbol + details::D +end + +""" + teststatistic(test::CopulaTest) + +Return the observed value of the test statistic. +""" +teststatistic(test::CopulaTest) = test.statistic_value + +""" + pvalue(test::CopulaTest) + +Return the p-value of `test`. +""" +pvalue(test::CopulaTest) = test.p + +StatsBase.nobs(test::CopulaTest) = test.n + +""" + testname(x) + +Return the display name for a copula hypothesis or test. + +This is an extension hook for new copula tests. It is intentionally not exported to avoid clashes with `HypothesisTests.testname`. +""" +testname(test::CopulaTest) = testname(test.hypothesis) + +""" + nullhypothesis(x) + +Return the textual null hypothesis for a copula hypothesis or test. This is an extension hook used by the generic display machinery. +""" +nullhypothesis(test::CopulaTest) = nullhypothesis(test.hypothesis) + +""" + _available_statistics(h::CopulaHypothesis) + +Return the statistic symbols available for `h`. The first entry is the default, mirroring `_available_fitting_methods`. +""" +function _available_statistics(h::CopulaHypothesis) + throw(ArgumentError("No statistics are implemented for $(nameof(typeof(h))).")) +end + +function default_statistic(h::CopulaHypothesis) + return _find_statistic(h, :default) +end + +""" + _available_calibrations(h::CopulaHypothesis, ::Val{statistic}) + +Return the calibration symbols available for a hypothesis/statistic pair. The first entry is the default calibration. +""" +function _available_calibrations(h::CopulaHypothesis, ::Val{statistic}) where {statistic} + throw(ArgumentError("Statistic :$statistic is not implemented for $(nameof(typeof(h))).")) +end + +function default_calibration(h::CopulaHypothesis, stat::Val) + return _find_calibration(h, stat, :default) +end + +_symbol_list(symbols) = join((":" * String(symbol) for symbol in symbols), ", ") + +function _find_statistic(h::CopulaHypothesis, statistic::Symbol) + statistics = _available_statistics(h) + isempty(statistics) && throw(ArgumentError("No statistics are available for $(nameof(typeof(h))).")) + statistic === :default && return first(statistics) + statistic in statistics || throw(ArgumentError("Statistic :$statistic is not available for $(nameof(typeof(h))). Available statistics: $(_symbol_list(statistics)).")) + return statistic +end + +function _find_calibration(h::CopulaHypothesis, stat::Val{statistic}, calibration::Symbol) where {statistic} + calibrations = _available_calibrations(h, stat) + isempty(calibrations) && throw(ArgumentError("No calibrations are available for statistic :$statistic under $(nameof(typeof(h))).")) + calibration === :default && return first(calibrations) + calibration in calibrations || throw(ArgumentError("Calibration :$calibration is not available for statistic :$statistic under $(nameof(typeof(h))). Available calibrations: $(_symbol_list(calibrations)).")) + return calibration +end + +function _teststatistic(h::CopulaHypothesis, ::Val{statistic}, U::AbstractMatrix; kwargs...) where {statistic} + throw(ArgumentError("Statistic :$statistic is not implemented for $(nameof(typeof(h))).")) +end + +function _calibrate(h::CopulaHypothesis, ::Val{calibration}, ::Val{statistic}, U::AbstractMatrix, observed::Real; kwargs...) where {calibration,statistic} + throw(ArgumentError("Calibration :$calibration is not implemented for statistic :$statistic under $(nameof(typeof(h))).")) +end + +""" + CopulaTest(hypothesis, U; statistic, calibration, N, pseudo_values, rng) + +Run a copula hypothesis test. + +`statistic` and `calibration` are public symbols and are converted internally to `Val` dispatch. New hypotheses, statistics, and calibrations extend the framework by adding methods, not by modifying this constructor. +""" +function CopulaTest(h::CopulaHypothesis, U::AbstractMatrix{<:Real}; statistic::Symbol=:default, calibration::Symbol=:default, + N::Integer=1000, pseudo_values::Bool=false, rng::Distributions.AbstractRNG=Random.default_rng(), kwargs...) + V, d, n = _test_pseudos(U, pseudo_values) + statistic = _find_statistic(h, statistic) + stat = Val(statistic) + calibration = _find_calibration(h, stat, calibration) + observed = _teststatistic(h, stat, V; kwargs...) + p, n_resamples, details = _calibrate(h, Val(calibration), stat, V, observed; N=Int(N), rng=rng, kwargs...) + return CopulaTest(h, n, d, observed, p, n_resamples, statistic, calibration, details) +end + +function _test_pseudos(U::AbstractMatrix{<:Real}, pseudo_values::Bool) + all(isfinite, U) || throw(ArgumentError("input data must be finite")) + d, n = size(U) + d >= 2 || throw(ArgumentError("at least two components are required")) + n >= 2 || throw(ArgumentError("at least two observations are required")) + + for j in 1:d + allunique(@view U[j, :]) || throw(ArgumentError( + "copula hypothesis tests currently require continuous, tie-free margins; " * + "ties were detected in margin $j. Tie-aware procedures are not yet implemented.")) + end + + V = pseudo_values ? Matrix{Float64}(U) : pseudos(U) + all(x -> 0 <= x <= 1, V) || throw(ArgumentError("pseudo-observations must lie in [0, 1]")) + return V, d, n +end + +function _empirical_copula_partial(Cn::EmpiricalCopula, u::AbstractVector, l::Integer, h::Real) + lo = Vector{Float64}(u) + hi = Vector{Float64}(u) + lo[l] = max(lo[l] - h, 0.0) + hi[l] = min(hi[l] + h, 1.0) + width = hi[l] - lo[l] + return width > 0 ? (Distributions.cdf(Cn, hi) - Distributions.cdf(Cn, lo)) / width : 0.0 +end + +function _exceedance_pvalue(exceedances::Integer, N::Integer; correction=0.5) + correction === nothing && return exceedances / N + return (correction + exceedances) / (N + 1) +end + +function _check_resamples(N::Integer) + N >= 1 || throw(ArgumentError("`N` must be positive.")) + return Int(N) +end + +function _simulation_sample(h::CopulaHypothesis, U::AbstractMatrix, rng::Distributions.AbstractRNG) + throw(ArgumentError("Simulation under the null is not implemented for $(nameof(typeof(h))).")) +end + +function _calibrate(h::CopulaHypothesis, ::Val{:simulation}, stat::Val, U::AbstractMatrix, observed::Real; N::Integer, rng::Distributions.AbstractRNG, kwargs...) + N = _check_resamples(N) + exceedances = 0 + for _ in 1:N + sample = pseudos(_simulation_sample(h, U, rng)) + exceedances += _teststatistic(h, stat, sample; kwargs...) >= observed + end + return _exceedance_pvalue(exceedances, N), N, (;) +end + +function _randomization_sample(h::CopulaHypothesis, U::AbstractMatrix, rng::Distributions.AbstractRNG) + throw(ArgumentError("Randomization under the null is not implemented for $(nameof(typeof(h))).")) +end + +_randomization_details(::CopulaHypothesis) = (;) +_randomization_pseudos(::CopulaHypothesis, sample::AbstractMatrix) = pseudos(sample) + +function _calibrate(h::CopulaHypothesis, ::Val{:randomization}, stat::Val, U::AbstractMatrix, observed::Real; N::Integer, rng::Distributions.AbstractRNG, kwargs...) + N = _check_resamples(N) + exceedances = 0 + for _ in 1:N + sample = _randomization_pseudos(h, _randomization_sample(h, U, rng),) + exceedances += _teststatistic(h, stat, sample; kwargs...) >= observed + end + return _exceedance_pvalue(exceedances, N), N, _randomization_details(h) +end + +function _multiplier_representation(h::CopulaHypothesis, ::Val{statistic}, U::AbstractMatrix) where {statistic} + throw(ArgumentError("Calibration :multiplier is not implemented for statistic :$statistic under $(nameof(typeof(h))).")) +end + +function _calibrate(h::CopulaHypothesis, ::Val{:multiplier}, stat::Val, U::AbstractMatrix, observed::Real; N::Integer, rng::Distributions.AbstractRNG, kwargs...) + N = _check_resamples(N) + rep = _multiplier_representation(h, stat, U) + p = _multiplier_pvalue(rep.matrices, observed, N, rng; + weights=get(rep, :weights, nothing), + scale=rep.scale, + strict=get(rep, :strict, false), + correction=get(rep, :correction, 0.5)) + return p, N, get(rep, :details, (;)) +end + +function _multiplier_pvalue(matrices, observed::Real, N::Integer, rng::Distributions.AbstractRNG; weights=nothing, scale::Real, strict::Bool=false, correction=0.5) + n = size(first(matrices), 2) + xi = Vector{Float64}(undef, n) + work = Vector{Float64}(undef, n) + inv_sqrt_n = inv(sqrt(n)) + exceedances = 0 + + for _ in 1:N + Random.randexp!(rng, xi) + xi .-= Statistics.mean(xi) + bootstrap_stat = 0.0 + + if weights === nothing + for Q in matrices + LinearAlgebra.mul!(work, Q, xi) + @inbounds for i in 1:n + bootstrap_stat += abs2(inv_sqrt_n * work[i]) + end + end + else + for (Q, w) in zip(matrices, weights) + LinearAlgebra.mul!(work, Q, xi) + @inbounds for i in 1:n + bootstrap_stat += abs2(inv_sqrt_n * work[i]) * w[i] + end + end + end + + value = scale * bootstrap_stat + exceedances += strict ? value > observed : value >= observed + end + + return _exceedance_pvalue(exceedances, N; correction) +end + +function _bootstrap_copula(h::CopulaHypothesis) + throw(ArgumentError("Parametric bootstrap is not implemented for $(nameof(typeof(h))).")) +end + +_bootstrap_hypothesis(h::CopulaHypothesis, ::AbstractMatrix) = h + +function _calibrate(h::CopulaHypothesis, ::Val{:parametric_bootstrap}, stat::Val, U::AbstractMatrix, observed::Real; N::Integer, rng::Distributions.AbstractRNG, kwargs...) + N = _check_resamples(N) + _, n = size(U) + exceedances = 0 + for _ in 1:N + sample = pseudos(rand(rng, _bootstrap_copula(h), n)) + bootstrap_hypothesis = _bootstrap_hypothesis(h, sample) + exceedances += _teststatistic(bootstrap_hypothesis, stat, sample; kwargs...) >= observed + end + return _exceedance_pvalue(exceedances, N), N, (;) +end + +################################################################################ +##### Independence +################################################################################ + +""" + IndependenceHypothesis() + +Hypothesis that the components of a copula are mutually independent. +""" +struct IndependenceHypothesis <: CopulaHypothesis end + +""" + IndependenceCopulaTest(U; statistic=:cvm, N=1000, calibration=:simulation, pseudo_values=false, rng=Random.default_rng()) + +Test mutual independence between the components of a random vector. +""" +const IndependenceCopulaTest = CopulaTest{IndependenceHypothesis} + +(::Type{CopulaTest{IndependenceHypothesis}})(U::AbstractMatrix{<:Real}; kwargs...) = CopulaTest(IndependenceHypothesis(), U; kwargs...) + +testname(::IndependenceHypothesis) = "Copula independence test" +nullhypothesis(::IndependenceHypothesis) = "The components are mutually independent." +_available_statistics(::IndependenceHypothesis) = (:cvm,) +_available_calibrations(::IndependenceHypothesis, ::Val{:cvm}) = (:simulation,) + +function _teststatistic(::IndependenceHypothesis, ::Val{:cvm}, U::AbstractMatrix; kwargs...) + Cn = EmpiricalCopula(U; pseudo_values=true) + s = 0.0 + @inbounds for u in eachcol(U) + s += abs2(Distributions.cdf(Cn, u) - prod(u)) + end + return s +end + +function _simulation_sample(::IndependenceHypothesis, U::AbstractMatrix, rng::Distributions.AbstractRNG) + sample = similar(U) + Random.rand!(rng, sample) + return sample +end + +################################################################################ +##### Exchangeability +################################################################################ + +""" + ExchangeabilityHypothesis(; permutations=:G2, weight=:wm2) + +Hypothesis that a copula is invariant under coordinate permutations. +""" +struct ExchangeabilityHypothesis{P} <: CopulaHypothesis + permutations::P + weight::Symbol +end + +ExchangeabilityHypothesis(; permutations=:G2, weight::Symbol=:wm2) = ExchangeabilityHypothesis(permutations, weight) + +""" + ExchangeabilityCopulaTest(U; statistic=:Sn, permutations=:G2, weight=:wm2, N=1000, calibration=:multiplier, pseudo_values=false, rng=Random.default_rng()) + +Test exchangeability of a copula in arbitrary dimension. +""" +const ExchangeabilityCopulaTest = CopulaTest{<:ExchangeabilityHypothesis} + +function (::Type{<:CopulaTest{<:ExchangeabilityHypothesis}})(U::AbstractMatrix{<:Real}; permutations=:G2, weight::Symbol=:wm2, kwargs...) + return CopulaTest(ExchangeabilityHypothesis(; permutations, weight), U; kwargs...) +end + +testname(::ExchangeabilityHypothesis) = "Copula exchangeability test" +nullhypothesis(::ExchangeabilityHypothesis) = "The copula is exchangeable." +_available_statistics(::ExchangeabilityHypothesis) = (:Sn,) +_available_calibrations(::ExchangeabilityHypothesis, ::Val{:Sn}) = (:multiplier,) + +function _teststatistic(h::ExchangeabilityHypothesis, ::Val{:Sn}, U::AbstractMatrix; kwargs...) + d, n = size(U) + + # `:all` must be checked before resolution because materializing d! + # permutations may itself be prohibitive. + _check_exchangeability_all_cost(h.permutations, d, n) + + permutations = _exchangeability_permutations(h.permutations, d) + + # Once resolved, guard the actual number of selected permutations. This also + # protects explicit custom collections. + _check_exchangeability_matrix_cost(length(permutations), n) + + return _exchangeability_sn_statistic(U, permutations, h.weight,) +end + +const _MAX_EXCHANGEABILITY_MATRIX_BYTES = 512 * 1024^2 + +function _check_exchangeability_matrix_cost(nperms::Integer, n::Integer; + label::AbstractString="the requested permutation collection") + matrix_bytes = big(nperms) * big(n)^2 * sizeof(Float64) + + matrix_bytes <= _MAX_EXCHANGEABILITY_MATRIX_BYTES && return nothing + + estimated_mib = Float64(matrix_bytes) / 1024^2 + limit_mib = _MAX_EXCHANGEABILITY_MATRIX_BYTES / 1024^2 + + throw(ArgumentError( + "$(label) would materialize $(nperms) dense $(n)×$(n) " * + "multiplier matrices (approximately $(round(estimated_mib; digits=1)) MiB), " * + "exceeding the current $(round(limit_mib; digits=0)) MiB safety limit. " * + "Use `permutations=:G1`, `:G2`, or provide a smaller custom collection." + )) +end + +function _check_exchangeability_all_cost(permutations, d::Integer, n::Integer) + permutations === :all || return nothing + + # Check the factorial-sized collection before materializing it. + nperms = factorial(big(d)) - 1 + return _check_exchangeability_matrix_cost( + nperms, + n; + label="`permutations=:all`", + ) +end + +function _exchangeability_permutations(permutations, d::Integer) + identity_perm = ntuple(i -> i, d) + raw = if permutations === :G2 + d == 2 ? ((2, 1),) : + ((2, 1, ntuple(i -> i + 2, d - 2)...), ntuple(i -> i == d ? 1 : i + 1, d)) + elseif permutations === :G1 + ntuple(i -> Tuple(j == 1 ? i + 1 : j == i + 1 ? 1 : j for j in 1:d), d - 1) + elseif permutations === :all + Combinatorics.permutations(1:d) + else + is_single = (permutations isa Tuple || permutations isa AbstractVector) && length(permutations) == d && all(x -> x isa Integer, permutations) + is_single ? (permutations,) : permutations + end + + result = NTuple{d,Int}[] + for perm in raw + p = Tuple(Int.(perm)) + length(p) == d || throw(ArgumentError("permutations must have length $d")) + sort(collect(p)) == collect(1:d) || throw(ArgumentError("invalid permutation `$perm`")) + p == identity_perm || push!(result, p) + end + isempty(result) && throw(ArgumentError("at least one non-identity permutation is required")) + return Tuple(result) +end + +function _exchangeability_weight(u::AbstractVector, perm::Tuple, weight::Symbol) + weight === :none && return 1.0 + weight === :wm2 || throw(ArgumentError("Only `weight=:wm2` and `weight=:none` are implemented.")) + + m = minimum(u) + omega = if count(i -> perm[i] != i, eachindex(perm)) == 2 && all(perm[perm[i]] == i for i in eachindex(perm)) + i = findfirst(k -> perm[k] != k, eachindex(perm)) + j = perm[i] + abs(u[i] - u[j]) + else + v = sort(collect(u)) + sum(v[i] - m for i in cld(length(v), 2) + 1:length(v)) + end + wm = min(m, omega, length(u) - 1 + m - sum(u)) + return abs2(max(wm, 0.0)) +end + +function _exchangeability_sn_statistic(U::AbstractMatrix, permutations, weight::Symbol) + d, n = size(U) + Cn = EmpiricalCopula(U; pseudo_values=true) + s = 0.0 + uperm = Vector{Float64}(undef, d) + + @inbounds for perm in permutations + for i in 1:n + u = @view U[:, i] + for k in 1:d + uperm[k] = u[perm[k]] + end + diff = Distributions.cdf(Cn, u) - Distributions.cdf(Cn, uperm) + s += abs2(diff) * _exchangeability_weight(u, perm, weight) + end + end + return s +end + +function _multiplier_representation(h::ExchangeabilityHypothesis, ::Val{:Sn}, U::AbstractMatrix) + d, n = size(U) + _check_exchangeability_all_cost(h.permutations, d, n) + permutations = _exchangeability_permutations(h.permutations, d) + _check_exchangeability_matrix_cost(length(permutations), n) + matrices, weights, bandwidth = _exchangeability_multiplier_matrices(U, permutations, h.weight) + return (;matrices, weights, scale=inv(n), strict=true, correction=nothing, + details=(; permutations=h.permutations, generator=permutations, weight=h.weight, multiplier=:exponential, derivative_bandwidth=bandwidth),) +end + +function _exchangeability_multiplier_matrices(U::AbstractMatrix, permutations, weight::Symbol) + d, n = size(U) + Cn = EmpiricalCopula(U; pseudo_values=true) + h = inv(sqrt(n)) + partials = Matrix{Float64}(undef, d, n) + q_matrices = Matrix{Float64}[] + weights = Vector{Float64}[] + + @inbounds for i in 1:n + u = @view U[:, i] + for l in 1:d + partials[l, i] = _empirical_copula_partial(Cn, u, l, h) + end + end + + @inbounds for perm in permutations + invperm = Vector{Int}(undef, d) + for k in 1:d + invperm[perm[k]] = k + end + + Q = Matrix{Float64}(undef, n, n) + w = Vector{Float64}(undef, n) + for i in 1:n + u = @view U[:, i] + w[i] = _exchangeability_weight(u, perm, weight) + for j in 1:n + le_u = true + le_up = true + for k in 1:d + U[k, j] <= u[k] || (le_u = false) + U[k, j] <= u[perm[k]] || (le_up = false) + end + + q = (le_u ? 1.0 : 0.0) - (le_up ? 1.0 : 0.0) + for l in 1:d + le_margin = U[l, j] <= u[l] + le_permuted_margin = U[invperm[l], j] <= u[l] + q -= partials[l, i] * + ((le_margin ? 1.0 : 0.0) - (le_permuted_margin ? 1.0 : 0.0)) + end + Q[i, j] = q + end + end + push!(q_matrices, Q) + push!(weights, w) + end + return q_matrices, weights, h +end + +################################################################################ +##### Radial Symmetry +################################################################################ + +""" + RadialSymmetryHypothesis() + +Hypothesis that a copula is radially symmetric. +""" +struct RadialSymmetryHypothesis <: CopulaHypothesis end + +""" + RadialSymmetryCopulaTest(U; statistic=:Sn, N=1000, calibration=:randomization, pseudo_values=false, rng=Random.default_rng()) + +Test radial symmetry of a copula. +""" +const RadialSymmetryCopulaTest = CopulaTest{RadialSymmetryHypothesis} + +(::Type{CopulaTest{RadialSymmetryHypothesis}})(U::AbstractMatrix{<:Real}; kwargs...) = CopulaTest(RadialSymmetryHypothesis(), U; kwargs...) + +testname(::RadialSymmetryHypothesis) = "Copula radial symmetry test" +nullhypothesis(::RadialSymmetryHypothesis) = "The copula is radially symmetric." +_available_statistics(::RadialSymmetryHypothesis) = (:Sn,) +_available_calibrations(::RadialSymmetryHypothesis, ::Val{:Sn}) = (:randomization,) + +function _teststatistic(::RadialSymmetryHypothesis, ::Val{:Sn}, U::AbstractMatrix; kwargs...) + Cn = EmpiricalCopula(U; pseudo_values=true) + Cbar = EmpiricalCopula(1 .- U; pseudo_values=true) + s = 0.0 + @inbounds for u in eachcol(U) + s += abs2(Distributions.cdf(Cn, u) - Distributions.cdf(Cbar, u)) + end + return s +end + +function _randomization_sample(::RadialSymmetryHypothesis, U::AbstractMatrix, rng::Distributions.AbstractRNG) + d, n = size(U) + sample = similar(U) + @inbounds for i in 1:n + reflected = rand(rng) < 0.5 + for j in 1:d + sample[j, i] = reflected ? 1 - U[j, i] : U[j, i] + end + end + return sample +end + +function _average_pseudos(sample::AbstractMatrix) + d, n = size(sample) + U = Matrix{Float64}(undef, d, n) + denom = n + 1 + + @inbounds for j in 1:d + U[j, :] .= StatsBase.tiedrank(@view sample[j, :]) ./ denom + end + + return U +end + +_randomization_pseudos(::RadialSymmetryHypothesis, sample::AbstractMatrix,) = _average_pseudos(sample) + +_randomization_details(::RadialSymmetryHypothesis) = (; reflection_probability=0.5,) + +################################################################################ +##### Extreme Value +################################################################################ + +""" + ExtremeValueHypothesis(; powers=3:5) + +Hypothesis that a copula belongs to the extreme-value class. +""" +struct ExtremeValueHypothesis{P} <: CopulaHypothesis + powers::P +end + +ExtremeValueHypothesis(; powers=3:5) = ExtremeValueHypothesis(powers) + +""" + ExtremeValueCopulaTest(U; statistic=:Sn, powers=3:5, N=1000, calibration=:multiplier, pseudo_values=false, rng=Random.default_rng()) + +Test whether a copula belongs to the extreme-value class. +""" +const ExtremeValueCopulaTest = CopulaTest{<:ExtremeValueHypothesis} + +function (::Type{<:CopulaTest{<:ExtremeValueHypothesis}})(U::AbstractMatrix{<:Real}; powers=3:5, kwargs...) + return CopulaTest(ExtremeValueHypothesis(; powers), U; kwargs...) +end + +testname(::ExtremeValueHypothesis) = "Extreme-value copula test" +nullhypothesis(::ExtremeValueHypothesis) = "The copula belongs to the extreme-value class." +_available_statistics(::ExtremeValueHypothesis) = (:Sn,) +_available_calibrations(::ExtremeValueHypothesis, ::Val{:Sn}) = (:multiplier,) + +function _teststatistic(h::ExtremeValueHypothesis, ::Val{:Sn}, U::AbstractMatrix; kwargs...) + return _extreme_value_sn_statistic(U, _max_stability_powers(h.powers)) +end + +function _max_stability_powers(powers) + raw = powers isa Real ? (powers,) : Tuple(powers) + result = Float64[] + for r in raw + isfinite(r) && r > 1 || + throw(ArgumentError("max-stability powers must be finite and greater than one")) + push!(result, Float64(r)) + end + isempty(result) && throw(ArgumentError("at least one max-stability power is required")) + return Tuple(result) +end + +function _extreme_value_sn_statistic(U::AbstractMatrix, powers) + d, n = size(U) + Cn = EmpiricalCopula(U; pseudo_values=true) + uroot = Vector{Float64}(undef, d) + s = 0.0 + + @inbounds for r in powers + invr = inv(r) + for u in eachcol(U) + for k in 1:d + uroot[k] = u[k]^invr + end + diff = Distributions.cdf(Cn, uroot)^r - Distributions.cdf(Cn, u) + s += abs2(diff) + end + end + return s +end + +function _multiplier_representation(h::ExtremeValueHypothesis, ::Val{:Sn}, U::AbstractMatrix) + powers = _max_stability_powers(h.powers) + matrices, bandwidth = _extreme_value_multiplier_matrices(U, powers) + _, n = size(U) + return (;matrices, scale=inv(n), strict=false, correction=0.5, + details=(; powers, multiplier=:exponential, derivative_bandwidth=bandwidth),) +end + +function _extreme_value_multiplier_matrices(U::AbstractMatrix, powers) + d, n = size(U) + Cn = EmpiricalCopula(U; pseudo_values=true) + h = inv(sqrt(n)) + uroot = Vector{Float64}(undef, d) + partials_u = Vector{Float64}(undef, d) + partials_root = Vector{Float64}(undef, d) + matrices = Matrix{Float64}[] + + @inbounds for r in powers + Q = Matrix{Float64}(undef, n, n) + invr = inv(r) + for i in 1:n + u = @view U[:, i] + for k in 1:d + uroot[k] = u[k]^invr + partials_u[k] = _empirical_copula_partial(Cn, u, k, h) + end + croot = Distributions.cdf(Cn, uroot) + factor = r * croot^(r - 1) + for k in 1:d + partials_root[k] = _empirical_copula_partial(Cn, uroot, k, h) + end + + for j in 1:n + le_u = true + le_root = true + for k in 1:d + U[k, j] <= u[k] || (le_u = false) + U[k, j] <= uroot[k] || (le_root = false) + end + + q_u = le_u ? 1.0 : 0.0 + q_root = le_root ? 1.0 : 0.0 + for k in 1:d + q_u -= partials_u[k] * (U[k, j] <= u[k] ? 1.0 : 0.0) + q_root -= partials_root[k] * (U[k, j] <= uroot[k] ? 1.0 : 0.0) + end + Q[i, j] = factor * q_root - q_u + end + end + push!(matrices, Q) + end + return matrices, h +end + +################################################################################ +##### Goodness of Fit +################################################################################ + +""" + GoodnessOfFitHypothesis(model) + +Hypothesis that data follow a specified copula or fitted copula model. The field `kind` distinguishes `:simple` and `:composite`. +""" +struct GoodnessOfFitHypothesis{M} <: CopulaHypothesis + model::M + kind::Symbol +end + +GoodnessOfFitHypothesis(C::Copula) = GoodnessOfFitHypothesis(C, :simple) +GoodnessOfFitHypothesis(M::CopulaModel) = GoodnessOfFitHypothesis(M, :composite) + +""" + GOFCopulaTest(C, U; statistic=:Sn, N=1000, calibration=:parametric_bootstrap, pseudo_values=false, rng=Random.default_rng()) + GOFCopulaTest(model, U; statistic=:Sn, N=1000, calibration=:parametric_bootstrap, pseudo_values=false, rng=Random.default_rng()) + GOFCopulaTest(model; kwargs...) + +Test goodness of fit for a copula or fitted copula model. + +`GOFCopulaTest(C, U)` treats `C` as a fixed specified copula and tests a simple +null hypothesis. + +`GOFCopulaTest(M)` tests the fitting sample stored by `M` under a composite null +hypothesis. The estimator specification that produced `M` is replayed in every +parametric-bootstrap replicate. + +`GOFCopulaTest(M, U)` first refits that estimator specification on `U`; the +resulting fitted model is used for the observed statistic, and the same fitting +procedure is repeated in every bootstrap replicate. If the fitting procedure is +not reproducibly specified, composite GOF throws an `ArgumentError`. +""" +const GOFCopulaTest = CopulaTest{<:GoodnessOfFitHypothesis} + +function (::Type{<:CopulaTest{<:GoodnessOfFitHypothesis}})(C::Copula, U::AbstractMatrix{<:Real}; kwargs...) + return CopulaTest(GoodnessOfFitHypothesis(C), U; kwargs...) +end + +function (::Type{<:CopulaTest{<:GoodnessOfFitHypothesis}})(M::CopulaModel, U::AbstractMatrix{<:Real}; pseudo_values::Bool=false, kwargs...) + # For a composite null, M describes the estimator specification. + # The observed statistic must use parameters estimated from the sample being + # tested, just as every bootstrap replicate is refitted. + V, _, _ = _test_pseudos(U, pseudo_values) + Mrefit = _refit(M, V) + return CopulaTest(GoodnessOfFitHypothesis(Mrefit), V; pseudo_values=true, kwargs...,) +end + +function (::Type{<:CopulaTest{<:GoodnessOfFitHypothesis}})(M::CopulaModel; kwargs...) + haskey(M.method_details, :U) || throw(ArgumentError("the fitted model does not store its fitting sample")) + + # Most copula fits receive pseudo-observations directly. Empirical fitting + # routines may instead have received raw data with `pseudo_values=false`; + # respect that metadata rather than silently treating the stored input as + # already ranked. + stored_pseudo_values = get(M.method_details, :pseudo_values, true) + + return CopulaTest(GoodnessOfFitHypothesis(M), M.method_details.U; pseudo_values=stored_pseudo_values, kwargs...,) +end + +testname(::GoodnessOfFitHypothesis) = "Copula goodness-of-fit test" +function nullhypothesis(h::GoodnessOfFitHypothesis) + h.kind === :simple && return "The data follow the specified copula." + return "The data belong to the specified copula family." +end + +_available_statistics(::GoodnessOfFitHypothesis) = (:Sn,) +_available_calibrations(::GoodnessOfFitHypothesis, ::Val{:Sn}) = (:parametric_bootstrap,) + +function _teststatistic(h::GoodnessOfFitHypothesis, ::Val{:Sn}, U::AbstractMatrix; kwargs...) + C = _gof_copula(h) + length(C) == size(U, 1) || throw(DimensionMismatch("model dimension does not match input data")) + return _gof_sn_statistic(U, C) +end + +_gof_copula(h::GoodnessOfFitHypothesis) = h.model isa CopulaModel ? _copula_of(h.model) : h.model + +function _gof_sn_statistic(U::AbstractMatrix, C::Copula) + Cn = EmpiricalCopula(U; pseudo_values=true) + s = 0.0 + @inbounds for u in eachcol(U) + s += abs2(Distributions.cdf(Cn, u) - Distributions.cdf(C, u)) + end + return s +end + +_bootstrap_copula(h::GoodnessOfFitHypothesis) = _gof_copula(h) + +function _bootstrap_hypothesis(h::GoodnessOfFitHypothesis{<:CopulaModel}, U::AbstractMatrix) + return GoodnessOfFitHypothesis(_refit(h.model, U)) +end diff --git a/src/Copulas.jl b/src/Copulas.jl index e60791ad6..7670997b7 100644 --- a/src/Copulas.jl +++ b/src/Copulas.jl @@ -21,6 +21,7 @@ module Copulas import Printf import TaylorSeries import ADTypes + import StatsAPI: HypothesisTest, pvalue # Main code include("utils.jl") @@ -136,10 +137,16 @@ module Copulas include("ArchimaxCopula.jl") + include("CopulaTest.jl") + include("show.jl") export pseudos, condition, subsetdims, rosenblatt, inverse_rosenblatt, Nataf export SklarDist, CopulaModel + export CopulaHypothesis, CopulaTest + export IndependenceCopulaTest, ExchangeabilityCopulaTest + export RadialSymmetryCopulaTest, ExtremeValueCopulaTest, GOFCopulaTest + export pvalue, teststatistic export WilliamsonGenerator, 𝒲, EmpiricalGenerator, DiscreteSpectralTail export ArchimedeanCopula, ExtremeValueCopula, LiouvilleCopula diff --git a/src/Fitting.jl b/src/Fitting.jl index 15c5d43eb..ac24682a3 100644 --- a/src/Fitting.jl +++ b/src/Fitting.jl @@ -62,6 +62,47 @@ struct CopulaModel{CT, TM<:Union{Nothing,AbstractMatrix}, TD<:NamedTuple} <: Sta end end +# Internal description of the estimator that produced a CopulaModel. +# +# `target` is the fitting target accepted by `fit(CopulaModel, target, U; ...)` +# (normally a copula type, but it can also be a runtime template such as a +# NestedArchimedeanCopula instance). `kwargs` contains estimator-defining +# keywords only: generic inference controls such as `vcov` and +# `derived_measures` are handled separately by `fit`. +struct _CopulaFitSpec{T,K<:NamedTuple} + target::T + method::Symbol + kwargs::K +end + +# Bootstrap/refit samples are already on the copula scale. If the original +# estimator accepted raw data through `pseudo_values=false`, replay the same +# estimator on the supplied pseudo-observations without ranking them again. +function _refit_kwargs(kwargs::NamedTuple) + haskey(kwargs, :pseudo_values) || return kwargs + return merge(kwargs, (; pseudo_values=true)) +end + +""" + _refit(M::CopulaModel, U) + +Refit the same estimator specification that produced `M` to pseudo-observations +`U`. + +This is an internal inference hook. A model is refittable only when its fitting +entry point recorded a reproducible `_CopulaFitSpec`. +""" +function _refit(M::CopulaModel, U::AbstractMatrix) + spec = get(M.method_details, :_fit_spec, nothing) + spec isa _CopulaFitSpec || throw(ArgumentError( + "this fitted model does not store a reproducible fitting specification; " * + "composite goodness-of-fit refitting is unavailable for this model")) + + kwargs = _refit_kwargs(spec.kwargs) + + return Distributions.fit(CopulaModel, spec.target, U; method=spec.method, derived_measures=false, vcov=false, kwargs...,) +end + # Fallbacks that throw if the interface is not implemented correctly. """ Distributions.params(C::Copula) @@ -207,6 +248,7 @@ function Distributions.fit(::Type{CopulaModel}, CT::Type{<:Copula}, U; throw(ArgumentError("unknown vcov method `$vcov_method`; expected one of $allowed_vcov")) d, n = size(U) method = _find_method(CT, d, method) + fit_spec = _CopulaFitSpec(CT, method, (; kwargs...)) t = @elapsed (rez = _fit(CT, U, Val{method}(); kwargs...)) C, meta = rez quick_fit && return (result=C,) # as soon as possible. @@ -230,8 +272,7 @@ function Distributions.fit(::Type{CopulaModel}, CT::Type{<:Copula}, U; meta = (; meta..., vcov, vmeta...) end - md = (; d, n, method, meta..., null_ll=0.0, - elapsed_sec=t, derived_measures, U=U) + md = (; d, n, method, meta..., null_ll=0.0, elapsed_sec=t, derived_measures, U=U, _fit_spec=fit_spec) return CopulaModel(C, n, ll, method; vcov = get(md, :vcov, nothing), diff --git a/src/NestedArchimedeanCopula.jl b/src/NestedArchimedeanCopula.jl index 41381db5e..3409da787 100644 --- a/src/NestedArchimedeanCopula.jl +++ b/src/NestedArchimedeanCopula.jl @@ -1088,7 +1088,8 @@ one of two ways: use `fit(CopulaModel, reparam, init, U).result`. """ # Shared optimiser + model assembly for a parametrisation `recon: α -> copula`. -function _fit_nested(recon, α₀::AbstractVector, U, d::Int, n::Int; quick_fit, derived_measures) +function _fit_nested(recon, α₀::AbstractVector, U, d::Int, n::Int; + quick_fit, derived_measures, fit_spec=nothing) loss(α) = -Distributions.loglikelihood(recon(α), U) t = @elapsed res = try Optim.optimize(loss, α₀, Optim.LBFGS(); autodiff = ADTypes.AutoForwardDiff()) @@ -1102,7 +1103,8 @@ function _fit_nested(recon, α₀::AbstractVector, U, d::Int, n::Int; quick_fit, # (type-positional reconstruction we do not have) is never reached. md = (; d, n, method = :mle, nparams = length(α₀), optimizer = Optim.summary(res), converged = Optim.converged(res), - iterations = Optim.iterations(res), elapsed_sec = t, derived_measures, U = U) + iterations = Optim.iterations(res), elapsed_sec = t, derived_measures, + U = U, _fit_spec = fit_spec) return CopulaModel(Chat, n, ll, :mle; vcov = nothing, converged = Optim.converged(res), iterations = Optim.iterations(res), elapsed_sec = t, method_details = md) @@ -1123,8 +1125,9 @@ function Distributions.fit(::Type{CopulaModel}, C0::NestedArchimedeanCopula{d}, method=:mle, quick_fit=false, vcov=false, derived_measures=true, kwargs...) where {d} method === :mle || throw(ArgumentError("NestedArchimedeanCopula supports only method=:mle (got $method).")) _validate_nested_fit_data(U, d) + fit_spec = _CopulaFitSpec(C0, :mle, (; kwargs...)) return _fit_nested(Base.Fix1(_nested_rebound, C0), _nested_unbound(C0), U, d, size(U, 2); - quick_fit, derived_measures) + quick_fit, derived_measures, fit_spec) end # Custom parametrisation: a map `reparam : α -> NestedArchimedeanCopula` and its diff --git a/src/show.jl b/src/show.jl index 1418f7b3d..f4c06b03d 100644 --- a/src/show.jl +++ b/src/show.jl @@ -281,3 +281,70 @@ function _print_marginals_section(io, S::SklarDist, Vm) end end end + +############################################################################### +##### Copula hypothesis tests +############################################################################### + +_show_test_model(::IO, ::CopulaHypothesis, ::Val, ::Val, ::NamedTuple) = nothing +function _show_test_model(io::IO, h::GoodnessOfFitHypothesis, ::Val, ::Val, ::NamedTuple) + label = + h.kind === :simple ? "Specified copula" : + "Fitted model" + + model = h.model isa CopulaModel ? _copula_of(h.model) : h.model + model_label = replace(string(typeof(model)), "Copulas." => "") + println(io, "Hypothesis: ", h.kind) + println(io, label, ": ", model_label) +end + +_show_test_details(::IO, ::CopulaHypothesis, ::Val, ::Val, ::NamedTuple) = nothing +function _show_test_details(io::IO, ::ExchangeabilityHypothesis, ::Val{:Sn}, + ::Val{:multiplier}, details::NamedTuple) + hasproperty(details, :permutations) || return nothing + println(io, "Permutations: ", details.permutations) + println(io, "Weight: ", details.weight) + if hasproperty(details, :multiplier) + println(io, "Multiplier: ", details.multiplier) + end + if hasproperty(details, :derivative_bandwidth) + println(io, "Derivative bandwidth: ", details.derivative_bandwidth) + end +end + +function _show_test_details(io::IO, ::RadialSymmetryHypothesis, ::Val{:Sn}, + ::Val{:randomization}, details::NamedTuple) + hasproperty(details, :reflection_probability) || return nothing + println(io, "Reflection probability: ", details.reflection_probability) +end + +function _show_test_details(io::IO, ::ExtremeValueHypothesis, ::Val{:Sn}, + ::Val{:multiplier}, details::NamedTuple) + hasproperty(details, :powers) || return nothing + println(io, "Powers: ", details.powers) + println(io, "Multiplier: ", details.multiplier) + println(io, "Derivative bandwidth: ", details.derivative_bandwidth) +end + +function Base.show(io::IO, ::MIME"text/plain", test::CopulaTest) + name = testname(test) + println(io, name) + println(io, repeat('-', length(name))) + h = test.hypothesis + statistic = Val(test.statistic) + calibration = Val(test.calibration) + _show_test_model(io, h, statistic, calibration, test.details) + println(io, "Number of observations: ", StatsBase.nobs(test)) + println(io, "Dimension: ", test.dimension) + println(io, "Statistic: ", replace(string(test.statistic), '_' => ' ')) + println(io, "Observed value: ", teststatistic(test)) + _show_test_details(io, h, statistic, calibration, test.details) + if test.n_resamples > 0 + println(io, "Number of resamples: ", test.n_resamples) + println(io, "Calibration: ", replace(string(test.calibration), '_' => ' ')) + end + println(io, "p-value: ", pvalue(test)) + println(io) + println(io, "Null hypothesis:") + print(io, nullhypothesis(test)) +end diff --git a/test/operations/hypothesis_testing.jl b/test/operations/hypothesis_testing.jl new file mode 100644 index 000000000..8b976ba05 --- /dev/null +++ b/test/operations/hypothesis_testing.jl @@ -0,0 +1,574 @@ +const COPULA_TEST_RESAMPLES = parse(Int, get(ENV, "COPULAS_TEST_RESAMPLES", "19")) +const COPULA_TEST_TINY_RESAMPLES = min(COPULA_TEST_RESAMPLES, 9) + +struct MockHypothesis <: CopulaHypothesis end + +Copulas.testname(::MockHypothesis) = "Mock copula hypothesis test" +Copulas.nullhypothesis(::MockHypothesis) = "The mock null hypothesis holds." +Copulas._available_statistics(::MockHypothesis) = (:mean, :sum) +Copulas._available_calibrations(::MockHypothesis, ::Val{:mean}) = (:simulation,) +Copulas._available_calibrations(::MockHypothesis, ::Val{:sum}) = (:simulation,) +Copulas._teststatistic(::MockHypothesis, ::Val{:mean}, U::AbstractMatrix; kwargs...) = + sum(U) / length(U) +Copulas._teststatistic(::MockHypothesis, ::Val{:sum}, U::AbstractMatrix; kwargs...) = sum(U) + +function Copulas._simulation_sample(::MockHypothesis, U::AbstractMatrix, rng::Distributions.AbstractRNG) + sample = similar(U) + Random.rand!(rng, sample) + return sample +end + +@testset "Copula hypothesis tests [copula_tests]" begin + @testset "Extensible framework" begin + U = rand(Xoshiro(123), 2, 40) + test = CopulaTest(MockHypothesis(), U; N=COPULA_TEST_TINY_RESAMPLES, rng=Xoshiro(1)) + + @test test isa CopulaTest{MockHypothesis} + @test Copulas.default_statistic(MockHypothesis()) === :mean + @test Copulas.default_calibration(MockHypothesis(), Val(:mean)) === :simulation + @test Copulas.testname(test) == "Mock copula hypothesis test" + @test Copulas.nullhypothesis(test) == "The mock null hypothesis holds." + @test test.statistic === :mean + @test test.calibration === :simulation + @test StatsBase.nobs(test) == 40 + @test isfinite(teststatistic(test)) + @test 0 < pvalue(test) < 1 + + io = IOBuffer() + show(io, MIME("text/plain"), test) + printed = String(take!(io)) + @test occursin("Mock copula hypothesis test", printed) + @test occursin("The mock null hypothesis holds.", printed) + + other = CopulaTest(MockHypothesis(), U; statistic=:sum, N=COPULA_TEST_TINY_RESAMPLES, rng=Xoshiro(1)) + @test other isa CopulaTest{MockHypothesis} + @test other.statistic === :sum + @test other.calibration === :simulation + @test isfinite(teststatistic(other)) + @test IndependenceCopulaTest(U; N=2, rng=Xoshiro(1)).statistic === :cvm + + stat_err = try + CopulaTest(MockHypothesis(), U; statistic=:missing, N=2) + catch err + err + end + @test stat_err isa ArgumentError + @test occursin("Statistic :missing", sprint(showerror, stat_err)) + @test occursin("Available statistics: :mean, :sum", sprint(showerror, stat_err)) + + cal_err = try + CopulaTest(MockHypothesis(), U; calibration=:multiplier, N=2) + catch err + err + end + @test cal_err isa ArgumentError + @test occursin("Calibration :multiplier", sprint(showerror, cal_err)) + @test occursin("Available calibrations: :simulation", sprint(showerror, cal_err)) + + @test !(:testname in names(Copulas)) + end + + @testset "Tied margins are rejected" begin + Utied = [ + 0.1 0.1 0.4 0.7 + 0.2 0.3 0.6 0.8 + ] + + err = try + IndependenceCopulaTest(Utied; N=2, rng=Xoshiro(10),) + catch e + e + end + + @test err isa ArgumentError + @test occursin("tie-free margins", sprint(showerror, err)) + @test occursin("margin 1", sprint(showerror, err)) + + # Supplying pseudo-observations must not bypass the tie check. + @test_throws ArgumentError IndependenceCopulaTest(Utied; pseudo_values=true, N=2, rng=Xoshiro(11),) + end + + @testset "RNG reproducibility" begin + # Simulation calibration + Uind = rand(Xoshiro(810), IndependentCopula(2), 30) + + tsim1 = IndependenceCopulaTest( + Uind; + N=COPULA_TEST_TINY_RESAMPLES, + rng=Xoshiro(811), + ) + tsim2 = IndependenceCopulaTest( + Uind; + N=COPULA_TEST_TINY_RESAMPLES, + rng=Xoshiro(811), + ) + + @test pvalue(tsim1) == pvalue(tsim2) + + # Multiplier calibration + Uex = rand(Xoshiro(812), GaussianCopula(3, 0.5), 30) + + tmul1 = ExchangeabilityCopulaTest( + Uex; + N=COPULA_TEST_TINY_RESAMPLES, + rng=Xoshiro(813), + ) + tmul2 = ExchangeabilityCopulaTest( + Uex; + N=COPULA_TEST_TINY_RESAMPLES, + rng=Xoshiro(813), + ) + + @test pvalue(tmul1) == pvalue(tmul2) + + # Randomization calibration + Urad = rand(Xoshiro(814), GaussianCopula(3, 0.5), 30) + + tran1 = RadialSymmetryCopulaTest( + Urad; + N=COPULA_TEST_TINY_RESAMPLES, + rng=Xoshiro(815), + ) + tran2 = RadialSymmetryCopulaTest( + Urad; + N=COPULA_TEST_TINY_RESAMPLES, + rng=Xoshiro(815), + ) + + @test pvalue(tran1) == pvalue(tran2) + + # Parametric-bootstrap calibration + Cgof = ClaytonCopula(2, 2.5) + Ugof = rand(Xoshiro(816), Cgof, 30) + + tboot1 = GOFCopulaTest( + Cgof, + Ugof; + N=COPULA_TEST_TINY_RESAMPLES, + rng=Xoshiro(817), + ) + tboot2 = GOFCopulaTest( + Cgof, + Ugof; + N=COPULA_TEST_TINY_RESAMPLES, + rng=Xoshiro(817), + ) + + @test pvalue(tboot1) == pvalue(tboot2) + end + + @testset "IndependenceCopulaTest" begin + U0 = rand(Xoshiro(123), IndependentCopula(2), 80) + t0 = IndependenceCopulaTest(U0; N=COPULA_TEST_RESAMPLES, rng=Xoshiro(1)) + + @test t0 isa IndependenceCopulaTest + @test t0.statistic === :cvm + @test t0.calibration === :simulation + @test Copulas.testname(t0) == "Copula independence test" + @test StatsBase.nobs(t0) == 80 + @test t0.dimension == 2 + @test isfinite(teststatistic(t0)) + @test 0 < pvalue(t0) < 1 + + U1 = rand(Xoshiro(456), ClaytonCopula(2, 8.0), 80) + t1 = IndependenceCopulaTest(U1; N=COPULA_TEST_RESAMPLES, rng=Xoshiro(1)) + + @test teststatistic(t1) > teststatistic(t0) + @test pvalue(t1) <= 0.05 + + io = IOBuffer() + show(io, MIME("text/plain"), t1) + printed = String(take!(io)) + @test occursin("Statistic:", printed) + @test occursin("Observed value:", printed) + + @test_throws ArgumentError IndependenceCopulaTest(U0; statistic=:ks, N=9, rng=Xoshiro(1)) + @test_throws ArgumentError IndependenceCopulaTest(U0; calibration=:bootstrap, N=9, rng=Xoshiro(1)) + @test_throws ArgumentError IndependenceCopulaTest(U0; N=0, rng=Xoshiro(1)) + end + + @testset "ExchangeabilityCopulaTest" begin + @testset "Published Sn normalization" begin + Udet = [ + 0.2 0.5 0.8 + 0.3 0.9 0.6 + ] + + permutations = ((2, 1),) + expected = 1 / 3 + + @test Copulas._exchangeability_sn_statistic(Udet, permutations, :none,) ≈ expected + + hdet = Copulas.ExchangeabilityHypothesis(permutations=(2, 1), weight=:none,) + rep = Copulas._multiplier_representation(hdet, Val(:Sn), Udet,) + @test rep.scale == inv(size(Udet, 2)) + Tdet = ExchangeabilityCopulaTest(Udet; permutations=(2, 1), weight=:none, pseudo_values=true, N=1, rng=Xoshiro(779),) + @test teststatistic(Tdet) ≈ expected + end + + U2 = rand(Xoshiro(123), ClaytonCopula(2, 3.0), 80) + t2 = ExchangeabilityCopulaTest(U2; N=COPULA_TEST_TINY_RESAMPLES, rng=Xoshiro(1)) + + @test t2 isa ExchangeabilityCopulaTest + @test t2.statistic === :Sn + @test t2.calibration === :multiplier + @test t2.details.permutations === :G2 + @test t2.details.generator == ((2, 1),) + @test t2.details.weight === :wm2 + @test StatsBase.nobs(t2) == 80 + @test t2.dimension == 2 + @test Copulas.testname(t2) == "Copula exchangeability test" + @test isfinite(teststatistic(t2)) + @test 0 <= pvalue(t2) <= 1 + + tall = ExchangeabilityCopulaTest(U2; permutations=:all, N=COPULA_TEST_TINY_RESAMPLES, rng=Xoshiro(1)) + @test tall.details.generator == ((2, 1),) + + Ue = rand(Xoshiro(234), ClaytonCopula(3, 3.0), 120) + te = ExchangeabilityCopulaTest(Ue; N=COPULA_TEST_RESAMPLES, rng=Xoshiro(1)) + @test te.details.generator == ((2, 1, 3), (2, 3, 1)) + @test pvalue(te) > 0.05 + + tc = ExchangeabilityCopulaTest(Ue; permutations=(2, 1, 3), N=COPULA_TEST_TINY_RESAMPLES, rng=Xoshiro(1)) + @test tc.details.generator == ((2, 1, 3),) + + x = rand(Xoshiro(2), 120) + y = x .+ 0.04 .* randn(Xoshiro(3), 120) + z = rand(Xoshiro(4), 120) + Ua = permutedims(hcat(x, y, z)) + ta = ExchangeabilityCopulaTest(Ua; N=COPULA_TEST_RESAMPLES, rng=Xoshiro(1)) + @test teststatistic(ta) > teststatistic(te) + @test pvalue(ta) <= 0.05 + + io = IOBuffer() + show(io, MIME("text/plain"), ta) + printed = String(take!(io)) + @test occursin("Permutations:", printed) + @test occursin("Weight:", printed) + + @testset "Cost guard for all permutations" begin + # Small problems may still use all permutations. + @test Copulas._check_exchangeability_all_cost(:all, 3, 20) === nothing + @test Copulas._check_exchangeability_all_cost(:G2, 20, 5000) === nothing + + # 5! - 1 = 119 dense 1000×1000 Float64 matrices would require + # roughly 908 MiB before accounting for auxiliary allocations. + @test_throws ArgumentError Copulas._check_exchangeability_all_cost(:all, 5, 1000,) + + Ularge = rand(Xoshiro(781), 5, 1000) + err = try + ExchangeabilityCopulaTest(Ularge; permutations=:all, N=1, rng=Xoshiro(782),) + catch e + e + end + + @test err isa ArgumentError + @test occursin("permutations=:all", sprint(showerror, err)) + @test occursin("safety limit", sprint(showerror, err)) + + # Explicit collections must be guarded by their resolved size as + # well; otherwise they bypass the preflight specific to `:all`. + explicit_permutations = fill((2, 1, 3), 70) + resolved = Copulas._exchangeability_permutations(explicit_permutations, 3) + + @test length(resolved) == 70 + @test Copulas._check_exchangeability_matrix_cost(64, 1000) === nothing + @test_throws ArgumentError Copulas._check_exchangeability_matrix_cost(70, 1000) + + Uexplicit = rand(Xoshiro(783), 3, 1000) + explicit_err = try + ExchangeabilityCopulaTest( + Uexplicit; + permutations=explicit_permutations, + N=1, + rng=Xoshiro(784), + ) + catch e + e + end + + @test explicit_err isa ArgumentError + @test occursin("requested permutation collection", sprint(showerror, explicit_err)) + @test occursin("70 dense", sprint(showerror, explicit_err)) + @test occursin("safety limit", sprint(showerror, explicit_err)) + end + + @test_throws ArgumentError ExchangeabilityCopulaTest(U2; statistic=:Rn, N=9, rng=Xoshiro(1)) + @test_throws ArgumentError ExchangeabilityCopulaTest(U2; calibration=:randomization, N=9, rng=Xoshiro(1)) + @test_throws ArgumentError ExchangeabilityCopulaTest(U2; weight=:wm, N=9, rng=Xoshiro(1)) + @test_throws ArgumentError ExchangeabilityCopulaTest(U2; permutations=(1, 1), N=9, rng=Xoshiro(1)) + @test_throws ArgumentError ExchangeabilityCopulaTest(U2; permutations=(1, 2), N=9, rng=Xoshiro(1)) + @test_throws ArgumentError ExchangeabilityCopulaTest(U2; N=0, rng=Xoshiro(1)) + end + + @testset "RadialSymmetryCopulaTest" begin + @testset "Randomization reranking of induced ties" begin + sample = [ + 0.25 0.25 0.75 0.75 + 0.20 0.40 0.60 0.80 + ] + + V = Copulas._randomization_pseudos(Copulas.RadialSymmetryHypothesis(), sample,) + + @test V[1, :] ≈ [0.3, 0.3, 0.7, 0.7] # Average ranks in the first margin: + # (1.5, 1.5, 3.5, 3.5) / (4 + 1) + @test V[2, :] ≈ [0.2, 0.4, 0.6, 0.8] # Without ties, tied ranks coincide with ordinary ranks. + @test V[1, 1] == V[1, 2] + @test V[1, 3] == V[1, 4] + end + + @testset "Published Sn normalization" begin + Udet = [ + 0.2 0.5 0.8 + 0.3 0.9 0.6 + ] + + expected = 1 / 9 + + hdet = Copulas.RadialSymmetryHypothesis() + @test Copulas._teststatistic(hdet, Val(:Sn), Udet,) ≈ expected + Tdet = RadialSymmetryCopulaTest(Udet; pseudo_values=true, N=1, rng=Xoshiro(780),) + @test teststatistic(Tdet) ≈ expected + end + + Us = rand(Xoshiro(123), GaussianCopula(3, 0.5), 100) + ts = RadialSymmetryCopulaTest(Us; N=COPULA_TEST_RESAMPLES, rng=Xoshiro(1)) + + @test ts isa RadialSymmetryCopulaTest + @test ts.statistic === :Sn + @test ts.calibration === :randomization + @test ts.details.reflection_probability == 0.5 + @test StatsBase.nobs(ts) == 100 + @test ts.dimension == 3 + @test Copulas.testname(ts) == "Copula radial symmetry test" + @test isfinite(teststatistic(ts)) + @test 0 < pvalue(ts) < 1 + + Ua = rand(Xoshiro(456), ClaytonCopula(3, 4.0), 100) + ta = RadialSymmetryCopulaTest(Ua; N=COPULA_TEST_RESAMPLES, rng=Xoshiro(1)) + + @test teststatistic(ta) > teststatistic(ts) + + io = IOBuffer() + show(io, MIME("text/plain"), ta) + printed = String(take!(io)) + @test occursin("Reflection probability:", printed) + @test occursin("The copula is radially symmetric.", printed) + + @test_throws ArgumentError RadialSymmetryCopulaTest(Us; statistic=:cvm, N=9, rng=Xoshiro(1)) + @test_throws ArgumentError RadialSymmetryCopulaTest(Us; calibration=:multiplier, N=9, rng=Xoshiro(1)) + @test_throws ArgumentError RadialSymmetryCopulaTest(Us; N=0, rng=Xoshiro(1)) + end + + @testset "ExtremeValueCopulaTest" begin + @testset "Published Sn normalization" begin + Udet = [ + 0.2 0.5 0.8 + 0.3 0.6 0.9 + ] + + expected = 8 / 81 + + @test Copulas._extreme_value_sn_statistic(Udet, (2.0,)) ≈ expected + + hdet = Copulas.ExtremeValueHypothesis(; powers=2) + rep = Copulas._multiplier_representation(hdet, Val(:Sn), Udet) + + @test rep.scale == inv(size(Udet, 2)) + Tdet = ExtremeValueCopulaTest(Udet; powers=2, pseudo_values=true, N=1, rng=Xoshiro(778),) + @test teststatistic(Tdet) ≈ expected + end + + Uev = rand(Xoshiro(123), GumbelCopula(3, 3.0), 100) + tev = ExtremeValueCopulaTest(Uev; N=COPULA_TEST_RESAMPLES, rng=Xoshiro(1)) + + @test tev isa ExtremeValueCopulaTest + @test tev.statistic === :Sn + @test tev.calibration === :multiplier + @test tev.details.powers == (3.0, 4.0, 5.0) + @test tev.details.multiplier === :exponential + @test tev.details.derivative_bandwidth == inv(sqrt(100)) + @test StatsBase.nobs(tev) == 100 + @test tev.dimension == 3 + @test Copulas.testname(tev) == "Extreme-value copula test" + @test isfinite(teststatistic(tev)) + @test 0 < pvalue(tev) < 1 + + Ucl = rand(Xoshiro(456), ClaytonCopula(3, 3.0), 100) + tcl = ExtremeValueCopulaTest(Ucl; N=COPULA_TEST_RESAMPLES, rng=Xoshiro(1)) + + @test teststatistic(tcl) > teststatistic(tev) + @test pvalue(tcl) <= 0.05 + + tp = ExtremeValueCopulaTest(Uev; powers=2, N=COPULA_TEST_TINY_RESAMPLES, rng=Xoshiro(1)) + @test tp.details.powers == (2.0,) + + io = IOBuffer() + show(io, MIME("text/plain"), tcl) + printed = String(take!(io)) + @test occursin("Powers:", printed) + @test occursin("The copula belongs to the extreme-value class.", printed) + + @test_throws ArgumentError ExtremeValueCopulaTest(Uev; statistic=:cvm, N=9, rng=Xoshiro(1)) + @test_throws ArgumentError ExtremeValueCopulaTest(Uev; calibration=:simulation, N=9, rng=Xoshiro(1)) + @test_throws ArgumentError ExtremeValueCopulaTest(Uev; powers=1, N=9, rng=Xoshiro(1)) + @test_throws ArgumentError ExtremeValueCopulaTest(Uev; powers=(), N=9, rng=Xoshiro(1)) + @test_throws ArgumentError ExtremeValueCopulaTest(Uev; N=0, rng=Xoshiro(1)) + end + + @testset "GOFCopulaTest" begin + @testset "Published Sn normalization" begin + Udet = [ + 0.2 0.5 0.8 + 0.3 0.6 0.9 + ] + Cdet = IndependentCopula(2) + + expected = 647 / 2250 + + @test Copulas._gof_sn_statistic(Udet, Cdet) ≈ expected + Tdet = GOFCopulaTest(Cdet, Udet; pseudo_values=true, N=1, rng=Xoshiro(777),) + @test teststatistic(Tdet) ≈ expected + end + + U = rand(Xoshiro(123), ClaytonCopula(2, 3.0), 60) + Ts = GOFCopulaTest(ClaytonCopula(2, 3.0), U; + N=COPULA_TEST_TINY_RESAMPLES, rng=Xoshiro(1)) + + @test Ts isa GOFCopulaTest + @test Ts.hypothesis.kind === :simple + @test Ts.statistic === :Sn + @test Ts.calibration === :parametric_bootstrap + @test StatsBase.nobs(Ts) == 60 + @test Ts.dimension == 2 + @test Copulas.testname(Ts) == "Copula goodness-of-fit test" + @test isfinite(teststatistic(Ts)) + @test 0 < pvalue(Ts) < 1 + + M = fit(CopulaModel, ClaytonCopula, U; vcov=false) + Tc = GOFCopulaTest(M; N=COPULA_TEST_TINY_RESAMPLES, rng=Xoshiro(1)) + @test Tc.hypothesis.kind === :composite + @test Tc.hypothesis.model === M + @test 0 < pvalue(Tc) < 1 + + @testset "Composite GOF refits the tested sample" begin + Unew = rand(Xoshiro(905), ClaytonCopula(2, 1.2), 40) + Vnew = pseudos(Unew) + + Mexpected = Copulas._refit(M, Vnew) + Tnew = GOFCopulaTest(M, Unew; N=1, rng=Xoshiro(906),) + + @test Tnew.hypothesis.kind === :composite + @test Tnew.hypothesis.model !== M + @test StatsBase.coef(Tnew.hypothesis.model) ≈ StatsBase.coef(Mexpected) + @test teststatistic(Tnew) ≈ Copulas._gof_sn_statistic(Vnew, Copulas._copula_of(Mexpected)) + end + + @testset "Composite GOF preserves runtime fitting structure" begin + # NestedArchimedeanCopula stores its tree structure and inner + # generator families in the instance, so refitting from typeof(C) + # is insufficient. + Cnested = NestedArchimedeanCopula(Copulas.ClaytonGenerator(1.5); leaves=[1], children=[ClaytonCopula(2, 3.0) => [2, 3]],) + Unested = rand(Xoshiro(901), Cnested, 30) + + Mnested = fit(CopulaModel, Cnested, Unested; vcov=false, derived_measures=false,) + + Mnested_refit = Copulas._refit(Mnested, Unested) + Cnested_refit = Copulas._copula_of(Mnested_refit) + + @test Cnested_refit isa NestedArchimedeanCopula + @test Cnested_refit.leafdims == Cnested.leafdims + @test length(Cnested_refit.children) == length(Cnested.children) + @test Cnested_refit.children[1][2] == Cnested.children[1][2] + @test typeof(Cnested_refit.G).name.wrapper === typeof(Cnested.G).name.wrapper + @test typeof(Cnested_refit.children[1][1].G).name.wrapper === typeof(Cnested.children[1][1].G).name.wrapper + + Tnested = GOFCopulaTest(Mnested; N=1, rng=Xoshiro(903),) + + @test Tnested.hypothesis.kind === :composite + @test isfinite(teststatistic(Tnested)) + @test 0 <= pvalue(Tnested) <= 1 + + # SurvivalCopula stores the flip pattern in the instance rather + # than in its concrete type. + Csurvival = SurvivalCopula(ClaytonCopula(3, 2.5), (1, 3)) + Usurvival = rand(Xoshiro(902), Csurvival, 40) + + Msurvival = fit(CopulaModel, typeof(Csurvival), Usurvival; method=:itau, flips=Csurvival.flipmask, vcov=false, derived_measures=false,) + + Msurvival_refit = Copulas._refit(Msurvival, Usurvival) + Csurvival_refit = Copulas._copula_of(Msurvival_refit) + + @test Csurvival_refit isa SurvivalCopula + @test Csurvival_refit.flipmask == Csurvival.flipmask + @test Csurvival_refit.flipmask == (true, false, true) + + Tsurvival = GOFCopulaTest(Msurvival; N=1, rng=Xoshiro(904),) + + @test Tsurvival.hypothesis.kind === :composite + @test isfinite(teststatistic(Tsurvival)) + @test 0 <= pvalue(Tsurvival) <= 1 + + # Estimator-defining runtime keywords must also survive refitting. + Uchecker = rand(Xoshiro(907), GaussianCopula(2, 0.4), 36) + Mchecker = fit(CopulaModel, CheckerboardCopula, Uchecker; m=3, vcov=false, derived_measures=false,) + Mchecker_refit = Copulas._refit(Mchecker, Uchecker) + + @test Tuple(Mchecker_refit.result.m) == (3, 3) + + # If an empirical fit ranked raw input internally, GOF(model) must + # apply the same preprocessing to the stored fitting sample. + Xraw = [ + 0.12 0.91 0.35 0.67 0.48 0.76 + 0.88 0.22 0.71 0.41 0.59 0.13 + ] + + Mraw = fit(CopulaModel, CheckerboardCopula, Xraw; m=2, pseudo_values=false, vcov=false, derived_measures=false,) + + Traw = GOFCopulaTest(Mraw; N=1, rng=Xoshiro(908),) + Vraw = pseudos(Xraw) + Mraw_refit = Copulas._refit(Mraw, Vraw) + + @test Mraw.method_details.pseudo_values === false + @test Mraw_refit.method_details.pseudo_values === true + @test Tuple(Mraw_refit.result.m) == (2, 2) + @test teststatistic(Traw) ≈ Copulas._gof_sn_statistic(Vraw, Copulas._copula_of(Mraw)) + + @testset "Unreproducible custom parametrisation is rejected" begin + # An arbitrary user-supplied parametrisation may capture external + # state and cannot be reconstructed safely from the fitted result. + # Composite GOF therefore rejects such models rather than silently + # changing the estimator used by the bootstrap. + reparam = α -> begin + θ = exp(α[1]) + NestedArchimedeanCopula(Copulas.ClaytonGenerator(θ); leaves=[1], children=[ClaytonCopula(2, θ + one(θ)) => [2, 3]],) + end + + Mcustom = fit(CopulaModel, reparam, [log(1.5)], Unested; vcov=false, derived_measures=false,) + @test Mcustom.method_details._fit_spec === nothing + + refit_err = try + Copulas._refit(Mcustom, Unested) + catch err + err + end + + @test refit_err isa ArgumentError + @test occursin("reproducible fitting specification", sprint(showerror, refit_err),) + @test_throws ArgumentError GOFCopulaTest(Mcustom; N=1, rng=Xoshiro(909),) + end + + end + + io = IOBuffer() + show(io, MIME("text/plain"), Tc) + printed = String(take!(io)) + @test occursin("Hypothesis:", printed) + @test occursin("Fitted model:", printed) + + @test_throws ArgumentError GOFCopulaTest(ClaytonCopula(2, 3.0), U; statistic=:ks, N=2) + @test_throws ArgumentError GOFCopulaTest(M, U; calibration=:multiplier, N=2) + @test_throws ArgumentError GOFCopulaTest(M, U; N=0) + @test_throws DimensionMismatch GOFCopulaTest(ClaytonCopula(3, 3.0), U; N=2) + end +end diff --git a/test/runtests.jl b/test/runtests.jl index b5720211d..116a17740 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -209,6 +209,7 @@ testfiles = ( "operations/rosenblatt.jl", "operations/dependence.jl", "operations/fitting.jl", + "operations/hypothesis_testing.jl", "operations/nataf.jl", "extensions/expectation_maximization.jl" )