# 11.1 Portfolio Optimization¶

In this section the Markowitz portfolio optimization problem and variants are implemented using the MOSEK optimizer API.

## 11.1.1 A Basic Portfolio Optimization Model¶

The classical Markowitz portfolio optimization problem considers investing in $$n$$ stocks or assets held over a period of time. Let $$x_j$$ denote the amount invested in asset $$j$$, and assume a stochastic model where the return of the assets is a random variable $$r$$ with known mean

$\mu = \mathbf{E} r$

and covariance

$\Sigma = \mathbf{E}(r-\mu)(r-\mu)^T.$

The return of the investment is also a random variable $$y = r^Tx$$ with mean (or expected return)

$\mathbf{E} y = \mu^T x$

and variance (or risk)

$(y - \mathbf{E} y)^2 = x^T \Sigma x.$

The problem facing the investor is to rebalance the portfolio to achieve a good compromise between risk and expected return, e.g., maximize the expected return subject to a budget constraint and an upper bound (denoted $$\gamma$$) on the tolerable risk. This leads to the optimization problem

(1)$\begin{split}\begin{array} {lrcl} \mbox{maximize} & \mu^T x & & \\ \st & e^T x & = & w + e^T x^0, \\ & x^T \Sigma x & \leq & \gamma^2, \\ & x & \geq & 0. \end{array}\end{split}$

The variables $$x$$ denote the investment i.e. $$x_j$$ is the amount invested in asset $$j$$ and $$x_j^0$$ is the initial holding of asset $$j$$. Finally, $$w$$ is the initial amount of cash available.

A popular choice is $$x^0=0$$ and $$w=1$$ because then $$x_j$$ may be interpreted as the relative amount of the total portfolio that is invested in asset $$j$$.

Since $$e$$ is the vector of all ones then

$e^T x = \sum_{j=1}^n x_j$

is the total investment. Clearly, the total amount invested must be equal to the initial wealth, which is

$w + e^T x^0.$

This leads to the first constraint

$e^T x = w + e^T x^0.$

The second constraint

$x^T \Sigma x \leq \gamma^2$

ensures that the variance, or the risk, is bounded by $$\gamma^2$$. Therefore, $$\gamma$$ specifies an upper bound of the standard deviation the investor is willing to undertake. Finally, the constraint

$x_j \geq 0$

excludes the possibility of short-selling. This constraint can of course be excluded if short-selling is allowed.

The covariance matrix $$\Sigma$$ is positive semidefinite by definition and therefore there exist a matrix $$G$$ such that

(2)$\Sigma = G G^T.$

In general the choice of $$G$$ is not unique and one possible choice of $$G$$ is the Cholesky factorization of $$\Sigma$$. However, in many cases another choice is better for efficiency reasons as discussed in Sec. 11.1.3 (Improving the Computational Efficiency).

For a given $$G$$ we have that

$\begin{split}\begin{array} {lcc} x^T \Sigma x & = & x^T G G^T x \\ & = & \left\|G^T x\right\|^2. \end{array}\end{split}$

Hence, we may write the risk constraint as

$\gamma \geq \left\|G^T x\right\|$

or equivalently

$[\gamma;G^T x] \in \mathcal{Q}^{n+1}.$

where $$\mathcal{Q}^{n+1}$$ is the $$n+1$$ dimensional quadratic cone. Therefore, problem (1) can be written as

(3)$\begin{split}\begin{array} {lrcl} \mbox{maximize} & \mu^T x & &\\ \mbox{subject to} & e^T x & = & w+e^T x^0,\\ & [\gamma;G^T x] & \in & \Q^{n+1},\\ & x & \geq & 0, \end{array}\end{split}$

which is a conic quadratic optimization problem that can easily be solved using MOSEK.

Example data

Subsequently we will use the following sample input taken from [CT07]. We set

$\begin{split}\mu = \left[ \begin{array} {c} 0.1073\\ 0.0737\\ 0.0627 \end{array} \right]\end{split}$

and

$\begin{split}\Sigma = 0.1 \left[ \begin{array} {ccc} 0.2778 & 0.0387 & 0.0021 \\ 0.0387 & 0.1112 & -0.0020 \\ 0.0021 & -0.0020 & 0.0115 \end{array} \right]\end{split}$

This implies

$\begin{split}G^T = \sqrt{0.1} \left[ \begin{array} {ccc} 0.5271 & 0.0734 & 0.0040\\ 0 & 0.3253 & -0.0070\\ 0 & 0 & 0.1069 \end{array} \right]\end{split}$

using 5 significant digits. Moreover, let

$\begin{split}x^0 = \left[ \begin{array} {c} 0.0\\ 0.0\\ 0.0 \end{array} \right]\end{split}$

and

$w = 1.0.$

### 11.1.1.1 Why a Conic Formulation?¶

Problem (1) is a convex quadratically constrained optimization problem that can be solved directly using MOSEK. Why then reformulate it as a conic quadratic optimization problem (3)? The main reason for choosing a conic model is that it is more robust and usually solves faster and more reliably. For instance it is not always easy to numerically validate that the matrix $$\Sigma$$ in (1) is positive semidefinite due to the presence of rounding errors. It is also very easy to make a mistake so $$\Sigma$$ becomes indefinite. These problems are completely eliminated in the conic formulation.

Moreover, observe the constraint

$\left\|G^T x\right\| \leq \gamma$

more numerically robust than

$x^T \Sigma x \leq \gamma^2$

