Skip to content

Instantly share code, notes, and snippets.

@tschm
Created September 3, 2022 19:27
Show Gist options
  • Save tschm/02c4f8629408c3088af8f00fe0c39590 to your computer and use it in GitHub Desktop.
Save tschm/02c4f8629408c3088af8f00fe0c39590 to your computer and use it in GitHub Desktop.
import time
import numpy as np
import cvxpy as cp
from mosek.fusion import Model, Domain, Expr, ObjectiveSense
def cvx_explicit(cov_factor, residual_var, F, mu, gamma, Lmax):
n, m = F.shape
w = cp.Variable(n)
f = cp.Variable(m)
ret = mu.T @ w
risk = cp.sum_squares(
np.transpose(np.linalg.cholesky(cov_factor)) @ f
) + cp.sum_squares(np.sqrt(residual_var) @ w)
problem = cp.Problem(
cp.Maximize(ret - gamma * risk),
[cp.sum(w) == 0, f == F.T @ w, cp.norm(w, 1) <= Lmax],
)
problem.solve()
return w.value
def cvx_implicit(cov_factor, residual_var, F, mu, gamma, Lmax):
n, m = F.shape
w = cp.Variable(n)
ret = mu.T @ w
risk = cp.sum_squares(
np.transpose(np.linalg.cholesky(cov_factor)) @ (F.T @ w)
) + cp.sum_squares(np.sqrt(residual_var) @ w)
problem = cp.Problem(
cp.Maximize(ret - gamma * risk),
[cp.sum(w) == 0, cp.norm(w, 1) <= Lmax],
)
problem.solve()
return w.value
# a case for monkey-patching?
def sum_squares(model, vector):
x = model.variable()
model.constraint(Expr.vstack(x, 0.5, vector), Domain.inRotatedQCone())
return x
# a case for monkey-patching
def one_norm(model, vector):
p = model.variable(n, Domain.unbounded())
model.constraint(Expr.hstack(p, vector), Domain.inQCone())
return Expr.sum(p)
def mosek_implicit(cov_factor, residual_var, F, mu, gamma, Lmax):
n, m = F.shape
with Model("cash-neutral-factor") as M:
w = M.variable("w", n, Domain.unbounded())
# position in factor space
f = Expr.mul(F.T, w)
expected_return = Expr.dot(mu, w)
# cash-neutral
M.constraint("cash-neutral", Expr.sum(w), Domain.equalsTo(0.0))
# Leverage
M.constraint("leverage", one_norm(M, w), Domain.lessThan(Lmax))
# Residual variances
res_var = sum_squares(M, Expr.dot(np.sqrt(residual_var), w))
# Variance in factor space
factor_var = sum_squares(
M, Expr.mul(np.transpose(np.linalg.cholesky(cov_factor)), f)
)
risk = Expr.add(factor_var, res_var)
M.objective(
"obj",
ObjectiveSense.Maximize,
Expr.sub(expected_return, Expr.mul(gamma, risk)),
)
M.solve()
return w.level()
if __name__ == "__main__":
n = 2000
m = 200
mu = np.random.randn(n)
# construct a factor covariance matrix
Sigma_tilde = np.random.randn(m, m)
Sigma_tilde = Sigma_tilde.T.dot(Sigma_tilde)
# construct idiosyncratic variances
D = np.random.uniform(0, 0.9, size=n)
# construct the factor loadings
F = np.random.randn(n, m)
Lmax = 2
gamma = 0.1
# solve the problem with the explicit construction of the f variable
t = time.time()
for i in range(5):
w = cvx_explicit(
cov_factor=Sigma_tilde, residual_var=D, F=F, mu=mu, Lmax=Lmax, gamma=gamma
)
print(mu.T @ w)
print(time.time() - t)
# solve the problem but without the explicit construction of the f variable
t = time.time()
for i in range(5):
w = cvx_implicit(
cov_factor=Sigma_tilde, residual_var=D, F=F, mu=mu, Lmax=Lmax, gamma=gamma
)
print(mu.T @ w)
print(time.time() - t)
# solve the same problem with Mosek fusion
t = time.time()
for i in range(5):
w = mosek_implicit(
cov_factor=Sigma_tilde, residual_var=D, F=F, mu=mu, Lmax=Lmax, gamma=gamma
)
print(mu.T @ w)
print(time.time() - t)
@tschm
Copy link
Author

tschm commented Sep 3, 2022

Careful please. This is work in progress. We are solving a standard Markowitz problem with realistic dimensions, e.g. 2000 assets and 200 factors. The core idea is to project the weights (in asset space) down into the lower dimensional factor space.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment