LaTeX macros (hidden cell) $ \newcommand{\Q}{\mathcal{Q}} \newcommand{\ECov}{\boldsymbol{\Sigma}} \newcommand{\EMean}{\boldsymbol{\mu}} \newcommand{\EAlpha}{\boldsymbol{\alpha}} \newcommand{\EBeta}{\boldsymbol{\beta}} $

Imports and configuration

In [14]:
import sys
import os
import re
import glob
import datetime as dt

import numpy as np
import pandas as pd
%matplotlib inline
import matplotlib
import matplotlib.pyplot as plt
from matplotlib.colors import LinearSegmentedColormap

from mosek.fusion import *

from notebook.services.config import ConfigManager

from portfolio_tools import data_download, DataReader
In [15]:
# Version checks
print(sys.version)
print('matplotlib: {}'.format(matplotlib.__version__))

# Jupyter configuration
c = ConfigManager()
c.update('notebook', {"CodeCell": {"cm_config": {"autoCloseBrackets": False}}})  

# Numpy options
np.set_printoptions(precision=5, linewidth=120, suppress=True)

# Pandas options
pd.set_option('display.max_rows', None)

# Matplotlib options
plt.rcParams['figure.figsize'] = [12, 8]
plt.rcParams['figure.dpi'] = 200
3.6.9 (default, Jan 26 2021, 15:33:00) 
[GCC 8.4.0]
matplotlib: 3.3.4

Prepare input data

In this example, the input data is given. It consists of the vector $\EMean$ of expected returns, and the covariance matrix $\ECov$.

In [16]:
# Linear return statistics on the investment horizon
mu = np.array([0.07197349, 0.15518171, 0.17535435, 0.0898094 , 0.42895777, 0.39291844, 0.32170722, 0.18378628])
Sigma = np.array([
        [0.09460323, 0.03735969, 0.03488376, 0.03483838, 0.05420885, 0.03682539, 0.03209623, 0.03271886],
        [0.03735969, 0.07746293, 0.03868215, 0.03670678, 0.03816653, 0.03634422, 0.0356449 , 0.03422235],
        [0.03488376, 0.03868215, 0.06241065, 0.03364444, 0.03949475, 0.03690811, 0.03383847, 0.02433733],
        [0.03483838, 0.03670678, 0.03364444, 0.06824955, 0.04017978, 0.03348263, 0.04360484, 0.03713009],
        [0.05420885, 0.03816653, 0.03949475, 0.04017978, 0.17243352, 0.07886889, 0.06999607, 0.05010711],
        [0.03682539, 0.03634422, 0.03690811, 0.03348263, 0.07886889, 0.09093307, 0.05364518, 0.04489357],
        [0.03209623, 0.0356449 , 0.03383847, 0.04360484, 0.06999607, 0.05364518, 0.09649728, 0.04419974],
        [0.03271886, 0.03422235, 0.02433733, 0.03713009, 0.05010711, 0.04489357, 0.04419974, 0.08159633]
      ])

Define the optimization model

The optimization problem we would like to solve is $$ \begin{array}{lrcl} \mbox{maximize} & \EMean^\mathsf{T}\mathbf{x} & &\\ \mbox{subject to} & \left(\gamma^2, \frac{1}{2}, \mathbf{G}^\mathsf{T}\mathbf{x}\right) & \in & \Q_\mathrm{r}^{N+2},\\ & \mathbf{1}^\mathsf{T}\mathbf{x} & = & 1,\\ & \mathbf{x} & \geq & 0.\\ \end{array} $$

Here we define this model in MOSEK Fusion.

In [23]:
# Define function solving the optimization model
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(gamma2, 0.5, Expr.mul(G.transpose(), x)), Domain.inRotatedQCone())

        # Solve optimization
        M.solve()
        
        returns = M.primalObjValue()
        portfolio = x.level()
        
    return returns, portfolio

Run the optimization

Define the parameters

The problem parameters are the number of securities $N$ and the risk limit $\gamma^2$.

In [20]:
N = mu.shape[0]  # Number of securities
gamma2 = 0.05   # Risk limit (variance)

Factorize the covariance matrix

Here we factorize $\ECov$ because the model is defined in conic form, and it expects a matrix $G$ such that $\ECov = GG^\mathsf{T}$.

In [21]:
G = np.linalg.cholesky(Sigma)  # Cholesky factor of S to use in conic risk constraint 

Solve the optimization problem

Next we call the function that defines the Fusion model and runs the optimization.

