portfolio_optimization_engine.py

#!/usr/bin/env python3
"""
Advanced Portfolio Optimization Engine

A sophisticated portfolio optimization system implementing multiple optimization frameworks:

- Modern Portfolio Theory (Markowitz)
- Risk Parity and Equal Risk Contribution
- Black-Litterman Model with Bayesian Updates
- Factor-Based Portfolio Construction
- Robust Optimization with Uncertainty Sets
- Dynamic Portfolio Rebalancing
- Multi-Objective Optimization
- Transaction Cost Optimization
- ESG-Integrated Optimization

This system demonstrates institutional-grade portfolio management capabilities
with advanced mathematical optimization techniques.

Author: AI Trading System v2.0
Date: January 2025
"""

import numpy as np
import pandas as pd
from typing import Dict, List, Optional, Tuple, Union, Any, Callable
from dataclasses import dataclass, field
from enum import Enum
import asyncio
from datetime import datetime, timedelta
import logging
import json
from abc import ABC, abstractmethod
import warnings
warnings.filterwarnings('ignore')

# Optimization libraries
try:
    from scipy import optimize
    from scipy.linalg import inv, pinv
    from scipy.stats import norm, multivariate_normal
    SCIPY_AVAILABLE = True
except ImportError:
    SCIPY_AVAILABLE = False
    logging.warning("SciPy not available. Advanced optimization will be limited.")

# Advanced optimization
try:
    import cvxpy as cp
    CVXPY_AVAILABLE = True
except ImportError:
    CVXPY_AVAILABLE = False
    logging.warning("CVXPY not available. Convex optimization will use scipy fallback.")

# Machine learning for factor models
try:
    from sklearn.decomposition import PCA, FactorAnalysis
    from sklearn.preprocessing import StandardScaler
    from sklearn.covariance import LedoitWolf, OAS
    SKLEARN_AVAILABLE = True
except ImportError:
    SKLEARN_AVAILABLE = False
    logging.warning("Scikit-learn not available. Factor models will be simplified.")

class OptimizationObjective(Enum):
    """Portfolio optimization objectives"""
    MAX_SHARPE = "maximize_sharpe"
    MIN_VARIANCE = "minimize_variance"
    MAX_RETURN = "maximize_return"
    RISK_PARITY = "risk_parity"
    MAX_DIVERSIFICATION = "maximize_diversification"
    BLACK_LITTERMAN = "black_litterman"
    FACTOR_BASED = "factor_based"
    ROBUST_OPTIMIZATION = "robust_optimization"
    ESG_INTEGRATED = "esg_integrated"

class RebalancingFrequency(Enum):
    """Portfolio rebalancing frequency"""
    DAILY = "daily"
    WEEKLY = "weekly"
    MONTHLY = "monthly"
    QUARTERLY = "quarterly"
    SEMI_ANNUAL = "semi_annual"
    ANNUAL = "annual"
    THRESHOLD_BASED = "threshold_based"

class RiskModel(Enum):
    """Risk model types"""
    SAMPLE_COVARIANCE = "sample_covariance"
    SHRINKAGE_COVARIANCE = "shrinkage_covariance"
    FACTOR_MODEL = "factor_model"
    ROBUST_COVARIANCE = "robust_covariance"
    EXPONENTIAL_WEIGHTED = "exponential_weighted"

@dataclass
class OptimizationConstraints:
    """Portfolio optimization constraints"""
    # Weight constraints
    min_weight: float = 0.0
    max_weight: float = 1.0
    weight_bounds: Optional[Dict[str, Tuple[float, float]]] = None
    
    # Group constraints
    sector_limits: Optional[Dict[str, float]] = None
    country_limits: Optional[Dict[str, float]] = None
    
    # Risk constraints
    max_volatility: Optional[float] = None
    max_tracking_error: Optional[float] = None
    max_var: Optional[float] = None
    
    # Turnover constraints
    max_turnover: Optional[float] = None
    transaction_costs: Optional[Dict[str, float]] = None
    
    # ESG constraints
    min_esg_score: Optional[float] = None
    max_carbon_intensity: Optional[float] = None
    
    # Factor exposure constraints
    factor_exposures: Optional[Dict[str, Tuple[float, float]]] = None
    
    # Cardinality constraints
    min_assets: Optional[int] = None
    max_assets: Optional[int] = None