for very small and very large values of $$\gamma$$. Indeed, if say $$\gamma \approx 10^4$$ then $$\gamma^2\approx 10^8$$, which introduces a scaling issue in the model. Hence, using conic formulation we work with the standard deviation instead of variance, which usually gives rise to a better scaled model.

### 11.1.1.2 Implementing the Portfolio Model¶

Creating a matrix formulation

The Optimizer API for Python requires that an optimization problem is entered in the following standard form:

(4)$\begin{split}\begin{array} {lrrcll} \mbox{maximize} & & & c^T \hat{x} & &\\ \mbox{subject to} & l^c & \leq & A \hat{x} & \leq & u^c,\\ & l^x & \leq & \hat{x} & \leq & u^x, \\ & & & \hat{x} \in \K. & & \end{array}\end{split}$

We refer to $$\hat{x}$$ as the API variable. It means we need to reformulate (3). The first step is to introduce auxiliary variables so that the conic constraint involves only unique variables:

(5)$\begin{split}\begin{array} {lrcl} \mbox{maximize} & \mu^T x & & \\ \mbox{subject to} & e^T x & = & w+e^T x^0, \\ & G^T x - t & = & 0, \\ & [s;t] & \in & \Q^{n+1},\\ & x & \geq & 0,\\ & s & = & \gamma. \end{array}\end{split}$

Here $$s$$ is an additional scalar variable and $$t$$ is a vector variable of dimension $$n$$. The next step is to concatenate all the variables into one long variable vector:

(6)$\begin{split}\hat{x} = \left[ x; s; t \right] = \left[ \begin{array} {c} x\\ s\\ t \end{array} \right]\end{split}$

The details of the concatenation are specified below.

Table 11 Storage layout of the $$\hat{x}$$ variable.
Variable Length Offset
$$x$$ $$n$$ $$0$$
$$s$$ $$1$$ $$n$$
$$t$$ $$n$$ $$n+1$$

The offset determines where the variable starts. (Note that all variables are indexed from 0). For instance

$\hat{x}_{n+1+i} = t_i.$

because the offset of the $$t$$ variable is $$n+1$$.

Given the ordering of the variables specified by (6) it is useful to visualize the linear constraints (4) in an explicit block matrix form:

(7)$\begin{split}\left[ \begin{array}{ccc|c|ccc} & 1 & & 0 & & 0 & \\ \hline & & & & -1 & & \\ & G^T & & 0 & & -1 & \\ & & & & & & -1 \end{array}\right] \cdot \left[ \begin{array}{c} \\ x \\ \\ \hline s \\ \hline \\ t \\ \\ \end{array}\right] = \left[ \begin{array}{c} w+e^Tx_0\\ \hline \\0\\ \\ \end{array}\right].\end{split}$

In other words, we should define the specific components of the problem description as follows:

(8)$\begin{split}\begin{array} {rcl} c & = & \left[ \begin{array} {ccc} \mu^T & 0 & 0_{n} \end{array}\right]^T, \\ A & = & \left[ \begin{array} {ccc} e^T & 0 & 0_{n} \\ G^T & 0_n & -I_n \end{array}\right], \\ l^c & = & {\left[\begin{array} {cc} w+e^T x^0 & 0_{n} \end{array}\right]}^T,\\ u^c & = & {\left[\begin{array} {cc} w+e^T x^0 & 0_{n} \end{array}\right]}^T,\\ l^x & = & {\left[\begin{array} {ccc} 0_{n} & \gamma & -\infty_{n} \end{array}\right]}^T,\\ u^x & = & {\left[\begin{array} {ccc} \infty_{n} & \gamma & \infty_{n} \end{array}\right]}^T. \end{array}\end{split}$

Source code example

From the block matrix form (7) and the explicit specification (8), using the offset information in Table 11 it is easy to calculate the index and value of each entry of the linear constraint matrix. The code below sets up the general optimization problem (3) and solves it for the example data. Of course it is only necessary to set non-zero entries of the linear constraint matrix.

Listing 27 Code implementing model (3). Click here to download.
import mosek
import sys

def streamprinter(text):
sys.stdout.write("%s" % text),

if __name__ == '__main__':

n = 3
gamma = 0.05
mu = [0.1073, 0.0737, 0.0627]
GT = [[0.1667, 0.0232, 0.0013],
[0.0000, 0.1033, -0.0022],
[0.0000, 0.0000, 0.0338]]
x0 = [0.0, 0.0, 0.0]
w = 1.0

inf = 0.0 # This value has no significance

with mosek.Env() as env:

# Constraints.

# Total budget constraint - set bounds l^c = u^c
rtemp = w + sum(x0)
task.putconbound(0, mosek.boundkey.fx, rtemp, rtemp)

# The remaining constraints GT * x - t = 0 - set bounds l^c = u^c
task.putconboundlist(range(1 + 0, 1 + n), [mosek.boundkey.fx] * n, [0.0] * n, [0.0] * n)
for j in range(1, 1 + n):
task.putconname(j, "GT[%d]" % j)

# Variables.
task.appendvars(1 + 2 * n)

# Offset of variables into the API variable.
offsetx = 0
offsets = n
offsett = n + 1

# x variables.
# Returns of assets in the objective
task.putclist(range(offsetx + 0, offsetx + n), mu)
# Coefficients in the first row of A
task.putaijlist([0] * n, range(offsetx + 0, offsetx + n), [1.0] * n)
# No short-selling - x^l = 0, x^u = inf
task.putvarboundslice(offsetx, offsetx + n, [mosek.boundkey.lo] * n, [0.0] * n, [inf] * n)
for j in range(0, n):
task.putvarname(offsetx + j, "x[%d]" % (1 + j))

