copula_core/elliptical/
gaussian.rs

1//! Gaussian (Normal) copula implementation.
2//!
3//! ## Bibliography
4//! - Embrechts, P., McNeil, A., & Straumann, D. (2002). Correlation and
5//!   dependence in risk management: properties and pitfalls. In *Risk
6//!   Management: Value at Risk and Beyond*.
7//! - McNeil, A. J., Frey, R., & Embrechts, P. (2015). *Quantitative Risk
8//!   Management: Concepts, Techniques and Tools*. Princeton University Press.
9//! - Nelsen, R. B. (2006). *An Introduction to Copulas*. Springer.
10
11#[cfg(feature = "estimation")]
12use crate::traits::FittableCopula;
13#[cfg(feature = "estimation")]
14use crate::utils::multivariate_kendall_tau;
15use crate::{utils::validate_correlation_matrix, Copula, CopulaError, Result};
16use mv_norm::tvpack::bvnd;
17use nalgebra::{DMatrix, DVector};
18use rand::Rng;
19use rand_distr::{Distribution, StandardNormal};
20use statrs::distribution::{ContinuousCDF, Normal};
21
22/// Gaussian copula placeholder
23#[derive(Debug, Clone)]
24pub struct GaussianCopula {
25    correlation: DMatrix<f64>,
26}
27
28validated_serde!("GaussianCopula", GaussianCopula { correlation: DMatrix<f64> } => GaussianCopula::new(correlation));
29
30impl GaussianCopula {
31    /// Create a Gaussian copula from a correlation matrix.
32    pub fn new(correlation: DMatrix<f64>) -> Result<Self> {
33        validate_correlation_matrix(&correlation)?;
34        Ok(Self { correlation })
35    }
36
37    /// Create a Gaussian copula with an identity correlation matrix of the given dimension.
38    pub fn new_identity(dim: usize) -> Result<Self> {
39        if dim < 2 {
40            return Err(CopulaError::invalid_parameter(
41                "Copula dimension must be at least 2",
42            ));
43        }
44        Ok(Self {
45            correlation: DMatrix::<f64>::identity(dim, dim),
46        })
47    }
48
49    /// The correlation matrix.
50    pub fn correlation(&self) -> &DMatrix<f64> {
51        &self.correlation
52    }
53
54    fn dim(&self) -> usize {
55        self.correlation.ncols()
56    }
57}
58
59impl Copula for GaussianCopula {
60    fn cdf(&self, u: &[f64]) -> Result<f64> {
61        if u.len() != self.dim() {
62            return Err(CopulaError::dimension_mismatch(self.dim(), u.len()));
63        }
64        crate::error::validate_unit_range(u)?;
65        // Quantile transforms are infinite on the boundary; the copula axioms
66        // give the exact value there.
67        if let Some(value) = crate::utils::copula_boundary_value(u) {
68            return Ok(value);
69        }
70        let normal = Normal::new(0.0, 1.0)
71            .map_err(|_| CopulaError::computation("failed to create Normal(0,1)"))?;
72        if self.dim() == 2 {
73            let x = normal.inverse_cdf(u[0]);
74            let y = normal.inverse_cdf(u[1]);
75            let r = self.correlation[(0, 1)];
76            return Ok(crate::utils::clamp_to_frechet_bounds(u, bvnd(-x, -y, r)));
77        }
78
79        // Monte Carlo approximation for higher dimensions
80        let dim = self.dim();
81        let quantiles: Vec<f64> = u.iter().map(|&ui| normal.inverse_cdf(ui)).collect();
82        let chol = self
83            .correlation
84            .clone()
85            .cholesky()
86            .ok_or_else(|| CopulaError::invalid_parameter("correlation not PD"))?;
87        let mut rng = rand::rng();
88        let normal = StandardNormal;
89
90        let mut count = 0usize;
91        let n_samples = 10_000usize;
92
93        for _ in 0..n_samples {
94            let z = DVector::from_iterator(dim, (0..dim).map(|_| normal.sample(&mut rng)));
95            let sample = chol.l() * z;
96            if sample.iter().zip(&quantiles).all(|(&s, &x)| s <= x) {
97                count += 1;
98            }
99        }
100
101        Ok(crate::utils::clamp_to_frechet_bounds(
102            u,
103            count as f64 / n_samples as f64,
104        ))
105    }
106
107    fn pdf(&self, u: &[f64]) -> Result<f64> {
108        if u.len() != self.dim() {
109            return Err(CopulaError::dimension_mismatch(self.dim(), u.len()));
110        }
111        crate::error::validate_unit_range(u)?;
112
113        if self.dim() == 1 {
114            return Ok(1.0);
115        }
116
117        let normal = Normal::new(0.0, 1.0)
118            .map_err(|_| CopulaError::computation("failed to create Normal(0,1)"))?;
119        let x = DVector::from_iterator(self.dim(), u.iter().map(|&ui| normal.inverse_cdf(ui)));
120        let inv = self
121            .correlation
122            .clone()
123            .try_inverse()
124            .ok_or_else(|| CopulaError::matrix_error("inverse", "singular"))?;
125        let det = self.correlation.determinant();
126        let quad = x.transpose() * (&inv * &x);
127        let norm_sq = x.dot(&x);
128        let exponent = -0.5 * (quad[(0, 0)] - norm_sq);
129        Ok(det.powf(-0.5) * exponent.exp())
130    }
131
132    fn sample<R: Rng + ?Sized>(&self, n: usize, rng: &mut R) -> Result<DMatrix<f64>> {
133        let dim = self.dim();
134        let chol = self
135            .correlation
136            .clone()
137            .cholesky()
138            .ok_or_else(|| CopulaError::invalid_parameter("correlation not PD"))?;
139        let normal = StandardNormal;
140        let std_normal = Normal::new(0.0, 1.0)
141            .map_err(|_| CopulaError::computation("failed to create Normal(0,1)"))?;
142        let mut samples = DMatrix::<f64>::zeros(n, dim);
143
144        for i in 0..n {
145            let z = DVector::from_iterator(dim, (0..dim).map(|_| normal.sample(rng)));
146            let x = chol.l() * z;
147            for j in 0..dim {
148                samples[(i, j)] = std_normal.cdf(x[j]);
149            }
150        }
151
152        Ok(samples)
153    }
154
155    fn dimension(&self) -> usize {
156        self.dim()
157    }
158}
159
160#[cfg(feature = "estimation")]
161impl FittableCopula for GaussianCopula {
162    type Parameters = DMatrix<f64>;
163
164    fn fit(&mut self, pseudo_obs: &DMatrix<f64>) -> Result<Self::Parameters> {
165        crate::utils::validate_pseudo_observations(pseudo_obs)?;
166        let n = pseudo_obs.nrows();
167        let dim = pseudo_obs.ncols();
168        let normal = Normal::new(0.0, 1.0).map_err(|_| {
169            CopulaError::computation("failed to create standard normal distribution")
170        })?;
171        let mut z = DMatrix::<f64>::zeros(n, dim);
172        for i in 0..n {
173            for j in 0..dim {
174                z[(i, j)] = normal.inverse_cdf(pseudo_obs[(i, j)]);
175            }
176        }
177
178        let mut corr = DMatrix::<f64>::identity(dim, dim);
179        for i in 0..dim {
180            for j in i + 1..dim {
181                let mut sum_i = 0.0;
182                let mut sum_j = 0.0;
183                for k in 0..n {
184                    sum_i += z[(k, i)];
185                    sum_j += z[(k, j)];
186                }
187                let mean_i = sum_i / n as f64;
188                let mean_j = sum_j / n as f64;
189                let mut cov = 0.0;
190                let mut var_i = 0.0;
191                let mut var_j = 0.0;
192                for k in 0..n {
193                    let xi = z[(k, i)] - mean_i;
194                    let xj = z[(k, j)] - mean_j;
195                    cov += xi * xj;
196                    var_i += xi * xi;
197                    var_j += xj * xj;
198                }
199                cov /= n as f64;
200                var_i /= n as f64;
201                var_j /= n as f64;
202                let r = cov / (var_i.sqrt() * var_j.sqrt());
203                corr[(i, j)] = r;
204                corr[(j, i)] = r;
205            }
206        }
207        validate_correlation_matrix(&corr)?;
208        self.correlation = corr.clone();
209        Ok(corr)
210    }
211
212    fn log_likelihood(&self, pseudo_obs: &DMatrix<f64>) -> Result<f64> {
213        crate::utils::validate_pseudo_observations(pseudo_obs)?;
214        if pseudo_obs.ncols() != self.dim() {
215            return Err(CopulaError::dimension_mismatch(
216                self.dim(),
217                pseudo_obs.ncols(),
218            ));
219        }
220        let n = pseudo_obs.nrows();
221        let normal = Normal::new(0.0, 1.0).map_err(|_| {
222            CopulaError::computation("failed to create standard normal distribution")
223        })?;
224        let mut ll = 0.0;
225        let inv = self
226            .correlation
227            .clone()
228            .try_inverse()
229            .ok_or_else(|| CopulaError::matrix_error("inverse", "singular"))?;
230        let det = self.correlation.determinant();
231        for i in 0..n {
232            let x = DVector::from_iterator(
233                self.dim(),
234                (0..self.dim()).map(|j| normal.inverse_cdf(pseudo_obs[(i, j)])),
235            );
236            let quad = x.transpose() * (&inv * &x);
237            let norm_sq = x.dot(&x);
238            ll += -0.5 * (det.ln() + quad[(0, 0)] - norm_sq);
239        }
240        Ok(ll)
241    }
242
243    fn fit_moments(&mut self, pseudo_obs: &DMatrix<f64>) -> Result<Self::Parameters> {
244        let tau = multivariate_kendall_tau(pseudo_obs)?;
245        let dim = tau.ncols();
246        let mut corr = DMatrix::<f64>::identity(dim, dim);
247        for i in 0..dim {
248            for j in (i + 1)..dim {
249                let val = (std::f64::consts::PI * 0.5 * tau[(i, j)]).sin();
250                corr[(i, j)] = val;
251                corr[(j, i)] = val;
252            }
253        }
254        validate_correlation_matrix(&corr)?;
255        self.correlation = corr.clone();
256        Ok(corr)
257    }
258
259    fn parameters(&self) -> Self::Parameters {
260        self.correlation.clone()
261    }
262
263    fn set_parameters(&mut self, params: Self::Parameters) -> Result<()> {
264        validate_correlation_matrix(&params)?;
265        self.correlation = params;
266        Ok(())
267    }
268}
269#[cfg(test)]
270mod tests {
271    use super::*;
272
273    #[test]
274    fn new_identity_sets_dimension() {
275        let cop = GaussianCopula::new_identity(3).unwrap();
276        assert_eq!(cop.dimension(), 3);
277    }
278
279    #[test]
280    fn cdf_identity_is_product() {
281        let cop = GaussianCopula::new_identity(2).unwrap();
282        let val = cop.cdf(&[0.1, 0.9]).unwrap();
283        assert!((val - 0.1 * 0.9).abs() < 1e-12);
284    }
285
286    #[test]
287    fn cdf_with_correlation() {
288        let corr = DMatrix::from_row_slice(2, 2, &[1.0, 0.5, 0.5, 1.0]);
289        let cop = GaussianCopula::new(corr).unwrap();
290        let normal = Normal::new(0.0, 1.0).unwrap();
291        let x = normal.inverse_cdf(0.4);
292        let y = normal.inverse_cdf(0.7);
293        let expected = bvnd(-x, -y, 0.5);
294        let val = cop.cdf(&[0.4, 0.7]).unwrap();
295        assert!((val - expected).abs() < 1e-12);
296    }
297
298    #[test]
299    fn cdf_higher_dimension_identity() {
300        let cop = GaussianCopula::new_identity(3).unwrap();
301        let val = cop.cdf(&[0.2, 0.3, 0.4]).unwrap();
302        let expected = 0.2 * 0.3 * 0.4;
303        assert!((val - expected).abs() < 0.02);
304    }
305
306    #[test]
307    fn pdf_identity_matches_one() {
308        let cop = GaussianCopula::new_identity(2).unwrap();
309        let pdf = cop.pdf(&[0.3, 0.7]).unwrap();
310        // Independence copula density is 1
311        assert!((pdf - 1.0).abs() < 1e-12);
312    }
313
314    #[test]
315    fn sample_dimensions() {
316        let mut rng = rand::rng();
317        let cop = GaussianCopula::new_identity(2).unwrap();
318        let samples = cop.sample(5, &mut rng).unwrap();
319        assert_eq!(samples.nrows(), 5);
320        assert_eq!(samples.ncols(), 2);
321        for i in 0..5 {
322            for j in 0..2 {
323                assert!(samples[(i, j)] > 0.0 && samples[(i, j)] < 1.0);
324            }
325        }
326    }
327}
copula_core/elliptical/gaussian.rs

copula_core/elliptical/
gaussian.rs
