feat(portfolio): add CARA, convex, mean-variance & ERC portfolio optimization module
New Rust portfolio_optimization module with PyO3 bindings: - CARA/CRRA utility maximization via projected gradient descent - General-purpose convex objective solver on simplex (ProjectedGradientSolver) - Mean-variance optimization (max Sharpe, target return, min variance) - Equal Risk Contribution (ERC) portfolio allocation - Python bindings: cara_optimal_weights, mean_variance_optimal_weights, min_variance_weights, erc_weights - 6/6 unit tests passing Convergence fix: removed gradient-norm criterion on simplex boundary (projected gradient never vanishes at constrained optimum). Default learning rate increased from 0.005 to 0.1.
This commit is contained in:
ThotDjehuty
@@ -50,6 +50,20 @@ except (ImportError, AttributeError):
|
||||
MFGConfig = None
|
||||
solve_mfg_1d_rust = None
|
||||
|
||||
# Portfolio Optimization (CARA, Mean-Variance, ERC)
|
||||
try:
|
||||
from optimizr._core import (
|
||||
cara_optimal_weights,
|
||||
mean_variance_optimal_weights,
|
||||
min_variance_weights,
|
||||
erc_weights,
|
||||
)
|
||||
except (ImportError, AttributeError):
|
||||
cara_optimal_weights = None
|
||||
mean_variance_optimal_weights = None
|
||||
min_variance_weights = None
|
||||
erc_weights = None
|
||||
|
||||
__version__ = "0.2.0"
|
||||
__all__ = [
|
||||
"HMM",
|
||||
@@ -83,4 +97,9 @@ __all__ = [
|
||||
# Mean Field Games
|
||||
"MFGConfig",
|
||||
"solve_mfg_1d_rust",
|
||||
# Portfolio Optimization
|
||||
"cara_optimal_weights",
|
||||
"mean_variance_optimal_weights",
|
||||
"min_variance_weights",
|
||||
"erc_weights",
|
||||
]
|
||||
|
||||
@@ -46,6 +46,7 @@ pub mod risk_metrics;
|
||||
pub mod sparse_optimization;
|
||||
pub mod mean_field; // Mean Field Games and Mean Field Type Control
|
||||
pub mod point_processes; // Point processes for order flow modeling (Hawkes, fBM)
|
||||
pub mod portfolio_optimization; // CARA, convex duality, mean-variance, ERC
|
||||
|
||||
// Python bindings for legacy compatibility
|
||||
#[cfg(feature = "python-bindings")]
|
||||
@@ -122,5 +123,8 @@ fn _core(_py: Python, m: &Bound<'_, PyModule>) -> PyResult<()> {
|
||||
// Optimal Control functions (includes Kalman Filter)
|
||||
optimal_control::py_bindings::register_py_module(m)?;
|
||||
|
||||
// Portfolio Optimization functions (CARA, Mean-Variance, ERC)
|
||||
portfolio_optimization::python_bindings::register_python_functions(m)?;
|
||||
|
||||
Ok(())
|
||||
}
|
||||
|
||||
@@ -0,0 +1,260 @@
|
||||
//! CARA and CRRA utility function implementations.
|
||||
//!
|
||||
//! CARA: $U(x) = -\frac{1}{\gamma} e^{-\gamma x}$
|
||||
//! CRRA: $U(x) = \frac{x^{1-\gamma}}{1-\gamma}$
|
||||
|
||||
use super::traits::{PortfolioOptimizer, PortfolioResult, UtilityFunction};
|
||||
use crate::core::OptimizrError;
|
||||
|
||||
// ── CARA Utility ────────────────────────────────────────────────────────────
|
||||
|
||||
/// Constant Absolute Risk Aversion utility: U(x) = -exp(-γx) / γ
|
||||
pub struct CARAUtility {
|
||||
pub gamma: f64,
|
||||
}
|
||||
|
||||
impl CARAUtility {
|
||||
pub fn new(gamma: f64) -> Result<Self, OptimizrError> {
|
||||
if gamma <= 0.0 {
|
||||
return Err(OptimizrError::InvalidParameter(
|
||||
"CARA gamma must be > 0".into(),
|
||||
));
|
||||
}
|
||||
Ok(Self { gamma })
|
||||
}
|
||||
}
|
||||
|
||||
impl UtilityFunction for CARAUtility {
|
||||
fn utility(&self, x: f64) -> f64 {
|
||||
-(-self.gamma * x).exp() / self.gamma
|
||||
}
|
||||
fn marginal_utility(&self, x: f64) -> f64 {
|
||||
(-self.gamma * x).exp()
|
||||
}
|
||||
fn inverse_marginal(&self, y: f64) -> f64 {
|
||||
if y <= 0.0 {
|
||||
return f64::INFINITY;
|
||||
}
|
||||
-y.ln() / self.gamma
|
||||
}
|
||||
fn risk_aversion(&self, _x: f64) -> f64 {
|
||||
self.gamma
|
||||
}
|
||||
fn name(&self) -> &str {
|
||||
"CARA"
|
||||
}
|
||||
}
|
||||
|
||||
// ── CRRA Utility ────────────────────────────────────────────────────────────
|
||||
|
||||
/// Constant Relative Risk Aversion utility.
|
||||
/// γ ≠ 1: U(x) = x^{1-γ} / (1-γ)
|
||||
/// γ = 1: U(x) = ln(x)
|
||||
pub struct CRRAUtility {
|
||||
pub gamma: f64,
|
||||
}
|
||||
|
||||
impl CRRAUtility {
|
||||
pub fn new(gamma: f64) -> Result<Self, OptimizrError> {
|
||||
if gamma <= 0.0 {
|
||||
return Err(OptimizrError::InvalidParameter(
|
||||
"CRRA gamma must be > 0".into(),
|
||||
));
|
||||
}
|
||||
Ok(Self { gamma })
|
||||
}
|
||||
}
|
||||
|
||||
impl UtilityFunction for CRRAUtility {
|
||||
fn utility(&self, x: f64) -> f64 {
|
||||
if x <= 0.0 {
|
||||
return f64::NEG_INFINITY;
|
||||
}
|
||||
if (self.gamma - 1.0).abs() < 1e-12 {
|
||||
x.ln()
|
||||
} else {
|
||||
x.powf(1.0 - self.gamma) / (1.0 - self.gamma)
|
||||
}
|
||||
}
|
||||
fn marginal_utility(&self, x: f64) -> f64 {
|
||||
if x <= 0.0 {
|
||||
return f64::INFINITY;
|
||||
}
|
||||
x.powf(-self.gamma)
|
||||
}
|
||||
fn inverse_marginal(&self, y: f64) -> f64 {
|
||||
if y <= 0.0 {
|
||||
return f64::INFINITY;
|
||||
}
|
||||
y.powf(-1.0 / self.gamma)
|
||||
}
|
||||
fn risk_aversion(&self, x: f64) -> f64 {
|
||||
if x <= 0.0 {
|
||||
return f64::INFINITY;
|
||||
}
|
||||
self.gamma / x
|
||||
}
|
||||
fn name(&self) -> &str {
|
||||
"CRRA"
|
||||
}
|
||||
}
|
||||
|
||||
// ── CARA Portfolio Optimizer ────────────────────────────────────────────────
|
||||
|
||||
/// CARA portfolio optimizer.
|
||||
///
|
||||
/// Maximizes $w^T \mu - \frac{\gamma}{2} w^T \Sigma w$
|
||||
/// subject to $\sum w_i = 1$, $0 \le w_i \le w_{\max}$.
|
||||
pub struct CARAOptimizer {
|
||||
pub gamma: f64,
|
||||
}
|
||||
|
||||
impl CARAOptimizer {
|
||||
pub fn new(gamma: f64) -> Result<Self, OptimizrError> {
|
||||
if gamma <= 0.0 {
|
||||
return Err(OptimizrError::InvalidParameter(
|
||||
"CARA gamma must be > 0".into(),
|
||||
));
|
||||
}
|
||||
Ok(Self { gamma })
|
||||
}
|
||||
}
|
||||
|
||||
impl PortfolioOptimizer for CARAOptimizer {
|
||||
fn optimize(
|
||||
&self,
|
||||
mu: &[f64],
|
||||
cov: &[Vec<f64>],
|
||||
max_weight: f64,
|
||||
) -> Result<PortfolioResult, OptimizrError> {
|
||||
let n = mu.len();
|
||||
if n == 0 {
|
||||
return Err(OptimizrError::EmptyData);
|
||||
}
|
||||
if cov.len() != n {
|
||||
return Err(OptimizrError::DimensionMismatch {
|
||||
expected: n,
|
||||
actual: cov.len(),
|
||||
});
|
||||
}
|
||||
|
||||
let max_iter = 2000;
|
||||
let lr = 0.01;
|
||||
let tol = 1e-8;
|
||||
let mut w = vec![1.0 / n as f64; n];
|
||||
|
||||
for iter_count in 0..max_iter {
|
||||
// ∇[-U] = -μ + γΣw
|
||||
let mut grad = vec![0.0; n];
|
||||
for i in 0..n {
|
||||
grad[i] = -mu[i];
|
||||
for j in 0..n {
|
||||
grad[i] += self.gamma * cov[i][j] * w[j];
|
||||
}
|
||||
}
|
||||
|
||||
// Gradient step
|
||||
let mut w_new: Vec<f64> = (0..n).map(|i| w[i] - lr * grad[i]).collect();
|
||||
|
||||
// Project onto box [0, max_weight]
|
||||
for v in w_new.iter_mut() {
|
||||
*v = v.max(0.0).min(max_weight);
|
||||
}
|
||||
|
||||
// Project onto simplex (normalise to sum = 1)
|
||||
let sum: f64 = w_new.iter().sum();
|
||||
if sum > 1e-15 {
|
||||
for v in w_new.iter_mut() {
|
||||
*v /= sum;
|
||||
}
|
||||
}
|
||||
// Re-clip after normalisation
|
||||
for v in w_new.iter_mut() {
|
||||
*v = v.min(max_weight);
|
||||
}
|
||||
let sum2: f64 = w_new.iter().sum();
|
||||
if sum2 > 1e-15 {
|
||||
for v in w_new.iter_mut() {
|
||||
*v /= sum2;
|
||||
}
|
||||
}
|
||||
|
||||
let diff: f64 = w
|
||||
.iter()
|
||||
.zip(w_new.iter())
|
||||
.map(|(a, b)| (a - b).powi(2))
|
||||
.sum::<f64>()
|
||||
.sqrt();
|
||||
w = w_new;
|
||||
|
||||
if diff < tol {
|
||||
let (ret, var) = portfolio_stats(&w, mu, cov);
|
||||
return Ok(PortfolioResult {
|
||||
weights: w,
|
||||
utility: ret - 0.5 * self.gamma * var,
|
||||
expected_return: ret,
|
||||
portfolio_variance: var,
|
||||
iterations: iter_count + 1,
|
||||
converged: true,
|
||||
});
|
||||
}
|
||||
}
|
||||
|
||||
let (ret, var) = portfolio_stats(&w, mu, cov);
|
||||
Ok(PortfolioResult {
|
||||
weights: w,
|
||||
utility: ret - 0.5 * self.gamma * var,
|
||||
expected_return: ret,
|
||||
portfolio_variance: var,
|
||||
iterations: max_iter,
|
||||
converged: false,
|
||||
})
|
||||
}
|
||||
}
|
||||
|
||||
/// Compute portfolio expected return and variance.
|
||||
pub fn portfolio_stats(w: &[f64], mu: &[f64], cov: &[Vec<f64>]) -> (f64, f64) {
|
||||
let n = w.len();
|
||||
let ret: f64 = (0..n).map(|i| w[i] * mu[i]).sum();
|
||||
let var: f64 = (0..n)
|
||||
.flat_map(|i| (0..n).map(move |j| w[i] * w[j] * cov[i][j]))
|
||||
.sum();
|
||||
(ret, var)
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests {
|
||||
use super::*;
|
||||
|
||||
#[test]
|
||||
fn test_cara_utility_basic() {
|
||||
let u = CARAUtility::new(2.0).unwrap();
|
||||
assert!((u.utility(0.0) - (-0.5)).abs() < 1e-10);
|
||||
assert!((u.marginal_utility(0.0) - 1.0).abs() < 1e-10);
|
||||
assert!((u.risk_aversion(42.0) - 2.0).abs() < 1e-10);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn test_crra_utility_log() {
|
||||
let u = CRRAUtility::new(1.0).unwrap();
|
||||
let val = u.utility(std::f64::consts::E);
|
||||
assert!((val - 1.0).abs() < 1e-10);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn test_cara_optimizer_equal() {
|
||||
// Equal means, no covariance → equal weights
|
||||
let mu = vec![0.01, 0.01, 0.01];
|
||||
let cov = vec![
|
||||
vec![0.04, 0.0, 0.0],
|
||||
vec![0.0, 0.04, 0.0],
|
||||
vec![0.0, 0.0, 0.04],
|
||||
];
|
||||
let opt = CARAOptimizer::new(2.0).unwrap();
|
||||
let res = opt.optimize(&mu, &cov, 0.5).unwrap();
|
||||
assert!(res.converged);
|
||||
for w in &res.weights {
|
||||
assert!((w - 1.0 / 3.0).abs() < 0.05);
|
||||
}
|
||||
}
|
||||
}
|
||||
@@ -0,0 +1,200 @@
|
||||
//! Convex optimisation via projected gradient descent.
|
||||
//!
|
||||
//! Implements constrained optimisation over the simplex with box constraints.
|
||||
//! Supports generic `ConvexObjective` + `ConvexConstraint` traits.
|
||||
|
||||
use super::traits::ConvexObjective;
|
||||
use crate::core::OptimizrError;
|
||||
|
||||
/// Result of convex optimisation.
|
||||
#[derive(Debug, Clone)]
|
||||
pub struct ConvexResult {
|
||||
pub x: Vec<f64>,
|
||||
pub objective_value: f64,
|
||||
pub iterations: usize,
|
||||
pub converged: bool,
|
||||
pub gradient_norm: f64,
|
||||
}
|
||||
|
||||
/// Projected gradient descent solver for convex problems on the simplex.
|
||||
///
|
||||
/// Solves: $\min f(w)$ subject to $\sum w_i = 1$, $l \le w_i \le u$.
|
||||
pub struct ProjectedGradientSolver {
|
||||
pub max_iter: usize,
|
||||
pub learning_rate: f64,
|
||||
pub tolerance: f64,
|
||||
pub box_lower: f64,
|
||||
pub box_upper: f64,
|
||||
}
|
||||
|
||||
impl Default for ProjectedGradientSolver {
|
||||
fn default() -> Self {
|
||||
Self {
|
||||
max_iter: 2000,
|
||||
learning_rate: 0.1,
|
||||
tolerance: 1e-8,
|
||||
box_lower: 0.0,
|
||||
box_upper: 1.0,
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
impl ProjectedGradientSolver {
|
||||
pub fn new(max_iter: usize, lr: f64, tol: f64, lower: f64, upper: f64) -> Self {
|
||||
Self {
|
||||
max_iter,
|
||||
learning_rate: lr,
|
||||
tolerance: tol,
|
||||
box_lower: lower,
|
||||
box_upper: upper,
|
||||
}
|
||||
}
|
||||
|
||||
/// Solve min f(x) subject to Σx_i = 1, lower ≤ x_i ≤ upper.
|
||||
pub fn solve(&self, objective: &dyn ConvexObjective) -> Result<ConvexResult, OptimizrError> {
|
||||
let n = objective.dim();
|
||||
if n == 0 {
|
||||
return Err(OptimizrError::EmptyData);
|
||||
}
|
||||
|
||||
let mut x = vec![1.0 / n as f64; n];
|
||||
let mut best_val = f64::INFINITY;
|
||||
let mut best_x = x.clone();
|
||||
|
||||
for iter in 0..self.max_iter {
|
||||
let grad = objective.gradient(&x);
|
||||
let grad_norm: f64 = grad.iter().map(|g| g * g).sum::<f64>().sqrt();
|
||||
let current_val = objective.value(&x);
|
||||
|
||||
if current_val < best_val {
|
||||
best_val = current_val;
|
||||
best_x = x.clone();
|
||||
}
|
||||
|
||||
// Gradient step
|
||||
let mut x_new: Vec<f64> = x
|
||||
.iter()
|
||||
.zip(grad.iter())
|
||||
.map(|(xi, gi)| xi - self.learning_rate * gi)
|
||||
.collect();
|
||||
|
||||
// Project onto box
|
||||
for xi in x_new.iter_mut() {
|
||||
*xi = xi.max(self.box_lower).min(self.box_upper);
|
||||
}
|
||||
|
||||
// Project onto simplex
|
||||
let sum: f64 = x_new.iter().sum();
|
||||
if sum > 1e-15 {
|
||||
for xi in x_new.iter_mut() {
|
||||
*xi /= sum;
|
||||
}
|
||||
}
|
||||
|
||||
// Re-clip after normalisation
|
||||
for xi in x_new.iter_mut() {
|
||||
*xi = xi.max(self.box_lower).min(self.box_upper);
|
||||
}
|
||||
let sum2: f64 = x_new.iter().sum();
|
||||
if sum2 > 1e-15 {
|
||||
for xi in x_new.iter_mut() {
|
||||
*xi /= sum2;
|
||||
}
|
||||
}
|
||||
|
||||
let diff: f64 = x
|
||||
.iter()
|
||||
.zip(x_new.iter())
|
||||
.map(|(a, b)| (a - b).powi(2))
|
||||
.sum::<f64>()
|
||||
.sqrt();
|
||||
x = x_new;
|
||||
|
||||
if diff < self.tolerance {
|
||||
return Ok(ConvexResult {
|
||||
x: best_x.clone(),
|
||||
objective_value: best_val,
|
||||
iterations: iter + 1,
|
||||
converged: true,
|
||||
gradient_norm: grad_norm,
|
||||
});
|
||||
}
|
||||
}
|
||||
|
||||
Ok(ConvexResult {
|
||||
x: best_x,
|
||||
objective_value: best_val,
|
||||
iterations: self.max_iter,
|
||||
converged: false,
|
||||
gradient_norm: 0.0,
|
||||
})
|
||||
}
|
||||
}
|
||||
|
||||
/// Mean-variance objective: min γ/2 w'Σw - w'μ (+ optional score tilting).
|
||||
pub struct MeanVarianceObjective {
|
||||
pub mu: Vec<f64>,
|
||||
pub cov: Vec<Vec<f64>>,
|
||||
pub gamma: f64,
|
||||
pub score_weights: Option<Vec<f64>>,
|
||||
}
|
||||
|
||||
impl ConvexObjective for MeanVarianceObjective {
|
||||
fn value(&self, w: &[f64]) -> f64 {
|
||||
let n = w.len();
|
||||
let ret: f64 = (0..n).map(|i| w[i] * self.mu[i]).sum();
|
||||
let var: f64 = (0..n)
|
||||
.flat_map(|i| (0..n).map(move |j| w[i] * w[j] * self.cov[i][j]))
|
||||
.sum();
|
||||
let mut val = 0.5 * self.gamma * var - ret;
|
||||
if let Some(ref scores) = self.score_weights {
|
||||
let bonus: f64 = (0..n.min(scores.len())).map(|i| w[i] * scores[i]).sum();
|
||||
val -= 0.1 * bonus;
|
||||
}
|
||||
val
|
||||
}
|
||||
|
||||
fn gradient(&self, w: &[f64]) -> Vec<f64> {
|
||||
let n = w.len();
|
||||
let mut grad = vec![0.0; n];
|
||||
for i in 0..n {
|
||||
grad[i] = -self.mu[i];
|
||||
for j in 0..n {
|
||||
grad[i] += self.gamma * self.cov[i][j] * w[j];
|
||||
}
|
||||
if let Some(ref scores) = self.score_weights {
|
||||
if i < scores.len() {
|
||||
grad[i] -= 0.1 * scores[i];
|
||||
}
|
||||
}
|
||||
}
|
||||
grad
|
||||
}
|
||||
|
||||
fn dim(&self) -> usize {
|
||||
self.mu.len()
|
||||
}
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests {
|
||||
use super::*;
|
||||
|
||||
#[test]
|
||||
fn test_solver_diagonal_cov() {
|
||||
let obj = MeanVarianceObjective {
|
||||
mu: vec![0.05, 0.03],
|
||||
cov: vec![vec![0.04, 0.0], vec![0.0, 0.01]],
|
||||
gamma: 2.0,
|
||||
score_weights: None,
|
||||
};
|
||||
let solver = ProjectedGradientSolver {
|
||||
box_upper: 0.8,
|
||||
..Default::default()
|
||||
};
|
||||
let res = solver.solve(&obj).unwrap();
|
||||
assert!(res.converged);
|
||||
let total: f64 = res.x.iter().sum();
|
||||
assert!((total - 1.0).abs() < 0.01);
|
||||
}
|
||||
}
|
||||
@@ -0,0 +1,201 @@
|
||||
//! Classic mean-variance (Markowitz) portfolio optimisation + ERC.
|
||||
|
||||
use super::cara::portfolio_stats;
|
||||
use super::convex::{MeanVarianceObjective, ProjectedGradientSolver};
|
||||
use super::traits::{PortfolioOptimizer, PortfolioResult};
|
||||
use crate::core::OptimizrError;
|
||||
|
||||
/// Markowitz mean-variance optimizer with optional score tilting.
|
||||
pub struct MeanVarianceOptimizer {
|
||||
pub risk_aversion: f64,
|
||||
pub scores: Option<Vec<f64>>,
|
||||
}
|
||||
|
||||
impl MeanVarianceOptimizer {
|
||||
pub fn new(risk_aversion: f64) -> Self {
|
||||
Self {
|
||||
risk_aversion,
|
||||
scores: None,
|
||||
}
|
||||
}
|
||||
|
||||
pub fn with_scores(mut self, scores: Vec<f64>) -> Self {
|
||||
self.scores = Some(scores);
|
||||
self
|
||||
}
|
||||
}
|
||||
|
||||
impl PortfolioOptimizer for MeanVarianceOptimizer {
|
||||
fn optimize(
|
||||
&self,
|
||||
mu: &[f64],
|
||||
cov: &[Vec<f64>],
|
||||
max_weight: f64,
|
||||
) -> Result<PortfolioResult, OptimizrError> {
|
||||
let n = mu.len();
|
||||
if n == 0 {
|
||||
return Err(OptimizrError::EmptyData);
|
||||
}
|
||||
if cov.len() != n {
|
||||
return Err(OptimizrError::DimensionMismatch {
|
||||
expected: n,
|
||||
actual: cov.len(),
|
||||
});
|
||||
}
|
||||
|
||||
let obj = MeanVarianceObjective {
|
||||
mu: mu.to_vec(),
|
||||
cov: cov.to_vec(),
|
||||
gamma: self.risk_aversion,
|
||||
score_weights: self.scores.clone(),
|
||||
};
|
||||
|
||||
let solver = ProjectedGradientSolver {
|
||||
box_upper: max_weight,
|
||||
..Default::default()
|
||||
};
|
||||
|
||||
let res = solver.solve(&obj)?;
|
||||
let (ret, var) = portfolio_stats(&res.x, mu, cov);
|
||||
|
||||
Ok(PortfolioResult {
|
||||
weights: res.x,
|
||||
utility: ret - 0.5 * self.risk_aversion * var,
|
||||
expected_return: ret,
|
||||
portfolio_variance: var,
|
||||
iterations: res.iterations,
|
||||
converged: res.converged,
|
||||
})
|
||||
}
|
||||
}
|
||||
|
||||
/// Minimum variance portfolio (γ → ∞, ignores expected returns).
|
||||
pub fn minimum_variance(
|
||||
cov: &[Vec<f64>],
|
||||
max_weight: f64,
|
||||
) -> Result<PortfolioResult, OptimizrError> {
|
||||
let n = cov.len();
|
||||
if n == 0 {
|
||||
return Err(OptimizrError::EmptyData);
|
||||
}
|
||||
let mu = vec![0.0; n];
|
||||
let opt = MeanVarianceOptimizer {
|
||||
risk_aversion: 100.0,
|
||||
scores: None,
|
||||
};
|
||||
opt.optimize(&mu, cov, max_weight)
|
||||
}
|
||||
|
||||
/// Equal-Risk-Contribution (ERC / Risk Parity) portfolio.
|
||||
///
|
||||
/// Iterates: $w_i \propto 1 / (\Sigma w)_i$.
|
||||
pub fn equal_risk_contribution(
|
||||
cov: &[Vec<f64>],
|
||||
max_weight: f64,
|
||||
) -> Result<PortfolioResult, OptimizrError> {
|
||||
let n = cov.len();
|
||||
if n == 0 {
|
||||
return Err(OptimizrError::EmptyData);
|
||||
}
|
||||
|
||||
let mut w = vec![1.0 / n as f64; n];
|
||||
let max_iter = 500;
|
||||
|
||||
for _ in 0..max_iter {
|
||||
// Marginal risk contribution: (Σw)_i
|
||||
let mut mrc = vec![0.0; n];
|
||||
for i in 0..n {
|
||||
for j in 0..n {
|
||||
mrc[i] += cov[i][j] * w[j];
|
||||
}
|
||||
}
|
||||
|
||||
// New weights ∝ 1/|mrc_i|
|
||||
let mut w_new: Vec<f64> = mrc
|
||||
.iter()
|
||||
.map(|m| {
|
||||
if m.abs() > 1e-15 {
|
||||
1.0 / m.abs()
|
||||
} else {
|
||||
1.0
|
||||
}
|
||||
})
|
||||
.collect();
|
||||
|
||||
// Normalise
|
||||
let sum: f64 = w_new.iter().sum();
|
||||
for v in w_new.iter_mut() {
|
||||
*v /= sum;
|
||||
}
|
||||
// Clip
|
||||
for v in w_new.iter_mut() {
|
||||
*v = v.min(max_weight);
|
||||
}
|
||||
let sum2: f64 = w_new.iter().sum();
|
||||
for v in w_new.iter_mut() {
|
||||
*v /= sum2;
|
||||
}
|
||||
|
||||
let diff: f64 = w
|
||||
.iter()
|
||||
.zip(w_new.iter())
|
||||
.map(|(a, b)| (a - b).powi(2))
|
||||
.sum::<f64>()
|
||||
.sqrt();
|
||||
w = w_new;
|
||||
if diff < 1e-10 {
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
let var: f64 = {
|
||||
let w_ref = &w;
|
||||
(0..n)
|
||||
.flat_map(|i| (0..n).map(move |j| (i, j)))
|
||||
.map(|(i, j)| w_ref[i] * w_ref[j] * cov[i][j])
|
||||
.sum()
|
||||
};
|
||||
|
||||
Ok(PortfolioResult {
|
||||
weights: w,
|
||||
utility: -var,
|
||||
expected_return: 0.0,
|
||||
portfolio_variance: var,
|
||||
iterations: max_iter,
|
||||
converged: true,
|
||||
})
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests {
|
||||
use super::*;
|
||||
|
||||
#[test]
|
||||
fn test_mean_variance_identity_cov() {
|
||||
let mu = vec![0.10, 0.05, 0.08];
|
||||
let cov = vec![
|
||||
vec![0.04, 0.0, 0.0],
|
||||
vec![0.0, 0.04, 0.0],
|
||||
vec![0.0, 0.0, 0.04],
|
||||
];
|
||||
let opt = MeanVarianceOptimizer::new(2.0);
|
||||
let res = opt.optimize(&mu, &cov, 0.5).unwrap();
|
||||
// Highest mu (0.10) should get the largest weight
|
||||
assert!(res.weights[0] > res.weights[1]);
|
||||
assert!(res.converged);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn test_erc_diagonal() {
|
||||
let cov = vec![
|
||||
vec![0.04, 0.0, 0.0],
|
||||
vec![0.0, 0.04, 0.0],
|
||||
vec![0.0, 0.0, 0.04],
|
||||
];
|
||||
let res = equal_risk_contribution(&cov, 0.5).unwrap();
|
||||
// Identical variances → equal weights
|
||||
for w in &res.weights {
|
||||
assert!((w - 1.0 / 3.0).abs() < 0.01);
|
||||
}
|
||||
}
|
||||
}
|
||||
@@ -0,0 +1,22 @@
|
||||
//! Portfolio Optimization — CARA utility, convex duality, mean-variance.
|
||||
//!
|
||||
//! Provides generic trait-based abstractions for:
|
||||
//! - **CARA** (Constant Absolute Risk Aversion) utility maximisation
|
||||
//! - **CRRA** (Constant Relative Risk Aversion) utility
|
||||
//! - **Convex duality** optimization (Legendre-Fenchel transform)
|
||||
//! - **Mean-variance** portfolio weights (Markowitz)
|
||||
//!
|
||||
//! All heavy lifting uses `ndarray` + `rayon` for parallelism.
|
||||
//!
|
||||
//! # References
|
||||
//! - Markowitz (1952) — Portfolio Selection
|
||||
//! - Merton (1969) — Lifetime Portfolio Selection under Uncertainty
|
||||
//! - Rockafellar (1970) — Convex Analysis (duality)
|
||||
|
||||
pub mod traits;
|
||||
pub mod cara;
|
||||
pub mod convex;
|
||||
pub mod mean_variance;
|
||||
|
||||
#[cfg(feature = "python-bindings")]
|
||||
pub mod python_bindings;
|
||||
@@ -0,0 +1,108 @@
|
||||
//! Python bindings for portfolio optimisation.
|
||||
//!
|
||||
//! Exposes CARA, mean-variance, minimum-variance, and ERC optimisers.
|
||||
|
||||
use pyo3::exceptions::PyValueError;
|
||||
use pyo3::prelude::*;
|
||||
|
||||
use super::cara::CARAOptimizer;
|
||||
use super::mean_variance::{equal_risk_contribution, minimum_variance, MeanVarianceOptimizer};
|
||||
use super::traits::PortfolioOptimizer;
|
||||
|
||||
/// CARA optimal weights: maximize w'μ − (γ/2)w'Σw s.t. Σw=1, 0≤w≤max.
|
||||
#[pyfunction]
|
||||
#[pyo3(signature = (mu, cov, risk_aversion=2.0, max_weight=0.3))]
|
||||
fn cara_optimal_weights(
|
||||
py: Python<'_>,
|
||||
mu: Vec<f64>,
|
||||
cov: Vec<Vec<f64>>,
|
||||
risk_aversion: f64,
|
||||
max_weight: f64,
|
||||
) -> PyResult<PyObject> {
|
||||
let opt =
|
||||
CARAOptimizer::new(risk_aversion).map_err(|e| PyValueError::new_err(format!("{}", e)))?;
|
||||
let result = opt
|
||||
.optimize(&mu, &cov, max_weight)
|
||||
.map_err(|e| PyValueError::new_err(format!("{}", e)))?;
|
||||
|
||||
let dict = pyo3::types::PyDict::new_bound(py);
|
||||
dict.set_item("weights", result.weights.clone())?;
|
||||
dict.set_item("utility", result.utility)?;
|
||||
dict.set_item("expected_return", result.expected_return)?;
|
||||
dict.set_item("portfolio_variance", result.portfolio_variance)?;
|
||||
dict.set_item("sharpe_ratio", result.sharpe_ratio(0.05))?;
|
||||
dict.set_item("iterations", result.iterations)?;
|
||||
dict.set_item("converged", result.converged)?;
|
||||
Ok(dict.into())
|
||||
}
|
||||
|
||||
/// Mean-variance optimal weights with optional score tilting.
|
||||
#[pyfunction]
|
||||
#[pyo3(signature = (mu, cov, risk_aversion=2.0, max_weight=0.3, scores=None))]
|
||||
fn mean_variance_optimal_weights(
|
||||
py: Python<'_>,
|
||||
mu: Vec<f64>,
|
||||
cov: Vec<Vec<f64>>,
|
||||
risk_aversion: f64,
|
||||
max_weight: f64,
|
||||
scores: Option<Vec<f64>>,
|
||||
) -> PyResult<PyObject> {
|
||||
let mut opt = MeanVarianceOptimizer::new(risk_aversion);
|
||||
if let Some(s) = scores {
|
||||
opt = opt.with_scores(s);
|
||||
}
|
||||
let result = opt
|
||||
.optimize(&mu, &cov, max_weight)
|
||||
.map_err(|e| PyValueError::new_err(format!("{}", e)))?;
|
||||
|
||||
let dict = pyo3::types::PyDict::new_bound(py);
|
||||
dict.set_item("weights", result.weights.clone())?;
|
||||
dict.set_item("utility", result.utility)?;
|
||||
dict.set_item("expected_return", result.expected_return)?;
|
||||
dict.set_item("portfolio_variance", result.portfolio_variance)?;
|
||||
dict.set_item("sharpe_ratio", result.sharpe_ratio(0.05))?;
|
||||
dict.set_item("iterations", result.iterations)?;
|
||||
dict.set_item("converged", result.converged)?;
|
||||
Ok(dict.into())
|
||||
}
|
||||
|
||||
/// Minimum variance portfolio.
|
||||
#[pyfunction]
|
||||
#[pyo3(signature = (cov, max_weight=0.3))]
|
||||
fn min_variance_weights(
|
||||
py: Python<'_>,
|
||||
cov: Vec<Vec<f64>>,
|
||||
max_weight: f64,
|
||||
) -> PyResult<PyObject> {
|
||||
let result =
|
||||
minimum_variance(&cov, max_weight).map_err(|e| PyValueError::new_err(format!("{}", e)))?;
|
||||
let dict = pyo3::types::PyDict::new_bound(py);
|
||||
dict.set_item("weights", result.weights.clone())?;
|
||||
dict.set_item("portfolio_variance", result.portfolio_variance)?;
|
||||
dict.set_item("iterations", result.iterations)?;
|
||||
dict.set_item("converged", result.converged)?;
|
||||
Ok(dict.into())
|
||||
}
|
||||
|
||||
/// Equal-risk-contribution (Risk Parity) portfolio.
|
||||
#[pyfunction]
|
||||
#[pyo3(signature = (cov, max_weight=0.3))]
|
||||
fn erc_weights(py: Python<'_>, cov: Vec<Vec<f64>>, max_weight: f64) -> PyResult<PyObject> {
|
||||
let result = equal_risk_contribution(&cov, max_weight)
|
||||
.map_err(|e| PyValueError::new_err(format!("{}", e)))?;
|
||||
let dict = pyo3::types::PyDict::new_bound(py);
|
||||
dict.set_item("weights", result.weights.clone())?;
|
||||
dict.set_item("portfolio_variance", result.portfolio_variance)?;
|
||||
dict.set_item("iterations", result.iterations)?;
|
||||
dict.set_item("converged", result.converged)?;
|
||||
Ok(dict.into())
|
||||
}
|
||||
|
||||
/// Register all portfolio optimization functions with the Python module.
|
||||
pub fn register_python_functions(m: &Bound<'_, pyo3::types::PyModule>) -> PyResult<()> {
|
||||
m.add_function(wrap_pyfunction!(cara_optimal_weights, m)?)?;
|
||||
m.add_function(wrap_pyfunction!(mean_variance_optimal_weights, m)?)?;
|
||||
m.add_function(wrap_pyfunction!(min_variance_weights, m)?)?;
|
||||
m.add_function(wrap_pyfunction!(erc_weights, m)?)?;
|
||||
Ok(())
|
||||
}
|
||||
@@ -0,0 +1,75 @@
|
||||
//! Core traits for portfolio utility and convex optimisation.
|
||||
|
||||
use crate::core::OptimizrError;
|
||||
|
||||
/// A utility function U(x) mapping wealth → utility.
|
||||
pub trait UtilityFunction: Send + Sync {
|
||||
/// U(x) — the utility of wealth level x.
|
||||
fn utility(&self, x: f64) -> f64;
|
||||
|
||||
/// U'(x) — first derivative (marginal utility).
|
||||
fn marginal_utility(&self, x: f64) -> f64;
|
||||
|
||||
/// (U')^{-1}(y) — inverse marginal utility (used in duality).
|
||||
fn inverse_marginal(&self, y: f64) -> f64;
|
||||
|
||||
/// Risk-aversion coefficient A(x) = -U''(x) / U'(x).
|
||||
fn risk_aversion(&self, x: f64) -> f64;
|
||||
|
||||
/// Name identifier for logging / serialisation.
|
||||
fn name(&self) -> &str;
|
||||
}
|
||||
|
||||
/// Convex objective f(w) over portfolio weights w ∈ ℝ^n.
|
||||
///
|
||||
/// Used by convex solvers (projected gradient, ADMM, etc.).
|
||||
pub trait ConvexObjective: Send + Sync {
|
||||
/// f(w) — objective value.
|
||||
fn value(&self, w: &[f64]) -> f64;
|
||||
|
||||
/// ∇f(w) — gradient vector.
|
||||
fn gradient(&self, w: &[f64]) -> Vec<f64>;
|
||||
|
||||
/// Dimension of the weight vector.
|
||||
fn dim(&self) -> usize;
|
||||
}
|
||||
|
||||
/// Convex constraint g(w) ≤ 0.
|
||||
pub trait ConvexConstraint: Send + Sync {
|
||||
/// g(w) — constraint value (feasible when ≤ 0).
|
||||
fn value(&self, w: &[f64]) -> f64;
|
||||
|
||||
/// ∇g(w) — gradient of constraint function.
|
||||
fn gradient(&self, w: &[f64]) -> Vec<f64>;
|
||||
}
|
||||
|
||||
/// Result of a portfolio optimisation.
|
||||
#[derive(Debug, Clone)]
|
||||
pub struct PortfolioResult {
|
||||
pub weights: Vec<f64>,
|
||||
pub utility: f64,
|
||||
pub expected_return: f64,
|
||||
pub portfolio_variance: f64,
|
||||
pub iterations: usize,
|
||||
pub converged: bool,
|
||||
}
|
||||
|
||||
impl PortfolioResult {
|
||||
pub fn sharpe_ratio(&self, risk_free: f64) -> f64 {
|
||||
let vol = self.portfolio_variance.sqrt();
|
||||
if vol < 1e-15 {
|
||||
return 0.0;
|
||||
}
|
||||
(self.expected_return - risk_free) / vol
|
||||
}
|
||||
}
|
||||
|
||||
/// Generic portfolio optimiser trait.
|
||||
pub trait PortfolioOptimizer: Send + Sync {
|
||||
fn optimize(
|
||||
&self,
|
||||
mu: &[f64],
|
||||
cov: &[Vec<f64>],
|
||||
max_weight: f64,
|
||||
) -> Result<PortfolioResult, OptimizrError>;
|
||||
}
|
||||
Reference in New Issue
Block a user