quantwave-core 0.9.0

//! Portfolio optimization engine (Rust/nalgebra port, quantwave-a465).
//!
//! Rust port of the Python-first v1 shipped at
//! `quantwave-py/python/quantwave/portfolio.py` (p2k0.5, PR #21). That module
//! is the numerical SPEC: sample/EWMA/Ledoit-Wolf covariance estimators;
//! min-variance / max-Sharpe (closed-form) / risk-parity (Spinu fixed-point)
//! / HRP (hand-rolled single-linkage) optimizers; capped-simplex projection
//! for box+budget constraints. This module mirrors that algorithmically,
//! operating on `nalgebra` `DMatrix<f64>` / `DVector<f64>` instead of numpy
//! arrays, so it can be embedded in the core engine and driven from Rust
//! callers (backtests, live sizing) without a Python round-trip. PyO3
//! bindings in `quantwave-py` expose it to Python as `quantwave._portfolio`
//! (see that crate for the numpy <-> nalgebra bridge).
//!
//! Conventions
//! -----------
//! * Returns matrices are `(rows = time, cols = assets)`.
//! * Every optimizer returns a weights vector that sums to 1.0.
//! * Degenerate inputs (singular covariance, a single asset, NaN returns)
//!   return [`PortfolioError`] rather than panicking or returning NaN
//!   weights (workspace lint: `unwrap_used`/`expect_used` = warn,
//!   `panic` = deny).

mod covariance;
mod hrp;
mod optimizers;
mod projection;

pub use covariance::{ewma_cov, ledoit_wolf_cov, sample_cov};
pub use hrp::{HrpInputKind, hrp};
pub use optimizers::{max_sharpe, min_variance, risk_parity};
pub use projection::{PortfolioConstraint, apply_bounds_and_budget};

use thiserror::Error;

/// Per-asset `(low, high)` weight bounds.
pub type Bounds = (f64, f64);

/// Errors from the portfolio-optimization engine.
///
/// Mirrors `quantwave.portfolio.PortfolioError` in the Python v1: raised for
/// degenerate or invalid inputs (singular/non-invertible covariance, NaN/Inf
/// returns, empty asset universe, or bounds that cannot satisfy the
/// full-investment budget constraint) rather than letting a raw linear
/// algebra failure or NaN weights escape.
#[derive(Error, Debug, Clone, PartialEq)]
pub enum PortfolioError {
    #[error("{0} is empty")]
    Empty(String),

    #[error("{0} contains NaN or Inf values")]
    NonFinite(String),

    #[error("returns must have at least 2 time observations")]
    InsufficientObservations,

    #[error("returns must have at least 1 asset")]
    NoAssets,

    #[error("covariance must be square 2-D, got shape=({0}, {1})")]
    NonSquareCov(usize, usize),

    #[error("halflife must be positive, got {0}")]
    InvalidHalflife(f64),

    #[error("single-asset covariance is non-positive or invalid")]
    InvalidSingleAssetCov,

    #[error(
        "covariance matrix is singular or ill-conditioned (condition number={0:.3e}); cannot invert"
    )]
    IllConditionedCov(f64),

    #[error("covariance matrix is degenerate: {0}")]
    DegenerateCov(String),

    #[error("solve produced non-finite weights")]
    NonFiniteSolve,

    #[error("mean length {0} does not match covariance size {1}")]
    MeanLengthMismatch(usize, usize),

    #[error("min_variance: Sigma^-1 . 1 sums to ~0, cannot normalize")]
    MinVarianceDegenerate,

    #[error("max_sharpe: Sigma^-1 . mu sums to ~0, cannot normalize")]
    MaxSharpeDegenerate,

    #[error(
        "risk_parity: covariance produced non-positive marginal risk; matrix may not be positive semi-definite"
    )]
    RiskParityNonPositiveMarginalRisk,

    #[error("risk_parity failed to converge to finite weights")]
    RiskParityDidNotConverge,

    #[error("hrp produced degenerate weights")]
    HrpDegenerate,

    #[error("bounds length {0} does not match number of assets {1}")]
    BoundsLengthMismatch(usize, usize),

    #[error("bounds have low > high for at least one asset")]
    BoundsLowGreaterThanHigh,

    #[error("bounds cannot satisfy full-investment budget: sum(upper bounds) < 1")]
    BoundsUpperSumTooLow,

    #[error("bounds cannot satisfy full-investment budget: sum(lower bounds) > 1")]
    BoundsLowerSumTooHigh,

    #[error("unable to project weights onto the feasible bounds set")]
    ProjectionFailed,
}

