Portfolio optimization

Portfolio optimization seeks to allocate assets in a way that maximizes the risk adjusted return,

\[\begin{split}\begin{array}{ll} \mbox{maximize} & \mu^T x - \gamma \left( x^T \Sigma x \right) \\ \mbox{subject to} & \boldsymbol{1}^T x = 1 \\ & x \ge 0 \end{array}\end{split}\]

where \(x \in \mathbf{R}^{n}\) represents the portfolio, \(\mu \in \mathbf{R}^{n}\) the vector of expected returns, \(\gamma > 0\) the risk aversion parameter, and \(\Sigma \in \mathbf{S}^{n}_{+}\) the risk model covariance matrix. The risk model is usually assumed to be the sum of a diagonal and a rank \(k < n\) matrix,

\[\Sigma = F F^T + D,\]

where \(F \in \mathbf{R}^{n \times k}\) is the factor loading matrix and \(D \in \mathbf{S}^{n}_{+}\) is a diagonal matrix describing the asset-specific risk. The resulting problem has the following equivalent form,

\[\begin{split}\begin{array}{ll} \mbox{minimize} & \frac{1}{2} x^T D x + \frac{1}{2} y^T y - \frac{1}{2\gamma}\mu^T x \\ \mbox{subject to} & y = F^T x \\ & \boldsymbol{1}^T x = 1 \\ & x \ge 0 \end{array}\end{split}\]

Python

import osqp
import numpy as np
import scipy as sp
from scipy import sparse

# Generate problem data
sp.random.seed(1)
n = 100
k = 10
F = sparse.random(n, k, density=0.7, format='csc')
D = sparse.diags(np.random.rand(n) * np.sqrt(k), format='csc')
mu = np.random.randn(n)
gamma = 1

# OSQP data
P = sparse.block_diag([D, sparse.eye(k)], format='csc')
q = np.hstack([-mu / (2*gamma), np.zeros(k)])
A = sparse.bmat([[F.T,             -sparse.eye(k)],
                 [np.ones((1, n)),  None],
                 [sparse.eye(n),    None]], format='csc')
l = np.hstack([np.zeros(k), 1., np.zeros(n)])
u = np.hstack([np.zeros(k), 1., np.ones(n)])

# Create an OSQP object
prob = osqp.OSQP()

# Setup workspace
prob.setup(P, q, A, l, u)

# Solve problem
res = prob.solve()

Matlab

% Generate problem data
rng(1)
n = 100;
k = 10;
F = sprandn(n, k, 0.7);
D = sparse(diag( sqrt(k)*rand(n,1) ));
mu = randn(n, 1);
gamma = 1;

% OSQP data
P = blkdiag(D, speye(k));
q = [-mu/(2*gamma); zeros(k, 1)];
A = [F', -speye(k);
     ones(1, n), zeros(1, k);
     speye(n), sparse(n, k)];
l = [zeros(k, 1); 1; zeros(n, 1)];
u = [zeros(k, 1); 1; ones(n, 1)];

% Create an OSQP object
prob = osqp;

% Setup workspace
prob.setup(P, q, A, l, u);

% Solve problem
res = prob.solve();

CVXPY

from cvxpy import *
import numpy as np
import scipy as sp
from scipy import sparse

# Generate problem data
sp.random.seed(1)
n = 100
k = 10
F = sparse.random(n, k, density=0.7, format='csc')
D = sparse.diags(np.random.rand(n) * np.sqrt(k), format='csc')
mu = np.random.randn(n)
gamma = 1
Sigma = F@F.T + D

# Define problem
x = Variable(n)
objective = mu.T@x - gamma*quad_form(x, Sigma)
constraints = [sum(x) == 1, x >= 0]

# Solve with OSQP
Problem(Maximize(objective), constraints).solve(solver=OSQP)

YALMIP

% Generate problem data
rng(1)
n = 100;
k = 10;
F = sprandn(n, k, 0.7);
D = sparse(diag( sqrt(k)*rand(n,1) ));
mu = randn(n, 1);
gamma = 1;
Sigma = F*F' + D;

% Define problem
x = sdpvar(n, 1);
objective = gamma * (x'*Sigma*x) - mu'*x;
constraints = [sum(x) == 1, x >= 0];

% Solve with OSQP
options = sdpsettings('solver', 'osqp');
optimize(constraints, objective, options);