fin_primitives/cross_asset/
mod.rs

1//! Cross-asset rolling correlations and PCA-based dimensionality reduction.
2//!
3//! ## Responsibility
4//! Rolling cross-asset correlation tracking and dimensionality reduction via
5//! Principal Component Analysis (PCA) on the correlation matrix.
6//!
7//! ## Guarantees
8//! - `CrossAssetCorrelation` returns `None` until `window` samples have been seen
9//! - Correlation values are clamped to `[-1, 1]` to absorb floating-point rounding
10//! - `PrincipalComponents` extracts up to 3 PCs using the power-iteration method;
11//!   returns `None` when fewer than 2 assets are tracked or matrix not ready
12//! - All arithmetic is in `f64` (sufficient for statistical computation)
13//! - Zero panics; validation errors returned as `FinError`
14//!
15//! ## NOT Responsible For
16//! - Causal inference or portfolio optimization
17//! - Persistence
18
19use crate::error::FinError;
20use std::collections::VecDeque;
21
22// ─────────────────────────────────────────
23//  CrossAssetCorrelation
24// ─────────────────────────────────────────
25
26/// Tracks rolling Pearson correlations between N instruments.
27///
28/// Feed one return observation per bar per instrument via [`update`](CrossAssetCorrelation::update).
29/// Once `window` samples have been seen, the full NxN [`CorrelationMatrix`] is available.
30///
31/// # Example
32/// ```rust
33/// use fin_primitives::cross_asset::CrossAssetCorrelation;
34///
35/// let mut cac = CrossAssetCorrelation::new(
36///     vec!["SPY".into(), "QQQ".into(), "IWM".into()],
37///     20,
38/// ).unwrap();
39///
40/// for i in 0..20 {
41///     let returns = vec![0.001 * i as f64, 0.0012 * i as f64, 0.0008 * i as f64];
42///     cac.update(&returns).unwrap();
43/// }
44///
45/// let matrix = cac.correlation_matrix();
46/// assert!(matrix.is_some());
47/// ```
48#[derive(Debug)]
49pub struct CrossAssetCorrelation {
50    /// Instrument names in column order.
51    names: Vec<String>,
52    /// Number of instruments.
53    n: usize,
54    /// Rolling window size.
55    window: usize,
56    /// Circular buffer: each entry is one bar's vector of `n` return values.
57    buf: VecDeque<Vec<f64>>,
58}
59
60impl CrossAssetCorrelation {
61    /// Constructs a `CrossAssetCorrelation` tracker.
62    ///
63    /// # Parameters
64    /// - `names`: instrument names (must have >= 2 elements, no duplicates).
65    /// - `window`: rolling window in bars (must be >= 2).
66    ///
67    /// # Errors
68    /// Returns [`FinError::InvalidInput`] if fewer than 2 instruments are provided or names
69    /// contain duplicates.
70    /// Returns [`FinError::InvalidPeriod`] if `window < 2`.
71    pub fn new(names: Vec<String>, window: usize) -> Result<Self, FinError> {
72        if names.len() < 2 {
73            return Err(FinError::InvalidInput(
74                "CrossAssetCorrelation requires at least 2 instruments".to_owned(),
75            ));
76        }
77        if window < 2 {
78            return Err(FinError::InvalidPeriod(window));
79        }
80        // Check for duplicate names
81        for (i, name) in names.iter().enumerate() {
82            for (j, other) in names.iter().enumerate() {
83                if i != j && name == other {
84                    return Err(FinError::InvalidInput(format!(
85                        "duplicate instrument name: '{name}'"
86                    )));
87                }
88            }
89        }
90        let n = names.len();
91        Ok(Self { names, n, window, buf: VecDeque::with_capacity(window) })
92    }
93
94    /// Returns the instrument names.
95    pub fn names(&self) -> &[String] {
96        &self.names
97    }
98
99    /// Returns the number of instruments tracked.
100    pub fn n_instruments(&self) -> usize {
101        self.n
102    }
103
104    /// Returns the rolling window size.
105    pub fn window(&self) -> usize {
106        self.window
107    }
108
109    /// Returns the number of samples currently buffered.
110    pub fn sample_count(&self) -> usize {
111        self.buf.len()
112    }
113
114    /// Returns `true` when enough samples have been collected.
115    pub fn is_ready(&self) -> bool {
116        self.buf.len() >= self.window
117    }
118
119    /// Records one bar's returns for all instruments.
120    ///
121    /// `returns.len()` must equal `self.n_instruments()`.
122    ///
123    /// # Errors
124    /// Returns [`FinError::InvalidInput`] if the length is wrong or any value is non-finite.
125    pub fn update(&mut self, returns: &[f64]) -> Result<(), FinError> {
126        if returns.len() != self.n {
127            return Err(FinError::InvalidInput(format!(
128                "expected {} returns, got {}",
129                self.n,
130                returns.len()
131            )));
132        }
133        for (i, r) in returns.iter().enumerate() {
134            if !r.is_finite() {
135                return Err(FinError::InvalidInput(format!(
136                    "return[{i}] is not finite: {r}"
137                )));
138            }
139        }
140        self.buf.push_back(returns.to_vec());
141        if self.buf.len() > self.window {
142            self.buf.pop_front();
143        }
144        Ok(())
145    }
146
147    /// Returns the Pearson correlation between instruments `i` and `j`.
148    ///
149    /// Returns `None` when fewer than `window` samples have been seen, indices are
150    /// out of bounds, or either instrument has zero variance.
151    pub fn get(&self, i: usize, j: usize) -> Option<f64> {
152        if !self.is_ready() || i >= self.n || j >= self.n {
153            return None;
154        }
155        if i == j {
156            return Some(1.0);
157        }
158        let n = self.buf.len() as f64;
159        let mut sx = 0.0_f64;
160        let mut sy = 0.0_f64;
161        let mut sxy = 0.0_f64;
162        let mut sx2 = 0.0_f64;
163        let mut sy2 = 0.0_f64;
164        for row in &self.buf {
165            let x = row[i];
166            let y = row[j];
167            sx += x;
168            sy += y;
169            sxy += x * y;
170            sx2 += x * x;
171            sy2 += y * y;
172        }
173        let num = n * sxy - sx * sy;
174        let den_sq = (n * sx2 - sx * sx) * (n * sy2 - sy * sy);
175        if den_sq <= 0.0 {
176            return None;
177        }
178        Some((num / den_sq.sqrt()).clamp(-1.0, 1.0))
179    }
180
181    /// Returns the index of an instrument by name, or `None` if not found.
182    pub fn index_of(&self, name: &str) -> Option<usize> {
183        self.names.iter().position(|n| n == name)
184    }
185
186    /// Returns the full NxN correlation matrix (row-major `Vec<f64>`).
187    ///
188    /// Element at `(i, j)` is at index `i * n + j`.
189    /// Returns `None` until ready.
190    pub fn correlation_matrix(&self) -> Option<CorrelationMatrix> {
191        if !self.is_ready() {
192            return None;
193        }
194        let mut mat = vec![0.0_f64; self.n * self.n];
195        for i in 0..self.n {
196            for j in 0..self.n {
197                mat[i * self.n + j] = self.get(i, j).unwrap_or(0.0);
198            }
199        }
200        Some(CorrelationMatrix {
201            n: self.n,
202            data: mat,
203            names: self.names.clone(),
204        })
205    }
206
207    /// Resets the tracker.
208    pub fn reset(&mut self) {
209        self.buf.clear();
210    }
211}
212
213// ─────────────────────────────────────────
214//  CorrelationMatrix
215// ─────────────────────────────────────────
216
217/// An NxN symmetric correlation matrix (row-major, `f64`).
218///
219/// Produced by [`CrossAssetCorrelation::correlation_matrix`].
220#[derive(Debug, Clone)]
221pub struct CorrelationMatrix {
222    /// Dimension.
223    n: usize,
224    /// Row-major data, length `n * n`.
225    data: Vec<f64>,
226    /// Names of the instruments (columns = rows).
227    names: Vec<String>,
228}
229
230impl CorrelationMatrix {
231    /// Returns the correlation between instruments `i` and `j`, or `None` if out of bounds.
232    pub fn get(&self, i: usize, j: usize) -> Option<f64> {
233        if i >= self.n || j >= self.n {
234            return None;
235        }
236        Some(self.data[i * self.n + j])
237    }
238
239    /// Returns the raw row-major data slice.
240    pub fn data(&self) -> &[f64] {
241        &self.data
242    }
243
244    /// Returns the matrix dimension N.
245    pub fn n(&self) -> usize {
246        self.n
247    }
248
249    /// Returns the instrument names.
250    pub fn names(&self) -> &[String] {
251        &self.names
252    }
253
254    /// Returns the most correlated instrument indices for instrument `i`,
255    /// sorted descending by absolute correlation.
256    pub fn most_correlated_with(&self, i: usize) -> Vec<(usize, f64)> {
257        if i >= self.n {
258            return vec![];
259        }
260        let mut pairs: Vec<(usize, f64)> = (0..self.n)
261            .filter(|&j| j != i)
262            .filter_map(|j| self.get(i, j).map(|r| (j, r)))
263            .collect();
264        pairs.sort_by(|a, b| {
265            b.1.abs()
266                .partial_cmp(&a.1.abs())
267                .unwrap_or(std::cmp::Ordering::Equal)
268        });
269        pairs
270    }
271}
272
273// ─────────────────────────────────────────
274//  PrincipalComponents
275// ─────────────────────────────────────────
276
277/// Principal Component Analysis on a correlation matrix.
278///
279/// Extracts the first `k` (up to 3) principal components via the power-iteration
280/// method (deflation). Useful for dimensionality reduction and regime identification.
281///
282/// # Example
283/// ```rust
284/// use fin_primitives::cross_asset::{CrossAssetCorrelation, PrincipalComponents};
285///
286/// let mut cac = CrossAssetCorrelation::new(
287///     vec!["A".into(), "B".into(), "C".into()],
288///     10,
289/// ).unwrap();
290/// for i in 0..10 {
291///     cac.update(&[i as f64 * 0.01, i as f64 * 0.012, i as f64 * 0.009]).unwrap();
292/// }
293/// if let Some(mat) = cac.correlation_matrix() {
294///     let pca = PrincipalComponents::from_matrix(&mat, 2).unwrap();
295///     assert!(pca.explained_variance_ratio()[0] >= 0.0);
296/// }
297/// ```
298#[derive(Debug, Clone)]
299pub struct PrincipalComponents {
300    /// Number of PCs extracted.
301    k: usize,
302    /// Principal component vectors, one per PC, length n each.
303    components: Vec<Vec<f64>>,
304    /// Eigenvalues (variance explained by each PC).
305    eigenvalues: Vec<f64>,
306    /// Total variance (sum of all diagonal entries of the correlation matrix = n).
307    total_variance: f64,
308}
309
310impl PrincipalComponents {
311    /// Extracts up to `k` principal components from a [`CorrelationMatrix`].
312    ///
313    /// Uses power iteration with deflation. `k` is clamped to `min(k, n - 1, 3)`.
314    ///
315    /// # Errors
316    /// Returns [`FinError::InvalidInput`] if `n < 2` or `k == 0`.
317    pub fn from_matrix(matrix: &CorrelationMatrix, k: usize) -> Result<Self, FinError> {
318        let n = matrix.n();
319        if n < 2 {
320            return Err(FinError::InvalidInput(
321                "PCA requires at least 2 instruments".to_owned(),
322            ));
323        }
324        if k == 0 {
325            return Err(FinError::InvalidInput("k must be at least 1".to_owned()));
326        }
327        let k_actual = k.min(n - 1).min(3);
328        let total_variance = n as f64; // trace of correlation matrix = n
329
330        // Work on a mutable copy (for deflation)
331        let mut work: Vec<f64> = matrix.data().to_vec();
332
333        let mut components = Vec::with_capacity(k_actual);
334        let mut eigenvalues = Vec::with_capacity(k_actual);
335
336        for _ in 0..k_actual {
337            let (eigval, eigvec) = power_iterate(&work, n, 200, 1e-8);
338            if eigval <= 0.0 {
339                break;
340            }
341            // Deflate: A = A - λ * v * v^T
342            for i in 0..n {
343                for j in 0..n {
344                    work[i * n + j] -= eigval * eigvec[i] * eigvec[j];
345                }
346            }
347            eigenvalues.push(eigval);
348            components.push(eigvec);
349        }
350
351        Ok(Self { k: components.len(), components, eigenvalues, total_variance })
352    }
353
354    /// Returns the number of PCs extracted.
355    pub fn k(&self) -> usize {
356        self.k
357    }
358
359    /// Returns the PC vectors (each of length N instruments).
360    pub fn components(&self) -> &[Vec<f64>] {
361        &self.components
362    }
363
364    /// Returns the eigenvalues (variance explained by each PC).
365    pub fn eigenvalues(&self) -> &[f64] {
366        &self.eigenvalues
367    }
368
369    /// Returns the proportion of variance explained by each PC (`eigenvalue / total_variance`).
370    pub fn explained_variance_ratio(&self) -> Vec<f64> {
371        if self.total_variance <= 0.0 {
372            return vec![0.0; self.k];
373        }
374        self.eigenvalues.iter().map(|e| e / self.total_variance).collect()
375    }
376
377    /// Projects a returns vector (length N) onto the first `m` PCs.
378    ///
379    /// Returns a `Vec<f64>` of length `min(m, self.k)`.
380    pub fn project(&self, returns: &[f64], m: usize) -> Vec<f64> {
381        let take = m.min(self.k);
382        self.components[..take]
383            .iter()
384            .map(|pc| pc.iter().zip(returns.iter()).map(|(a, b)| a * b).sum())
385            .collect()
386    }
387}
388
389/// Power-iteration algorithm to find the dominant eigenvector and eigenvalue.
390///
391/// Returns `(eigenvalue, eigenvector)`. The eigenvector is L2-normalized.
392/// If the matrix is zero or degenerate, returns `(0.0, vec![0; n])`.
393fn power_iterate(matrix: &[f64], n: usize, max_iter: usize, tol: f64) -> (f64, Vec<f64>) {
394    // Initialize with ones (uniform start)
395    let mut v: Vec<f64> = vec![1.0 / (n as f64).sqrt(); n];
396
397    for _ in 0..max_iter {
398        // w = A * v
399        let mut w = vec![0.0_f64; n];
400        for i in 0..n {
401            for j in 0..n {
402                w[i] += matrix[i * n + j] * v[j];
403            }
404        }
405        // Compute eigenvalue estimate (Rayleigh quotient numerator)
406        let norm: f64 = w.iter().map(|x| x * x).sum::<f64>().sqrt();
407        if norm == 0.0 {
408            return (0.0, vec![0.0; n]);
409        }
410        let new_v: Vec<f64> = w.iter().map(|x| x / norm).collect();
411        // Check convergence
412        let diff: f64 = v.iter().zip(&new_v).map(|(a, b)| (a - b).powi(2)).sum::<f64>().sqrt();
413        v = new_v;
414        if diff < tol {
415            break;
416        }
417    }
418
419    // Eigenvalue = v^T A v
420    let mut eigenvalue = 0.0_f64;
421    for i in 0..n {
422        let mut av_i = 0.0_f64;
423        for j in 0..n {
424            av_i += matrix[i * n + j] * v[j];
425        }
426        eigenvalue += v[i] * av_i;
427    }
428
429    (eigenvalue, v)
430}
431
432#[cfg(test)]
433mod tests {
434    use super::*;
435
436    fn make_cac(n: usize, window: usize) -> CrossAssetCorrelation {
437        let names: Vec<String> = (0..n).map(|i| format!("asset{i}")).collect();
438        CrossAssetCorrelation::new(names, window).unwrap()
439    }
440
441    // ── CrossAssetCorrelation ─────────────────────────────────────────────
442
443    #[test]
444    fn test_too_few_assets_fails() {
445        assert!(CrossAssetCorrelation::new(vec!["A".into()], 10).is_err());
446    }
447
448    #[test]
449    fn test_duplicate_name_fails() {
450        assert!(CrossAssetCorrelation::new(
451            vec!["A".into(), "A".into()],
452            10
453        )
454        .is_err());
455    }
456
457    #[test]
458    fn test_window_1_fails() {
459        assert!(CrossAssetCorrelation::new(vec!["A".into(), "B".into()], 1).is_err());
460    }
461
462    #[test]
463    fn test_not_ready_before_window() {
464        let mut cac = make_cac(2, 5);
465        cac.update(&[0.01, 0.02]).unwrap();
466        assert!(!cac.is_ready());
467        assert!(cac.correlation_matrix().is_none());
468    }
469
470    #[test]
471    fn test_perfect_correlation_detected() {
472        let mut cac = make_cac(2, 5);
473        for i in 1..=5 {
474            // asset0 and asset1 are identical → r = 1.0
475            cac.update(&[i as f64 * 0.01, i as f64 * 0.01]).unwrap();
476        }
477        let r = cac.get(0, 1).unwrap();
478        assert!((r - 1.0).abs() < 1e-9, "r={r}");
479    }
480
481    #[test]
482    fn test_self_correlation_is_one() {
483        let mut cac = make_cac(2, 5);
484        for i in 1..=5 {
485            cac.update(&[i as f64 * 0.01, i as f64 * 0.02]).unwrap();
486        }
487        let r = cac.get(0, 0).unwrap();
488        assert_eq!(r, 1.0);
489    }
490
491    #[test]
492    fn test_correlation_matrix_shape() {
493        let mut cac = make_cac(3, 4);
494        for i in 1..=4 {
495            cac.update(&[i as f64, i as f64 * 2.0, i as f64 * 0.5]).unwrap();
496        }
497        let mat = cac.correlation_matrix().unwrap();
498        assert_eq!(mat.data().len(), 9);
499        assert_eq!(mat.n(), 3);
500        // Diagonal should be 1.0
501        assert!((mat.get(0, 0).unwrap() - 1.0).abs() < 1e-9);
502        assert!((mat.get(1, 1).unwrap() - 1.0).abs() < 1e-9);
503        assert!((mat.get(2, 2).unwrap() - 1.0).abs() < 1e-9);
504    }
505
506    #[test]
507    fn test_wrong_length_update_fails() {
508        let mut cac = make_cac(3, 5);
509        assert!(cac.update(&[0.01, 0.02]).is_err());
510    }
511
512    #[test]
513    fn test_nan_update_fails() {
514        let mut cac = make_cac(2, 5);
515        assert!(cac.update(&[f64::NAN, 0.01]).is_err());
516    }
517
518    #[test]
519    fn test_index_of() {
520        let cac = CrossAssetCorrelation::new(
521            vec!["SPY".into(), "QQQ".into()],
522            5,
523        )
524        .unwrap();
525        assert_eq!(cac.index_of("SPY"), Some(0));
526        assert_eq!(cac.index_of("QQQ"), Some(1));
527        assert_eq!(cac.index_of("MISSING"), None);
528    }
529
530    #[test]
531    fn test_reset_clears_buffer() {
532        let mut cac = make_cac(2, 3);
533        for i in 0..3 {
534            cac.update(&[i as f64, i as f64 * 2.0]).unwrap();
535        }
536        assert!(cac.is_ready());
537        cac.reset();
538        assert!(!cac.is_ready());
539        assert_eq!(cac.sample_count(), 0);
540    }
541
542    #[test]
543    fn test_most_correlated_with_sorted() {
544        let mut cac = make_cac(3, 5);
545        for i in 1..=5 {
546            let v = i as f64;
547            cac.update(&[v, v * 2.0, -v]).unwrap();
548        }
549        let mat = cac.correlation_matrix().unwrap();
550        let corrs = mat.most_correlated_with(0);
551        assert_eq!(corrs.len(), 2);
552        assert!(corrs[0].1.abs() >= corrs[1].1.abs());
553    }
554
555    // ── PrincipalComponents ───────────────────────────────────────────────
556
557    #[test]
558    fn test_pca_explained_variance_sums_to_at_most_one() {
559        let mut cac = make_cac(3, 10);
560        for i in 1..=10 {
561            let v = i as f64;
562            cac.update(&[v, v * 1.1, v * 0.9]).unwrap();
563        }
564        let mat = cac.correlation_matrix().unwrap();
565        let pca = PrincipalComponents::from_matrix(&mat, 3).unwrap();
566        let total: f64 = pca.explained_variance_ratio().iter().sum();
567        assert!(total <= 1.0 + 1e-9, "total explained variance ratio={total}");
568        assert!(total >= 0.0);
569    }
570
571    #[test]
572    fn test_pca_k_zero_fails() {
573        let mut cac = make_cac(2, 5);
574        for i in 1..=5 {
575            cac.update(&[i as f64, i as f64 * 2.0]).unwrap();
576        }
577        let mat = cac.correlation_matrix().unwrap();
578        assert!(PrincipalComponents::from_matrix(&mat, 0).is_err());
579    }
580
581    #[test]
582    fn test_pca_project_length() {
583        let mut cac = make_cac(3, 10);
584        for i in 1..=10 {
585            let v = i as f64;
586            cac.update(&[v, v * 1.5, -v]).unwrap();
587        }
588        let mat = cac.correlation_matrix().unwrap();
589        let pca = PrincipalComponents::from_matrix(&mat, 2).unwrap();
590        let proj = pca.project(&[0.01, 0.02, -0.01], 2);
591        assert!(proj.len() <= 2);
592    }
593
594    #[test]
595    fn test_pca_first_eigenvalue_largest() {
596        let mut cac = make_cac(3, 15);
597        for i in 1..=15 {
598            let v = i as f64;
599            cac.update(&[v, v * 1.2, v * 0.8]).unwrap();
600        }
601        let mat = cac.correlation_matrix().unwrap();
602        let pca = PrincipalComponents::from_matrix(&mat, 3).unwrap();
603        let evs = pca.eigenvalues();
604        if evs.len() >= 2 {
605            assert!(evs[0] >= evs[1], "first eigenvalue should be largest");
606        }
607    }
608}
fin_primitives/cross_asset/mod.rs

fin_primitives/cross_asset/
mod.rs