# s variable is a constant equal to gamma
task.putvarbound(offsets + 0, mosek.boundkey.fx, gamma, gamma)
task.putvarname(offsets + 0, "s")

# t variables (t = GT*x).
# Copying the GT matrix in the appropriate block of A
for j in range(0, n):
[1 + j] * n, range(offsetx + 0, offsetx + n), GT[j])
# Diagonal -1 entries in a block of A
task.putaijlist(range(1, n + 1), range(offsett + 0, offsett + n), [-1.0] * n)
# Free - no bounds
task.putvarboundslice(offsett + 0, offsett + n, [mosek.boundkey.fr] * n, [-inf] * n, [inf] * n)
for j in range(0, n):
task.putvarname(offsett + j, "t[%d]" % (1 + j))

# Define the cone spanned by variables (s, t), i.e. dimension = n + 1
task.appendcone(mosek.conetype.quad, 0.0, [offsets] + list(range(offsett, offsett + n)))

# Dump the problem to a human readable OPF file.

# Display solution summary for quick inspection of results.

# Retrieve results
xx = [0.] * (n + 1)
task.getxxslice(mosek.soltype.itr, offsetx + 0, offsets + 1, xx)
expret = sum(mu[j] * xx[j] for j in range(offsetx, offsetx + n))
stddev = xx[offsets]

print("\nExpected return %e for gamma %e\n" % (expret, stddev))


The above code produces the result:

Listing 28 Output from the solver.
Interior-point solution summary
Problem status  : PRIMAL_AND_DUAL_FEASIBLE
Solution status : OPTIMAL
Primal.  obj: 7.4766507287e-02    nrm: 1e+00    Viol.  con: 2e-08    var: 0e+00    cones: 2e-08
Dual.    obj: 7.4766554102e-02    nrm: 3e-01    Viol.  con: 0e+00    var: 3e-08    cones: 0e+00

Expected return 7.476651e-02 for gamma 5.000000e-02


The source code is a direct translation of the model (5) using the explicit block matrix specification (8) but a few comments are nevertheless in place.

In the lines

            # Offset of variables into the API variable.
offsetx = 0
offsets = n
offsett = n + 1


offsets into the MOSEK API variable are stored as in Table 11. The code

            # Returns of assets in the objective
task.putclist(range(offsetx + 0, offsetx + n), mu)
# Coefficients in the first row of A
task.putaijlist([0] * n, range(offsetx + 0, offsetx + n), [1.0] * n)
# No short-selling - x^l = 0, x^u = inf
task.putvarboundslice(offsetx, offsetx + n, [mosek.boundkey.lo] * n, [0.0] * n, [inf] * n)
for j in range(0, n):
task.putvarname(offsetx + j, "x[%d]" % (1 + j))


sets up the data for $$x$$ variables. For instance

            # Returns of assets in the objective
task.putclist(range(offsetx + 0, offsetx + n), mu)


inputs the objective coefficients for the $$x$$ variables. Moreover, the code

            for j in range(0, n):
task.putvarname(offsetx + j, "x[%d]" % (1 + j))


assigns meaningful names to the API variables. This is not needed but it makes debugging easier.

Note that the solution values are only accessed for the interesting variables; for instance the auxiliary variable $$t$$ is omitted from this process.

### 11.1.1.3 Debugging Tips¶

Implementing an optimization model in Optimizer API for Python can be error-prone. In order to check the code for accidental errors it is very useful to dump the problem to a file in a human readable form for visual inspection. The line

            # Dump the problem to a human readable OPF file.


does that and it produces a file with the content:

Listing 29 Problem (5) stored in OPF format.
[comment]
Written by MOSEK version 8.1.0.24
Date 11-09-17
Time 14:34:24
[/comment]

[hints]
[hint NUMVAR] 7 [/hint]
[hint NUMCON] 4 [/hint]
[hint NUMANZ] 12 [/hint]
[hint NUMQNZ] 0 [/hint]
[hint NUMCONE] 1 [/hint]
[/hints]

[variables disallow_new_variables]
'x[1]' 'x[2]' 'x[3]' s 't[1]'
't[2]' 't[3]'
[/variables]

[objective maximize]
1.073e-01 'x[1]' + 7.37e-02 'x[2]' + 6.270000000000001e-02 'x[3]'
[/objective]

[constraints]
[con 'budget']  'x[1]' + 'x[2]' + 'x[3]' = 1e+00 [/con]
[con 'GT[1]']  1.667e-01 'x[1]' + 2.32e-02 'x[2]' + 1.3e-03 'x[3]' - 't[1]' = 0e+00 [/con]
[con 'GT[2]']  1.033e-01 'x[2]' - 2.2e-03 'x[3]' - 't[2]' = 0e+00 [/con]
[con 'GT[3]']  3.38e-02 'x[3]' - 't[3]' = 0e+00 [/con]
[/constraints]

[bounds]
[b] 0e+00      <= 'x[1]','x[2]','x[3]' [/b]
[b]               s =  5e-02 [/b]
[b]               't[1]','t[2]','t[3]' free [/b]
[cone quad 'stddev'] s, 't[1]', 't[2]', 't[3]' [/cone]
[/bounds]


Since the API variables have been given meaningful names it is easy to verify by hand that the model is correct.

## 11.1.2 The efficient Frontier¶