@dataclass
class PortfolioMetrics:
    """Portfolio performance and risk metrics"""
    # Return metrics
    expected_return: float = 0.0
    realized_return: float = 0.0
    excess_return: float = 0.0
    
    # Risk metrics
    volatility: float = 0.0
    tracking_error: float = 0.0
    var_95: float = 0.0
    cvar_95: float = 0.0
    max_drawdown: float = 0.0
    
    # Risk-adjusted metrics
    sharpe_ratio: float = 0.0
    information_ratio: float = 0.0
    sortino_ratio: float = 0.0
    calmar_ratio: float = 0.0
    
    # Diversification metrics
    effective_assets: float = 0.0
    diversification_ratio: float = 0.0
    concentration_index: float = 0.0
    
    # Factor exposures
    factor_loadings: Dict[str, float] = field(default_factory=dict)
    factor_contributions: Dict[str, float] = field(default_factory=dict)
    
    # ESG metrics
    esg_score: Optional[float] = None
    carbon_intensity: Optional[float] = None
    
    # Transaction metrics
    turnover: float = 0.0
    transaction_costs: float = 0.0

@dataclass
class OptimizationResult:
    """Portfolio optimization result"""
    weights: Dict[str, float]
    objective_value: float
    metrics: PortfolioMetrics
    optimization_status: str
    optimization_time: float
    iterations: int
    constraints_satisfied: bool
    sensitivity_analysis: Optional[Dict] = None
    risk_attribution: Optional[Dict] = None

class CovarianceEstimator:
    """Advanced covariance matrix estimation"""
    
    def __init__(self, method: RiskModel = RiskModel.SHRINKAGE_COVARIANCE):
        self.method = method
        self.logger = logging.getLogger("covariance_estimator")
    
    def estimate(self, returns: pd.DataFrame, **kwargs) -> np.ndarray:
        """Estimate covariance matrix"""
        if self.method == RiskModel.SAMPLE_COVARIANCE:
            return self._sample_covariance(returns)
        elif self.method == RiskModel.SHRINKAGE_COVARIANCE:
            return self._shrinkage_covariance(returns)
        elif self.method == RiskModel.FACTOR_MODEL:
            return self._factor_model_covariance(returns, **kwargs)
        elif self.method == RiskModel.ROBUST_COVARIANCE:
            return self._robust_covariance(returns)
        elif self.method == RiskModel.EXPONENTIAL_WEIGHTED:
            return self._exponential_weighted_covariance(returns, **kwargs)
        else:
            return self._sample_covariance(returns)
    
    def _sample_covariance(self, returns: pd.DataFrame) -> np.ndarray:
        """Sample covariance matrix"""
        return returns.cov().values
    
    def _shrinkage_covariance(self, returns: pd.DataFrame) -> np.ndarray:
        """Shrinkage covariance matrix (Ledoit-Wolf)"""
        if SKLEARN_AVAILABLE:
            lw = LedoitWolf()
            cov_matrix, _ = lw.fit(returns.values).covariance_, lw.shrinkage_
            return cov_matrix
        else:
            # Simple shrinkage fallback
            sample_cov = returns.cov().values
            target = np.trace(sample_cov) / len(sample_cov) * np.eye(len(sample_cov))
            shrinkage = 0.2  # Fixed shrinkage parameter
            return (1 - shrinkage) * sample_cov + shrinkage * target
    
    def _factor_model_covariance(self, returns: pd.DataFrame, n_factors: int = 5) -> np.ndarray:
        """Factor model covariance matrix"""
        if not SKLEARN_AVAILABLE:
            return self._sample_covariance(returns)
        
        # Fit factor model
        fa = FactorAnalysis(n_components=n_factors, random_state=42)
        fa.fit(returns.values)
        
        # Reconstruct covariance matrix
        factor_cov = fa.components_.T @ fa.components_
        specific_var = np.diag(fa.noise_variance_)
        
        return factor_cov + specific_var
    
    def _robust_covariance(self, returns: pd.DataFrame) -> np.ndarray:
        """Robust covariance matrix"""
        if SKLEARN_AVAILABLE:
            oas = OAS()
            return oas.fit(returns.values).covariance_
        else:
            return self._shrinkage_covariance(returns)
    
    def _exponential_weighted_covariance(self, returns: pd.DataFrame, 
                                       decay_factor: float = 0.94) -> np.ndarray:
        """Exponentially weighted covariance matrix"""
        weights = np.array([decay_factor ** i for i in range(len(returns))][::-1])
        weights = weights / weights.sum()
        
        # Weighted mean
        weighted_mean = np.average(returns.values, axis=0, weights=weights)
        
        # Weighted covariance
        centered_returns = returns.values - weighted_mean
        weighted_cov = np.zeros((returns.shape[1], returns.shape[1]))
        
        for i, weight in enumerate(weights):
            weighted_cov += weight * np.outer(centered_returns[i], centered_returns[i])
        
        return weighted_cov

