copula_core/archimedean/
amh.rs

1//! Ali-Mikhail-Haq (AMH) copula implementation.
2//!
3//! ## Bibliography
4//! - Ali, M. M., Mikhail, N. N., & Haq, M. S. (1978). A class of bivariate
5//!   distribution functions. *Journal of Multivariate Analysis*, 8(3), 405-412.
6//! - Nelsen, R. B. (2006). *An Introduction to Copulas*. Springer.
7//! - Joe, H. (2014). *Dependence Modeling with Copulas*. CRC Press.
8
9use crate::{ArchimedeanCopula, Copula, CopulaError, Result};
10use nalgebra::DMatrix;
11use rand::{Rng, RngExt};
12
13/// Ali-Mikhail-Haq copula with parameter `theta` in (-1, 1).
14#[derive(Debug, Clone)]
15pub struct AMHCopula {
16    /// Copula parameter θ ∈ (-1, 1)
17    theta: f64,
18}
19
20validated_serde!("AMHCopula", AMHCopula { theta: f64 } => AMHCopula::new(theta));
21
22impl AMHCopula {
23    /// Create a new AMH copula with parameter `theta`.
24    pub fn new(theta: f64) -> Result<Self> {
25        if !theta.is_finite() || theta.abs() >= 1.0 {
26            return Err(CopulaError::invalid_parameter(
27                "theta must be finite and in (-1, 1)",
28            ));
29        }
30        Ok(Self { theta })
31    }
32}
33
34impl Copula for AMHCopula {
35    fn cdf(&self, u: &[f64]) -> Result<f64> {
36        if u.len() != 2 {
37            return Err(CopulaError::dimension_mismatch(2, u.len()));
38        }
39        crate::error::validate_unit_range(u)?;
40        let denom = 1.0 - self.theta * (1.0 - u[0]) * (1.0 - u[1]);
41        Ok(crate::utils::clamp_to_frechet_bounds(
42            u,
43            u[0] * u[1] / denom,
44        ))
45    }
46
47    fn pdf(&self, u: &[f64]) -> Result<f64> {
48        if u.len() != 2 {
49            return Err(CopulaError::dimension_mismatch(2, u.len()));
50        }
51        crate::error::validate_unit_range(u)?;
52
53        let theta = self.theta;
54        let u1_bar = 1.0 - u[0];
55        let u2_bar = 1.0 - u[1];
56        let denom = 1.0 - theta * u1_bar * u2_bar;
57
58        // AMH copula PDF: c(u,v) = [1 - θ + 2θ(1-u)(1-v)] / [1 - θ(1-u)(1-v)]^2
59        let numerator = 1.0 - theta + 2.0 * theta * u1_bar * u2_bar;
60        let denominator = denom.powi(2);
61
62        Ok(numerator / denominator)
63    }
64
65    fn sample<R: Rng + ?Sized>(&self, n: usize, rng: &mut R) -> Result<DMatrix<f64>> {
66        let mut samples = DMatrix::<f64>::zeros(n, 2);
67
68        for i in 0..n {
69            let u1: f64 = rng.random::<f64>();
70            let v: f64 = rng.random::<f64>();
71
72            // Binary search for u2 using conditional CDF
73            let mut u2_low: f64 = 1e-10;
74            let mut u2_high: f64 = 1.0 - 1e-10;
75            let mut u2: f64 = 0.5;
76
77            for _ in 0..50 {
78                u2 = (u2_low + u2_high) / 2.0;
79
80                let denom = 1.0 - self.theta * (1.0 - u1) * (1.0 - u2);
81                let denom2 = denom.powi(2);
82
83                // Conditional CDF: ∂C/∂u1
84                let cond_cdf = u2 * (1.0 - self.theta * (1.0 - u2)) / denom2;
85
86                if (cond_cdf - v).abs() < 1e-10 {
87                    break;
88                }
89
90                if cond_cdf < v {
91                    u2_low = u2;
92                } else {
93                    u2_high = u2;
94                }
95            }
96
97            samples[(i, 0)] = u1;
98            samples[(i, 1)] = u2;
99        }
100
101        Ok(samples)
102    }
103
104    fn dimension(&self) -> usize {
105        2
106    }
107}
108
109impl ArchimedeanCopula for AMHCopula {
110    fn phi(&self, t: f64) -> Result<f64> {
111        if t <= 0.0 || t > 1.0 {
112            return Err(CopulaError::invalid_range(vec![t]));
113        }
114        // φ(t) = ln[(1-θ(1-t))/t]
115        let numerator = 1.0 - self.theta * (1.0 - t);
116        if numerator <= 0.0 || t <= 0.0 {
117            return Err(CopulaError::numerical("phi argument out of valid range"));
118        }
119        Ok((numerator / t).ln())
120    }
121
122    fn phi_inv(&self, s: f64) -> Result<f64> {
123        if s < 0.0 {
124            return Err(CopulaError::invalid_range(vec![s]));
125        }
126        // φ^(-1)(s) = (1-θ)/[exp(s) - θ]
127        let exp_s = s.exp();
128        let denom = exp_s - self.theta;
129        if denom.abs() < 1e-15 {
130            return Err(CopulaError::numerical("phi_inv denominator too small"));
131        }
132        Ok((1.0 - self.theta) / denom)
133    }
134
135    fn phi_inv_deriv(&self, s: f64, k: usize) -> Result<f64> {
136        if s < 0.0 {
137            return Err(CopulaError::invalid_range(vec![s]));
138        }
139
140        let exp_s = s.exp();
141        let denom = exp_s - self.theta;
142
143        match k {
144            1 => {
145                // First derivative: -(1-θ)exp(s) / [exp(s) - θ]^2
146                Ok(-(1.0 - self.theta) * exp_s / denom.powi(2))
147            }
148            2 => {
149                // Second derivative
150                let num = (1.0 - self.theta) * exp_s * (2.0 * exp_s - self.theta);
151                Ok(num / denom.powi(3))
152            }
153            _ => Err(CopulaError::not_implemented("phi_inv_deriv k>2")),
154        }
155    }
156}
157
158#[cfg(test)]
159mod tests {
160    use super::*;
161
162    #[test]
163    fn new_rejects_out_of_bounds_theta() {
164        assert!(AMHCopula::new(1.0).is_err());
165        assert!(AMHCopula::new(-1.0).is_err());
166    }
167
168    #[test]
169    fn valid_new_returns_copula() {
170        let cop = AMHCopula::new(0.5).unwrap();
171        assert_eq!(cop.dimension(), 2);
172    }
173
174    #[test]
175    fn cdf_matches_formula() {
176        let cop = AMHCopula::new(0.2).unwrap();
177        let cdf = cop.cdf(&[0.5, 0.5]).unwrap();
178        let expected = 0.5 * 0.5 / (1.0 - 0.2 * 0.5 * 0.5);
179        assert!((cdf - expected).abs() < 1e-12);
180    }
181}
copula_core/archimedean/amh.rs

copula_core/archimedean/
amh.rs