The portfolio computed by the Markowitz model is efficient in the sense that there is no other portfolio giving a strictly higher return for the same amount of risk. An efficient portfolio is also sometimes called a Pareto optimal portfolio. Clearly, an investor should only invest in efficient portfolios and therefore it may be relevant to present the investor with all efficient portfolios so the investor can choose the portfolio that has the desired tradeoff between return and risk. This leads to the concept of efficient frontier.

Given a nonnegative $$\alpha$$ the optimization problem

(9)$\begin{split}\begin{array} {lccc} \mbox{maximize} & \mu^T x - \alpha s & & \\ \mbox{subject to} & e^T x & = & w + e^T x^0,\\ & [s;G^T x] & \in & \Q^{n+1},\\ & x & \geq & 0. \end{array}\end{split}$

computes an efficient portfolio which maximizes expected return while minimizing risk, where the tradeoff between the two is controlled by $$\alpha$$. Ideally the problem (9) should be solved for all values $$\alpha \geq 0$$ but in practice that is impossible.

For the example data from Sec. 11.1.1 (A Basic Portfolio Optimization Model), the optimal values of return and risk for a range of $$\alpha$$s are listed below:

Listing 30 Results obtained solving problem (9) for different values of $$\alpha$$.
   alpha         exp ret       std dev
0.000e+000    1.073e-001    7.261e-001
2.500e-001    1.033e-001    1.499e-001
5.000e-001    6.976e-002    3.735e-002
7.500e-001    6.766e-002    3.383e-002
1.000e+000    6.679e-002    3.281e-002
1.500e+000    6.599e-002    3.214e-002
2.000e+000    6.560e-002    3.192e-002
2.500e+000    6.537e-002    3.181e-002
3.000e+000    6.522e-002    3.176e-002
3.500e+000    6.512e-002    3.173e-002
4.000e+000    6.503e-002    3.170e-002
4.500e+000    6.497e-002    3.169e-002


Source code example

The example code in Listing 31 demonstrates how to compute the efficient portfolios for several values of $$\alpha$$. The code is mostly similar to the one in Sec. 11.1.1 (A Basic Portfolio Optimization Model), except the problem is re-optimized in a loop for varying $$\alpha$$.

Listing 31 Code implementing model (9). Click here to download.
import mosek

def streamprinter(text):
print("%s" % text),

if __name__ == '__main__':

n = 3
gamma = 0.05
mu = [0.1073, 0.0737, 0.0627]
GT = [[0.1667, 0.0232, 0.0013],
[0.0000, 0.1033, -0.0022],
[0.0000, 0.0000, 0.0338]]
x0 = [0.0, 0.0, 0.0]
w = 1.0
alphas = [0.0, 0.25, 0.5, 0.75, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5]

inf = 0.0 # This value has no significance

with mosek.Env() as env:

rtemp = w
for j in range(0, n):
rtemp += x0[j]

# Constraints.
task.putconbound(0, mosek.boundkey.fx, rtemp, rtemp)

task.putconboundlist(range(1 + 0, 1 + n), n *
[mosek.boundkey.fx], n * [0.0], n * [0.0])
for j in range(1, 1 + n):
task.putconname(j, "GT[%d]" % j)

# Variables.
task.appendvars(1 + 2 * n)

offsetx = 0   # Offset of variable x into the API variable.
offsets = n   # Offset of variable x into the API variable.
offsett = n + 1 # Offset of variable t into the API variable.

# x variables.
task.putclist(range(offsetx + 0, offsetx + n), mu)
n * [0], range(offsetx + 0, offsetx + n), n * [1.0])
for j in range(0, n):
n * [1 + j], range(offsetx + 0, offsetx + n), GT[j])

range(offsetx + 0, offsetx + n), n * [mosek.boundkey.lo], n * [0.0], n * [inf])
for j in range(0, n):
task.putvarname(offsetx + j, "x[%d]" % (1 + j))

# s variable.
task.putvarbound(offsets + 0, mosek.boundkey.fr, gamma, gamma)
task.putvarname(offsets + 0, "s")

# t variables.
task.putaijlist(range(1, n + 1), range(offsett +
0, offsett + n), n * [-1.0])
task.putvarboundlist(range(offsett + 0, offsett + n),
n * [mosek.boundkey.fr], n * [-inf], n * [inf])
for j in range(0, n):
task.putvarname(offsett + j, "t[%d]" % (1 + j))

offsets] + list(range(offsett, offsett + n)))

# Turn all log output off.

for alpha in alphas:
# Dump the problem to a human readable OPF file.

task.putcj(offsets + 0, -alpha)

# Display the solution summary for quick inspection of results.

if solsta in [mosek.solsta.optimal, mosek.solsta.near_optimal]:
expret = 0.0
x = [0.] * n
offsetx + 0, offsetx + n, x)
for j in range(0, n):
expret += mu[j] * x[j]

stddev = [0.]
offsets + 0, offsets + 1, stddev)

print("\nExpected return %e for gamma %e" %
(expret, stddev[0])),
else:
print("An error occurred when solving for alpha=%e\n" % alpha)


## 11.1.3 Improving the Computational Efficiency¶

In practice it is often important to solve the portfolio problem very quickly. Therefore, in this section we discuss how to improve computational efficiency at the modelling stage.

The computational cost is of course to some extent dependent on the number of constraints and variables in the optimization problem. However, in practice a more important factor is the sparsity: the number of nonzeros used to represent the problem. Indeed it is often better to focus on the number of nonzeros in $$G$$ see (2) and try to reduce that number by for instance changing the choice of $$G$$.

