diff --git a/Project.toml b/Project.toml index 4ec85e6b..5cd45aa8 100644 --- a/Project.toml +++ b/Project.toml @@ -18,6 +18,7 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" +StatsFuns = "4c63d2b9-4356-54db-8cca-17b64c39e42c" TaylorSeries = "6aa5eb33-94cf-58f4-a9d0-e4b2c4fc25ea" WilliamsonTransforms = "48feb556-9bdd-43a2-8e10-96100ec25e22" @@ -39,6 +40,7 @@ Roots = "1, 2" SpecialFunctions = "2" StableRNGs = "1" StatsBase = "0.33, 0.34" +StatsFuns = "0.9, 1.3" TaylorSeries = "0.12, 0.13, 0.14, 0.15, 0.16, 0.17" Test = "1.6" TestItemRunner = "v0.2" diff --git a/README.md b/README.md index 66baa03d..f8c36d1c 100644 --- a/README.md +++ b/README.md @@ -77,6 +77,7 @@ There are competing packages in Julia, such as [`BivariateCopulas.jl`](https://g | - Classic Bivariate | ✔️ | ✔️ | ✔️ | | - Classic Multivariate | ✔️ | ✔️ | ❌ | | - Archimedeans | ✔️ (All of them) | ⚠️ Selected ones | ⚠️Selected ones | +| - Bivariate Extreme Value| ✔️ | ❌ | ❌ | | - Obscure Bivariate | ✔️ | ❌ | ❌ | | - Archimedean Chains | ❌ | ✔️ | ❌ | diff --git a/docs/make.jl b/docs/make.jl index f6906a04..0f13293e 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -28,6 +28,7 @@ makedocs(; "Elliptical Copulas" => "elliptical/generalities.md", "Archimedean Copulas" => "archimedean/generalities.md", "Liouville Copulas" => "Liouville.md", + "Extreme Value Copulas" => "extremevalue/generalities.md", "Empirical Copulas" => "empirical/generalities.md", "Dependence measures" => "dependence_measures.md", ], @@ -35,6 +36,7 @@ makedocs(; "Bestiary" => [ "elliptical/available_models.md", "archimedean/available_models.md", + "extremevalue/available_models.md", "empirical/available_models.md", "Other Copulas" => "miscellaneous.md", "Transformed Copulas" => "transformations.md", diff --git a/docs/src/assets/references.bib b/docs/src/assets/references.bib index 20ba245f..e2d4e96a 100644 --- a/docs/src/assets/references.bib +++ b/docs/src/assets/references.bib @@ -737,3 +737,112 @@ @article{Raftery2023 year={2023}, publisher={MDPI} } + +@article{tawn1988bivariate, + title={Bivariate extreme value theory: models and estimation}, + author={Tawn, Jonathan A}, + journal={Biometrika}, + volume={75}, + number={3}, + pages={397--415}, + year={1988}, + publisher={Oxford University Press} +} + +@article{mai2011bivariate, + title={Bivariate extreme-value copulas with discrete Pickands dependence measure}, + author={Mai, Jan-Frederik and Scherer, Matthias}, + journal={Extremes}, + volume={14}, + pages={311--324}, + year={2011}, + publisher={Springer} +} + +@article{nikoloulopoulos2009extreme, + title={Extreme value properties of multivariate t copulas}, + author={Nikoloulopoulos, Aristidis K and Joe, Harry and Li, Haijun}, + journal={Extremes}, + volume={12}, + pages={129--148}, + year={2009}, + publisher={Springer} +} + +@book{mai2012simulating, + title={Simulating copulas: stochastic models, sampling algorithms, and applications}, + author={Mai, Jan-Frederik and Scherer, Matthias}, + volume={4}, + year={2012}, + publisher={World Scientific} +} + +@article{husler1989maxima, + title={Maxima of normal random vectors: between independence and complete dependence}, + author={H{\"u}sler, J{\"u}rg and Reiss, Rolf-Dieter}, + journal={Statistics \& Probability Letters}, + volume={7}, + number={4}, + pages={283--286}, + year={1989}, + publisher={Elsevier} +} + +@article{galambos1975order, + title={Order statistics of samples from multivariate distributions}, + author={Galambos, Janos}, + journal={Journal of the American Statistical Association}, + volume={70}, + number={351a}, + pages={674--680}, + year={1975}, + publisher={Taylor \& Francis} +} + +@article{ghoudi1998proprietes, + title={Propri{\'e}t{\'e}s statistiques des copules de valeurs extr{\^e}mes bidimensionnelles}, + author={Ghoudi, Kilani and Khoudraji, Abdelhaq and Rivest, Et Louis-Paul}, + journal={Canadian Journal of Statistics}, + volume={26}, + number={1}, + pages={187--197}, + year={1998}, + publisher={Wiley Online Library} +} + +@inproceedings{gudendorf2010extreme, + title={Extreme-value copulas}, + author={Gudendorf, Gordon and Segers, Johan}, + booktitle={Copula Theory and Its Applications: Proceedings of the Workshop Held in Warsaw, 25-26 September 2009}, + pages={127--145}, + year={2010}, + organization={Springer} +} + +@article{Joe1990, + title={Families of min-stable multivariate exponential and multivariate extreme value distributions}, + author={Joe, Harry}, + journal={Statistics \& probability letters}, + volume={9}, + number={1}, + pages={75--81}, + year={1990}, + publisher={Elsevier} +} + +@article{deheuvels1991limiting, + title={On the limiting behavior of the Pickands estimator for bivariate extreme-value distributions}, + author={Deheuvels, Paul}, + journal={Statistics \& Probability Letters}, + volume={12}, + number={5}, + pages={429--439}, + year={1991}, + publisher={Elsevier} +} +@book{mai2014financial, + title={Financial engineering with copulas explained}, + author={Mai, Jan-Frederik and Scherer, Matthias}, + year={2014}, + publisher={Springer} +} \ No newline at end of file diff --git a/docs/src/extremevalue/available_models.md b/docs/src/extremevalue/available_models.md new file mode 100644 index 00000000..c549c86b --- /dev/null +++ b/docs/src/extremevalue/available_models.md @@ -0,0 +1,60 @@ +```@meta +CurrentModule = Copulas +``` + +# [Extreme Values Copulas](@id available_extreme_models) + +## `AsymGalambosCopula` +```@docs +AsymGalambosCopula +``` + +## `AsymLogCopula` +```@docs +AsymLogCopula +``` + +## `AsymMixedCopula` +```@docs +AsymMixedCopula +``` + +## `BC2Copula` +```@docs +BC2Copula +``` + +## `CuadrasAugeCopula` +```@docs +CuadrasAugeCopula +``` + +## `GalambosCopula` +```@docs +GalambosCopula +``` + +## `HuslerReissCopula` +```@docs +HuslerReissCopula +``` + +## `LogCopula` +```@docs +LogCopula +``` + +## `MixedCopula` +```@docs +MixedCopula +``` + +## `MOCopula` +```@docs +MOCopula +``` + +## `tEVCopula` +```@docs +tEVCopula +``` \ No newline at end of file diff --git a/docs/src/extremevalue/generalities.md b/docs/src/extremevalue/generalities.md new file mode 100644 index 00000000..1824c4bf --- /dev/null +++ b/docs/src/extremevalue/generalities.md @@ -0,0 +1,99 @@ +```@meta +CurrentModule = Copulas +``` + +# [Extreme Value Copulas](@id Extreme_theory) + +*Extreme value copulas* are fundamental in the study of rare and extreme events due to their ability to model dependency in situations of extreme risk. This package proposes a wide selection of bivariate extreme values copulas, while multivariate cases are not implemented yet. Feel free to open an issue and/or propose pull requests if you want an implementation of a multivariate case. + + +A bivariate extreme value copula [gudendorf2010extreme](@cite) $C$ is a copula that has the following caracteristic property: + +$$C(u_1^t, u_2^t)=C(u_1,u_2)^t, t > 0.$$ + +It can be represented through its stable tail dependence function $\ell(\cdot)$: + +$$C(u_1, u_2)=\exp\{-\ell(\log(u_1),\log(u_2))\},$$ + +or through a convex function $A: [0,1] \to [1/2, 1]$ satisfying $\max(t, t-1)\leq A(t) \leq 1,$ called its Pickands dependence function: + +$$C(u_1,u_2)=\exp\left\{\log(u_1u_2)A\left(\frac{\log(u_1)}{\log(u_1u_2)}\right)\right\},$$ + +In the context of bivariate extreme value copulas, the functions $\ell$ and $A$ are related as follows: + +$$\ell(u_1,u_2)=(u_1+u_2)A\left(\frac{u_1}{u_1+u_2}\right).$$ + +Therefore, in this implementation, it is sufficient to provide the Pickands dependence function $A$ to construct the implementation structure of an extreme value copula. + +In this package, there is an abstract type [`ExtremeValueCopula`](@ref) that provides a foundation for defining extreme value copulas. Many extreme value copulas are already implemented for you! See [this list](@ref available_extreme_models) to get an overview. + +If you do not find the one you need, you may define it yourself by subtyping `ExtremeValueCopula`. The API does not require much information, which is really convenient. Only a method for the pickand function `A(C::ExtremeValueCopula) = ...` is required. By providing this functions, you can easily create a new extreme value copula that fits your specific needs: + +```julia +struct MyExtremeValueCopula{P} <: ExtremeValueCopula{P} + θ::P +end + +A(C::ExtremeValueCopula, t) = (t^C.θ + (1 - t)^C.θ)^(1/C.θ) # This is the Pickands function of the Logistic (Gumbel) Copula +``` + +!!! info "Nomenclature information" + We have called `A()` the Pickands function, which is necessary for constructing the Extreme Value Copula. This binding is very generic and thus not exported from the package, you can use it through `Copulas.A()` and/or by importing it. + +# Advanced Concepts + +Here, we present some important concepts from the theory of extreme value copulas that are useful for the development of this package. + +Let $(X,Y) \sim C$ where $C$ is a bivariate extreme value copula. We have the following results: + +> **Proposition 1 (Ghoudi 1998 [ghoudi1998proprietes](@cite)):** Let $(X, Y) \sim C$, where $C$ is an extreme value copula. The joint distribution of $X$ and $Z = \frac{\log(X)}{\log(XY)}$ is given by: +> +> $$P(Z \leq z, X \leq x) =G(z,x)=\left(z + z(1-z)\frac{A'(z)}{A(z)}\right)x^{A(z)/z}, \quad 0\leq x,z \leq 1$$ +> +> where $A'(z)$ denotes the derivate of function $A(z)$ at point $z.$ + +Since $A$ is a convex function defined on $[0, 1]$ and satisfies $-1 \leq A'(z) \leq 1$, by extension, we define $A'(1)$ as the supremum of $A'(z)$ over $(0, 1)$. By setting $x = 1$ in the previous result, we obtain the marginal distribution of $Z$: +$$P(Z \leq z) = G_Z(z) = z + z(1 - z) \frac{A'(z)}{A(z)}, \quad 0 \leq z \leq 1.$$ + +This result was demonstrated by Deheuvels (1991) [deheuvels1991limiting](@cite) in the case where $A$ admits a second derivative. + + +## Simulation of Bivariate Extreme Value Distributions + +To simulate a bivariate extreme value distribution $C(x, y)$, remark that if $F_1$ and $F_2$ are univariate extreme value distributions, then the pair $( F_1^{-1}(X), F_2^{-1}(Y) )$ is distributed according to a bivariate extreme value distribution. The proposed algorithm in Ghoudi,1998, [ghoudi1998proprietes](@cite) allows simulating such a distribution. + +Assume $A$ has a second derivative, making the distribution absolutely continuous. In this case, $Z$ is also absolutely continuous and has a density $g_Z(z)$ given by: + +$$g_Z(z) = \frac{d}{dz} G_Z(z) = 1 + (1 - z)^{-1} \left(A(z) - z A'(z)\right)$$ + +The conditional distribution of $W$ given $Z$ is: + +$$F(w|z) = \frac{1}{g_Z(z)} \frac{d}{dz} F(z, w),$$ + +which simplifies to: + +$$F(w|z) = w \frac{z(1 - z) A'(z)}{A(z) g_Z(z)} + (w - w \log w) \left(1 - \frac{z(1 - z) A''(z)}{A(z) g_Z(z)} \right)$$ + +Given $Z$, the distribution of $W$ is uniform on $(0, 1)$ with probability $p(Z)$ and equals the product of two independent uniforms on $(0, 1)$ with probability $1 - p(Z)$, where: + +$$p(z) = \frac{z(1 - z) A'(z)}{A(z) g_Z(z)}$$ + +Since $g_Z(z)$ is the derivative of the cumulative distribution function of $Z$, it holds that $0 \leq p(z) \leq 1$. + +For the class of Extreme Value Copulas, We follow the methodology proposed by Ghoudi,1998. page 191. [ghoudi1998proprietes](@cite). Here, is a detailed algorithm for sampling from bivariate Extreme Value Copulas: + +> **Algotithm 1: Bivariate Extreme Value Copulas** +> +> *(1)* Simulate $U_1,U_2 \sim \mathcal{U}[0,1]$ +> +> *(2)* Simulate $Z \sim G_Z(z)$ +> +> *(3)* Select $W=U_1$ with probability $p(Z)$ and $W=U_1U_2$ with probability $1-p(Z)$ +> +> *(4)* Return $X=W^{Z/A(Z)}$ and $Y=W^{(1-Z)/A(Z)}$ + +Note that all functions present in the algorithm were previously defined in previous sections to ensure that the implemented methodology has a solid theoretical basis. + +```@docs +ExtremeValueCopula +``` diff --git a/src/Copulas.jl b/src/Copulas.jl index 1d1571a4..8f3fd192 100644 --- a/src/Copulas.jl +++ b/src/Copulas.jl @@ -7,6 +7,7 @@ module Copulas import Roots import Distributions import StatsBase + import StatsFuns import TaylorSeries import ForwardDiff import Cubature @@ -49,6 +50,7 @@ module Copulas include("UnivariateDistribution/AlphaStable.jl") include("UnivariateDistribution/ClaytonWilliamsonDistribution.jl") include("UnivariateDistribution/WilliamsonFromFrailty.jl") + include("UnivariateDistribution/ExtremeDist.jl") # Archimedean generators include("Generator.jl") @@ -80,6 +82,33 @@ module Copulas InvGaussianCopula, JoeCopula + # bivariate Extreme Value Copulas + include("ExtremeValueCopula.jl") + include("ExtremeValueCopulas/AsymGalambosCopula.jl") + include("ExtremeValueCopulas/AsymLogCopula.jl") + include("ExtremeValueCopulas/AsymMixedCopula.jl") + include("ExtremeValueCopulas/BC2Copula.jl") + include("ExtremeValueCopulas/CuadrasAugeCopula.jl") + include("ExtremeValueCopulas/GalambosCopula.jl") + include("ExtremeValueCopulas/HuslerReissCopula.jl") + include("ExtremeValueCopulas/LogCopula.jl") + include("ExtremeValueCopulas/MixedCopula.jl") + include("ExtremeValueCopulas/MOCopula.jl") + include("ExtremeValueCopulas/tEVCopula.jl") + + export AsymGalambosCopula, + AsymLogCopula, + AsymMixedCopula, + BC2Copula, + CuadrasAugeCopula, + GalambosCopula, + HuslerReissCopula, + LogCopula, + MixedCopula, + MOCopula, + tEVCopula + + # Subsetting include("SubsetCopula.jl") # not exported yet. diff --git a/src/ExtremeValueCopula.jl b/src/ExtremeValueCopula.jl new file mode 100644 index 00000000..2ddf19a8 --- /dev/null +++ b/src/ExtremeValueCopula.jl @@ -0,0 +1,167 @@ +""" + ExtremeValueCopula{P} + +Fields: + - P::Parameters: Parameters that define the copula. + +Constructor: + ExtremeValueCopula(P) + +Represents a bivariate extreme value copula parameterized by `P`. Extreme value copulas are used to model the dependence structure between two random variables in the tails of their distribution, making them particularly useful in risk management, environmental studies, and finance. + +In the bivariate case, an extreme value copula can be expressed as: + +``math +C(u, v) = \\exp(-\\ell(\\log(u), \\log(v))). +`` +where ``\\ell(\\cdot)`` is a tail dependence function associated with the bivariate extreme value copula. Furthermore, ``A(t)`` is a function ``A: [0, 1] \\to [0.5, 1] `` that is convex on the interval [0,1] and satisfies the boundary conditions ``A(0) = A(1) = 1``. This is denominated Pickands representation or Pickands function. + +It is possible to relate these functions in the following way + +``math +\\ell(u, v) = \\frac{u}{u+v}A\\left(\\frac{u}{u+v}\\right). +`` + + +In this way, in order to define a bivariate copula of extreme values, it is only necessary to introduce the function ``A``. + +A generic bivariate Extreme Values copula can be constructed as follows: + +```julia +struct GalambosCopula{P} <: ExtremeValueCopula{P} +A(C::GalambosCopula, t::Real) = 1 - (t^(-C.θ) + (1 - t)^(-C.θ))^(-1/C.θ) # You can define your own Pickands representation +param = 2.5 +C = GalambosCopula(param) +``` + +The obtained model can be used as follows: + +```julia +samples = rand(C,1000) # sampling +cdf(C,samples) # cdf +pdf(C,samples) # pdf +``` + +References: + +* [gudendorf2010extreme](@cite) G., & Segers, J. (2010). Extreme-value copulas. In Copula Theory and Its Applications (pp. 127-145). Springer. +* [joe2014](@cite) Joe, H. (2014). Dependence Modeling with Copulas. CRC Press. +* [mai2014financial](@cite) Mai, J. F., & Scherer, M. (2014). Financial engineering with copulas explained (p. 168). London: Palgrave Macmillan. +""" + +abstract type ExtremeValueCopula{P} <: Copula{2} end + +# Función genérica para A +function A(C::ExtremeValueCopula, t::Real) + throw(ArgumentError("Function A must be defined for specific copula")) +end + +function dA(C::ExtremeValueCopula, t::Real) + ForwardDiff.derivative(t -> A(C, t), t) +end + +function d²A(C::ExtremeValueCopula, t::Real) + ForwardDiff.derivative(t -> dA(C, t), t) +end + +function ℓ(C::ExtremeValueCopula, t::Vector) + sumu = sum(t) + vectw = t[1] / sumu + return sumu * A(C, vectw) +end + +# Función CDF para ExtremeValueCopula +function _cdf(C::ExtremeValueCopula, u::AbstractArray{<:Real}) + t = abs.(log.(u)) # 0 <= u <= 1 so abs == neg, but return corectly 0 instead of -0 when u = 1. + return exp(-ℓ(C, t)) +end + +# Función genérica para calcular derivadas parciales de ℓ +function D_B_ℓ(C::ExtremeValueCopula, t::Vector{Float64}, B::Vector{Int}) + f = x -> ℓ(C, x) + + if length(B) == 1 + return ForwardDiff.gradient(f, t)[B[1]] + elseif length(B) == 2 + return ForwardDiff.hessian(f, t)[B[1], B[2]] + else + throw(ArgumentError("Higher order partial derivatives are not required for bivariate case")) + end +end + +# Función PDF para ExtremeValueCopula usando ℓ +function _pdf(C::ExtremeValueCopula, u::AbstractArray{<:Real}) + t = -log.(u) + c = exp(-ℓ(C, t)) + D1 = D_B_ℓ(C, t, [1]) + D2 = D_B_ℓ(C, t, [2]) + D12 = D_B_ℓ(C, t, [1, 2]) + return c * (-D12 + D1 * D2) / (u[1] * u[2]) +end +function Distributions._logpdf(C::ExtremeValueCopula, u::AbstractArray{<:Real}) + t = -log.(u) + c = exp(-ℓ(C, t)) + D1 = D_B_ℓ(C, t, [1]) + D2 = D_B_ℓ(C, t, [2]) + D12 = D_B_ℓ(C, t, [1, 2]) + return log(c) + log(-D12 + D1 * D2) - log(u[1] * u[2]) +end +# Definir la función para calcular τ +function τ(C::ExtremeValueCopula) + integrand(x) = begin + a = A(C, x) + da = dA(C, x) + return (x * (1 - x) / a) * da + end + + integrate, _ = QuadGK.quadgk(integrand, 0.0, 1.0) + return integrate +end + +function ρₛ(C::ExtremeValueCopula) + integrand(x) = 1 / (1 + A(C, x))^2 + + integral, _ = QuadGK.quadgk(integrand, 0, 1) + + ρs = 12 * integral - 3 + return ρs +end +# Función para calcular el coeficiente de dependencia en el límite superior +function λᵤ(C::ExtremeValueCopula) + return 2(1 - A(C, 0.5)) +end + +function λₗ(C::ExtremeValueCopula) + if A(C, 0.5) > 0.5 + return 0 + else + return 1 + end +end + +function probability_z(C::ExtremeValueCopula, z) + num = z*(1 - z)*d²A(C, z) + dem = A(C, z)*_pdf(ExtremeDist(C), z) + p = num / dem + return clamp(p, 0.0, 1.0) +end + +function Distributions._rand!(rng::Distributions.AbstractRNG, C::ExtremeValueCopula, x::AbstractVector{T}) where {T<:Real} + u1, u2 = rand(rng, Distributions.Uniform(0,1), 2) + z = rand(rng, ExtremeDist(C)) + p = probability_z(C, z) + if p < -eps() || p > eps() + p = 0 + end + c = rand(rng, Distributions.Bernoulli(p)) + w = 0 + if c == 1 + w = u1 + else + w = u1*u2 + end + a = A(C, z) + x[1] = w^(z/a) + x[2] = w^((1-z)/a) + return x +end \ No newline at end of file diff --git a/src/ExtremeValueCopulas/AsymGalambosCopula.jl b/src/ExtremeValueCopulas/AsymGalambosCopula.jl new file mode 100644 index 00000000..986825d0 --- /dev/null +++ b/src/ExtremeValueCopulas/AsymGalambosCopula.jl @@ -0,0 +1,59 @@ +""" + AsymGalambosCopula{P} + +Fields: + + - α::Real - Dependency parameter + - θ::Vector - Asymmetry parameters (size 2) + +Constructor + + AsymGalambosCopula(α, θ) + +The Asymmetric bivariate Galambos copula is parameterized by one dependence parameter ``\\alpha \\in [0, \\infty]`` and two asymmetry parameters ``\\theta_{i} \\in [0,1], i=1,2``. It is an Extreme value copula with Pickands function: + +```math +\\A(t) = 1 - ((\\theta_1t)^{-\\alpha}+(\\theta_2(1-t))^{-\\alpha})^{-\\frac{1}{\\alpha}} +``` + +It has a few special cases: + +- When α = 0, it is the Independent Copula +- When θ₁ = θ₂ = 0, it is the Independent Copula +- When θ₁ = θ₂ = 1, it is the Galambos Copula + +References: +* [Joe1990](@cite) Families of min-stable multivariate exponential and multivariate extreme value distributions. Statist. Probab, 1990. +""" +struct AsymGalambosCopula{P} <: ExtremeValueCopula{P} + α::P # Dependence parameter + θ::Vector{P} # Asymmetry parameters (size 2) + function AsymGalambosCopula(α::P, θ::Vector{P}) where {P} + if length(θ) != 2 + throw(ArgumentError("The vector θ must have 2 elements for the bivariate case")) + elseif !(0 <= α) + throw(ArgumentError("The parameter α must be greater than or equal to 0")) + elseif !(0 <= θ[1] <= 1) || !(0 <= θ[2] <= 1) + throw(ArgumentError("All parameters θ must be in the interval [0, 1]")) + elseif α == 0 || (θ[1] == 0 && θ[2] == 0) + return IndependentCopula(2) + elseif θ[1] == 1 && θ[2] == 1 + return GalambosCopula(α) + else + return new{P}(α, θ) + end + end +end + +function A(C::AsymGalambosCopula, t::Real) + α = C.α + θ = C.θ + + term1 = (θ[1] * t)^(-α) + term2 = (θ[2] * (1 - t))^(-α) + + inner_term = term1 + term2 + + result = 1 - inner_term^(-1 / α) + return result +end \ No newline at end of file diff --git a/src/ExtremeValueCopulas/AsymLogCopula.jl b/src/ExtremeValueCopulas/AsymLogCopula.jl new file mode 100644 index 00000000..329f75fc --- /dev/null +++ b/src/ExtremeValueCopulas/AsymLogCopula.jl @@ -0,0 +1,44 @@ +""" + AsymLogCopula{P} + +Fields: + + - α::Real - Dependency parameter + - θ::Vector - Asymmetry parameters (size 2) + +Constructor + + AsymLogCopula(α, θ) + +The Asymmetric bivariate Logistic copula is parameterized by one dependence parameter ``\\alpha \\in [1, \\infty]`` and two asymmetry parameters ``\\theta_{i} \\in [0,1], i=1,2``. It is an Extreme value copula with Pickands dependence function: + +```math +A(t) = (\\theta_1^{\\alpha}(1-t)^{\\alpha} + \\theta_2^{\\alpha}t^{\\alpha})^{\\frac{1}{\\alpha}} + (\\theta_1 - \\theta_2)t + 1 - \\theta_1 +``` + +References: +* [tawn1988bivariate](@cite) : Tawn, Jonathan A. "Bivariate extreme value theory: models and estimation." Biometrika 75.3 (1988): 397-415. +""" +struct AsymLogCopula{P} <: ExtremeValueCopula{P} + α::P # Dependence Parameter + θ::Vector{P} # Asymmetry parameters (size 2) + function AsymLogCopula(α::P, θ::Vector{P}) where {P} + if length(θ) != 2 + throw(ArgumentError("The vector θ must have 2 elements for the bivariate case")) + elseif !(1 <= α) + throw(ArgumentError("The parameter α must be greater than or equal to 1")) + elseif !(0 <= θ[1] <= 1) || !(0 <= θ[2] <= 1) + throw(ArgumentError("All parameters θ must be in the interval [0, 1]")) + else + return new{P}(α, θ) + end + end +end + +function A(C::AsymLogCopula, t::Real) + α = C.α + θ = C.θ + + A = ((θ[1]^α)*(1-t)^α + (θ[2]^α)*(t^α))^(1/α)+(θ[1]- θ[2])*t + 1 -θ[1] + return A +end \ No newline at end of file diff --git a/src/ExtremeValueCopulas/AsymMixedCopula.jl b/src/ExtremeValueCopulas/AsymMixedCopula.jl new file mode 100644 index 00000000..e3b524e2 --- /dev/null +++ b/src/ExtremeValueCopulas/AsymMixedCopula.jl @@ -0,0 +1,59 @@ +""" + AsymMixedCopula{P} + +Fields: + + - θ::Vector - parameters (size 2) + +Constructor + + AsymMixedCopula(θ) + +The Asymmetric bivariate Mixed copula is parameterized by two parameters ``\\theta_{i}, i=1,2`` that must meet the following conditions: +* θ₁ ≥ 0 +* θ₁+θ₂ ≤ 1 +* θ₁+2θ₂ ≤ 1 +* θ₁+3θ₂ ≥ 0 + +It is an Extreme value copula with Pickands dependence function: + +```math +A(t) = \\theta_{2}t^3 + \\theta_{1}t^2-(\\theta_1+\\theta_2)t+1 +``` + +It has a few special cases: + +- When θ₁ = θ₂ = 0, it is the Independent Copula +- When θ₂ = 0, it is the Mixed Copula + +References: +* [tawn1988bivariate](@cite) : Tawn, Jonathan A. "Bivariate extreme value theory: models and estimation." Biometrika 75.3 (1988): 397-415. +""" +struct AsymMixedCopula{P} <: ExtremeValueCopula{P} + θ::Vector{P} # Asymmetry parameters + + function AsymMixedCopula(θ::Vector{P}) where {P} + if length(θ) != 2 + throw(ArgumentError("The vector θ must have 2 elements for the bivariate case")) + elseif !(0 <= θ[1]) + throw(ArgumentError("The parameter θ₁ must be greater than or equal to 0")) + elseif !(θ[1]+θ[2] <= 1) + throw(ArgumentError("the sum of θ₁+θ₂ ≤ 1")) + elseif !(θ[1]+2*θ[2] <= 1) + throw(ArgumentError("the sum of θ₁+2θ₂ ≤ 1")) + elseif !(0 <= θ[1]+3*θ[2]) + throw(ArgumentError("the sum of 0 ≤ θ₁+3θ₂")) + elseif θ[1] == 0 && θ[2] == 0 + return IndependentCopula(2) + else + return new{P}(θ) + end + end +end + +function A(C::AsymMixedCopula, t::Real) + θ = C.θ + + a = θ[2]*t^3 + θ[1]*t^2-(θ[1]+θ[2])*t+1 + return a +end \ No newline at end of file diff --git a/src/ExtremeValueCopulas/BC2Copula.jl b/src/ExtremeValueCopulas/BC2Copula.jl new file mode 100644 index 00000000..9ffdb7b1 --- /dev/null +++ b/src/ExtremeValueCopulas/BC2Copula.jl @@ -0,0 +1,77 @@ +""" + BC2Copula{P} + +Fields: + + - θ1::Real - parameter + - θ1::Real - parameter + +Constructor + + BC2Copula(θ1, θ2) + +The bivariate BC₂ copula is parameterized by two parameters ``\\theta_{i} \\in [0,1], i=1,2``. It is an Extreme value copula with Pickands dependence function: + +```math +A(t) = \\max\\{\\theta_1 t, \\theta_2(1-t) \\} + \\max\\{(1-\\theta_1)t, (1-\\theta_2)(1-t)\\} +``` + +References: +* [mai2011bivariate](@cite) Mai, J. F., & Scherer, M. (2011). Bivariate extreme-value copulas with discrete Pickands dependence measure. Extremes, 14, 311-324. Springer, 2011. +""" +struct BC2Copula{P} <: ExtremeValueCopula{P} + θ1::P + θ2::P + + function BC2Copula(θ::Vararg{Real}) + if length(θ) !== 2 + throw(ArgumentError("BC2Copula requires only 2 arguments.")) + end + T = promote_type(typeof(θ[1]), typeof(θ[2])) + θ1, θ2 = T(θ[1]), T(θ[2]) + + if !(0 <= θ1 <= 1) || !(0 <= θ2 <= 1) + throw(ArgumentError("All θ parameters must be in [0,1]")) + end + return new{T}(θ1, θ2) + end +end + +function ℓ(C::BC2Copula, t::Vector) + θ1, θ2 = C.θ1, C.θ2 + t₁, t₂ = t + return max(θ1*t₁, θ2*t₂) + max((1-θ1)*t₁, (1-θ2)*t₂) +end + +function A(C::BC2Copula, t::Real) + θ1, θ2 = C.θ1, C.θ2 + return max(θ1*t, θ2*(1-t)) + max((1-θ1)*t, (1-θ2)*(1-t)) +end + +function dA(C::BC2Copula, t::Float64) + θ1, θ2 = C.θ1, C.θ2 + + # Conditions for the derivative of the first part + if θ1*t >= θ2*(1-t) + f1_derivative = θ1 + else + f1_derivative = -θ2 + end + + # Conditions for the derivative of the second part + if (1-θ1)*t >= (1-θ2)*(1-t) + f2_derivative = 1-θ1 + else + f2_derivative = -(1-θ2) + end + + return f1_derivative + f2_derivative +end + +function Distributions._rand!(rng::Distributions.AbstractRNG, C::BC2Copula, u::AbstractVector{T}) where {T<:Real} + θ1, θ2 = C.θ1, C.θ2 + v1, v2 = rand(rng, Distributions.Uniform(0,1), 2) + u[1] = max(v1^(1/θ1),v2^(1/(1-θ1))) + u[2] = max(v1^(1/θ2),v2^(1/(1-θ2))) + return u +end \ No newline at end of file diff --git a/src/ExtremeValueCopulas/CuadrasAugeCopula.jl b/src/ExtremeValueCopulas/CuadrasAugeCopula.jl new file mode 100644 index 00000000..7e8f1419 --- /dev/null +++ b/src/ExtremeValueCopulas/CuadrasAugeCopula.jl @@ -0,0 +1,58 @@ +""" + CuadrasAugeCopula{P} + +Fields: + + - α::Real - parameter + +Constructor + + CuadrasAugeCopula(α) + +The bivariate Cuadras Auge copula is parameterized by ``\\alpha \\in [0,1]``. It is an Extreme value copula with Pickands dependence function: + +```math +A(t) = \\max\\{t, 1-t \\} + (1-\\theta)\\max\\{t, 1-t\\} +``` + +References: +* [mai2012simulating](@cite) Mai, J. F., & Scherer, M. (2012). Simulating copulas: stochastic models, sampling algorithms, and applications (Vol. 4). World Scientific. +""" +struct CuadrasAugeCopula{P} <: ExtremeValueCopula{P} + θ::P # Copula parameter + function CuadrasAugeCopula(θ) + if !(0 <= θ <= 1) + throw(ArgumentError("Theta must be in [0,1]")) + elseif θ == 0 + return IndependentCopula(2) + elseif θ == 1 + return MCopula(2) + else + return new{typeof(θ)}(θ) + end + end +end + +A(C::CuadrasAugeCopula, t::Real) = max(t, 1-t) + (1-C.θ) * min(t, 1-t) + +dA(C::CuadrasAugeCopula, t::Real) = t <= 0.5 ? -C.θ : C.θ + +τ(C::CuadrasAugeCopula) = C.θ/(2-C.θ) + +ρₛ(C::CuadrasAugeCopula) = (3*C.θ)/(4-C.θ) + +# specific ℓ específica of Cuadras-Augé Copula +function ℓ(C::CuadrasAugeCopula, t::Vector) + θ = C.θ + t₁, t₂ = t + return max(t₁, t₂) + (1-θ) * min(t₁, t₂) +end + +function Distributions._rand!(rng::Distributions.AbstractRNG, C::CuadrasAugeCopula, x::AbstractVector{T}) where {T<:Real} + θ = C.θ + E₁, E₂ = rand(rng, Distributions.Exponential(θ/(1-θ)),2) + E₁₂ = rand(rng, Distributions.Exponential()) + x[1] = exp(-(1/θ)*min(E₁,E₁₂)) + x[2] = exp(-(1/θ)*min(E₂,E₁₂)) + return x +end \ No newline at end of file diff --git a/src/ExtremeValueCopulas/GalambosCopula.jl b/src/ExtremeValueCopulas/GalambosCopula.jl new file mode 100644 index 00000000..b58e70ac --- /dev/null +++ b/src/ExtremeValueCopulas/GalambosCopula.jl @@ -0,0 +1,41 @@ +""" + GalambosCopula{P} + +Fields: + + - θ::Real - parameter + +Constructor + + GalambosCopula(θ) + +The bivariate Galambos copula is parameterized by ``\\alpha \\in [0,\\infty)``. It is an Extreme value copula with Pickands dependence function: + +```math +A(t) = 1 - (t^{-\\theta}+(1-t)^{-\\theta})^{-\\frac{1}{\\theta}} +``` + +It has a few special cases: + +- When θ = 0, it is the Independent Copula +- When θ = ∞, it is the M Copula (Upper Frechet-Hoeffding bound) + +References: +* [galambos1975order](@cite) Galambos, J. (1975). Order statistics of samples from multivariate distributions. Journal of the American Statistical Association, 70(351a), 674-680. +""" +struct GalambosCopula{P} <: ExtremeValueCopula{P} + θ::P # Copula parameter + function GalambosCopula(θ) + if θ < 0 + throw(ArgumentError("Theta must be >= 0")) + elseif θ == 0 + return IndependentCopula(2) + elseif θ == Inf + return MCopula(2) + else + return new{typeof(θ)}(θ) + end + end +end + +A(C::GalambosCopula, t::Real) = 1 - (t^(-C.θ) + (1 - t)^(-C.θ))^(-1/C.θ) \ No newline at end of file diff --git a/src/ExtremeValueCopulas/HuslerReissCopula.jl b/src/ExtremeValueCopulas/HuslerReissCopula.jl new file mode 100644 index 00000000..30f5a5cb --- /dev/null +++ b/src/ExtremeValueCopulas/HuslerReissCopula.jl @@ -0,0 +1,72 @@ +""" + HuslerReissCopula{P} + +Fields: + + - θ::Real - parameter + +Constructor + + HuslerReissCopula(θ) + +The bivariate Husler-Reiss copula is parameterized by ``\\theta \\in [0,\\infty)``. It is an Extreme value copula with Pickands dependence function: + +```math +A(t) = t\\Phi(\\theta^{-1}+\\frac{1}{2}\\theta\\log(\\frac{t}{1-t})) +(1-t)\\Phi(\\theta^{-1}+\\frac{1}{2}\\theta\\log(\\frac{1-t}{t})) +``` +Where ``\\Phi``is the cumulative distribution function (CDF) of the standard normal distribution. + +It has a few special cases: + +- When θ = 0, it is the Independent Copula +- When θ = ∞, it is the M Copula (Upper Frechet-Hoeffding bound) + +References: +* [husler1989maxima](@cite) Hüsler, J., & Reiss, R. D. (1989). Maxima of normal random vectors: between independence and complete dependence. Statistics & Probability Letters, 7(4), 283-286. +""" +struct HuslerReissCopula{P} <: ExtremeValueCopula{P} + θ::P # Copula parameter + function HuslerReissCopula(θ) + if θ < 0 + throw(ArgumentError("Theta must be ≥ 0")) + elseif θ == 0 + return IndependentCopula(2) + elseif θ == Inf + return MCopula(2) + else + return new{typeof(θ)}(θ) + end + end +end +# specific ℓ funcion of HuslerReissCopula +function ℓ(H::HuslerReissCopula, t::Vector) + θ = H.θ + t₁, t₂ = t + return t₁*Distributions.cdf(Distributions.Normal(),θ^(-1)+0.5*θ*log(t₁/t₂))+t₂*Distributions.cdf(Distributions.Normal(),θ^(-1)+0.5*θ*log(t₂/t₁)) +end + +# specific A funcion of HuslerReissCopula +function A(H::HuslerReissCopula, t::Real) + θ = H.θ + term1 = t * Distributions.cdf(Distributions.Normal(), θ^(-1) + 0.5 * θ * log(t / (1 - t))) + term2 = (1 - t) * Distributions.cdf(Distributions.Normal(), θ^(-1) + 0.5 * θ * log((1 - t) / t)) + + a = term1 + term2 + + + return a +end + +function dA(H::HuslerReissCopula, t::Real) + θ = H.θ + # Derivada of A(x) respected to t + dA_term1 = Distributions.cdf(Distributions.Normal(), θ^(-1) + 0.5 * θ * log(t / (1 - t))) + + t * Distributions.pdf(Distributions.Normal(), θ^(-1) + 0.5 * θ * log(t / (1 - t))) * (0.5 * θ * (1 / t + 1 / (1 - t))) + + dA_term2 = -Distributions.cdf(Distributions.Normal(), θ^(-1) + 0.5 * θ * log((1 - t) / t)) + + (1 - t) * Distributions.pdf(Distributions.Normal(), θ^(-1) + 0.5 * θ * log((1 - t) / t)) * (0.5 * θ * (-1 / t - 1 / (1 - t))) + + da = dA_term1 + dA_term2 + + return da +end \ No newline at end of file diff --git a/src/ExtremeValueCopulas/LogCopula.jl b/src/ExtremeValueCopulas/LogCopula.jl new file mode 100644 index 00000000..f9483792 --- /dev/null +++ b/src/ExtremeValueCopulas/LogCopula.jl @@ -0,0 +1,49 @@ +""" + LogCopula{P} + +Fields: + + - θ::Real - parameter + +Constructor + + LogCopula(θ) + +The bivariate Logistic copula (or Gumbel Copula) is parameterized by ``\\theta \\in [1,\\infty)``. It is an Extreme value copula with Pickands dependence function: + +```math +A(t) = (t^{\\theta}+(1-t)^{\\theta})^{\\frac{1}{\\theta}} +``` + +It has a few special cases: +- When θ = 1, it is the IndependentCopula +- When θ = ∞, is is the MCopula (Upper Frechet-Hoeffding bound) + +References: +* [tawn1988bivariate](@cite) : Tawn, Jonathan A. "Bivariate extreme value theory: models and estimation." Biometrika 75.3 (1988): 397-415. +""" +struct LogCopula{P} <: ExtremeValueCopula{P} + θ::P # Copula parameter + function LogCopula(θ) + if !(1 <= θ) + throw(ArgumentError(" The param θ must be in [1, ∞)")) + elseif θ == 1 + return IndependentCopula(2) + elseif θ == Inf + return MCopula(2) + else + return new{typeof(θ)}(θ) + end + end +end +# # specific ℓ funcion of LogCopula +function ℓ(G::LogCopula, t::Vector) + θ = G.θ + t₁, t₂ = t + return (t₁^θ + t₂^θ)^(1/θ) +end +# # specific A funcion of HuslerReissCopula +function A(C::LogCopula, t::Real) + θ = C.θ + return (t^θ + (1 - t)^θ)^(1/θ) +end \ No newline at end of file diff --git a/src/ExtremeValueCopulas/MOCopula.jl b/src/ExtremeValueCopulas/MOCopula.jl new file mode 100644 index 00000000..0b74ac19 --- /dev/null +++ b/src/ExtremeValueCopulas/MOCopula.jl @@ -0,0 +1,69 @@ +""" + MOCopula{P} + +Fields: + + - λ1::Real - parameter + - λ2::Real - parameter + - λ12::Real - parameter + +Constructor + + MOCopula(θ) + +The bivariate Marshall-Olkin copula is parameterized by ``\\lambda_i \\in [0,\\infty), i = 1, 2, \\{1,2\\}``. It is an Extreme value copula with Pickands dependence function: + +```math +A(t) = \\frac{\\lambda_1 (1-t)}{\\lambda_1 + \\lambda_{1,2}} + \\frac{\\lambda_2 t}{\\lambda_2 + \\lambda_{1,2}} + \\lambda_{1,2}\\max\\left \\{\\frac{1-t}{\\lambda_1 + \\lambda_{1,2}}, \\frac{t}{\\lambda_2 + \\lambda_{1,2}} \\right \\} +``` + +References: +* [mai2012simulating](@cite) Mai, J. F., & Scherer, M. (2012). Simulating copulas: stochastic models, sampling algorithms, and applications (Vol. 4). World Scientific. +""" +struct MOCopula{P} <: ExtremeValueCopula{P} + λ1::P + λ2::P + λ12::P + + function MOCopula(λ::Vararg{Real}) + if length(λ) !== 3 + throw(ArgumentError("MOCopula requires only 3 arguments.")) + end + + T = promote_type(typeof(λ[1]), typeof(λ[2]), typeof(λ[3])) + λ1, λ2, λ12 = T(λ[1]), T(λ[2]), T(λ[3]) + + if λ1 < 0 || λ2 < 0 || λ12 < 0 + throw(ArgumentError("All λ parameters must be >= 0")) + end + + return new{T}(λ1, λ2, λ12) + end +end + +function A(C::MOCopula, t::Real) + λ1, λ2, λ12 = C.λ1, C.λ2, C.λ12 + return ((λ1 * (1 - t))/(λ1 + λ12)) + ((λ2 * t)/(λ2 + λ12)) + λ12 * max((1-t)/(λ1 + λ12), t/(λ2 + λ12)) +end + +function ℓ(C::MOCopula, t::Vector) + λ1, λ2, λ12 = C.λ1, C.λ2, C.λ12 + t₁, t₂ = t + return ((λ1 * t₂)/(λ1 + λ12)) + ((λ2 * t₁)/(λ2 + λ12)) + λ12 * max(t₂/(λ1 + λ12), t₁/(λ2 + λ12)) +end + +function _cdf(C::MOCopula, u::AbstractArray{<:Real}) + λ1, λ2, λ12 = C.λ1, C.λ2, C.λ12 + u₁, u₂ = u + return min(u₂*u₁^(1 - λ12/(λ1 + λ12)), u₁ * u₂^(1 - λ12/(λ2 + λ12))) +end + +function Distributions._rand!(rng::Distributions.AbstractRNG, C::MOCopula, u::AbstractVector{T}) where {T<:Real} + λ1, λ2, λ12 = C.λ1, C.λ2, C.λ12 + r, s, t = rand(rng, Distributions.Uniform(0,1),3) + x = min(-log(r)/λ1, -log(t)/λ12) + y = min(-log(s)/λ2, -log(t)/λ12) + u[1] = exp(-(λ1+λ12)*x) + u[2] = exp(-(λ2+λ12)*y) + return u +end diff --git a/src/ExtremeValueCopulas/MixedCopula.jl b/src/ExtremeValueCopulas/MixedCopula.jl new file mode 100644 index 00000000..56dd3a44 --- /dev/null +++ b/src/ExtremeValueCopulas/MixedCopula.jl @@ -0,0 +1,37 @@ +""" + MixedCopula{P} + +Fields: + + - θ::Real - parameter + +Constructor + + MixedCopula(θ) + +The bivariate Mixed copula is parameterized by ``\\alpha \\in [0,1]``. It is an Extreme value copula with Pickands dependence function: + +```math +A(t) = \\theta t^2 - \\theta t + 1 +``` + +It has a few special cases: +- When θ = 0, it is the IndependentCopula + +References: +* [tawn1988bivariate](@cite) : Tawn, Jonathan A. "Bivariate extreme value theory: models and estimation." Biometrika 75.3 (1988): 397-415. +""" +struct MixedCopula{P} <: ExtremeValueCopula{P} + θ::P # Copula paremeter + function MixedCopula(θ) + if !(0 <= θ <= 1) + throw(ArgumentError("The parameter θ must be in the interval [0, 1]")) + elseif θ == 0 + return IndependentCopula(2) + else + return new{typeof(θ)}(θ) + end + end +end + +A(C::MixedCopula, t::Real) = C.θ*t^2 - C.θ*t + 1 \ No newline at end of file diff --git a/src/ExtremeValueCopulas/tEVCopula.jl b/src/ExtremeValueCopulas/tEVCopula.jl new file mode 100644 index 00000000..84b3f56e --- /dev/null +++ b/src/ExtremeValueCopulas/tEVCopula.jl @@ -0,0 +1,167 @@ +""" + tEVCopula{P} + +Fields: + - ν::Real - paremeter + - θ::Real - Parameter + +Constructor + + tEVCopula(ν, θ) + +The bivariate extreme t copula is parameterized by ``\\nu \\in [0,\\infty)`` and \\theta \\in (-1,1]. It is an Extreme value copula with Pickands dependence function: + +```math +A(x) = xt_{\\nu+1}(Z_x) +(1-x)t_{\\nu+1}(Z_{1-x}) +``` +Where ``t_{\\nu + 1}``is the cumulative distribution function (CDF) of the standard t distribution with \\nu + 1 degrees of freedom and + +```math +Z_x = \\frac{(1+\\nu)^{1/2}{\\sqrt{1-\\theta^2}}\\left [ \\left (\\frac{x}{1-x} \\right )^{1/\\nu} - \\theta \\right ] +``` + +It has a few special cases: + +- When θ = 0, it is the Independent Copula +- When θ = ∞, it is the M Copula (Upper Frechet-Hoeffding bound) + +References: +* [nikoloulopoulos2009extreme](@cite) Nikoloulopoulos, A. K., Joe, H., & Li, H. (2009). Extreme value properties of multivariate t copulas. Extremes, 12, 129-148. +""" +struct tEVCopula{df, P} <: ExtremeValueCopula{P} + ρ::P # correlation paremeter + ν::df # degree of freedom + + function tEVCopula(ν::df, ρ::P) where {df<:Real, P<:Real} + if ν <= 0 + throw(ArgumentError("The degrees of freedom ν must be positive real")) + end + if !(-1 < ρ <= 1) + throw(ArgumentError("The correlation parameter ρ must be in (-1, 1]")) + elseif ρ == 0 + return IndependentCopula(2) + elseif ρ == 1 + return MCopula(2) + end + return new{df, typeof(ρ)}(ρ, ν) + end +end +# # specific ℓ funcion of Extreme t Copula +function ℓ(T::tEVCopula{P}, t::Vector) where P + ρ = T.ρ + ν = T.ν + t₁, t₂ = t + b = sqrt(ν + 1) / sqrt(1 - ρ^2) + term1 = t₁ * StatsFuns.tdistcdf(ν + 1, b * ((t₁ / t₂)^(1 / ν) - ρ)) + term2 = t₂ * StatsFuns.tdistcdf(ν + 1, b * ((t₂ / t₁)^(1 / ν) - ρ)) + return term1 + term2 +end +function z(T::tEVCopula, t) + ρ = T.ρ + ν = T.ν + return ((1+ν)^(1/2))*((t/(1-t))^(1/ν) - ρ)*(1-ρ^2)^(-1/2) +end +# specific A funcion of Extreme t Copula +function A(T::tEVCopula, t::Real) + ρ = T.ρ + ν = T.ν + t = clamp(t, 0, 1) + zt = z(T,t) + tt_minus = z(T,1-t) + term1 = t * StatsFuns.tdistcdf(ν + 1, zt) + term2 = (1-t) * StatsFuns.tdistcdf(ν + 1, tt_minus) + return term1 + term2 +end + +function dA(C::tEVCopula, t::Real) + h = 1e-5 + t_h_clamped = clamp(t - h, 0, 1) + t_h_clamped_plus = clamp(t + h, 0, 1) + dA_minus = A(C, t_h_clamped) + dA_plus = A(C, t_h_clamped_plus) + da = (dA_plus - dA_minus) / (2 * h) + return da +end + +# Approximation of the second derivative of A because the t distribution is not compatible with ForwarDiff +function d²A(C::tEVCopula, t::Real) + h = 1e-5 + t_h_clamped = clamp(t - h, 0, 1) + t_h_clamped_plus = clamp(t + h, 0, 1) + dA_minus = dA(C, t_h_clamped) + dA_plus = dA(C, t_h_clamped_plus) + d2A = (dA_plus - dA_minus) / (2 * h) + return d2A +end + +# PDF function for ExtremeValueCopula using ℓ +function _pdf(C::tEVCopula, u::AbstractArray{<:Real}) + t = -log.(u) + c = exp(-ℓ(C, t)) + D1 = D_B_ℓ(C, t, [1]) + D2 = D_B_ℓ(C, t, [2]) + D12 = D_B_ℓ(C, t, [1, 2]) + return c * (-D12 + D1 * D2) / (u[1] * u[2]) +end +function Distributions._logpdf(C::tEVCopula, u::AbstractArray{<:Real}) + t = -log.(u) + c = exp(-ℓ(C, t)) + D1 = D_B_ℓ(C, t, [1]) + D2 = D_B_ℓ(C, t, [2]) + D12 = D_B_ℓ(C, t, [1, 2]) + return log(c) + log(-D12 + D1 * D2) - log(u[1] * u[2]) +end + +function D_B_ℓ(C::tEVCopula, t::Vector{Float64}, B::Vector{Int}) + h = 1e-5 + if length(B) == 1 + # First partial derivative + return partial_derivative_1(C, t, B[1], h) + elseif length(B) == 2 + # Second partial derivative or mixed derivative + return partial_derivative_2(C, t, B[1], B[2], h) + else + throw(ArgumentError("Higher order partial derivatives are not required for bivariate case")) + end +end + +function partial_derivative_1(C::tEVCopula, t::Vector{Float64}, i::Int, h::Float64) + t_plus = copy(t) + t_minus = copy(t) + t_plus[i] += h + t_minus[i] -= h + + return (ℓ(C, t_plus) - ℓ(C, t_minus)) / (2 * h) +end + +function partial_derivative_2(C::tEVCopula, t::Vector{Float64}, i::Int, j::Int, h::Float64) + if i == j + # Second partial derivative + t_plus = copy(t) + t_minus = copy(t) + t_plus[i] += h + t_minus[i] -= h + + d_plus = partial_derivative_1(C, t_plus, i, h) + d_minus = partial_derivative_1(C, t_minus, i, h) + + return (d_plus - d_minus) / (2 * h) + else + # Mixed derivative + t_plus_plus = copy(t) + t_plus_minus = copy(t) + t_minus_plus = copy(t) + t_minus_minus = copy(t) + + t_plus_plus[i] += h + t_plus_plus[j] += h + t_plus_minus[i] += h + t_plus_minus[j] -= h + t_minus_plus[i] -= h + t_minus_plus[j] += h + t_minus_minus[i] -= h + t_minus_minus[j] -= h + + return (ℓ(C, t_plus_plus) - ℓ(C, t_plus_minus) - ℓ(C, t_minus_plus) + ℓ(C, t_minus_minus)) / (4 * h^2) + end +end \ No newline at end of file diff --git a/src/UnivariateDistribution/ExtremeDist.jl b/src/UnivariateDistribution/ExtremeDist.jl new file mode 100644 index 00000000..8e272482 --- /dev/null +++ b/src/UnivariateDistribution/ExtremeDist.jl @@ -0,0 +1,43 @@ +struct ExtremeDist{C} <: Distributions.ContinuousUnivariateDistribution + G::C + function ExtremeDist(G) + return new{typeof(G)}(G) + end +end + +function Distributions.cdf(d::ExtremeDist, z) + if z < 0 + return 0.0 + elseif z > 1 + return 1.0 + else + copula = d.G + return z + z * (1 - z) * (dA(copula, z) / A(copula, z)) + end +end + +function _pdf(d::ExtremeDist, z) + if z < 0 || z > 1 + return 0.0 + else + copula = d.G + _A = A(copula, z) + A_prime = dA(copula, z) + A_double_prime = d²A(copula, z) + return 1 + (1 - 2z) * A_prime / _A + z * (1 - z) * (A_double_prime * _A - A_prime^2) / _A^2 + end +end + +function Distributions.quantile(d::ExtremeDist, p) + if p < 0 || p > 1 + throw(ArgumentError("p must be between 0 and 1")) + end + cdf_func(x) = Distributions.cdf(d, x) - p + return Roots.find_zero(cdf_func, (eps(), 1-eps()), Roots.Brent()) +end + +# Generate random samples from the radial distribution using the quantile function +function Distributions.rand(rng::Distributions.AbstractRNG, d::ExtremeDist) + u = rand(rng, Distributions.Uniform(0,1)) + return Distributions.quantile(d, u) +end \ No newline at end of file diff --git a/test/Extreme_value_test.jl b/test/Extreme_value_test.jl new file mode 100644 index 00000000..561e9702 --- /dev/null +++ b/test/Extreme_value_test.jl @@ -0,0 +1,326 @@ +@testitem "Extreme Value Copulas Tests" begin + using InteractiveUtils + using Copulas, Distributions + using Random + using StableRNGs + + rng = StableRNG(123) + d = 100 + @testset "AsymGalambosCopula - sampling, pdf, cdf" begin + for α in [0.0, rand(rng, Uniform()), rand(rng, Uniform(5.0, 9.0)), rand(rng, Uniform(10.0, 15.0))] + for θ in [[rand(rng, Uniform(0, 1)), rand(rng, Uniform(0, 1))], [0.0, 0.0], [1.0, 1.0]] + C = AsymGalambosCopula(α, θ) + data = rand(rng, C, 100) + + @test size(data) == (2, 100) + + for i in 1:d + u = data[:,i] + cdf_value = cdf(C, u) + pdf_value = pdf(C, u) + @test 0.0 <= cdf_value <= 1.0 + @test pdf_value >= 0.0 + end + end + end + end + + @testset "AsymLogCopula - sampling, pdf, cdf" begin + for α in [1.0, rand(rng, Uniform(1.0, 5.0)), rand(rng, Uniform(10.0, 15.0))] + for θ in [[rand(rng, Uniform(0, 1)), rand(rng, Uniform(0, 1))], [0.0, 0.0], [1.0, 1.0]] + C = AsymLogCopula(α, θ) + data = rand(rng, C, 100) + + @test size(data) == (2, 100) + + for i in 1:d + u = data[:,i] + cdf_value = cdf(C, u) + pdf_value = pdf(C, u) + @test 0.0 <= cdf_value <= 1.0 + @test pdf_value >= 0.0 + end + end + end + end + + @testset "AsymMixedCopula - sampling, pdf, cdf" begin + for θ in [[0.1, 0.2], [0.0, 0.0], [0.3, 0.4], [0.2, 0.4]] + try + C = AsymMixedCopula(θ) + data = rand(rng, C, 100) + + @test size(data) == (2, 100) + + for i in 1:d + u = data[:,i] + cdf_value = cdf(C, u) + pdf_value = pdf(C, u) + + @test 0.0 <= cdf_value <= 1.0 || error("CDF failure: θ=$θ, u=$u, cdf_value=$cdf_value") + @test pdf_value >= 0.0 || error("PDF failure: θ=$θ, u=$u, cdf_value=$cdf_value") + end + catch e + @test e isa ArgumentError + println("Could not construct AsymMixedCopula with θ=$θ: ", e) + end + end + end + + @testset "BC2Copula - sampling, pdf, cdf" begin + for θ in [[rand(rng), rand(rng)], [1.0, 0.0], [0.5, 0.5], [-0.1, 0.2], [1.1, 0.5], [0.5, -0.2]] + try + C = BC2Copula(θ...) + data = rand(rng, C, 100) + + @test size(data) == (2, 100) + + for i in 1:d + u = data[:,i] + cdf_value = cdf(C, u) + pdf_value = pdf(C, u) + + @test 0.0 <= cdf_value <= 1.0 || error("CDF failure: θ=$θ, u=$u, cdf_value=$cdf_value") + @test pdf_value >= 0.0 || error("PDF failure: θ=$θ, u=$u, cdf_value==$pdf_value") + end + catch e + @test e isa ArgumentError + println("Could not build BC2Copula with θ=$θ: ", e) + end + end + end + + @testset "CuadrasAugeCopula - sampling, pdf, cdf" begin + for θ in [rand(rng), 0.0, 1.0, -0.1, 1.1, rand(rng)] + try + C = CuadrasAugeCopula(θ) + data = rand(rng, C, 100) + + @test size(data) == (2, 100) + + for i in 1:d + u = data[:,i] + cdf_value = cdf(C, u) + #if θ == 1.0 + # @test_throws ArgumentError pdf(C, u) + #else + pdf_value = pdf(C, u) + @test 0.0 <= cdf_value <= 1.0 || error("CDF failure: θ=$θ, u=$u, cdf_value=$cdf_value") + @test pdf_value >= 0.0 || error("PDF failure: θ=$θ, u=$u, cdf_value=$pdf_value") + #end + end + catch e + @test e isa ArgumentError + println("Could not construct CuadrasAugeCopula with θ=$θ: ", e) + end + end + end + + @testset "GalambosCopula - sampling, pdf, cdf" begin + for θ in [rand(rng), 0.0, Inf, -1.0, rand(rng, Uniform(1.0, 5.0)), rand(rng, Uniform(5.0, 10.0))] + try + C = GalambosCopula(θ) + data = rand(rng, C, 100) + + @test size(data) == (2, 100) + + for i in 1:d + u = data[:,i] + cdf_value = cdf(C, u) + #if θ == Inf + # @test_throws ArgumentError pdf(C, u) + #else + pdf_value = pdf(C, u) + @test 0.0 <= cdf_value <= 1.0 || error("CDF failure: θ=$θ, u=$u, cdf_value=$cdf_value") + @test pdf_value >= 0.0 || error("PDF failure: θ=$θ, u=$u, cdf_value=$pdf_value") + #end + end + catch e + @test e isa ArgumentError + println("Could not construct GalambosCopula with θ=$θ: ", e) + end + end + end + + @testset "HuslerReissCopula - sampling, pdf, cdf" begin + for θ in [rand(rng), 0.0, Inf, -1.0, rand(rng, Uniform(1.0, 5.0)), rand(rng, Uniform(5.0, 10.0))] + try + C = HuslerReissCopula(θ) + data = rand(rng, C, 100) + + @test size(data) == (2, 100) + + for i in 1:d + u = data[:,i] + cdf_value = cdf(C, u) + #if θ == Inf + # @test_throws ArgumentError pdf(C, u) + #else + pdf_value = pdf(C, u) + @test 0.0 <= cdf_value <= 1.0 || error("CDF failure: θ=$θ, u=$u, cdf_value=$cdf_value") + @test pdf_value >= 0.0 || error("PDF failure: θ=$θ, u=$u, cdf_value=$pdf_value") + #end + end + catch e + @test e isa ArgumentError + println("Could not construct HuslerReissCopula with θ=$θ: ", e) + end + end + end + + @testset "LogCopula - sampling, pdf, cdf" begin + for θ in [1.0, Inf, 0.5, rand(rng, Uniform(1.0, 10.0))] + try + C1 = LogCopula(θ) + C2 = GumbelCopula(2, θ) + data = rand(rng, C1, 100) + + @test size(data) == (2, 100) + + for i in 1:d + u = data[:,i] + cdf_value_C1 = cdf(C1, u) + cdf_value_C2 = cdf(C2, u) + #if θ == Inf + # @test_throws ArgumentError pdf(C1, u) + # @test_throws ArgumentError pdf(C2, u) + #else + pdf_value_C1 = pdf(C1, u) + pdf_value_C2 = pdf(C2, u) + @test 0.0 <= cdf_value_C1 <= 1.0 || error("CDF failure LogCopula: θ=$θ, u=$u, cdf_value_C1=$cdf_value_C1") + @test 0.0 <= cdf_value_C2 <= 1.0 || error("CDF failure de GumbelCopula: θ=$θ, u=$u, cdf_value_C2=$cdf_value_C2") + @test pdf_value_C1 >= 0.0 || error("PDF failure LogCopula: θ=$θ, u=$u, pdf_value_C1=$pdf_value_C1") + @test pdf_value_C2 >= 0.0 || error("PDF failure GumbelCopula: θ=$θ, u=$u, pdf_value_C2=$pdf_value_C2") + + @test isapprox(cdf_value_C1, cdf_value_C2, atol=1e-6) || error("CDF LogCopula and GumbelCopula do not match: θ=$θ, u=$u, cdf_value_C1=$cdf_value_C1, cdf_value_C2=$cdf_value_C2") + @test isapprox(pdf_value_C1, pdf_value_C2, atol=1e-6) || error("PDF LogCopula and GumbelCopula do not match: θ=$θ, u=$u, pdf_value_C2=$pdf_value_C1, pdf_value_C2=$pdf_value_C2") + #end + end + catch e + @test e isa ArgumentError + println("Could not construct LogCopula with θ=$θ: ", e) + end + end + end + + @testset "MixedCopula - sampling, pdf, cdf" begin + for θ in [0.0, 1.0, 0.2, -rand(rng), 0.5] + try + C = MixedCopula(θ) + data = rand(rng, C, 100) + + @test size(data) == (2, 100) + + for i in 1:d + u = data[:,i] + cdf_value = cdf(C, u) + pdf_value = pdf(C, u) + @test 0.0 <= cdf_value <= 1.0 || error("CDF failure: θ=$θ, u=$u, cdf_value=$cdf_value") + @test pdf_value >= 0.0 || error("PDF failure: θ=$θ, u=$u, cdf_value=$pdf_value") + end + catch e + @test e isa ArgumentError + println("Could not construct MixedCopula with θ=$θ: ", e) + end + end + end + + @testset "MOCopula - sampling, pdf, cdf" begin + param_sets = [ + (0.1, 0.2, 0.3), + (1.0, 1.0, 1.0), + (0.5, 0.5, 0.5), + (-0.1, 0.2, 0.3), # Invalid case: negative parameter + (0.1, -0.2, 0.3), # Invalid case: negative parameter + (0.1, 0.2, -0.3), # Invalid case: negative parameter + (rand(rng), rand(rng), rand(rng)) # Random Params + ] + + for λ in param_sets + try + C = MOCopula(λ...) + data = rand(rng, C, 100) + + @test size(data) == (2, 100) + + for i in 1:d + u = data[:,i] + cdf_value = cdf(C, u) + pdf_value = pdf(C, u) + + @test 0.0 <= cdf_value <= 1.0 || error("cdf failure: λ=$λ, u=$u, cdf_value=$cdf_value") + @test pdf_value >= 0.0 || error("pdf failure: λ=$λ, u=$u, pdf_value=$pdf_value") + end + catch e + @test e isa ArgumentError + println("Could not construct MOCopula with λ=$λ: ", e) + end + end + end + + @testset "tEVCopula - sampling, pdf, cdf" begin + param_sets = [ + (2.0, 0.5), # ν > 0 y -1 < ρ <= 1 + (5.0, -0.5), # ν > 0 y -1 < ρ <= 1 + (10.0, 1.0), # ν > 0 y ρ == 1, MCopula + (3.0, 0.0), # ν > 0 y ρ == 0, Independent Copula + (-2.0, 0.5), # Invalid case: ν <= 0 + (3.0, 1.1), # Invalid case: ρ out of range + (3.0, -1.1), # Invalid case: ρ out of range + (rand(rng, Uniform(4.0, 10.0)), rand(rng, Uniform(-0.9, 1.0))) # random params inside of range + ] + + for (ν, ρ) in param_sets + try + C = tEVCopula(ν, ρ) + data = rand(rng, C, 100) + + @test size(data) == (2, 100) + + for i in 1:d + u = data[:,i] + + cdf_value = cdf(C, u) + #if ρ == 1 + # @test_throws ArgumentError pdf(C, u) + #else + pdf_value = pdf(C, u) + @test 0.0 <= cdf_value <= 1.0 || error("CDF Failure: ν=$ν, ρ=$ρ, u=$u, cdf_value=$cdf_value") + @test pdf_value >= 0.0 || error("PDF Failure: ν=$ν, ρ=$ρ, u=$u, pdf_value=$pdf_value") + #end + end + catch e + @test e isa ArgumentError + println("Could not construct tEVCopula with ν=$ν, ρ=$ρ: ", e) + end + end + end + + @testset "Pickands Function test" begin + cops = ( + AsymGalambosCopula(5.0, [0.8, 0.3]), + AsymLogCopula(1.5, [0.5, 0.2]), + AsymMixedCopula([0.1, 0.2]), + BC2Copula(0.5, 0.3), + CuadrasAugeCopula(0.8), + GalambosCopula(4.3), + HuslerReissCopula(3.5), + LogCopula(5.5), + MixedCopula(0.5), + MOCopula(0.1, 0.2, 0.3), + tEVCopula(4.0, 0.5)) + + for C in cops + for C in cops + @test Copulas.A(C, 0.0) == 1 + @test Copulas.A(C, 1.0) == 1 + t = rand() + A_value = Copulas.A(C, t) + @test 0.0 <= A_value <= 1.0 + @test isapprox(A_value, max(t, 1-t); atol=1e-6) || A_value >= max(t, 1-t) + @test A_value <= 1.0 + end + end + + end +end \ No newline at end of file diff --git a/test/margins_uniformity.jl b/test/margins_uniformity.jl index 6b36f710..e2634325 100644 --- a/test/margins_uniformity.jl +++ b/test/margins_uniformity.jl @@ -29,15 +29,38 @@ WCopula(2), ArchimedeanCopula(2,i𝒲(LogNormal(),2)), PlackettCopula(2.0), - EmpiricalCopula(randn(2,100),pseudo_values=false), + EmpiricalCopula(randn(2,2000),pseudo_values=false), SurvivalCopula(ClaytonCopula(2,-0.7),(1,2)), RafteryCopula(2, 0.2), RafteryCopula(3, 0.5), + AsymGalambosCopula(5.0, [0.8, 0.3]), + AsymLogCopula(1.5, [0.5, 0.2]), + AsymMixedCopula([0.1, 0.2]), + BC2Copula(0.5, 0.3), + CuadrasAugeCopula(0.8), + GalambosCopula(4.3), + HuslerReissCopula(3.5), + LogCopula(5.5), + MixedCopula(0.5), + MOCopula(0.1, 0.2, 0.3), + tEVCopula(4.0, 0.5), Copulas.SubsetCopula(RafteryCopula(3, 0.5), (2,1)), # Others ? Yes probably others too ! - ) - + # I added your copulas here + # You need to provide at least one example per type, but maybe you want to provide more ? Especially if some parameter changes yield very different behaviors. + AsymGalambosCopula(0.1, [0.2,0.6]), + AsymLogCopula(1.2, [0.3,0.6]), + AsymMixedCopula([0.1,0.2]), + BC2Copula(0.7,0.3), + CuadrasAugeCopula(0.1), + GalambosCopula(0.3), + HuslerReissCopula(0.1), + LogCopula(1.5), + MixedCopula(0.2), + MOCopula(0.1,0.2,0.3), + tEVCopula(2.0,0.5), + ) #### Try to ensure that every copula in the package is indeed in this list, to remmember contributors to add their model here: function _subtypes(type::Type) @@ -55,7 +78,7 @@ for CT in _subtypes(Copulas.Copula) # Check that every copula type has been used @test any(isa(C,CT) for C in cops) end - for TG in _subtypes(Copulas.Generator) # Check that every generator has been used + for TG in _subtypes(Copulas.Generator) # Check that every archimedean generator has been used @test any(isa(C.G,TG) for C in cops if typeof(C)<:Copulas.ArchimedeanCopula) end @@ -90,37 +113,35 @@ end return false end - n = 1000 - U = Uniform(0,1) for C in cops d = length(C) CT = typeof(C) rng = StableRNG(123) spl = rand(rng,C,n) - if !(CT<:TCopula) - # Check that the cdf has special values at the bounds: - @test cdf(C,zeros(d)) == 0 - @test cdf(C,ones(d)) == 1 - - # Check that the cdf values are in [0,1] - @test all(0 .<= cdf(C,spl) .<= 1) - end - # Check that samples are in [0,1]: - @test all(0 <= x <= 1 for x in spl) + @test iszero(cdf(C,zeros(d))) + @test isone(cdf(C,ones(d))) + @test 0 <= cdf(C,rand(rng,d)) <= 1 + @test all(0 .<= spl .<= 1) # Check uniformity of each marginal : if !(CT<:EmpiricalCopula) # this one is not a true copula :) for i in 1:d - # On the samples - @test pvalue(ApproximateOneSampleKSTest(spl[i,:], U),tail=:right) > 0.009 # this is weak but enough to catch mistakes. - - # On the cdf: - u = ones(d) - for val in [0,1,rand(rng,10)...] + # Check Samples uniformity + # this is weak but enough to catch impementation mistakes. + + # ks_test = pvalue(ExactOneSampleKSTest(spl[i,:], Uniform())) > 1e-5 + # @test ks_test + + # Check CDF marginals uniformity. + for val in rand(rng,5) + u = ones(d) u[i] = val + if !(cdf(C,u) ≈ val) + @show C, u + end @test cdf(C,u) ≈ val atol=1e-5 end # extra check for zeros: @@ -130,32 +151,7 @@ end end - # Conditionally on the applicability of the pdf method... - # Finally we do not check pdf, as it is too broken in a lot of cases... - - - # Something should be made to revamp this test - # if applicable(pdf,C,spl) - - # # if archimedean, check also that monotonicity is good: - # if !(CT<:ArchimedeanCopula) || ((Copulas.max_monotony(C.G) > d) && !(typeof(Copulas.williamson_dist(C.G,d))<:WilliamsonTransforms.𝒲₋₁)) - - # # check that pdf values are positives: - # @test all(pdf(C,spl) .>= 0) - - # # also check that pdf values are indeed derivatives of the cdf values: - # begin - # for _ in 1:10 - # u = rand(rng,d) - # @test isapprox(get_numerical_pdf(C,u),pdf(C,u),atol=1e-5) - # end - # end - # end - # end - - - # only check archimedeans for tau ∘ tau_inv - + # Extra checks, only for archimedeans. if is_archimedean_with_agenerator(CT) if applicable(Copulas.τ,C.G) @@ -180,16 +176,7 @@ @test Copulas.ρ(Copulas.generatorof(CT)(Copulas.ρ⁻¹(CT,rho))) ≈ rho end end - fit(CT,spl) - end - - # Check that fitting works: - # if additional_condition(CT) - # fit(CT,spl) - # end - # @test true - end end