In [24]:
# Run optimization 
f, x = Markowitz(N, mu, G, gamma2)
print("========================\n")
print("RESULTS:")
print(f"Optimal expected portfolio return: {f*100:.4f}%")
print(f"Optimal portfolio weights: {x}")
print(f"Sum of weights: {np.sum(x)}")
Problem
  Name                   : markowitz       
  Objective sense        : max             
  Type                   : CONIC (conic optimization problem)
  Constraints            : 11              
  Cones                  : 1               
  Scalar variables       : 19              
  Matrix variables       : 0               
  Integer variables      : 0               

Optimizer started.
Presolve started.
Linear dependency checker started.
Linear dependency checker terminated.
Eliminator started.
Freed constraints in eliminator : 0
Eliminator terminated.
Eliminator - tries                  : 1                 time                   : 0.00            
Lin. dep.  - tries                  : 1                 time                   : 0.00            
Lin. dep.  - number                 : 0               
Presolve terminated. Time: 0.00    
Problem
  Name                   : markowitz       
  Objective sense        : max             
  Type                   : CONIC (conic optimization problem)
  Constraints            : 11              
  Cones                  : 1               
  Scalar variables       : 19              
  Matrix variables       : 0               
  Integer variables      : 0               

Optimizer  - threads                : 20              
Optimizer  - solved problem         : the primal      
Optimizer  - Constraints            : 8
Optimizer  - Cones                  : 1
Optimizer  - Scalar variables       : 17                conic                  : 10              
Optimizer  - Semi-definite variables: 0                 scalarized             : 0               
Factor     - setup time             : 0.00              dense det. time        : 0.00            
Factor     - ML order time          : 0.00              GP order time          : 0.00            
Factor     - nonzeros before factor : 36                after factor           : 36              
Factor     - dense dim.             : 0                 flops                  : 6.00e+02        
ITE PFEAS    DFEAS    GFEAS    PRSTATUS   POBJ              DOBJ              MU       TIME  
0   1.0e+00  1.6e+00  1.4e+00  0.00e+00   0.000000000e+00   3.889087297e-01   1.0e+00  0.02  
1   1.6e-01  2.5e-01  4.9e-02  7.08e-01   4.275713587e-01   5.950950777e-01   1.6e-01  0.03  
2   5.5e-02  8.6e-02  8.7e-03  1.69e+00   3.610061700e-01   4.006907873e-01   5.5e-02  0.03  
3   3.1e-02  4.8e-02  4.0e-03  1.36e+00   3.098345093e-01   3.290205233e-01   3.1e-02  0.04  
4   8.9e-03  1.4e-02  6.4e-04  1.27e+00   2.910409783e-01   2.958703450e-01   8.9e-03  0.04  
5   3.5e-03  5.5e-03  1.7e-04  1.05e+00   2.845165201e-01   2.862826893e-01   3.5e-03  0.04  
6   1.3e-03  2.1e-03  4.4e-05  8.02e-01   2.789487823e-01   2.796238224e-01   1.3e-03  0.04  
7   1.6e-04  2.5e-04  1.8e-06  1.02e+00   2.770260963e-01   2.771103927e-01   1.6e-04  0.04  
8   2.1e-05  3.3e-05  8.5e-08  1.00e+00   2.767556228e-01   2.767667414e-01   2.1e-05  0.04  
9   1.4e-06  2.2e-06  1.4e-09  9.99e-01   2.767193888e-01   2.767201074e-01   1.4e-06  0.04  
10  7.3e-09  1.1e-08  5.4e-13  1.00e+00   2.767173176e-01   2.767173215e-01   7.3e-09  0.04  
Optimizer terminated. Time: 0.05    


Interior-point solution summary
  Problem status  : PRIMAL_AND_DUAL_FEASIBLE
  Solution status : OPTIMAL
  Primal.  obj: 2.7671731760e-01    nrm: 1e+00    Viol.  con: 5e-09    var: 0e+00    cones: 2e-09  
  Dual.    obj: 2.7671732147e-01    nrm: 6e+00    Viol.  con: 0e+00    var: 4e-09    cones: 0e+00  
========================

RESULTS:
Optimal expected portfolio return: 27.6717%
Optimal portfolio weights: [0.      0.09126 0.26911 0.      0.02531 0.32162 0.17652 0.11618]
Sum of weights: 0.9999999951391261

Test result

In [32]:
expected_x = np.array([0., 0.09126, 0.26911, 0., 0.02531, 0.32162, 0.17652, 0.11618])
diff = np.sum(np.abs(expected_x - x))
assert diff < 1e-4, f"Resulting portfolio does not match expected one. Difference is {diff}"