In [21]:
import sys
print(sys.version)

import numpy as np
np.set_printoptions(precision=5, linewidth=120, suppress=True)
%matplotlib inline
import matplotlib
import matplotlib.pyplot as plt
print('matplotlib: {}'.format(matplotlib.__version__))
from matplotlib.colors import LinearSegmentedColormap

from mosek.fusion import *

plt.rcParams['figure.figsize'] = [12, 8]
plt.rcParams['figure.dpi'] = 200

from notebook.services.config import ConfigManager
c = ConfigManager()
c.update('notebook', {"CodeCell": {"cm_config": {"autoCloseBrackets": False}}})

3.6.9 (default, Jan 26 2021, 15:33:00)
[GCC 8.4.0]
matplotlib: 3.3.4

Out[21]:
{'CodeCell': {'cm_config': {'autoCloseBrackets': False}}}
In [22]:
# Generate a large random covariance matrix with few significant eigenvalues
def random_factor_model(N, K, T):
# Generate K uncorrelated, zero mean factors, with weighted variance
S_F = np.diag(range(1, K + 1))
Z_F = np.random.default_rng(seed=1).multivariate_normal(np.zeros(K), S_F, T).T

# Generate random factor model parameters
B = np.random.default_rng(seed=2).normal(size=(N, K))
a = np.random.default_rng(seed=3).normal(loc=1, size=(N, 1))
e = np.random.default_rng(seed=4).multivariate_normal(np.zeros(N), np.eye(N), T).T

# Generate N time-series from the factors
Z = a + B @ Z_F + e

# Residual covariance
S_theta = np.cov(e)
S_theta = np.diag(np.diag(S_theta))

# Optimization parameters
m = np.mean(Z, axis=1)
S = np.cov(Z)
#print(np.linalg.eigvalsh(np.corrcoef(Z))[-20:])

return m, S, B, S_F, S_theta

In [ ]:
# Solve optimization
def Markowitz(N, m, G, gamma2):
with Model("markowitz") as M:
# Settings
#M.setLogHandler(sys.stdout)

# Decision variable (fraction of holdings in each security)
# The variable x is restricted to be positive, which imposes the constraint of no short-selling.
x = M.variable("x", N, Domain.greaterThan(0.0))

# Budget constraint
M.constraint('budget', Expr.sum(x), Domain.equalsTo(1))

# Objective
M.objective('obj', ObjectiveSense.Maximize, Expr.dot(m, x))

# Imposes a bound on the risk
M.constraint('risk', Expr.vstack(np.sqrt(gamma2), Expr.mul(G.T, x)), Domain.inQCone())

# Solve optimization
M.solve()

returns = M.primalObjValue()
portfolio = x.level()
time = M.getSolverDoubleInfo("optimizerTime")

return returns, time

In [ ]:
# Risk limit
gamma2 = 0.1

# Number of factors
K = 10

# Generate runtime data
# NOTE: This can have a long runtime, depending on the range given for n below!
list_runtimes_orig = []
list_runtimes_factor = []
for n in range(5, 13):
N = 2**n
T = N + 2**(n-1)
m, S, B, S_F, S_theta = random_factor_model(N, K, T)

F = np.linalg.cholesky(S_F)
G_factor = np.block([[B @ F, np.sqrt(S_theta)]])
G_orig = np.linalg.cholesky(S)

optimum_orig, runtime_orig = Markowitz(N, m, G_orig, gamma2)
optimum_factor, runtime_factor = Markowitz(N, m, G_factor, gamma2)
list_runtimes_orig.append((N, runtime_orig))
list_runtimes_factor.append((N, runtime_factor))

tup_N_orig, tup_time_orig = list(zip(*list_runtimes_orig))
tup_N_factor, tup_time_factor = list(zip(*list_runtimes_factor))

In [ ]:
# Runtime plot
plt.plot(tup_N_orig, tup_time_orig, "-o")
plt.plot(tup_N_factor, tup_time_factor, "-o")
plt.xlabel("N")
plt.ylabel("runtime (s)")
ax = plt.gca()
ax.set_xscale('log', base=2)
ax.set_yscale('log')
ax.grid()
legend = ["Cholesky", "factor model"]
plt.legend(legend)