In other words if the computational efficiency should be improved then it is always good idea to start with focusing at the covariance matrix. As an example assume that

$\Sigma = D + VV^T$

where $$D$$ is a positive definite diagonal matrix. Moreover, $$V$$ is a matrix with $$n$$ rows and $$p$$ columns. Such a model for the covariance matrix is called a factor model and usually $$p$$ is much smaller than $$n$$. In practice $$p$$ tends to be a small number independent of $$n$$, say less than 100.

One possible choice for $$G$$ is the Cholesky factorization of $$\Sigma$$ which requires storage proportional to $$n(n+1)/2$$. However, another choice is

$\begin{split}G^T = \left[ \begin{array}{c} D^{1/2}\\ V^T \end{array} \right]\end{split}$

because then

$GG^T = D + VV^T.$

This choice requires storage proportional to $$n+pn$$ which is much less than for the Cholesky choice of $$G$$. Indeed assuming $$p$$ is a constant storage requirements are reduced by a factor of $$n$$.

The example above exploits the so-called factor structure and demonstrates that an alternative choice of $$G$$ may lead to a significant reduction in the amount of storage used to represent the problem. This will in most cases also lead to a significant reduction in the solution time.

The lesson to be learned is that it is important to investigate how the covariance matrix is formed. Given this knowledge it might be possible to make a special choice for $$G$$ that helps reducing the storage requirements and enhance the computational efficiency. More details about this process can be found in [And13].

## 11.1.4 Slippage Cost¶

The basic Markowitz model assumes that there are no costs associated with trading the assets and that the returns of the assets are independent of the amount traded. Neither of those assumptions is usually valid in practice. Therefore, a more realistic model is

(10)$\begin{split}\begin{array}{lrcl} \mbox{maximize} & \mu^T x & & \\ \mbox{subject to} & e^T x + \sum_{j=1}^{n} C_j(x_j - x_j^0) & = & w + e^T x^0,\\ & x^T \Sigma x & \leq & \gamma^2,\\ & x & \geq & 0, \end{array}\end{split}$

where the function

$C_j( x_j - x_j^0)$

specifies the transaction costs when the holding of asset $$j$$ is changed from its initial value.

## 11.1.5 Market Impact Costs¶

If the initial wealth is fairly small and no short selling is allowed, then the holdings will be small and the traded amount of each asset must also be small. Therefore, it is reasonable to assume that the prices of the assets are independent of the amount traded. However, if a large volume of an asset is sold or purchased, the price, and hence return, can be expected to change. This effect is called market impact costs. It is common to assume that the market impact cost for asset $$j$$ can be modelled by

$C_j= m_j \sqrt{|x_j-x_j^0|}$

where $$m_j$$ is a constant that is estimated in some way by the trader. See [GK00] [p. 452] for details. Hence, we have

$C_j(x_j-x_j^0) = m_j |x_j - x_j^0| \sqrt{|x_j-x_j^0|} = m_j |x_j - x_j^0|^{3/2}.$

From [MOSEKApS12] it is known that

$\{(c,z): c \geq z^{3/2}, z\geq 0\} = \{ (c,z): (v,c,z),(z,1/8,v) \in \Q_r^3 \}$

where $$\Q_r^3$$ is the 3-dimensional rotated quadratic cone. Hence, it follows

$\begin{split}\begin{array}{rl} z_j & = |x_j - x_j^0|, \\ (v_j,c_j,z_j),(z_j,1/8,v_j) & \in \Q_r^3,\\ \sum_{j=1}^n C_j(x_j - x_j^0) & = \sum_{j=1}^n c_j. \end{array}\end{split}$

Unfortunately this set of constraints is nonconvex due to the constraint

(11)$z_j = |x_j - x_j^0|$

but in many cases the constraint may be replaced by the relaxed constraint

(12)$z_j \geq |x_j - x_j^0|,$

which is equivalent to

(13)$\begin{split}\begin{array}{l} z_j \geq x_j - x_j^0, \\ z_j \geq -(x_j - x_j^0). \end{array}\end{split}$

For instance if the universe of assets contains a risk free asset then

(14)$z_j > |x_j-x_j^0|$

cannot hold for an optimal solution.

If the optimal solution has the property (14) then the market impact cost within the model is larger than the true market impact cost and hence money are essentially considered garbage and removed by generating transaction costs. This may happen if a portfolio with very small risk is requested because the only way to obtain a small risk is to get rid of some of the assets by generating transaction costs. We generally assume that this is not the case and hence the models (11) and (12) are equivalent.

The above observations lead to

(15)$\begin{split}\begin{array}{lrcll} \mbox{maximize} & \mu^T x & & &\\ \mbox{subject to} & e^T x + m^T c & = & w + e^T x^0, &\\ & [\gamma;G^T x] & \in & \Q^{n+1}, & \\ & z_j & \geq & x_j - x_j^0, & j=1,\ldots,n,\\ & z_j & \geq & x_j^0 - x_j, & j=1,\ldots,n, \\ & [v_j;c_j;z_j],[z_j;1/8;v_j] & \in & \Q_r^3, & j=1,\ldots,n,\\ & x & \geq & 0. & \end{array}\end{split}$

The revised budget constraint

$e^T x + m^Tc = w + e^T x^0$

specifies that the initial wealth covers the investment and the transaction costs. Moreover, $$v$$ and $$z$$ are auxiliary variables that model the market impact cost so that $$z_j \geq | x_j - x_j^0|$$ and $$c_j \geq z_j^{3/2}$$.

It should be mentioned that transaction costs of the form

