QP Portfolio › Implementation

import numpy as np
from scipy.optimize import minimize


def optimize_portfolio(expected_returns, cov_matrix, target_return=None):
    """
    Optimize portfolio using quadratic programming

    Parameters:
    - expected_returns: array of expected returns for each asset
    - cov_matrix: covariance matrix of returns
    - target_return: target expected return (if None, finds minimum variance)

    Returns:
    - optimal weights as numpy array
    """
    n_assets = len(expected_returns)

    # Objective function: minimize portfolio variance
    def objective(weights):
        return np.dot(weights.T, np.dot(cov_matrix, weights))

    # Budget constraint: weights must sum to 1
    constraints = [{"type": "eq", "fun": lambda w: np.sum(w) - 1}]

    # Return constraint (if target return is specified)
    if target_return is not None:
        constraints.append(
            {
                "type": "eq",
                "fun": lambda w: np.sum(w * expected_returns) - target_return,
            }
        )

    # Bounds: no short selling (weights between 0 and 1)
    bounds = tuple((0, 1) for _ in range(n_assets))

    # Initial guess: equal weights
    x0 = np.array([1 / n_assets] * n_assets)

    # Solve using Sequential Least Squares Programming (SLSQP)
    result = minimize(
        objective, x0, method="SLSQP", bounds=bounds, constraints=constraints
    )

    return result.x if result.success else None

# Optimize for minimum variance portfolio
optimal_weights = optimize_portfolio(expected_returns.values, cov_matrix.values)

# Calculate portfolio statistics
portfolio_return = np.sum(optimal_weights * expected_returns.values)
portfolio_variance = np.dot(
    optimal_weights, np.dot(cov_matrix.values, optimal_weights)
)
portfolio_risk = np.sqrt(portfolio_variance)
sharpe_ratio = portfolio_return / portfolio_risk

# Optimize for 8% target return
target_return = 0.08
weights_target = optimize_portfolio(
    expected_returns.values, cov_matrix.values, target_return=target_return
)

# Calculate statistics for target return portfolio
ret_target = np.sum(weights_target * expected_returns.values)
var_target = np.dot(weights_target, np.dot(cov_matrix.values, weights_target))
risk_target = np.sqrt(var_target)
sharpe_target = ret_target / risk_target

# Create comparison data
comparison_data = {
    "Metric": ["Expected Return", "Portfolio Risk", "Sharpe Ratio"],
    "Minimum Variance": [
        f"{portfolio_return * 100:.2f}%",
        f"{portfolio_risk * 100:.2f}%",
        f"{sharpe_ratio:.3f}",
    ],
    "Target Return (8%)": [
        f"{ret_target * 100:.2f}%",
        f"{risk_target * 100:.2f}%",
        f"{sharpe_target:.3f}",
    ],
}

def solve_qp_custom(cov_matrix, expected_returns, target_return=None):
    """
    Custom QP solver using the method of Lagrange multipliers
    Solves the minimum variance portfolio problem analytically
    """
    n = len(expected_returns)

    if target_return is None:
        # Minimum variance portfolio
        # Solve: min (1/2) * w^T * Sigma * w subject to sum(w) = 1
        # Using method of Lagrange multipliers

        # The optimality conditions are:
        #   Sigma * w - lambda * 1 = 0  (gradient condition)
        #   1^T * w = 1                  (constraint)
        #
        # Rearranging: Sigma * w = lambda * 1
        # In matrix form: [2*Sigma, -1; 1^T, 0] * [w; lambda] = [0; 1]
        # The factor of 2 comes from the Hessian (second derivative) of (1/2)w^T*Sigma*w
        A = np.zeros((n + 1, n + 1))
        A[:n, :n] = 2 * cov_matrix  # 2 * Sigma (Hessian of the objective)
        A[:n, n] = -1  # negative constraint coefficients for lambda
        A[n, :n] = 1  # constraint coefficients (1^T)
        A[n, n] = 0  # zero in bottom-right

        b = np.zeros(n + 1)
        b[n] = 1  # constraint value

        solution = np.linalg.solve(A, b)
        weights = solution[:n]

    else:
        # Target return portfolio
        # Solve: min w^T * Sigma * w subject to sum(w) = 1 and sum(w*mu) = target
        A = np.zeros((n + 2, n + 2))
        A[:n, :n] = 2 * cov_matrix
        A[:n, n] = -1  # budget constraint coefficient for lambda1
        A[
            :n, n + 1
        ] = -expected_returns  # return constraint coefficient for lambda2
        A[n, :n] = 1
        A[n + 1, :n] = expected_returns

        b = np.zeros(n + 2)
        b[n] = 1  # budget constraint
        b[n + 1] = target_return  # return constraint

        solution = np.linalg.solve(A, b)
        weights = solution[:n]

    return weights


# Test the custom solver
weights_custom = solve_qp_custom(cov_matrix.values, expected_returns.values)
ret_custom = np.sum(weights_custom * expected_returns.values)
var_custom = np.dot(weights_custom, np.dot(cov_matrix.values, weights_custom))
risk_custom = np.sqrt(var_custom)