copula_core/archimedean/
clayton.rs

1//! Clayton copula implementation.
2//!
3//! ## Bibliography
4//! - Clayton, D. G. (1978). A model for association in bivariate life tables
5//!   and its application in epidemiological studies of familial tendency in
6//!   chronic disease incidence. *Biometrika*, 65(1), 141-151.
7//! - Nelsen, R. B. (2006). *An Introduction to Copulas*. Springer.
8//! - Joe, H. (2014). *Dependence Modeling with Copulas*. CRC Press.
9
10use crate::traits::BoundedParameters;
11#[cfg(feature = "estimation")]
12use crate::traits::FittableCopula;
13#[cfg(feature = "estimation")]
14use crate::utils::kendall_tau;
15use crate::{ArchimedeanCopula, Copula, CopulaError, Result};
16use nalgebra::DMatrix;
17use rand::Rng;
18#[cfg(feature = "estimation")]
19use statrs::distribution::ContinuousCDF;
20
21/// Clayton copula with positive parameter `theta`.
22#[derive(Debug, Clone)]
23pub struct ClaytonCopula {
24    /// Copula parameter θ > 0
25    theta: f64,
26}
27
28validated_serde!("ClaytonCopula", ClaytonCopula { theta: f64 } => ClaytonCopula::new(theta));
29
30impl ClaytonCopula {
31    /// Create a new Clayton copula with parameter `theta`.
32    pub fn new(theta: f64) -> Result<Self> {
33        if theta <= 0.0 || !theta.is_finite() {
34            return Err(CopulaError::invalid_parameter("theta must be positive"));
35        }
36        Ok(Self { theta })
37    }
38}
39
40impl BoundedParameters for ClaytonCopula {
41    fn parameter_bounds() -> Vec<(f64, f64)> {
42        vec![(1e-4, 10.0)]
43    }
44
45    fn check_bounds(&self) -> Result<()> {
46        let (min, max) = Self::parameter_bounds()[0];
47        if self.theta < min || self.theta > max {
48            Err(CopulaError::invalid_parameter("theta out of bounds"))
49        } else {
50            Ok(())
51        }
52    }
53}
54
55impl Copula for ClaytonCopula {
56    fn cdf(&self, u: &[f64]) -> Result<f64> {
57        if u.len() != 2 {
58            return Err(CopulaError::dimension_mismatch(2, u.len()));
59        }
60        crate::error::validate_unit_range(u)?;
61        let sum = u[0].powf(-self.theta) + u[1].powf(-self.theta) - 1.0;
62        if sum <= 0.0 {
63            Ok(0.0)
64        } else {
65            Ok(crate::utils::clamp_to_frechet_bounds(
66                u,
67                sum.powf(-1.0 / self.theta),
68            ))
69        }
70    }
71
72    fn pdf(&self, u: &[f64]) -> Result<f64> {
73        if u.len() != 2 {
74            return Err(CopulaError::dimension_mismatch(2, u.len()));
75        }
76        crate::error::validate_unit_range(u)?;
77        let (u1, u2) = (u[0], u[1]);
78        let theta = self.theta;
79        let sum = u1.powf(-theta) + u2.powf(-theta) - 1.0;
80        if sum <= 0.0 {
81            return Ok(0.0);
82        }
83        let base = sum.powf(-2.0 - 1.0 / theta);
84        let density = (1.0 + theta) * u1.powf(-1.0 - theta) * u2.powf(-1.0 - theta) * base;
85        Ok(density)
86    }
87
88    fn sample<R: Rng + ?Sized>(&self, n: usize, rng: &mut R) -> Result<DMatrix<f64>> {
89        use rand_distr::{Distribution, Exp1, Gamma};
90
91        let gamma = Gamma::<f64>::new(1.0 / self.theta, 1.0)
92            .map_err(|e| CopulaError::invalid_parameter(format!("gamma distribution: {}", e)))?;
93        let mut samples = DMatrix::<f64>::zeros(n, 2);
94        for i in 0..n {
95            let w: f64 = gamma.sample(rng);
96            for j in 0..2 {
97                let e: f64 = Exp1.sample(rng);
98                samples[(i, j)] = f64::powf(1.0 + e / w, -1.0 / self.theta);
99            }
100        }
101        Ok(samples)
102    }
103
104    fn dimension(&self) -> usize {
105        2
106    }
107
108    fn tail_dependence(&self) -> Result<(f64, f64)> {
109        let lambda_l = 2f64.powf(-1.0 / self.theta);
110        Ok((lambda_l, 0.0))
111    }
112}
113
114impl ArchimedeanCopula for ClaytonCopula {
115    fn phi(&self, t: f64) -> Result<f64> {
116        if t <= 0.0 || t > 1.0 {
117            return Err(CopulaError::invalid_range(vec![t]));
118        }
119        Ok((t.powf(-self.theta) - 1.0) / self.theta)
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        Ok((1.0 + self.theta * s).powf(-1.0 / self.theta))
127    }
128
129    fn phi_inv_deriv(&self, s: f64, k: usize) -> Result<f64> {
130        if s < 0.0 {
131            return Err(CopulaError::invalid_range(vec![s]));
132        }
133        let base = 1.0 + self.theta * s;
134        let pow = -1.0 / self.theta;
135        match k {
136            1 => Ok(-base.powf(pow - 1.0)),
137            2 => Ok((1.0 + self.theta) * base.powf(pow - 2.0)),
138            _ => Err(CopulaError::not_implemented("phi_inv_deriv k>2")),
139        }
140    }
141}
142
143#[cfg(feature = "estimation")]
144impl FittableCopula for ClaytonCopula {
145    type Parameters = f64;
146
147    fn fit(&mut self, pseudo_obs: &DMatrix<f64>) -> Result<Self::Parameters> {
148        use argmin::core::{CostFunction, Error as ArgminError, Executor, State};
149        use argmin::solver::brent::BrentOpt;
150
151        if pseudo_obs.ncols() != 2 {
152            return Err(CopulaError::dimension_mismatch(2, pseudo_obs.ncols()));
153        }
154        crate::utils::validate_pseudo_observations(pseudo_obs)?;
155
156        struct Nll<'a> {
157            data: &'a DMatrix<f64>,
158        }
159
160        impl<'a> CostFunction for Nll<'a> {
161            type Param = f64;
162            type Output = f64;
163
164            fn cost(&self, theta: &Self::Param) -> std::result::Result<f64, ArgminError> {
165                if *theta <= 0.0 {
166                    return Ok(f64::INFINITY);
167                }
168                let mut nll = 0.0;
169                for i in 0..self.data.nrows() {
170                    let u1 = self.data[(i, 0)];
171                    let u2 = self.data[(i, 1)];
172                    let sum = u1.powf(-*theta) + u2.powf(-*theta) - 1.0;
173                    if sum <= 0.0 {
174                        return Ok(f64::INFINITY);
175                    }
176                    let lp = (1.0 + theta).ln()
177                        + (-1.0 - theta) * (u1.ln() + u2.ln())
178                        + (-2.0 - 1.0 / theta) * sum.ln();
179                    nll -= lp;
180                }
181                Ok(nll)
182            }
183        }
184
185        let op = Nll { data: pseudo_obs };
186        let (min, max) = Self::parameter_bounds()[0];
187        let solver = BrentOpt::new(min, max);
188        let res = Executor::new(op, solver)
189            .configure(|state| state.max_iters(100))
190            .run()
191            .map_err(|e| CopulaError::optimization(e.to_string()))?;
192        let theta = *res
193            .state()
194            .get_best_param()
195            .ok_or_else(|| CopulaError::optimization("optimizer returned no best parameter"))?;
196        self.theta = theta;
197        Ok(theta)
198    }
199
200    fn log_likelihood(&self, pseudo_obs: &DMatrix<f64>) -> Result<f64> {
201        if pseudo_obs.ncols() != 2 {
202            return Err(CopulaError::dimension_mismatch(2, pseudo_obs.ncols()));
203        }
204        crate::utils::validate_pseudo_observations(pseudo_obs)?;
205        let mut ll = 0.0;
206        for i in 0..pseudo_obs.nrows() {
207            let u1 = pseudo_obs[(i, 0)];
208            let u2 = pseudo_obs[(i, 1)];
209            let sum = u1.powf(-self.theta) + u2.powf(-self.theta) - 1.0;
210            if sum <= 0.0 {
211                return Err(CopulaError::numerical("log_likelihood invalid sum"));
212            }
213            ll += (1.0 + self.theta).ln()
214                + (-1.0 - self.theta) * (u1.ln() + u2.ln())
215                + (-2.0 - 1.0 / self.theta) * sum.ln();
216        }
217        Ok(ll)
218    }
219
220    fn fit_moments(&mut self, pseudo_obs: &DMatrix<f64>) -> Result<Self::Parameters> {
221        if pseudo_obs.ncols() != 2 {
222            return Err(CopulaError::dimension_mismatch(2, pseudo_obs.ncols()));
223        }
224        for j in 0..2 {
225            crate::error::validate_unit_range(pseudo_obs.column(j).as_slice())?;
226        }
227        let u: Vec<f64> = pseudo_obs.column(0).iter().copied().collect();
228        let v: Vec<f64> = pseudo_obs.column(1).iter().copied().collect();
229        let tau = kendall_tau(&u, &v)?;
230        if (1.0 - tau).abs() < f64::EPSILON {
231            return Err(CopulaError::invalid_parameter("tau must be < 1"));
232        }
233        let theta = 2.0 * tau / (1.0 - tau);
234        self.theta = theta;
235        Ok(theta)
236    }
237
238    fn parameters(&self) -> Self::Parameters {
239        self.theta
240    }
241
242    fn set_parameters(&mut self, params: Self::Parameters) -> Result<()> {
243        *self = Self::new(params)?;
244        Ok(())
245    }
246
247    fn standard_errors(&self, pseudo_obs: &DMatrix<f64>) -> Result<Self::Parameters> {
248        if pseudo_obs.ncols() != 2 {
249            return Err(CopulaError::dimension_mismatch(2, pseudo_obs.ncols()));
250        }
251        crate::utils::validate_pseudo_observations(pseudo_obs)?;
252
253        let h = 1e-5;
254        let theta = self.theta;
255        let ll = |t: f64| {
256            let cop = Self { theta: t };
257            cop.log_likelihood(pseudo_obs)
258        };
259        let ll_p = ll(theta + h)?;
260        let ll_m = ll(theta - h)?;
261        let ll_0 = ll(theta)?;
262        let second = (ll_p - 2.0 * ll_0 + ll_m) / (h * h);
263        if second >= 0.0 || !second.is_finite() {
264            return Err(CopulaError::numerical("invalid hessian"));
265        }
266        let var = -1.0 / second;
267        Ok(var.sqrt())
268    }
269
270    fn confidence_intervals(
271        &self,
272        pseudo_obs: &DMatrix<f64>,
273        confidence_level: f64,
274    ) -> Result<(Self::Parameters, Self::Parameters)> {
275        if !(0.0 < confidence_level && confidence_level < 1.0) {
276            return Err(CopulaError::invalid_parameter("confidence_level"));
277        }
278        let se = self.standard_errors(pseudo_obs)?;
279        let z = statrs::distribution::Normal::new(0.0, 1.0)
280            .map_err(|_| CopulaError::computation("failed to create standard normal distribution"))?
281            .inverse_cdf(0.5 + confidence_level / 2.0);
282        let lower = self.theta - z * se;
283        let upper = self.theta + z * se;
284        Ok((lower, upper))
285    }
286}
287
288#[cfg(test)]
289mod tests {
290    use super::*;
291
292    #[test]
293    fn cdf_matches_formula() {
294        let cop = ClaytonCopula::new(2.0).unwrap();
295        let cdf = cop.cdf(&[0.5, 0.5]).unwrap();
296        // Formula: (u^{-θ} + v^{-θ} - 1)^{-1/θ}
297        let expected = (0.5_f64.powf(-2.0) + 0.5_f64.powf(-2.0) - 1.0).powf(-0.5);
298        assert!((cdf - expected).abs() < 1e-10);
299    }
300
301    #[test]
302    fn phi_inverse_derivatives() {
303        let cop = ClaytonCopula::new(1.5).unwrap();
304        let s = 0.2;
305        let phi_inv = cop.phi_inv(s).unwrap();
306        // numerical first derivative via finite difference
307        let h = 1e-6;
308        let fd = (cop.phi_inv(s + h).unwrap() - phi_inv) / h;
309        let analytic = cop.phi_inv_deriv(s, 1).unwrap();
310        assert!((fd - analytic).abs() < 1e-4);
311    }
312
313    #[test]
314    fn new_rejects_invalid_theta() {
315        assert!(ClaytonCopula::new(0.0).is_err());
316        assert!(ClaytonCopula::new(-1.0).is_err());
317    }
318
319    #[test]
320    fn cdf_validates_input() {
321        let cop = ClaytonCopula::new(1.5).unwrap();
322        assert!(cop.cdf(&[0.5]).is_err());
323        assert!(cop.cdf(&[0.5, 1.2]).is_err());
324    }
325
326    #[test]
327    fn pdf_matches_formula() {
328        let cop = ClaytonCopula::new(2.0).unwrap();
329        let pdf = cop.pdf(&[0.4, 0.6]).unwrap();
330        let theta = 2.0;
331        let sum = 0.4_f64.powf(-theta) + 0.6_f64.powf(-theta) - 1.0;
332        let expected = (1.0 + theta)
333            * 0.4_f64.powf(-1.0 - theta)
334            * 0.6_f64.powf(-1.0 - theta)
335            * sum.powf(-2.0 - 1.0 / theta);
336        assert!((pdf - expected).abs() < 1e-10);
337    }
338
339    #[test]
340    fn sample_produces_valid_data() {
341        let mut rng = rand::rng();
342        let cop = ClaytonCopula::new(1.2).unwrap();
343        let samples = cop.sample(10, &mut rng).unwrap();
344        assert_eq!(samples.ncols(), 2);
345        assert_eq!(samples.nrows(), 10);
346        for i in 0..10 {
347            for j in 0..2 {
348                assert!(samples[(i, j)] > 0.0 && samples[(i, j)] < 1.0);
349            }
350        }
351    }
352
353    #[test]
354    fn tail_dependence_coefficients() {
355        let cop = ClaytonCopula::new(2.0).unwrap();
356        let (lower, upper) = cop.tail_dependence().unwrap();
357        assert!((lower - 2f64.powf(-0.5)).abs() < 1e-12);
358        assert_eq!(upper, 0.0);
359    }
360
361    #[test]
362    fn parameter_bounds_are_enforced() {
363        let bad = ClaytonCopula { theta: -1.0 };
364        assert!(bad.check_bounds().is_err());
365        let good = ClaytonCopula::new(2.0).unwrap();
366        assert!(good.check_bounds().is_ok());
367    }
368}
copula_core/archimedean/clayton.rs

copula_core/archimedean/
clayton.rs