$c_j \geq z_j^{p/q}$

where $$p$$ and $$q$$ are both integers and $$p\geq q$$ can be modelled using quadratic cones. See [MOSEKApS12] for details.

Creating a matrix formulation

One more reformulation of (15) is needed to bring it to the standard form (4).

(16)$\begin{split}\begin{array} {lcccc} \mbox{maximize} & \mu^T x & & & \\ \mbox{subject to} & e^T x + m^T c & = & w + e^T x^0, &\\ & G^T x - t & = & 0, &\\ & z_j - x_j & \geq & - x_j^0, & j=1,\ldots,n,\\ & z_j + x_j & \geq & x_j^0, & j=1,\ldots,n, \\ & [v_j;c_j;z_j] - [f_{j,1};f_{j,2};f_{j,3}] & = & 0, & j=1,\ldots,n,\\ & [z_j;0;v_j] - [g_{j,1};g_{j,2};g_{j,3}] & = & [0;-1/8;0], & j=1,\ldots,n,\\ & [s;t] & \in & \Q^{n+1}, & \\ & [f_{j,1};f_{j,2};f_{j,3}] & \in & \Q_r^3, & j=1,\ldots,n,\\ & [g_{j,1};g_{j,2};g_{j,3}] & \in & \Q_r^3, & j=1,\ldots,n,\\ & x & \geq & 0, & \\ & s & = & \gamma, & \end{array}\end{split}$

where $$f,g \in \real^{n \times 3}$$. The additional variables $$f$$ and $$g$$ are introduced to ensure that each variable appears at most once in any cone.

The formulation (16) is not the most compact possible, but it is easy to implement. MOSEK presolve will automatically simplify it.

The first step in developing the implementation is to chose an ordering of the variables. We will choose the following ordering:

$\hat{x} = [x; s; t; c; v; z; f; g]$

Table 12 shows the mapping between the $$\hat{x}$$ vector and the model variables.

Table 12 Storage layout for the $$\hat{x}$$
Variable Length Offset
$$x$$ $$n$$ $$0$$
$$s$$ $$1$$ $$n$$
$$t$$ $$n$$ $$n+1$$
$$c$$ $$n$$ $$2n+1$$
$$v$$ $$n$$ $$3n+1$$
$$z$$ $$n$$ $$4n+1$$
$$f(:)^T$$ $$3n$$ $$5n+1$$
$$g(:)^T$$ $$3n$$ $$8n+1$$

The next step is to consider how the linear constraint matrix $$A$$ and the remaining data vectors are laid out. Reusing the idea in Sec. 11.1.1 (A Basic Portfolio Optimization Model) we can write the data in block matrix form and read off all the required coordinates. This extension of the code setting up the constraint $$G^Tx-t=0$$ from Sec. 11.1.1 (A Basic Portfolio Optimization Model) is shown below.

Source code example

The example code in Listing 32 demonstrates how to implement the model (16).

Listing 32 Code implementing model (16). Click here to download.
import mosek

def streamprinter(text):
print("%s" % text),

if __name__ == '__main__':

n = 3
gamma = 0.05
mu = [0.1073, 0.0737, 0.0627]
GT = [[0.1667, 0.0232, 0.0013],
[0.0000, 0.1033, -0.0022],
[0.0000, 0.0000, 0.0338]]
x0 = [0.0, 0.0, 0.0]
w = 1.0
m = [0.01, 0.01, 0.01]

# This value has no significance.
inf = 0.0

with mosek.Env() as env:

rtemp = w
for j in range(0, n):
rtemp += x0[j]

# Constraints.
task.appendcons(1 + 9 * n)
task.putconbound(0, mosek.boundkey.fx, rtemp, rtemp)

task.putconboundlist(range(1 + 0, 1 + n), n *
[mosek.boundkey.fx], n * [0.0], n * [0.0])
for j in range(1, 1 + n):
task.putconname(j, "GT[%d]" % j)

1 + n, 1 + 2 * n), n * [mosek.boundkey.lo], [-x0[j] for j in range(0, n)], n * [inf])
for i in range(0, n):
task.putconname(1 + n + i, "zabs1[%d]" % (1 + i))

task.putconboundlist(range(1 + 2 * n, 1 + 3 * n),
n * [mosek.boundkey.lo], x0, n * [inf])
for i in range(0, n):
task.putconname(1 + 2 * n + i, "zabs2[%d]" % (1 + i))

task.putconboundlist(range(1 + 3 * n, 1 + 3 * n + 3 * n),
3 * n * [mosek.boundkey.fx], 3 * n * [0.], 3 * n * [0.0])
for i in range(0, n):
for k in range(0, n):
task.putconname(1 + 3 * n + 3 * i + k,
"f[%d,%d]" % (1 + i, 1 + k))

task.putconboundlist(range(1 + 6 * n, 1 + 9 * n), 3 * n * [mosek.boundkey.fx],
3 * [0.0, -1.0 / 8.0, 0.0], 3 * [0.0, -1.0 / 8.0, 0.0])
for i in range(0, n):
for k in range(0, n):
task.putconname(1 + 6 * n + 3 * i + k,
"g[%d,%d]" % (1 + i, 1 + k))

# Offset of variables into the API variable.
offsetx = 0
offsets = n
offsett = n + 1
offsetc = 2 * n + 1
offsetv = 3 * n + 1
offsetz = 4 * n + 1
offsetf = 5 * n + 1
offsetg = 8 * n + 1

# Variables.
task.appendvars(1 + 11 * n)

