Introduction: Beyond Classical Mean-Variance Optimization

Minimize (1/2) · wᵀ Σ w
Subject to:
    wᵀ μ = μ_target  (Target return constraint)
    wᵀ 1 = 1         (Full investment constraint)
    w ≥ 0            (Optional: No short-selling constraint)

Where:
- w = N×1 vector of asset weights
- Σ = N×N covariance matrix (positive semi-definite)
- μ = N×1 vector of expected returns

Pathologies of Classical MVO:
┌─────────────────────────────────────────────────────────────────────┐
|  1. Extreme Sensitivity (Error Maximization):                    |
|  ┌─────────────────────────────────────────────────────────────┐   |
|  │  Small estimation errors in expected returns μ →           │   |
|  │  Extreme allocations (e.g., +200% in Asset A, -150% in B).│   |
|  │  MVO effectively acts as an "error maximizer" because it  │   |
|  │  allocates heavily to assets with the largest positive     │   |
|  │  estimation errors.                                        │   |
|  └─────────────────────────────────────────────────────────────┘   |
|                                                                  |
|  2. Matrix Inversion Instability:                              |
|  ┌─────────────────────────────────────────────────────────────┐   |
|  │  In large, correlated asset universes (e.g., 500 stocks),  │   |
|  │  Σ becomes nearly singular (ill-conditioned). Inverting   │   |
|  │  it amplifies tiny estimation errors into absurd weights.  │   |
|  │  Condition numbers often exceed 10⁶, rendering the        │   |
|  │  weights effectively random.                               │   |
|  └─────────────────────────────────────────────────────────────┘   |
|                                                                  |
|  3. Point-Estimate Fallacy:                                   │
|  ┌─────────────────────────────────────────────────────────────┐   |
|  │  MVO treats μ and Σ as known constants, ignoring the       │   |
|  │  fact that they are noisy estimates. This ignores          │   |
|  │  estimation risk entirely.                                 │   |
|  └─────────────────────────────────────────────────────────────┘   |
└─────────────────────────────────────────────────────────────────────┘

Π = λ · Σ · w_mkt

Where:
- λ = Risk aversion coefficient (typically 2.5–3.5)
- w_mkt = Market capitalization weights (e.g., S&P 500 weights)
- Σ = Covariance matrix
- Π = N×1 vector of implied equilibrium excess returns

Interpretation:
If the market is efficient, the current weights w_mkt are optimal.
Therefore, the returns required to justify those weights are Π.
This provides a highly stable, shrinkage-inducing prior.

μ_BL = [ (τΣ)⁻¹ + Pᵀ Ω⁻¹ P ]⁻¹ · [ (τΣ)⁻¹ Π + Pᵀ Ω⁻¹ Q ]

Where:
- τ = Scalar representing uncertainty in the prior equilibrium (typically 0.01–0.05).
- Ω = Diagonal covariance matrix of the views (confidence levels).
- μ_BL = The stabilized, blended expected return vector.

