#![allow(dead_code)]
use scirs2_core::ndarray::{Array1, Array2};
use std::collections::HashMap;
use std::f64::consts::PI;
pub struct PortfolioOptimizer {
returns: Array1<f64>,
covariance: Array2<f64>,
risk_aversion: f64,
constraints: PortfolioConstraints,
method: OptimizationMethod,
}
#[derive(Debug, Clone, Default)]
pub struct PortfolioConstraints {
min_investment: Option<Array1<f64>>,
max_investment: Option<Array1<f64>>,
target_return: Option<f64>,
max_risk: Option<f64>,
sector_constraints: Option<SectorConstraints>,
cardinality: Option<usize>,
transaction_costs: Option<Array1<f64>>,
current_portfolio: Option<Array1<f64>>,
}
#[derive(Debug, Clone)]
pub struct SectorConstraints {
asset_sectors: Vec<usize>,
min_sector_allocation: Vec<f64>,
max_sector_allocation: Vec<f64>,
}
#[derive(Debug, Clone)]
pub enum OptimizationMethod {
MeanVariance,
CVaR { confidence_level: f64 },
MaxSharpe { risk_free_rate: f64 },
RiskParity,
BlackLitterman { views: MarketViews },
KellyCriterion,
}
#[derive(Debug, Clone)]
pub struct MarketViews {
view_matrix: Array2<f64>,
view_expectations: Array1<f64>,
view_confidence: Array2<f64>,
}
impl PortfolioOptimizer {
pub fn new(
returns: Array1<f64>,
covariance: Array2<f64>,
risk_aversion: f64,
) -> Result<Self, String> {
if returns.len() != covariance.shape()[0] {
return Err("Returns and covariance dimensions mismatch".to_string());
}
Ok(Self {
returns,
covariance,
risk_aversion,
constraints: PortfolioConstraints::default(),
method: OptimizationMethod::MeanVariance,
})
}
pub fn with_method(mut self, method: OptimizationMethod) -> Self {
self.method = method;
self
}
pub fn with_constraints(mut self, constraints: PortfolioConstraints) -> Self {
self.constraints = constraints;
self
}
pub fn build_qubo(
&self,
num_bits_per_asset: usize,
) -> Result<(Array2<f64>, HashMap<String, usize>), String> {
let n_assets = self.returns.len();
let total_vars = n_assets * num_bits_per_asset;
let mut qubo = Array2::zeros((total_vars, total_vars));
let mut var_map = HashMap::new();
for i in 0..n_assets {
for j in 0..num_bits_per_asset {
let var_name = format!("asset_{i}_bit_{j}");
var_map.insert(var_name, i * num_bits_per_asset + j);
}
}
match &self.method {
OptimizationMethod::MeanVariance => {
self.add_mean_variance_objective(&mut qubo, num_bits_per_asset)?;
}
OptimizationMethod::CVaR { confidence_level } => {
self.add_cvar_objective(&mut qubo, num_bits_per_asset, *confidence_level)?;
}
OptimizationMethod::MaxSharpe { risk_free_rate } => {
self.add_sharpe_objective(&mut qubo, num_bits_per_asset, *risk_free_rate)?;
}
_ => {
return Err("Optimization method not yet implemented".to_string());
}
}
self.add_constraints(&mut qubo, num_bits_per_asset)?;
Ok((qubo, var_map))
}
fn add_mean_variance_objective(
&self,
qubo: &mut Array2<f64>,
bits_per_asset: usize,
) -> Result<(), String> {
let n_assets = self.returns.len();
for i in 0..n_assets {
for k in 0..bits_per_asset {
let weight = (1 << k) as f64 / ((1 << bits_per_asset) - 1) as f64;
let var_idx = i * bits_per_asset + k;
qubo[[var_idx, var_idx]] -= self.returns[i] * weight;
}
}
for i in 0..n_assets {
for j in 0..n_assets {
for k1 in 0..bits_per_asset {
for k2 in 0..bits_per_asset {
let w1 = (1 << k1) as f64 / ((1 << bits_per_asset) - 1) as f64;
let w2 = (1 << k2) as f64 / ((1 << bits_per_asset) - 1) as f64;
let var_idx1 = i * bits_per_asset + k1;
let var_idx2 = j * bits_per_asset + k2;
if var_idx1 <= var_idx2 {
let coef = self.risk_aversion * self.covariance[[i, j]] * w1 * w2;
if var_idx1 == var_idx2 {
qubo[[var_idx1, var_idx1]] += coef;
} else {
qubo[[var_idx1, var_idx2]] += coef / 2.0;
qubo[[var_idx2, var_idx1]] += coef / 2.0;
}
}
}
}
}
}
Ok(())
}
fn add_cvar_objective(
&self,
qubo: &mut Array2<f64>,
bits_per_asset: usize,
confidence_level: f64,
) -> Result<(), String> {
self.add_mean_variance_objective(qubo, bits_per_asset)?;
let tail_weight = 1.0 / (1.0 - confidence_level);
for i in 0..qubo.shape()[0] {
for j in 0..qubo.shape()[1] {
qubo[[i, j]] *= tail_weight;
}
}
Ok(())
}
fn add_sharpe_objective(
&self,
qubo: &mut Array2<f64>,
bits_per_asset: usize,
risk_free_rate: f64,
) -> Result<(), String> {
let adjusted_returns = &self.returns - risk_free_rate;
let n_assets = adjusted_returns.len();
for i in 0..n_assets {
for k in 0..bits_per_asset {
let weight = (1 << k) as f64 / ((1 << bits_per_asset) - 1) as f64;
let var_idx = i * bits_per_asset + k;
let mean_var = self
.covariance
.diag()
.mean()
.ok_or_else(|| "Empty covariance diagonal".to_string())?;
let vol_scale = 1.0 / mean_var.sqrt();
qubo[[var_idx, var_idx]] -= adjusted_returns[i] * weight * vol_scale;
}
}
self.add_mean_variance_objective(qubo, bits_per_asset)?;
Ok(())
}
fn add_constraints(&self, qubo: &mut Array2<f64>, bits_per_asset: usize) -> Result<(), String> {
let penalty = 100.0;
self.add_budget_constraint(qubo, bits_per_asset, penalty)?;
if let Some(max_assets) = self.constraints.cardinality {
self.add_cardinality_constraint(qubo, bits_per_asset, max_assets, penalty)?;
}
if let Some(target) = self.constraints.target_return {
self.add_return_constraint(qubo, bits_per_asset, target, penalty)?;
}
if let Some(max_risk) = self.constraints.max_risk {
self.add_risk_constraint(qubo, bits_per_asset, max_risk, penalty)?;
}
Ok(())
}
fn add_budget_constraint(
&self,
qubo: &mut Array2<f64>,
bits_per_asset: usize,
penalty: f64,
) -> Result<(), String> {
let n_assets = self.returns.len();
for i in 0..n_assets {
for j in 0..n_assets {
for k1 in 0..bits_per_asset {
for k2 in 0..bits_per_asset {
let w1 = (1 << k1) as f64 / ((1 << bits_per_asset) - 1) as f64;
let w2 = (1 << k2) as f64 / ((1 << bits_per_asset) - 1) as f64;
let var_idx1 = i * bits_per_asset + k1;
let var_idx2 = j * bits_per_asset + k2;
qubo[[var_idx1, var_idx2]] += penalty * w1 * w2;
}
}
}
}
for i in 0..n_assets {
for k in 0..bits_per_asset {
let weight = (1 << k) as f64 / ((1 << bits_per_asset) - 1) as f64;
let var_idx = i * bits_per_asset + k;
qubo[[var_idx, var_idx]] -= 2.0 * penalty * weight;
}
}
Ok(())
}
fn add_cardinality_constraint(
&self,
qubo: &mut Array2<f64>,
bits_per_asset: usize,
_max_assets: usize,
penalty: f64,
) -> Result<(), String> {
let n_assets = self.returns.len();
for i in 0..n_assets {
for j in i + 1..n_assets {
for k1 in 0..bits_per_asset {
for k2 in 0..bits_per_asset {
let var_idx1 = i * bits_per_asset + k1;
let var_idx2 = j * bits_per_asset + k2;
if i != j {
qubo[[var_idx1, var_idx2]] += penalty / (n_assets * n_assets) as f64;
}
}
}
}
}
Ok(())
}
fn add_return_constraint(
&self,
qubo: &mut Array2<f64>,
bits_per_asset: usize,
target_return: f64,
penalty: f64,
) -> Result<(), String> {
let n_assets = self.returns.len();
for i in 0..n_assets {
for j in 0..n_assets {
for k1 in 0..bits_per_asset {
for k2 in 0..bits_per_asset {
let w1 = (1 << k1) as f64 / ((1 << bits_per_asset) - 1) as f64;
let w2 = (1 << k2) as f64 / ((1 << bits_per_asset) - 1) as f64;
let var_idx1 = i * bits_per_asset + k1;
let var_idx2 = j * bits_per_asset + k2;
qubo[[var_idx1, var_idx2]] +=
penalty * self.returns[i] * self.returns[j] * w1 * w2;
}
}
}
}
for i in 0..n_assets {
for k in 0..bits_per_asset {
let weight = (1 << k) as f64 / ((1 << bits_per_asset) - 1) as f64;
let var_idx = i * bits_per_asset + k;
qubo[[var_idx, var_idx]] -=
2.0 * penalty * target_return * self.returns[i] * weight;
}
}
Ok(())
}
fn add_risk_constraint(
&self,
qubo: &mut Array2<f64>,
bits_per_asset: usize,
max_risk: f64,
penalty: f64,
) -> Result<(), String> {
let risk_scale = penalty / max_risk;
let n_assets = self.returns.len();
for i in 0..n_assets {
for j in 0..n_assets {
for k1 in 0..bits_per_asset {
for k2 in 0..bits_per_asset {
let w1 = (1 << k1) as f64 / ((1 << bits_per_asset) - 1) as f64;
let w2 = (1 << k2) as f64 / ((1 << bits_per_asset) - 1) as f64;
let var_idx1 = i * bits_per_asset + k1;
let var_idx2 = j * bits_per_asset + k2;
qubo[[var_idx1, var_idx2]] +=
risk_scale * self.covariance[[i, j]] * w1 * w2;
}
}
}
}
Ok(())
}
pub fn decode_solution(
&self,
solution: &HashMap<String, bool>,
bits_per_asset: usize,
) -> Array1<f64> {
let n_assets = self.returns.len();
let mut weights = Array1::zeros(n_assets);
for i in 0..n_assets {
let mut weight = 0.0;
for j in 0..bits_per_asset {
let var_name = format!("asset_{i}_bit_{j}");
if *solution.get(&var_name).unwrap_or(&false) {
weight += (1 << j) as f64 / ((1 << bits_per_asset) - 1) as f64;
}
}
weights[i] = weight;
}
let sum: f64 = weights.sum();
if sum > 0.0 {
weights /= sum;
}
weights
}
pub fn calculate_metrics(&self, weights: &Array1<f64>) -> PortfolioMetrics {
let expected_return = weights.dot(&self.returns);
let variance = weights.dot(&self.covariance.dot(weights));
let volatility = variance.sqrt();
let sharpe_ratio = if volatility > 0.0 {
expected_return / volatility
} else {
0.0
};
let weighted_vol: f64 = weights
.iter()
.zip(self.covariance.diag().iter())
.map(|(&w, &var)| w * var.sqrt())
.sum();
let diversification_ratio = if volatility > 0.0 {
weighted_vol / volatility
} else {
1.0
};
PortfolioMetrics {
expected_return,
volatility,
sharpe_ratio,
diversification_ratio,
max_drawdown: self.estimate_max_drawdown(weights),
value_at_risk: self.calculate_var(weights, 0.95),
conditional_value_at_risk: self.calculate_cvar(weights, 0.95),
}
}
fn estimate_max_drawdown(&self, weights: &Array1<f64>) -> f64 {
let variance = weights.dot(&self.covariance.dot(weights));
let volatility = variance.sqrt();
2.0 * volatility
}
fn calculate_var(&self, weights: &Array1<f64>, _confidence: f64) -> f64 {
let expected_return = weights.dot(&self.returns);
let variance = weights.dot(&self.covariance.dot(weights));
let volatility = variance.sqrt();
let z_score = 1.645; expected_return - z_score * volatility
}
fn calculate_cvar(&self, weights: &Array1<f64>, confidence: f64) -> f64 {
let var = self.calculate_var(weights, confidence);
let volatility = weights.dot(&self.covariance.dot(weights)).sqrt();
let z_score = 1.645;
let phi = (-z_score * z_score / 2.0_f64).exp() / (2.0_f64 * PI).sqrt();
var - volatility * phi / (1.0 - confidence)
}
}
#[derive(Debug, Clone)]
pub struct PortfolioMetrics {
pub expected_return: f64,
pub volatility: f64,
pub sharpe_ratio: f64,
pub diversification_ratio: f64,
pub max_drawdown: f64,
pub value_at_risk: f64,
pub conditional_value_at_risk: f64,
}
pub struct RiskParityOptimizer {
covariance: Array2<f64>,
target_contributions: Option<Array1<f64>>,
tolerance: f64,
max_iterations: usize,
}
impl RiskParityOptimizer {
pub const fn new(covariance: Array2<f64>) -> Self {
Self {
covariance,
target_contributions: None,
tolerance: 1e-6,
max_iterations: 1000,
}
}
pub fn with_target_contributions(mut self, targets: Array1<f64>) -> Self {
self.target_contributions = Some(targets);
self
}
pub fn build_qubo(
&self,
num_bits_per_asset: usize,
) -> Result<(Array2<f64>, HashMap<String, usize>), String> {
let n_assets = self.covariance.shape()[0];
let total_vars = n_assets * num_bits_per_asset;
let mut qubo = Array2::zeros((total_vars, total_vars));
let mut var_map = HashMap::new();
for i in 0..n_assets {
for j in 0..num_bits_per_asset {
let var_name = format!("rp_asset_{i}_bit_{j}");
var_map.insert(var_name, i * num_bits_per_asset + j);
}
}
let initial_weights = Array1::from_elem(n_assets, 1.0 / n_assets as f64);
let _risk_contribs = self.calculate_risk_contributions(&initial_weights);
let default_targets = Array1::from_elem(n_assets, 1.0 / n_assets as f64);
let _targets = self
.target_contributions
.as_ref()
.unwrap_or(&default_targets);
for i in 0..n_assets {
for k in 0..num_bits_per_asset {
let weight = (1 << k) as f64 / ((1 << num_bits_per_asset) - 1) as f64;
let var_idx = i * num_bits_per_asset + k;
let gradient = self.risk_contribution_gradient(&initial_weights, i);
qubo[[var_idx, var_idx]] += gradient * weight;
}
}
self.add_budget_constraint_rp(&mut qubo, num_bits_per_asset, 100.0)?;
Ok((qubo, var_map))
}
fn calculate_risk_contributions(&self, weights: &Array1<f64>) -> Array1<f64> {
let portfolio_risk = weights.dot(&self.covariance.dot(weights)).sqrt();
let marginal_risks = self.covariance.dot(weights);
weights * &marginal_risks / portfolio_risk
}
fn risk_contribution_gradient(&self, weights: &Array1<f64>, asset: usize) -> f64 {
let eps = 1e-6;
let mut weights_plus = weights.clone();
weights_plus[asset] += eps;
let rc_plus = self.calculate_risk_contributions(&weights_plus)[asset];
let rc = self.calculate_risk_contributions(weights)[asset];
(rc_plus - rc) / eps
}
fn add_budget_constraint_rp(
&self,
qubo: &mut Array2<f64>,
bits_per_asset: usize,
penalty: f64,
) -> Result<(), String> {
let n_assets = self.covariance.shape()[0];
for i in 0..n_assets {
for j in 0..n_assets {
for k1 in 0..bits_per_asset {
for k2 in 0..bits_per_asset {
let w1 = (1 << k1) as f64 / ((1 << bits_per_asset) - 1) as f64;
let w2 = (1 << k2) as f64 / ((1 << bits_per_asset) - 1) as f64;
let var_idx1 = i * bits_per_asset + k1;
let var_idx2 = j * bits_per_asset + k2;
qubo[[var_idx1, var_idx2]] += penalty * w1 * w2;
}
}
}
}
Ok(())
}
}
pub struct BlackLittermanOptimizer {
equilibrium_returns: Array1<f64>,
covariance: Array2<f64>,
risk_aversion: f64,
tau: f64,
views: MarketViews,
}
impl BlackLittermanOptimizer {
pub fn new(market_weights: Array1<f64>, covariance: Array2<f64>, risk_aversion: f64) -> Self {
let equilibrium_returns = risk_aversion * covariance.dot(&market_weights);
Self {
equilibrium_returns,
covariance,
risk_aversion,
tau: 0.05,
views: MarketViews {
view_matrix: Array2::zeros((0, 0)),
view_expectations: Array1::zeros(0),
view_confidence: Array2::zeros((0, 0)),
},
}
}
pub fn with_views(mut self, views: MarketViews) -> Self {
self.views = views;
self
}
pub fn posterior_returns(&self) -> Array1<f64> {
let _tau_sigma = self.tau * &self.covariance;
if self.views.view_matrix.shape()[0] == 0 {
return self.equilibrium_returns.clone();
}
let p = &self.views.view_matrix;
let q = &self.views.view_expectations;
let omega = &self.views.view_confidence;
let mut posterior = self.equilibrium_returns.clone();
for i in 0..p.shape()[0] {
let view_assets: Vec<_> = (0..p.shape()[1])
.filter(|&j| p[[i, j]].abs() > 1e-10)
.collect();
if !view_assets.is_empty() {
let view_return = q[i];
let confidence = 1.0 / omega[[i, i]];
for &asset in &view_assets {
posterior[asset] +=
confidence * (view_return - posterior[asset]) * p[[i, asset]];
}
}
}
posterior
}
pub fn build_qubo(
&self,
num_bits_per_asset: usize,
) -> Result<(Array2<f64>, HashMap<String, usize>), String> {
let posterior_returns = self.posterior_returns();
let optimizer = PortfolioOptimizer::new(
posterior_returns,
self.covariance.clone(),
self.risk_aversion,
)?;
optimizer.build_qubo(num_bits_per_asset)
}
}
pub struct TransactionCostOptimizer {
current_weights: Array1<f64>,
target_optimizer: PortfolioOptimizer,
cost_model: TransactionCostModel,
threshold: f64,
}
#[derive(Debug, Clone)]
pub enum TransactionCostModel {
Linear { rates: Array1<f64> },
Quadratic { impact_coefficients: Array1<f64> },
FixedPlusProportional {
fixed: f64,
proportional: Array1<f64>,
},
Custom,
}
impl TransactionCostOptimizer {
pub const fn new(
current_weights: Array1<f64>,
target_optimizer: PortfolioOptimizer,
cost_model: TransactionCostModel,
) -> Self {
Self {
current_weights,
target_optimizer,
cost_model,
threshold: 0.01,
}
}
pub fn build_qubo(
&self,
num_bits_per_asset: usize,
) -> Result<(Array2<f64>, HashMap<String, usize>), String> {
let (mut qubo, var_map) = self.target_optimizer.build_qubo(num_bits_per_asset)?;
self.add_transaction_costs(&mut qubo, &var_map, num_bits_per_asset)?;
Ok((qubo, var_map))
}
fn add_transaction_costs(
&self,
qubo: &mut Array2<f64>,
var_map: &HashMap<String, usize>,
bits_per_asset: usize,
) -> Result<(), String> {
let n_assets = self.current_weights.len();
if let TransactionCostModel::Linear { rates } = &self.cost_model {
for i in 0..n_assets {
let current_w = self.current_weights[i];
let rate = rates[i];
for k in 0..bits_per_asset {
let weight = (1 << k) as f64 / ((1 << bits_per_asset) - 1) as f64;
let var_name = format!("asset_{i}_bit_{k}");
if let Some(&var_idx) = var_map.get(&var_name) {
qubo[[var_idx, var_idx]] += rate * (weight - current_w).powi(2);
}
}
}
} else {
}
Ok(())
}
}
pub struct RebalancingScheduler {
frequency: RebalancingFrequency,
triggers: Vec<RebalancingTrigger>,
cost_threshold: f64,
}
#[derive(Debug, Clone)]
pub enum RebalancingFrequency {
Fixed { days: usize },
Calendar { schedule: String },
Adaptive,
}
#[derive(Debug, Clone)]
pub enum RebalancingTrigger {
Deviation { threshold: f64 },
VolatilitySpike { threshold: f64 },
CorrelationBreakdown { threshold: f64 },
Custom { name: String },
}
impl RebalancingScheduler {
pub const fn new(frequency: RebalancingFrequency) -> Self {
Self {
frequency,
triggers: Vec::new(),
cost_threshold: 0.001,
}
}
pub fn add_trigger(mut self, trigger: RebalancingTrigger) -> Self {
self.triggers.push(trigger);
self
}
pub fn should_rebalance(
&self,
current_weights: &Array1<f64>,
target_weights: &Array1<f64>,
market_data: &MarketData,
) -> bool {
for trigger in &self.triggers {
match trigger {
RebalancingTrigger::Deviation { threshold } => {
let max_deviation = current_weights
.iter()
.zip(target_weights.iter())
.map(|(&c, &t)| (c - t).abs())
.fold(0.0, f64::max);
if max_deviation > *threshold {
return true;
}
}
RebalancingTrigger::VolatilitySpike { threshold } => {
if market_data.volatility > *threshold {
return true;
}
}
_ => {}
}
}
false
}
}
#[derive(Debug, Clone)]
pub struct MarketData {
pub returns: Array1<f64>,
pub volatility: f64,
pub correlation: Array2<f64>,
pub volume: Array1<f64>,
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn test_portfolio_optimizer() {
let returns = Array1::from_vec(vec![0.10, 0.12, 0.08]);
let covariance = Array2::from_shape_vec(
(3, 3),
vec![0.01, 0.002, 0.001, 0.002, 0.015, 0.002, 0.001, 0.002, 0.008],
)
.expect("Failed to create covariance matrix from shape vector");
let optimizer = PortfolioOptimizer::new(returns, covariance, 2.0)
.expect("Failed to create portfolio optimizer");
let (qubo, var_map) = optimizer.build_qubo(3).expect("Failed to build QUBO");
assert_eq!(var_map.len(), 9); }
#[test]
fn test_portfolio_metrics() {
let returns = Array1::from_vec(vec![0.10, 0.12, 0.08]);
let covariance = Array2::eye(3) * 0.01;
let optimizer = PortfolioOptimizer::new(returns, covariance, 2.0)
.expect("Failed to create portfolio optimizer");
let mut weights = Array1::from_vec(vec![0.3, 0.5, 0.2]);
let mut metrics = optimizer.calculate_metrics(&weights);
assert!((metrics.expected_return - 0.106).abs() < 0.001);
assert!(metrics.volatility > 0.0);
assert!(metrics.sharpe_ratio > 0.0);
}
#[test]
fn test_risk_parity() {
let covariance = Array2::from_shape_vec(
(3, 3),
vec![0.01, 0.002, 0.001, 0.002, 0.015, 0.002, 0.001, 0.002, 0.008],
)
.expect("Failed to create covariance matrix from shape vector");
let rp_optimizer = RiskParityOptimizer::new(covariance);
let (qubo, var_map) = rp_optimizer.build_qubo(3).expect("Failed to build QUBO");
assert!(!var_map.is_empty());
}
}