# x variables.
task.putclist(range(offsetx + 0, offsetx + n), mu)
n * [0], range(offsetx + 0, offsetx + n), n * [1.0])
for j in range(0, n):
n * [1 + j], range(offsetx + 0, offsetx + n), GT[j])
task.putaij(1 + n + j, offsetx + j, -1.0)
task.putaij(1 + 2 * n + j, offsetx + j, 1.0)

range(offsetx + 0, offsetx + n), n * [mosek.boundkey.lo], n * [0.0], n * [inf])
for j in range(0, n):
task.putvarname(offsetx + j, "x[%d]" % (1 + j))

# s variable.
task.putvarbound(offsets + 0, mosek.boundkey.fx, gamma, gamma)
task.putvarname(offsets + 0, "s")

# t variables.
task.putaijlist(range(1, n + 1), range(offsett +
0, offsett + n), n * [-1.0])
task.putvarboundlist(range(offsett + 0, offsett + n),
n * [mosek.boundkey.fr], n * [-inf], n * [inf])
for j in range(0, n):
task.putvarname(offsett + j, "t[%d]" % (1 + j))

# c variables.
task.putaijlist(n * [0], range(offsetc, offsetc + n), m)
task.putaijlist(range(1 + 3 * n + 1, 1 + 6 * n + 1, 3),
range(offsetc, offsetc + n), n * [1.0])
task.putvarboundlist(range(offsetc, offsetc + n),
n * [mosek.boundkey.fr], n * [-inf], n * [inf])
for j in range(0, n):
task.putvarname(offsetc + j, "c[%d]" % (1 + j))

# v variables.
task.putaijlist(range(1 + 3 * n + 0, 1 + 6 * n + 0, 3),
range(offsetv, offsetv + n), n * [1.0])
task.putaijlist(range(1 + 6 * n + 2, 1 + 9 * n + 2, 3),
range(offsetv, offsetv + n), n * [1.0])
task.putvarboundlist(range(offsetv, offsetv + n),
n * [mosek.boundkey.fr], n * [-inf], n * [inf])
for j in range(0, n):
task.putvarname(offsetv + j, "v[%d]" % (1 + j))

# z variables.
task.putaijlist(range(1 + 1 * n, 1 + 2 * n),
range(offsetz, offsetz + n), n * [1.0])
task.putaijlist(range(1 + 2 * n, 1 + 3 * n),
range(offsetz, offsetz + n), n * [1.0])
task.putaijlist(range(1 + 3 * n + 2, 1 + 6 * n + 2, 3),
range(offsetz, offsetz + n), n * [1.0])
task.putaijlist(range(1 + 6 * n + 0, 1 + 9 * n + 0, 3),
range(offsetz, offsetz + n), n * [1.0])
task.putvarboundlist(range(offsetz, offsetz + n),
n * [mosek.boundkey.fr], n * [-inf], n * [inf])
for j in range(0, n):
task.putvarname(offsetz + j, "z[%d]" % (1 + j))

# f variables.
for j in range(0, n):
for k in range(0, n):
task.putaij(1 + 3 * n + 3 * j + k,
offsetf + 3 * j + k, -1.0)
task.putvarbound(offsetf + 3 * j + k,
mosek.boundkey.fr, -inf, inf)
task.putvarname(offsetf + 3 * j + k,
"f[%d,%d]" % (1 + j, 1 + k))

# g variables.
for j in range(0, n):
for k in range(0, n):
task.putaij(1 + 6 * n + 3 * j + k,
offsetg + 3 * j + k, -1.0)
task.putvarbound(offsetg + 3 * j + k,
mosek.boundkey.fr, -inf, inf)
task.putvarname(offsetg + 3 * j + k,
"g[%d,%d]" % (1 + j, 1 + k))

offsets] + list(range(offsett, offsett + n)))

for k in range(0, n):
0.0, 3, offsetf + 3 * k)
task.putconename(1 + k, "f[%d]" % (1 + k))

for k in range(0, n):
0.0, 3, offsetg + 3 * k)
task.putconename(1 + n + k, "g[%d]" % (1 + k))

# Turn all log output off.

# Dump the problem to a human readable OPF file.

# Display the solution summary for quick inspection of results.

expret = 0.0
x = [0.] * n
task.getxxslice(mosek.soltype.itr, offsetx + 0, offsetx + n, x)
for j in range(0, n):
expret += mu[j] * x[j]

stddev = [0.]
0, offsets + 1, stddev)

print("\nExpected return %e for gamma %e\n" % (expret, stddev[0]))


The example code above produces the result

Interior-point solution summary
Problem status  : PRIMAL_AND_DUAL_FEASIBLE
Solution status : OPTIMAL
Primal.  obj: 7.4390660847e-02    nrm: 1e+00    Viol.  con: 6e-09    var: 0e+00    cones: 4e-09
Dual.    obj: 7.4390675795e-02    nrm: 3e-01    Viol.  con: 1e-19    var: 8e-09    cones: 0e+00

Expected return 7.439066e-02 for gamma 5.000000e-02


If the problem is dumped to an OPF file, it has the following content.

Listing 33 OPF file for problem (16).
[comment]
Written by MOSEK version 8.1.0.24
Date 12-09-17
Time 12:34:27
[/comment]

[hints]
[hint NUMVAR] 34 [/hint]
[hint NUMCON] 28 [/hint]
[hint NUMANZ] 60 [/hint]
[hint NUMQNZ] 0 [/hint]
[hint NUMCONE] 7 [/hint]
[/hints]