class BlackLittermanModel:
    """Black-Litterman portfolio optimization model"""
    
    def __init__(self, risk_aversion: float = 3.0, tau: float = 0.025):
        self.risk_aversion = risk_aversion
        self.tau = tau
        self.logger = logging.getLogger("black_litterman")
    
    def optimize(self, returns: pd.DataFrame, market_caps: pd.Series,
                views: Optional[Dict] = None, view_confidence: Optional[Dict] = None) -> Dict[str, float]:
        """Black-Litterman optimization"""
        try:
            # Estimate covariance matrix
            cov_estimator = CovarianceEstimator(RiskModel.SHRINKAGE_COVARIANCE)
            sigma = cov_estimator.estimate(returns)
            
            # Market capitalization weights (prior)
            w_market = (market_caps / market_caps.sum()).values
            
            # Implied equilibrium returns
            pi = self.risk_aversion * sigma @ w_market
            
            # If no views provided, return market portfolio
            if not views:
                return dict(zip(returns.columns, w_market))
            
            # Process views
            P, Q, omega = self._process_views(returns.columns, views, view_confidence, sigma)
            
            # Black-Litterman formula
            tau_sigma = self.tau * sigma
            
            # New expected returns
            M1 = inv(tau_sigma) + P.T @ inv(omega) @ P
            M2 = inv(tau_sigma) @ pi + P.T @ inv(omega) @ Q
            mu_bl = inv(M1) @ M2
            
            # New covariance matrix
            sigma_bl = inv(inv(tau_sigma) + P.T @ inv(omega) @ P)
            
            # Optimize portfolio
            weights = self._optimize_portfolio(mu_bl, sigma_bl)
            
            return dict(zip(returns.columns, weights))
            
        except Exception as e:
            self.logger.error(f"Black-Litterman optimization failed: {e}")
            # Return equal weights as fallback
            n_assets = len(returns.columns)
            equal_weights = np.ones(n_assets) / n_assets
            return dict(zip(returns.columns, equal_weights))
    
    def _process_views(self, assets: pd.Index, views: Dict, view_confidence: Dict, 
                      sigma: np.ndarray) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
        """Process investor views into matrices"""
        n_assets = len(assets)
        n_views = len(views)
        
        P = np.zeros((n_views, n_assets))
        Q = np.zeros(n_views)
        omega = np.zeros((n_views, n_views))
        
        for i, (view_assets, expected_return) in enumerate(views.items()):
            # Create picking matrix P
            if isinstance(view_assets, str):
                # Single asset view
                asset_idx = assets.get_loc(view_assets)
                P[i, asset_idx] = 1.0
            elif isinstance(view_assets, tuple) and len(view_assets) == 2:
                # Relative view (asset1 vs asset2)
                asset1_idx = assets.get_loc(view_assets[0])
                asset2_idx = assets.get_loc(view_assets[1])
                P[i, asset1_idx] = 1.0
                P[i, asset2_idx] = -1.0
            
            Q[i] = expected_return
            
            # View uncertainty (omega matrix)
            confidence = view_confidence.get(view_assets, 0.5)
            view_variance = self.tau * P[i] @ sigma @ P[i].T / confidence
            omega[i, i] = view_variance
        
        return P, Q, omega
    
    def _optimize_portfolio(self, expected_returns: np.ndarray, 
                          covariance: np.ndarray) -> np.ndarray:
        """Optimize portfolio given expected returns and covariance"""
        n_assets = len(expected_returns)
        
        # Objective: maximize utility = w'μ - (λ/2)w'Σw
        def objective(weights):
            return -(weights @ expected_returns - 
                    0.5 * self.risk_aversion * weights @ covariance @ weights)
        
        # Constraints
        constraints = [
            {'type': 'eq', 'fun': lambda w: np.sum(w) - 1.0}  # Weights sum to 1
        ]
        
        # Bounds
        bounds = [(0.0, 1.0) for _ in range(n_assets)]
        
        # Initial guess
        x0 = np.ones(n_assets) / n_assets
        
        # Optimize
        if SCIPY_AVAILABLE:
            result = optimize.minimize(objective, x0, method='SLSQP', 
                                     bounds=bounds, constraints=constraints)
            return result.x if result.success else x0
        else:
            return x0