pub(crate) type PortfolioResult<T> = Result<T, PortfolioError>;

/// Numerical parity tests against the numpy v1
/// (`quantwave-py/python/quantwave/portfolio.py`). Expected values below
/// were generated by running that module directly (bypassing the PyO3
/// package `__init__`, which needs a built native extension) via
/// `.venv/bin/python` against the fixed inputs reproduced here; see the
/// task report (bead quantwave-a465) for the exact script used.
#[cfg(test)]
mod parity_tests {
    use super::*;
    use approx::assert_relative_eq;
    use nalgebra::{DMatrix, DVector};

    fn toy_returns() -> DMatrix<f64> {
        DMatrix::from_row_slice(
            5,
            3,
            &[
                0.01, 0.02, -0.01, 0.015, -0.005, 0.02, -0.02, 0.01, 0.005, 0.005, 0.03, -0.01,
                0.0, 0.01, 0.015,
            ],
        )
    }

    fn cov4() -> DMatrix<f64> {
        DMatrix::from_row_slice(
            4,
            4,
            &[
                0.04, 0.03, 0.001, 0.0, 0.03, 0.05, 0.0, 0.001, 0.001, 0.0, 0.02, 0.015, 0.0,
                0.001, 0.015, 0.03,
            ],
        )
    }

    #[test]
    fn min_variance_matches_numpy_v1() {
        let r = toy_returns();
        let cov = sample_cov(&r).unwrap();
        let w = min_variance(&cov, None).unwrap();
        let expected = [0.06852607456263189, 0.4824772231848401, 0.4489967022525281];
        for i in 0..3 {
            assert_relative_eq!(w[i], expected[i], epsilon = 1e-8);
        }
    }

    #[test]
    fn max_sharpe_matches_numpy_v1() {
        let r = toy_returns();
        let cov = sample_cov(&r).unwrap();
        let mean = DVector::from_vec(vec![
            r.column(0).sum() / 5.0,
            r.column(1).sum() / 5.0,
            r.column(2).sum() / 5.0,
        ]);
        let w = max_sharpe(&mean, &cov, None).unwrap();
        let expected = [0.04134877384196476, 0.5069482288828321, 0.4517029972752033];
        for i in 0..3 {
            assert_relative_eq!(w[i], expected[i], epsilon = 1e-8);
        }
    }

    #[test]
    fn hrp_from_cov_matches_numpy_v1() {
        let r = toy_returns();
        let cov = sample_cov(&r).unwrap();
        let w = hrp(&cov, None, HrpInputKind::Covariance).unwrap();
        let expected = [0.32939214692886076, 0.35832698650221917, 0.3122808665689201];
        for i in 0..3 {
            assert_relative_eq!(w[i], expected[i], epsilon = 1e-8);
        }
    }

    #[test]
    fn hrp_from_returns_matches_numpy_v1() {
        let r = toy_returns();
        let w = hrp(&r, None, HrpInputKind::Returns).unwrap();
        let expected = [0.32939214692886076, 0.35832698650221917, 0.3122808665689201];
        for i in 0..3 {
            assert_relative_eq!(w[i], expected[i], epsilon = 1e-8);
        }
    }