[variables disallow_new_variables]
'x[1]' 'x[2]' 'x[3]' s 't[1]'
't[2]' 't[3]' 'c[1]' 'c[2]' 'c[3]'
'v[1]' 'v[2]' 'v[3]' 'z[1]' 'z[2]'
'z[3]' 'f[1,1]' 'f[1,2]' 'f[1,3]' 'f[2,1]'
'f[2,2]' 'f[2,3]' 'f[3,1]' 'f[3,2]' 'f[3,3]'
'g[1,1]' 'g[1,2]' 'g[1,3]' 'g[2,1]' 'g[2,2]'
'g[2,3]' 'g[3,1]' 'g[3,2]' 'g[3,3]'
[/variables]

[objective maximize]
1.073e-01 'x[1]' + 7.37e-02 'x[2]' + 6.270000000000001e-02 'x[3]'
[/objective]

[constraints]
[con 'budget']  'x[1]' + 'x[2]' + 'x[3]' + 1e-02 'c[1]' + 1e-02 'c[2]'
+ 1e-02 'c[3]' = 1e+00 [/con]
[con 'GT[1]']  1.667e-01 'x[1]' + 2.32e-02 'x[2]' + 1.3e-03 'x[3]' - 't[1]' = 0e+00 [/con]
[con 'GT[2]']  1.033e-01 'x[2]' - 2.2e-03 'x[3]' - 't[2]' = 0e+00 [/con]
[con 'GT[3]']  3.38e-02 'x[3]' - 't[3]' = 0e+00 [/con]
[con 'zabs1[1]'] 0e+00 <= - 'x[1]' + 'z[1]' [/con]
[con 'zabs1[2]'] 0e+00 <= - 'x[2]' + 'z[2]' [/con]
[con 'zabs1[3]'] 0e+00 <= - 'x[3]' + 'z[3]' [/con]
[con 'zabs2[1]'] 0e+00 <= 'x[1]' + 'z[1]' [/con]
[con 'zabs2[2]'] 0e+00 <= 'x[2]' + 'z[2]' [/con]
[con 'zabs2[3]'] 0e+00 <= 'x[3]' + 'z[3]' [/con]
[con 'f[1,1]']  'v[1]' - 'f[1,1]' = 0e+00 [/con]
[con 'f[1,2]']  'c[1]' - 'f[1,2]' = 0e+00 [/con]
[con 'f[1,3]']  'z[1]' - 'f[1,3]' = 0e+00 [/con]
[con 'f[2,1]']  'v[2]' - 'f[2,1]' = 0e+00 [/con]
[con 'f[2,2]']  'c[2]' - 'f[2,2]' = 0e+00 [/con]
[con 'f[2,3]']  'z[2]' - 'f[2,3]' = 0e+00 [/con]
[con 'f[3,1]']  'v[3]' - 'f[3,1]' = 0e+00 [/con]
[con 'f[3,2]']  'c[3]' - 'f[3,2]' = 0e+00 [/con]
[con 'f[3,3]']  'z[3]' - 'f[3,3]' = 0e+00 [/con]
[con 'g[1,1]']  'z[1]' - 'g[1,1]' = 0e+00 [/con]
[con 'g[1,2]']  - 'g[1,2]' = -1.25e-01 [/con]
[con 'g[1,3]']  'v[1]' - 'g[1,3]' = 0e+00 [/con]
[con 'g[2,1]']  'z[2]' - 'g[2,1]' = 0e+00 [/con]
[con 'g[2,2]']  - 'g[2,2]' = -1.25e-01 [/con]
[con 'g[2,3]']  'v[2]' - 'g[2,3]' = 0e+00 [/con]
[con 'g[3,1]']  'z[3]' - 'g[3,1]' = 0e+00 [/con]
[con 'g[3,2]']  - 'g[3,2]' = -1.25e-01 [/con]
[con 'g[3,3]']  'v[3]' - 'g[3,3]' = 0e+00 [/con]
[/constraints]

[bounds]
[b] 0e+00      <= 'x[1]','x[2]','x[3]' [/b]
[b]               s =  5e-02 [/b]
[b]               't[1]','t[2]','t[3]','c[1]','c[2]','c[3]' free [/b]
[b]               'v[1]','v[2]','v[3]','z[1]','z[2]','z[3]' free [/b]
[b]               'f[1,1]','f[1,2]','f[1,3]','f[2,1]','f[2,2]','f[2,3]' free [/b]
[b]               'f[3,1]','f[3,2]','f[3,3]','g[1,1]','g[1,2]','g[1,3]' free [/b]
[b]               'g[2,1]','g[2,2]','g[2,3]','g[3,1]','g[3,2]','g[3,3]' free [/b]
[cone quad 'stddev'] s, 't[1]', 't[2]', 't[3]' [/cone]
[cone rquad 'f[1]'] 'f[1,1]', 'f[1,2]', 'f[1,3]' [/cone]
[cone rquad 'f[2]'] 'f[2,1]', 'f[2,2]', 'f[2,3]' [/cone]
[cone rquad 'f[3]'] 'f[3,1]', 'f[3,2]', 'f[3,3]' [/cone]
[cone rquad 'g[1]'] 'g[1,1]', 'g[1,2]', 'g[1,3]' [/cone]
[cone rquad 'g[2]'] 'g[2,1]', 'g[2,2]', 'g[2,3]' [/cone]
[cone rquad 'g[3]'] 'g[3,1]', 'g[3,2]', 'g[3,3]' [/cone]
[/bounds]


The file verifies that the correct problem has been set up.