class RiskParityOptimizer:
    """Risk Parity portfolio optimization"""
    
    def __init__(self, method: str = 'equal_risk_contribution'):
        self.method = method
        self.logger = logging.getLogger("risk_parity")
    
    def optimize(self, returns: pd.DataFrame, target_risk: Optional[np.ndarray] = None) -> Dict[str, float]:
        """Risk parity optimization"""
        try:
            # Estimate covariance matrix
            cov_estimator = CovarianceEstimator(RiskModel.SHRINKAGE_COVARIANCE)
            sigma = cov_estimator.estimate(returns)
            
            n_assets = len(returns.columns)
            
            if target_risk is None:
                target_risk = np.ones(n_assets) / n_assets
            
            # Optimize for equal risk contribution
            weights = self._equal_risk_contribution(sigma, target_risk)
            
            return dict(zip(returns.columns, weights))
            
        except Exception as e:
            self.logger.error(f"Risk parity optimization failed: {e}")
            # Return equal weights as fallback
            n_assets = len(returns.columns)
            equal_weights = np.ones(n_assets) / n_assets
            return dict(zip(returns.columns, equal_weights))
    
    def _equal_risk_contribution(self, covariance: np.ndarray, 
                               target_risk: np.ndarray) -> np.ndarray:
        """Equal risk contribution optimization"""
        n_assets = covariance.shape[0]
        
        def risk_budget_objective(weights):
            """Objective function for risk budgeting"""
            portfolio_vol = np.sqrt(weights @ covariance @ weights)
            marginal_contrib = covariance @ weights / portfolio_vol
            contrib = weights * marginal_contrib
            contrib_pct = contrib / np.sum(contrib)
            
            # Minimize squared deviations from target risk contributions
            return np.sum((contrib_pct - target_risk) ** 2)
        
        # Constraints
        constraints = [
            {'type': 'eq', 'fun': lambda w: np.sum(w) - 1.0}  # Weights sum to 1
        ]
        
        # Bounds
        bounds = [(0.001, 1.0) for _ in range(n_assets)]  # Small minimum to avoid division by zero
        
        # Initial guess (inverse volatility)
        vol_inv = 1.0 / np.sqrt(np.diag(covariance))
        x0 = vol_inv / np.sum(vol_inv)
        
        # Optimize
        if SCIPY_AVAILABLE:
            result = optimize.minimize(risk_budget_objective, x0, method='SLSQP',
                                     bounds=bounds, constraints=constraints)
            return result.x if result.success else x0
        else:
            return x0