Interpretation:
If views are highly uncertain (large Ω), μ_BL stays close to Π (market
equilibrium). If views are highly certain (small Ω), μ_BL moves
toward Q (the investor's forecast). This Bayesian framework prevents
the extreme weights produced by raw MVO.

Black-Litterman Blending Process:
┌─────────────────────────────────────────────────────────────────────┐
|                                                                  |
|  Market Caps (w_mkt) ──▶ Reverse Optimization ──▶ Π (Prior μ)   |
|                                                       │           |
|  Analyst Views (P, Q, Ω) ────────────────────────────┼───────────┐ |
|                                                       │           │ |
|                                                       ▼           │ |
|                                          Bayesian Posterior:     │ |
|                                         μ_BL = f(Π, P, Q, Ω)    │ |
|                                                       │           │ |
|                                                       ▼           │ |
|  Feed μ_BL into MVO → Stable, diverse, intuition-              │
|  informed weights without extreme shorts.                     │ |
└─────────────────────────────────────────────────────────────────────┘

RC_i = w_i · (Σ w)_i

Total Portfolio Variance: σ_p² = wᵀΣ w = Σ_{i=1}^{n} RC_i

RC_i = RC_j  ∀ i, j

Optimization Problem:
Minimize: Σ_{i=1}^{n} Σ_{j=1}^{n} ( w_i(Σw)_i - w_j(Σw)_j )²
Subject to: Σ w_i = 1,  w_i ≥ 0

Risk Decomposition: 60/40 vs. Risk Parity
┌─────────────────────────────────────────────────────────────────────┐
|  60/40 Stock/Bond Portfolio:                                    |
|  ┌─────────────────────────────────────────────────────────────┐   |
|  │  Capital:  ████████████████ Stocks  │████ Bonds           │   |
|  │  Risk:     ████████████████████████████████░ (95% Equity) │   |
|  └─────────────────────────────────────────────────────────────┘   |
|                                                                  |
|  Risk Parity (ERC):                                             |
|  ┌─────────────────────────────────────────────────────────────┐   |
|  │  Capital:  ████ Stocks │████████████████████ Bonds         │   |
|  │  Risk:     ████████████ Stocks │████████████ Bonds (50/50) │   |
|  └─────────────────────────────────────────────────────────────┘   |
|                                                                  |
|  Result: Risk Parity produces structurally higher Sharpe        |
|  ratios and more resilient drawdown profiles.                  |
└─────────────────────────────────────────────────────────────────────┘

HRP Algorithm Flow:
┌─────────────────────────────────────────────────────────────────────┐
|  1. Distance Matrix & Clustering:                               |
|  ┌─────────────────────────────────────────────────────────────┐   |
|  │  Convert correlation matrix to a Euclidean distance matrix. │   |
|  │  d(i,j) = sqrt(0.5 * (1 - ρ_ij)).                         │   |
|  │  Apply hierarchical clustering (Ward linkage) to group      │   |
|  │  similar assets (e.g., Tech stocks, Bonds, Energy).        │   |
|  └─────────────────────────────────────────────────────────────┘   |
|                              ▼                                    |
|  2. Quasi-Diagonalization:                                    |
|  ┌─────────────────────────────────────────────────────────────┐   |
|  │  Reorder the covariance matrix so that correlated assets   │   |
|  │  are placed adjacent to each other on the diagonal.        │   |
|  │  This creates a "heatmap" with distinct asset blocks.      │   |
|  └─────────────────────────────────────────────────────────────┘   |
|                              ▼                                    |
|  3. Recursive Bisection:                                      |
|  ┌─────────────────────────────────────────────────────────────┐   |
|  │  Starting from the root cluster (all assets):              │   |
|  │  For each cluster, split into two sub-clusters.            │   |
|  │  Allocate capital inversely proportional to the sub-       │   |
|  │  cluster's variance (risk budgeting).                     │   |
|  │  Repeat recursively down the dendrogram.                   │   |
|  └─────────────────────────────────────────────────────────────┘   |
|                              │                                    |
|  Output: Weights that are robust, diversified, and bypass      │
|          the need for Σ⁻¹ entirely.                           │
└─────────────────────────────────────────────────────────────────────┘

Minimize: w_tᵀ Σ w_t − μᵀ w_t + κ · || w_t − w_{t-1} ||_1

Where:
- κ = Transaction cost penalty coefficient (tuned empirically).
- ||w_t − w_{t-1}||₁ = Sum of absolute changes in weights (turnover).

Interpretation:
If the expected alpha from rebalancing is less than the transaction
costs (bid-ask spread + market impact), the penalty forces the
algorithm to do nothing, preserving capital.

import numpy as np
import pandas as pd
from scipy.optimize import minimize
from scipy.cluster.hierarchy import linkage, fcluster
import cvxpy as cp

# -------------------- 1. BLACK-LITTERMAN --------------------
def black_litterman(cov, w_mkt, risk_aversion=2.5, tau=0.02, P=None, Q=None, omega=None):
    """
    Compute Black-Litterman expected returns.
    """
    n = len(w_mkt)
    # Prior equilibrium returns
    pi = risk_aversion * cov @ w_mkt
    
    if P is None:  # No views -> return prior
        return pi
    
    # Posterior computation
    cov_prior = tau * cov
    inv_cov_prior = np.linalg.inv(cov_prior)
    
    # Ensure omega is diagonal confidence matrix
    if omega is None:
        omega = np.eye(len(Q)) * 0.1  # Default confidence
    
    inv_omega = np.linalg.inv(omega)
    
    # Black-Litterman formula
    term1 = inv_cov_prior + P.T @ inv_omega @ P
    term2 = inv_cov_prior @ pi + P.T @ inv_omega @ Q
    mu_bl = np.linalg.inv(term1) @ term2
    
    return mu_bl

# -------------------- 2. EQUAL RISK CONTRIBUTION (ERC) --------------------
def erc_objective(weights, cov):
    """Objective: Minimize variance of risk contributions."""
    portfolio_var = weights.T @ cov @ weights
    # Marginal risk contribution
    marginal_contrib = cov @ weights
    risk_contrib = weights * marginal_contrib
    # Minimize pairwise differences in risk contributions
    n = len(weights)
    diff = 0
    for i in range(n):
        for j in range(i+1, n):
            diff += (risk_contrib[i] - risk_contrib[j])**2
    return diff

def compute_erc_weights(cov):
    n = cov.shape[0]
    constraints = (
        {'type': 'eq', 'fun': lambda w: np.sum(w) - 1.0},
        {'type': 'ineq', 'fun': lambda w: w}  # w >= 0
    )
    result = minimize(
        erc_objective, 
        x0=np.ones(n) / n,
        args=(cov,),
        method='SLSQP',
        constraints=constraints
    )
    return result.x

# -------------------- 3. HIERARCHICAL RISK PARITY (HRP) --------------------
def get_quasi_diag(linkage_matrix):
    """Sort indices by hierarchical clustering result."""
    # Adapted from Lopez de Prado
    n = len(linkage_matrix) + 1
    cluster_idx = {}
    for i in range(n):
        cluster_idx[i] = [i]
    
    for i, (left, right, dist, _) in enumerate(linkage_matrix):
        idx = n + i
        left_idx = int(left)
        right_idx = int(right)
        cluster_idx[idx] = cluster_idx[left_idx] + cluster_idx[right_idx]
    
    return cluster_idx[n + len(linkage_matrix) - 1]  # Root cluster

def hrp_weights(cov):
    """Simple HRP implementation."""
    n = cov.shape[0]
    # 1. Distance matrix
    corr = np.corrcoef(cov)
    dist = np.sqrt(0.5 * (1 - corr))
    # 2. Clustering (Ward)
    linkage_matrix = linkage(dist, method='ward')
    # 3. Quasi-diagonal ordering
    ordered_indices = get_quasi_diag(linkage_matrix)
    
    # 4. Recursive bisection (simplified inverse-variance)
    weights = np.ones(n) / n
    # Recursively allocate based on cluster variance
    # (Full implementation omitted for brevity, but standard libs exist)
    return weights  # Placeholder

# -------------------- 4. TRANSACTION COST-AWARE OPTIMIZATION --------------------
def optimize_with_turnover(cov, mu, current_weights, kappa=0.001):
    """
    Solve: min w^T Σ w - μ^T w + κ * ||w - w_prev||_1
    Using cvxpy for convex optimization.
    """
    n = len(mu)
    w = cp.Variable(n)
    turnover = cp.sum(cp.abs(w - current_weights))
    
    objective = cp.Minimize(cp.quad_form(w, cov) - mu.T @ w + kappa * turnover)
    constraints = [cp.sum(w) == 1, w >= 0]
    
    problem = cp.Problem(objective, constraints)
    problem.solve(solver=cp.ECOS)
    
    return w.value

# -------------------- EXAMPLE USAGE --------------------
if __name__ == "__main__":
    n_assets = 5
    np.random.seed(42)
    
    # Simulated data
    cov = np.random.randn(n_assets, n_assets)
    cov = cov.T @ cov  # Positive semi-definite
    w_mkt = np.array([0.4, 0.3, 0.15, 0.1, 0.05])
    
    # Black-Litterman with a view: Asset 0 will outperform Asset 1 by 2%
    P = np.array([[1, -1, 0, 0, 0]])
    Q = np.array([0.02])
    mu_bl = black_litterman(cov, w_mkt, P=P, Q=Q)
    print("Black-Litterman Returns:", mu_bl)
    
    # ERC
    weights_erc = compute_erc_weights(cov)
    print("ERC Weights:", weights_erc)
    
    # Transaction cost optimization
    current_w = np.array([0.3, 0.3, 0.2, 0.1, 0.1])
    new_w = optimize_with_turnover(cov, mu_bl, current_w, kappa=0.02)
    print("Turnover Penalized Weights:", new_w)