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)