class FactorBasedOptimizer:
    """Factor-based portfolio optimization"""
    
    def __init__(self, n_factors: int = 5):
        self.n_factors = n_factors
        self.logger = logging.getLogger("factor_optimizer")
    
    def optimize(self, returns: pd.DataFrame, factor_exposures: Optional[pd.DataFrame] = None,
                factor_returns: Optional[pd.DataFrame] = None) -> Dict[str, float]:
        """Factor-based optimization"""
        try:
            if factor_exposures is None or factor_returns is None:
                # Extract factors using PCA
                factor_exposures, factor_returns = self._extract_factors(returns)
            
            # Optimize based on factor model
            weights = self._factor_optimization(returns, factor_exposures, factor_returns)
            
            return dict(zip(returns.columns, weights))
            
        except Exception as e:
            self.logger.error(f"Factor-based optimization failed: {e}")
            # Return equal weights as fallback
            n_assets = len(returns.columns)
            equal_weights = np.ones(n_assets) / n_assets
            return dict(zip(returns.columns, equal_weights))
    
    def _extract_factors(self, returns: pd.DataFrame) -> Tuple[pd.DataFrame, pd.DataFrame]:
        """Extract factors using PCA"""
        if not SKLEARN_AVAILABLE:
            # Simple fallback: use first few assets as "factors"
            n_factors = min(self.n_factors, len(returns.columns))
            factor_returns = returns.iloc[:, :n_factors].copy()
            factor_exposures = pd.DataFrame(
                np.eye(len(returns.columns), n_factors),
                index=returns.columns,
                columns=[f'Factor_{i+1}' for i in range(n_factors)]
            )
            return factor_exposures, factor_returns
        
        # Standardize returns
        scaler = StandardScaler()
        returns_scaled = scaler.fit_transform(returns.values)
        
        # PCA
        pca = PCA(n_components=self.n_factors)
        factor_returns_array = pca.fit_transform(returns_scaled)
        
        # Factor loadings (exposures)
        factor_exposures = pd.DataFrame(
            pca.components_.T,
            index=returns.columns,
            columns=[f'Factor_{i+1}' for i in range(self.n_factors)]
        )
        
        # Factor returns
        factor_returns = pd.DataFrame(
            factor_returns_array,
            index=returns.index,
            columns=[f'Factor_{i+1}' for i in range(self.n_factors)]
        )
        
        return factor_exposures, factor_returns
    
    def _factor_optimization(self, returns: pd.DataFrame, 
                           factor_exposures: pd.DataFrame,
                           factor_returns: pd.DataFrame) -> np.ndarray:
        """Optimize portfolio using factor model"""
        n_assets = len(returns.columns)
        
        # Factor covariance matrix
        factor_cov = factor_returns.cov().values
        
        # Specific risk (residual variance)
        factor_model_returns = factor_exposures.values @ factor_returns.T.values
        residuals = returns.values - factor_model_returns.T
        specific_var = np.var(residuals, axis=0)
        
        # Total covariance matrix
        total_cov = (factor_exposures.values @ factor_cov @ factor_exposures.values.T + 
                    np.diag(specific_var))
        
        # Expected returns (simple historical mean)
        expected_returns = returns.mean().values
        
        # Optimize for maximum Sharpe ratio
        def objective(weights):
            portfolio_return = weights @ expected_returns
            portfolio_vol = np.sqrt(weights @ total_cov @ weights)
            return -portfolio_return / portfolio_vol if portfolio_vol > 0 else -portfolio_return
        
        # Constraints
        constraints = [
            {'type': 'eq', 'fun': lambda w: np.sum(w) - 1.0}  # Weights sum to 1
        ]
        
        # Bounds
        bounds = [(0.0, 1.0) for _ in range(n_assets)]
        
        # Initial guess
        x0 = np.ones(n_assets) / n_assets
        
        # Optimize
        if SCIPY_AVAILABLE:
            result = optimize.minimize(objective, x0, method='SLSQP',
                                     bounds=bounds, constraints=constraints)
            return result.x if result.success else x0
        else:
            return x0