    #[test]
    fn hrp4_matches_numpy_v1() {
        let cov = cov4();
        let w = hrp(&cov, None, HrpInputKind::Covariance).unwrap();
        let expected = [
            0.18967334035827263,
            0.15173867228661872,
            0.39515279241306456,
            0.26343519494204404,
        ];
        for i in 0..4 {
            assert_relative_eq!(w[i], expected[i], epsilon = 1e-8);
        }
    }

    #[test]
    fn min_variance4_bounded_matches_numpy_v1() {
        let cov = cov4();
        let bounds = vec![(0.0, 0.3); 4];
        let w = min_variance(&cov, Some(&bounds)).unwrap();
        let expected = [
            0.27721197871550607,
            0.18718599183269521,
            0.29999999999999744,
            0.2356020294518014,
        ];
        for i in 0..4 {
            assert_relative_eq!(w[i], expected[i], epsilon = 1e-8);
        }
    }

    #[test]
    fn risk_parity4_matches_numpy_v1() {
        let cov = cov4();
        let w = risk_parity(&cov, None, 200, 1e-10).unwrap();
        let expected = [
            0.22127417224640838,
            0.19884787598270773,
            0.3185765076370055,
            0.2613014441338784,
        ];
        for i in 0..4 {
            assert_relative_eq!(w[i], expected[i], epsilon = 1e-6);
        }
    }

    #[test]
    fn risk_parity_diag_matches_numpy_v1() {
        let cov = DMatrix::from_diagonal(&DVector::from_row_slice(&[0.04, 0.09, 0.01]));
        let w = risk_parity(&cov, None, 200, 1e-10).unwrap();
        let expected = [0.2727272727272733, 0.18181818181818327, 0.5454545454545433];
        for i in 0..3 {
            assert_relative_eq!(w[i], expected[i], epsilon = 1e-6);
        }
    }

    #[test]
    fn sample_cov_matches_numpy_v1() {
        let r = toy_returns();
        let cov = sample_cov(&r).unwrap();
        let expected = DMatrix::from_row_slice(
            3,
            3,
            &[
                0.00018250000000000002,
                -1.3749999999999999e-05,
                2.5e-06,
                -1.3749999999999999e-05,
                0.00017,
                -0.000165,
                2.5e-06,
                -0.000165,
                0.0001925,
            ],
        );
        for i in 0..3 {
            for j in 0..3 {
                assert_relative_eq!(cov[(i, j)], expected[(i, j)], epsilon = 1e-12);
            }
        }
    }

    #[test]
    fn ewma_cov_matches_numpy_v1() {
        let r = toy_returns();
        let cov = ewma_cov(&r, 3.0).unwrap();
        let expected = DMatrix::from_row_slice(
            3,
            3,
            &[
                0.0001629382102880254,
                -1.5433060639633344e-06,
                -2.4151258161425324e-06,
                -1.5433060639633435e-06,
                0.00016378886197253209,
                -0.00015799870649033202,
                -2.415125816142528e-06,
                -0.00015799870649033202,
                0.00018564391979524118,
            ],
        );
        for i in 0..3 {
            for j in 0..3 {
                assert_relative_eq!(cov[(i, j)], expected[(i, j)], epsilon = 1e-12);
            }
        }
    }

    #[test]
    fn ledoit_wolf_cov_matches_numpy_v1() {
        let r = toy_returns();
        let cov = ledoit_wolf_cov(&r).unwrap();
        let expected = DMatrix::from_row_slice(
            3,
            3,
            &[
                0.0001542774782099034,
                -3.3089325216010282e-06,
                6.016240948365506e-07,
                -3.3089325216010282e-06,
                0.00015126935773572065,
                -3.970719025921234e-05,
                6.016240948365506e-07,
                -3.970719025921234e-05,
                0.0001566839745892496,
            ],
        );
        for i in 0..3 {
            for j in 0..3 {
                assert_relative_eq!(cov[(i, j)], expected[(i, j)], epsilon = 1e-12);
            }
        }
    }
}