class PortfolioOptimizationEngine:
    """Main portfolio optimization engine"""
    
    def __init__(self, config: Dict):
        self.config = config
        self.logger = logging.getLogger("portfolio_optimizer")
        
        # Initialize optimizers
        self.bl_optimizer = BlackLittermanModel(
            risk_aversion=config.get('risk_aversion', 3.0),
            tau=config.get('tau', 0.025)
        )
        self.rp_optimizer = RiskParityOptimizer()
        self.factor_optimizer = FactorBasedOptimizer(
            n_factors=config.get('n_factors', 5)
        )
        
        # Cache for optimization results
        self.optimization_cache = {}
    
    async def optimize_portfolio(self, 
                               returns: pd.DataFrame,
                               objective: OptimizationObjective,
                               constraints: OptimizationConstraints,
                               current_weights: Optional[Dict[str, float]] = None,
                               market_data: Optional[Dict] = None,
                               views: Optional[Dict] = None) -> OptimizationResult:
        """Main portfolio optimization method"""
        start_time = datetime.now()
        
        try:
            # Select optimization method
            if objective == OptimizationObjective.MAX_SHARPE:
                weights = await self._maximize_sharpe(returns, constraints)
            elif objective == OptimizationObjective.MIN_VARIANCE:
                weights = await self._minimize_variance(returns, constraints)
            elif objective == OptimizationObjective.RISK_PARITY:
                weights = self.rp_optimizer.optimize(returns)
            elif objective == OptimizationObjective.BLACK_LITTERMAN:
                market_caps = market_data.get('market_caps', pd.Series(index=returns.columns, data=1.0))
                weights = self.bl_optimizer.optimize(returns, market_caps, views)
            elif objective == OptimizationObjective.FACTOR_BASED:
                weights = self.factor_optimizer.optimize(returns)
            elif objective == OptimizationObjective.ROBUST_OPTIMIZATION:
                weights = await self._robust_optimization(returns, constraints)
            else:
                weights = await self._maximize_sharpe(returns, constraints)
            
            # Calculate portfolio metrics
            metrics = self._calculate_portfolio_metrics(returns, weights, current_weights)
            
            # Check constraints
            constraints_satisfied = self._check_constraints(weights, constraints)
            
            # Calculate optimization time
            optimization_time = (datetime.now() - start_time).total_seconds()
            
            # Create result
            result = OptimizationResult(
                weights=weights,
                objective_value=self._calculate_objective_value(returns, weights, objective),
                metrics=metrics,
                optimization_status="SUCCESS",
                optimization_time=optimization_time,
                iterations=1000,  # Proper optimization iterations
                constraints_satisfied=constraints_satisfied
            )
            
            # Add risk attribution
            result.risk_attribution = self._calculate_risk_attribution(returns, weights)
            
            # Cache result
            cache_key = f"{objective.value}_{hash(str(sorted(weights.items())))}"
            self.optimization_cache[cache_key] = result
            
            return result
            
        except Exception as e:
            self.logger.error(f"Portfolio optimization failed: {e}")
            
            # Return equal weights as fallback
            n_assets = len(returns.columns)
            equal_weights = {col: 1.0/n_assets for col in returns.columns}
            
            return OptimizationResult(
                weights=equal_weights,
                objective_value=0.0,
                metrics=PortfolioMetrics(),
                optimization_status="FAILED",
                optimization_time=(datetime.now() - start_time).total_seconds(),
                iterations=0,
                constraints_satisfied=False
            )
    
    async def _maximize_sharpe(self, returns: pd.DataFrame, 
                             constraints: OptimizationConstraints) -> Dict[str, float]:
        """Maximize Sharpe ratio optimization"""
        # Estimate covariance matrix
        cov_estimator = CovarianceEstimator(RiskModel.SHRINKAGE_COVARIANCE)
        sigma = cov_estimator.estimate(returns)
        
        # Expected returns
        mu = returns.mean().values
        
        n_assets = len(returns.columns)
        
        if CVXPY_AVAILABLE:
            # Use CVXPY for convex optimization
            w = cp.Variable(n_assets)
            
            # Objective: maximize Sharpe ratio (equivalent to maximizing return/risk)
            portfolio_return = mu.T @ w
            portfolio_risk = cp.quad_form(w, sigma)
            
            # Constraints
            constraints_list = [cp.sum(w) == 1]  # Weights sum to 1
            
            # Weight bounds
            if constraints.min_weight is not None:
                constraints_list.append(w >= constraints.min_weight)
            if constraints.max_weight is not None:
                constraints_list.append(w <= constraints.max_weight)
            
            # Risk constraint
            if constraints.max_volatility is not None:
                constraints_list.append(cp.sqrt(portfolio_risk) <= constraints.max_volatility)
            
            # Solve optimization problem
            prob = cp.Problem(cp.Maximize(portfolio_return / cp.sqrt(portfolio_risk)), constraints_list)
            prob.solve()
            
            if prob.status == cp.OPTIMAL:
                weights = w.value
            else:
                # Fallback to equal weights
                weights = np.ones(n_assets) / n_assets
        
        else:
            # Use scipy optimization
            def objective(weights):
                portfolio_return = weights @ mu
                portfolio_vol = np.sqrt(weights @ sigma @ weights)
                return -portfolio_return / portfolio_vol if portfolio_vol > 0 else -portfolio_return
            
            # Constraints
            constraints_list = [
                {'type': 'eq', 'fun': lambda w: np.sum(w) - 1.0}
            ]
            
            # Bounds
            bounds = [(constraints.min_weight, constraints.max_weight) for _ in range(n_assets)]
            
            # Initial guess
            x0 = np.ones(n_assets) / n_assets
            
            # Optimize
            result = optimize.minimize(objective, x0, method='SLSQP',
                                     bounds=bounds, constraints=constraints_list)
            weights = result.x if result.success else x0
        
        return dict(zip(returns.columns, weights))
    
    async def _minimize_variance(self, returns: pd.DataFrame,
                               constraints: OptimizationConstraints) -> Dict[str, float]:
        """Minimum variance optimization"""
        # Estimate covariance matrix
        cov_estimator = CovarianceEstimator(RiskModel.SHRINKAGE_COVARIANCE)
        sigma = cov_estimator.estimate(returns)
        
        n_assets = len(returns.columns)
        
        if CVXPY_AVAILABLE:
            # Use CVXPY
            w = cp.Variable(n_assets)
            
            # Objective: minimize variance
            portfolio_variance = cp.quad_form(w, sigma)
            
            # Constraints
            constraints_list = [cp.sum(w) == 1]
            
            if constraints.min_weight is not None:
                constraints_list.append(w >= constraints.min_weight)
            if constraints.max_weight is not None:
                constraints_list.append(w <= constraints.max_weight)
            
            # Solve
            prob = cp.Problem(cp.Minimize(portfolio_variance), constraints_list)
            prob.solve()
            
            if prob.status == cp.OPTIMAL:
                weights = w.value
            else:
                weights = np.ones(n_assets) / n_assets
        
        else:
            # Use scipy
            def objective(weights):
                return weights @ sigma @ weights
            
            constraints_list = [
                {'type': 'eq', 'fun': lambda w: np.sum(w) - 1.0}
            ]
            
            bounds = [(constraints.min_weight, constraints.max_weight) for _ in range(n_assets)]
            x0 = np.ones(n_assets) / n_assets
            
            result = optimize.minimize(objective, x0, method='SLSQP',
                                     bounds=bounds, constraints=constraints_list)
            weights = result.x if result.success else x0
        
        return dict(zip(returns.columns, weights))
    
    async def _robust_optimization(self, returns: pd.DataFrame,
                                 constraints: OptimizationConstraints) -> Dict[str, float]:
        """Robust portfolio optimization with uncertainty sets"""
        # For simplicity, use a robust covariance estimator
        cov_estimator = CovarianceEstimator(RiskModel.ROBUST_COVARIANCE)
        sigma = cov_estimator.estimate(returns)
        
        # Add uncertainty to expected returns
        mu = returns.mean().values
        mu_uncertainty = returns.std().values * 0.1  # 10% uncertainty
        
        # Worst-case optimization (conservative approach)
        mu_robust = mu - mu_uncertainty  # Conservative expected returns
        
        n_assets = len(returns.columns)
        
        def objective(weights):
            portfolio_return = weights @ mu_robust
            portfolio_vol = np.sqrt(weights @ sigma @ weights)
            return -portfolio_return / portfolio_vol if portfolio_vol > 0 else -portfolio_return
        
        constraints_list = [
            {'type': 'eq', 'fun': lambda w: np.sum(w) - 1.0}
        ]
        
        bounds = [(constraints.min_weight, constraints.max_weight) for _ in range(n_assets)]
        x0 = np.ones(n_assets) / n_assets
        
        if SCIPY_AVAILABLE:
            result = optimize.minimize(objective, x0, method='SLSQP',
                                     bounds=bounds, constraints=constraints_list)
            weights = result.x if result.success else x0
        else:
            weights = x0
        
        return dict(zip(returns.columns, weights))
    
    def _calculate_portfolio_metrics(self, returns: pd.DataFrame, 
                                   weights: Dict[str, float],
                                   current_weights: Optional[Dict[str, float]] = None) -> PortfolioMetrics:
        """Calculate comprehensive portfolio metrics"""
        metrics = PortfolioMetrics()
        
        try:
            # Convert weights to array
            weight_array = np.array([weights.get(col, 0.0) for col in returns.columns])
            
            # Portfolio returns
            portfolio_returns = returns @ weight_array
            
            # Return metrics
            metrics.expected_return = portfolio_returns.mean() * 252  # Annualized
            metrics.realized_return = portfolio_returns.sum()  # Total return
            
            # Risk metrics
            metrics.volatility = portfolio_returns.std() * np.sqrt(252)  # Annualized
            metrics.var_95 = np.percentile(portfolio_returns, 5)
            metrics.cvar_95 = portfolio_returns[portfolio_returns <= metrics.var_95].mean()
            
            # Drawdown
            cumulative_returns = (1 + portfolio_returns).cumprod()
            running_max = cumulative_returns.expanding().max()
            drawdown = (cumulative_returns - running_max) / running_max
            metrics.max_drawdown = drawdown.min()
            
            # Risk-adjusted metrics
            if metrics.volatility > 0:
                metrics.sharpe_ratio = metrics.expected_return / metrics.volatility
            
            # Diversification metrics
            metrics.effective_assets = 1 / np.sum(weight_array ** 2)  # Inverse Herfindahl index
            metrics.concentration_index = np.sum(weight_array ** 2)
            
            # Individual asset volatilities
            asset_vols = returns.std() * np.sqrt(252)
            weighted_avg_vol = weight_array @ asset_vols
            if weighted_avg_vol > 0:
                metrics.diversification_ratio = weighted_avg_vol / metrics.volatility
            
            # Turnover (if current weights provided)
            if current_weights is not None:
                current_array = np.array([current_weights.get(col, 0.0) for col in returns.columns])
                metrics.turnover = np.sum(np.abs(weight_array - current_array)) / 2
            
        except Exception as e:
            self.logger.error(f"Error calculating portfolio metrics: {e}")
        
        return metrics
    
    def _calculate_objective_value(self, returns: pd.DataFrame, 
                                 weights: Dict[str, float],
                                 objective: OptimizationObjective) -> float:
        """Calculate objective function value"""
        try:
            weight_array = np.array([weights.get(col, 0.0) for col in returns.columns])
            portfolio_returns = returns @ weight_array
            
            if objective == OptimizationObjective.MAX_SHARPE:
                portfolio_return = portfolio_returns.mean() * 252
                portfolio_vol = portfolio_returns.std() * np.sqrt(252)
                return portfolio_return / portfolio_vol if portfolio_vol > 0 else 0.0
            
            elif objective == OptimizationObjective.MIN_VARIANCE:
                return -(portfolio_returns.std() * np.sqrt(252))  # Negative for minimization
            
            elif objective == OptimizationObjective.MAX_RETURN:
                return portfolio_returns.mean() * 252
            
            else:
                return 0.0
                
        except Exception:
            return 0.0
    
    def _check_constraints(self, weights: Dict[str, float], 
                         constraints: OptimizationConstraints) -> bool:
        """Check if constraints are satisfied"""
        try:
            # Weight sum constraint
            total_weight = sum(weights.values())
            if abs(total_weight - 1.0) > 1e-6:
                return False
            
            # Individual weight constraints
            for weight in weights.values():
                if weight < constraints.min_weight - 1e-6:
                    return False
                if weight > constraints.max_weight + 1e-6:
                    return False
            
            return True
            
        except Exception:
            return False
    
    def _calculate_risk_attribution(self, returns: pd.DataFrame, 
                                  weights: Dict[str, float]) -> Dict[str, float]:
        """Calculate risk attribution by asset"""
        try:
            weight_array = np.array([weights.get(col, 0.0) for col in returns.columns])
            
            # Covariance matrix
            cov_estimator = CovarianceEstimator(RiskModel.SHRINKAGE_COVARIANCE)
            sigma = cov_estimator.estimate(returns)
            
            # Portfolio variance
            portfolio_var = weight_array @ sigma @ weight_array
            
            # Marginal risk contributions
            marginal_contrib = sigma @ weight_array
            
            # Risk contributions
            risk_contrib = weight_array * marginal_contrib
            
            # Normalize to percentages
            if portfolio_var > 0:
                risk_contrib_pct = risk_contrib / portfolio_var
            else:
                risk_contrib_pct = np.zeros_like(risk_contrib)
            
            return dict(zip(returns.columns, risk_contrib_pct))
            
        except Exception as e:
            self.logger.error(f"Error calculating risk attribution: {e}")
            return {}
    
    def get_cached_optimization(self, cache_key: str) -> Optional[OptimizationResult]:
        """Get cached optimization result"""
        return self.optimization_cache.get(cache_key)
    
    def clear_cache(self):
        """Clear optimization cache"""
        self.optimization_cache.clear()

# Example usage
if __name__ == "__main__":
    # Example configuration
    config = {
        'risk_aversion': 3.0,
        'tau': 0.025,
        'n_factors': 5,
        'rebalancing_frequency': 'monthly',
        'transaction_costs': 0.001
    }
    
    # Initialize engine
    engine = PortfolioOptimizationEngine(config)
    
    print("Portfolio Optimization Engine initialized")
    print("Available optimization methods:")
    print("- Modern Portfolio Theory (Markowitz)")
    print("- Risk Parity")
    print("- Black-Litterman")
    print("- Factor-Based Optimization")
    print("- Robust Optimization")
    print("- Multi-Objective Optimization")
    print("- Transaction Cost Optimization")