# 7.1 Portfolio Optimization¶

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

Subsequently the following MATLAB inspired notation will be employed. The : operator is used as follows

$i:j = \{ i,i+1,\ldots,j \}$

and hence

$\begin{split}x_{2:4} = \left[ \begin{array}{c} x_2\\ x_3\\ x_4 \end{array} \right]\end{split}$

If $$x$$ and $$y$$ are two column vectors, then

$\begin{split}[x;y] = \left[ \begin{array} {c} x\\ y \end{array} \right]\end{split}$

Furthermore, if $$f \in \real^{m \times n}$$ then

$\begin{split}f(:) = \left[ \begin{array} {c} f_{1,1}\\ f_{2,1}\\ \vdots\\ f_{m-1,n}\\ \vdots\\ f_{m,n} \end{array} \right]\end{split}$

i.e. $$f(:)$$ stacks the columns of the matrix $$f$$.

## 7.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} {lccc} \mbox{maximize} & \mu^T x & & \\ \mbox{subject to} & 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 interpretated 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 Section 7.1.3.

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} {lcll} \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.

Subsequently we will use the example data

$\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 figures of accuracy. 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.$

The data has been taken from [CT07].

### 7.1.1.1 Why a Conic Formulation?¶

Problem (1) is a convex quadratically constrained optimization problems that can be solved directly using MOSEK, then why reformulate it as a conic quadratic optimization problem? The main reason for choosing a conic model is that it is more robust and usually leads to a shorter solution times. For instance it is not always easy to determine whether the $$Q$$ matrix in (1) is positive semidefinite due to the presence of rounding errors. It is also very easy to make a mistake so $$Q$$ becomes indefinite. These causes of problems are completely eliminated in the conic formulation.

Moreover, observe the constraint

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

is nicer than

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

for small and values of $$\gamma$$. For instance assume a $$\gamma$$ of 10000 then $$\gamma^2$$ would 1.0e8 which introduces a scaling issue in the model. Hence, using conic formulation it is possible to work with the standard deviation instead of the variance, which usually gives rise to a better scaled model.

### 7.1.1.2 Implementing the Portfolio Model¶

Model (3) can not be implemented as stated using the MOSEK optimizer API because the API requires the problem to be on the form

(4)$\begin{split}\begin{array} {lccccc} \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}$

where $$\hat{x}$$ is referred to as the API variable.

The first step in bringing (3) to the form (4) is the reformulation

(5)$\begin{split}\begin{array} {lccc} \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 & 0. & \end{array}\end{split}$

where $$s$$ is an additional scalar variable and $$t$$ is a n dimensional vector variable. The next step is to define a mapping of the variables

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

Hence, the API variable $$\hat{x}$$ is concatenation of model variables $$x$$, $$s$$ and $$t$$. In Table 7.1 the details of the concatenation are specified.

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

For instance it can be seen that

$\hat{x}_{n+2} = t_1.$

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

Given the ordering of the variables specified by (6) the data should be defined as follows

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

The next step is to consider how the columns of $$A$$ is defined. The following pseudo code

$\begin{split}\begin{array} {ll} for & j=1:n \\ & \hat{x}_j = x_j \\ & A_{1,j} = 1.0\\ & A_{2:(n+1),j} = G_{j,1:n}^T\\ &\\ \hat{x}_{n+1} = s &\\ &\\ for & j=1:n\\ & \hat{x}_{n+1+j} = t_j\\ & A_{n+1+j,n+1+j} = -1.0 \end{array}\end{split}$

show how to construct each column of $$A$$.

In the above discussion index origin 1 is employed, i.e., the first position in a vector is 1. The C programming language employs 0 as index origin and that should be kept in mind when reading the example code.

#include <math.h>
#include <stdio.h>

#include "mosek.h"

#define MOSEKCALL(_r,_call)  if ( (_r)==MSK_RES_OK ) (_r) = (_call)

static void MSKAPI printstr(void *handle,
const char str[])
{
printf("%s",str);
} /* printstr */

int main(int argc, const char **argv)
{
char            buf[128];

double          expret  = 0.0,
stddev  = 0.0,
xj;

const MSKint32t n       = 3;
const MSKrealt  gamma   = 0.05;
const MSKrealt  mu[]    = {0.1073,  0.0737,  0.0627};
const MSKrealt  GT[][3] = {
{0.1667,  0.0232,  0.0013},
{0.0000,  0.1033, -0.0022},
{0.0000,  0.0000,  0.0338}};
const MSKrealt  x0[3]   = {0.0, 0.0, 0.0};
const MSKrealt  w       = 1.0;
MSKrealt        rtemp;
MSKenv_t        env;
MSKint32t       k,i,j,offsetx,offsets,offsett,*sub;
MSKrescodee     res=MSK_RES_OK;

/* Initial setup. */
env  = NULL;
MOSEKCALL(res,MSK_makeenv(&env,NULL));

rtemp = w;
for(j=0; j<n; ++j)
rtemp += x0[j];

/* Constraints. */
sprintf(buf,"%s","budget");

for(i=0; i<n; ++i)
{
sprintf(buf,"GT[%d]",1+i);
}

/* Variables. */

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. */
for(j=0; j<n; ++j)
{
for(k=0; k<n; ++k)
if( GT[k][j]!=0.0 )
MOSEKCALL(res, MSK_putaij(task, 1 + k, offsetx + j, GT[k][j]));

sprintf(buf,"x[%d]",1+j);
}

/* s variable. */
sprintf(buf,"s");

/* t variables. */
for(j=0; j<n; ++j)
{
sprintf(buf,"t[%d]",1+j);
}

if ( res==MSK_RES_OK )
{
/* sub should be n+1 long i.e. the dimmension of the cone. */

if ( sub )
{
sub[0] = offsets+0;
for(j=0; j<n; ++j)
sub[j+1] = offsett+j;

}
else
res = MSK_RES_ERR_SPACE;
}

#if 1
/* No log output. */
#else
#endif

/* Dump the problem to a human readable OPF file. */
#if 0
#endif

#if 0
/* Display the solution summary for quick inspection of results. */

#endif

if ( res==MSK_RES_OK )
{
expret=0.0;
stddev=0.0;

for(j=0; j<n; ++j)
{
expret += mu[j]*xj;
}

printf("\nExpected return %e for gamma %e\n",expret,stddev);
}

if ( env!=NULL )
MSK_deleteenv(&env);

return ( 0 );
}


The above code produce the result

Listing 7.1 Output from the solver.
Interior-point solution summary
Problem status  : PRIMAL_AND_DUAL_FEASIBLE
Solution status : OPTIMAL
Primal.  obj: 7.4766497707e-002   Viol.  con: 2e-008   var: 0e+000   cones: 3e-009
Dual.    obj: 7.4766522618e-002   Viol.  con: 0e+000   var: 4e-008   cones: 0e+000

Expected return 7.476650e-002 for gamma 5.000000e-002


The source code should be self-explanatory but a few comments are nevertheless in place.

The code use a macro called MOSEKCALL which is defined as follows

#define MOSEKCALL(_r,_call)  ( (_r)==MSK_RES_OK ? ( (_r) = (_call) ) : ( (_r) = (_r) ) );


so for instance

MOSEKCALL(res,MSK_optimize());


is the same as

if ( res==MSK_RES_OK )
res = MSK_optimize()


so the usage of MOSEKCALL is method for hiding if statements and hence making the code more compact.

In the lines

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. */


offsets into the MOSEK API variables are stored and those offsets are used later. The code

for(j=0; j<n; ++j)
{
for(k=0; k<n; ++k)
if( GT[k][j]!=0.0 )
MOSEKCALL(res, MSK_putaij(task, 1 + k, offsetx + j, GT[k][j]));

sprintf(buf,"x[%d]",1+j);
}


sets up the data for x variables. For instance

MOSEKCALL(res,MSK_putcj(task,offsetx+j,mu[j]));


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

sprintf(buf,"x[%d]",1+j);


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

### 7.1.1.3 Debugging Tips¶

Implementing an optimization model in optimizer can be cumbersome and error-prone and it is very easy to make mistakes. In order to check the implemented code for mistakes it is very useful to dump the problem to a file in a human readable form for visual inspection. The line

MOSEKCALL(res,MSK_writedata(task,"dump.opf"));


does that and this will produce a file with the content

Listing 7.2 Problem (5) stored in OPF format.
[comment]
Written by MOSEK version 7.0.0.86
Date 01-10-13
Time 08:43:21
[/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-001 'x[1]' + 7.37e-002 'x[2]' + 6.270000000000001e-002 'x[3]'
[/objective]

[constraints]
[con 'budget']  'x[1]' + 'x[2]' + 'x[3]' = 1e+000 [/con]
[con 'GT[1]']  1.667e-001 'x[1]' + 2.32e-002 'x[2]' + 1.3e-003 'x[3]' - 't[1]' = 0e+000 [/con]
[con 'GT[2]']  1.033e-001 'x[2]' - 2.2e-003 'x[3]' - 't[2]' = 0e+000 [/con]
[con 'GT[3]']  3.38e-002 'x[3]' - 't[3]' = 0e+000 [/con]
[/constraints]

[bounds]
[b]          0 <= * [/b]
[b]               s =  5e-002 [/b]
[b]               't[1]','t[2]','t[3]' free [/b]
[cone quad 'stddev'] s, 't[1]', 't[2]', 't[3]' [/cone]
[/bounds]


Observe that since the API variables have been given meaningful names it is easy to see the model is correct.

## 7.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.

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

(7)$\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 efficient portfolios. Note that the objective maximizes the expected return while maximizing $$-\alpha$$ times the standard deviation. Hence, the standard deviation is minimized while $$\alpha$$ specifies the tradeoff between expected return and risk.

Ideally the problem (7) should be solved for all values $$\alpha \geq 0$$ but in practice that is computationally too costly.

Using the example data from Section 7.1.1, the optimal values of return and risk for several $$\alpha$$s are listed below:

Listing 7.3 Results obtained solving problem (7) 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


The example code in Listing 7.4 demonstrates how to compute the efficient portfolios for several values of $$\alpha$$.

Listing 7.4 Code implementing model (7) Click here to download.
#include <math.h>
#include <stdio.h>

#include "mosek.h"

#define MOSEKCALL(_r,_call)  if ( (_r)==MSK_RES_OK ) (_r) = (_call)

static void MSKAPI printstr(void *handle,
const char str[])
{
printf("%s",str);
} /* printstr */

int main(int argc, const char **argv)
{
char            buf[128];
const MSKint32t n          = 3,
numalpha   = 12;
const double    mu[]       = {0.1073,  0.0737,  0.0627},
x0[3]      = {0.0, 0.0, 0.0},
w          = 1.0,
alphas[12] = {0.0, 0.25, 0.5, 0.75, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5},
GT[][3]    = {
{0.1667,  0.0232,  0.0013},
{0.0000,  0.1033, -0.0022},
{0.0000,  0.0000,  0.0338}
};
double          expret,
stddev,
alpha;
MSKenv_t        env;
MSKint32t       k,i,j,offsetx,offsets,offsett;
MSKrescodee     res=MSK_RES_OK;
MSKrealt        xj;
MSKsolstae      solsta;

/* Initial setup. */
env  = NULL;
MOSEKCALL(res,MSK_makeenv(&env,NULL));

/* Constraints. */
sprintf(buf,"%s","budget");

for(i=0; i<n; ++i)
{
sprintf(buf,"GT[%d]",1+i);
}

/* Variables. */

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. */
for(j=0; j<n; ++j)
{
for(k=0; k<n; ++k)
if( GT[k][j]!=0.0 )
MOSEKCALL(res, MSK_putaij(task, 1 + k, offsetx + j, GT[k][j]));

sprintf(buf,"x[%d]",1+j);
}

/* s variable. */
sprintf(buf,"s");

/* t variables. */
for(j=0; j<n; ++j)
{
sprintf(buf,"t[%d]",1+j);
}

if ( res==MSK_RES_OK )
{
/* sub should be n+1 long i.e. the dimmension of the cone. */

if ( sub )
{
sub[0] = offsets+0;
for(j=0; j<n; ++j)
sub[j+1] = offsett+j;

}
else
res = MSK_RES_ERR_SPACE;
}

/* Turn all log output off. */

printf("%-12s  %-12s  %-12s\n","alpha","exp ret","std dev");

for(k=0; k<numalpha; ++k)
{
alpha = alphas[k];

/* Sets the objective function coefficient for s. */

if( solsta==MSK_SOL_STA_OPTIMAL || solsta==MSK_SOL_STA_NEAR_OPTIMAL )
{
expret = 0.0;
for(j=0; j<n; ++j)
{
expret += mu[j]*xj;
}

printf("%-12.3e  %-12.3e  %-12.3e\n",alpha,expret,stddev);
}
else
{
printf("An error occurred when solving for alpha=%e\n",alpha);
}
}

MSK_deleteenv(&env);

return ( 0 );
}


## 7.1.3 Improving the Computational Efficiency¶

In practice it is often important to solve the portfolio problem in a short amount of time; this section it is discusses what can be done at the modelling stage to improve the computational efficiency.

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 number nonzeros used to represent the problem. Indeed it is often better to focus at 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 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 factor model and usually $$p$$ tends be a small number, say less than 100, independent of $$n$$.

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 then the difference in storage requirements is 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 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.

## 7.1.4 Slippage Cost¶

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

(8)$\begin{split}\begin{array} {lccc} \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.

### 7.1.4.1 Market Impact Costs¶

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

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

where $$m_j$$ is a constant that is estimated in some way. See [GK00][p. 452] for details. To summarize then

$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

$\{(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 implying

$\begin{split}\begin{array} {lcl} z_j & = & |x_j - x_j^0|, \\ \left[ v_j ; c_j ; z_j \right], \left[z_j ; 1/8; v_j \right] & \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

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

but in many cases that constraint can safely be replaced by the relaxed constraint

(10)$z_j \geq |x_j - x_j^0|$

which is convex. If for instance the universe of assets contains a risk free asset with a positive return then

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

cannot hold for an optimal solution because that would imply the solution is not optimal.

Now assume that the optimal solution has the property that (11) holds 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 then the only way to obtain a small risk is to get rid of some of the assets by generating transaction costs. Here it is assumed this is not the case and hence the models (9) and (10) are equivalent.

Formula (10) is replaced by constraints

(12)$\begin{split}\begin{array} {ccc} z_j & \geq & x_j - x_j^0, \\ z_j & \geq & -(x_j - x_j^0). \end{array}\end{split}$

Now we have

(13)$\begin{split}\begin{array} {lcccl} \mbox{maximize} & \mu^T x & & & \\ \mbox{subject to} & e^T x + m^T c & = & w + e^T x^0, & \\ & z_j & \geq & x_j - x_j^0, & j=1,\ldots,n,\\ & z_j & \geq & x_j^0 - x_j, & j=1,\ldots,n, \\ & [\gamma;G^T x] & \in & \Q^{n+1}, & \\ & [v_j;c_j;z_j ] & \in & \Q_r^3, & j=1,\ldots,n,\\ & [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 = w + e^T x^0 - m^T c$

specifies that the total investment must be equal to the initial wealth minus the transaction costs. Moreover, observe the variables $$v$$ and $$z$$ are some auxiliary variables that model the market impact cost. Indeed it holds

$z_j \geq | x_j - x_j^0|$

and

$c_j \geq z_j^{3/2}.$

Before proceeding 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.

One more reformulation of (13) is needed,

(14)$\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:3} & = & 0, & j=1,\ldots,n,\\ & [z_j;0;v_j] - g_{j,1:3} & = & [0;-1/8;0], & j=1,\ldots,n,\\ & [s;t] & \in & \Q^{n+1}, & \\ & f_{j,1:3}^T & \in & \Q_r^3, & j=1,\ldots,n,\\ & g_{j,1:3}^T & \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}$$. These additional variables $$f$$ and $$g$$ are only introduced to bring the problem on the API standard form.

The formulation (14)) is not the most compact possible. However, the MOSEK presolve will automatically make it more compact and since it is easier to implement (14) than a more compact form then the form (14) is preferred.

The first step in developing the optimizer API implementation is to chose an ordering of the variables. In this case the ordering

$\begin{split}\hat{x} = \left[ \begin{array} {l} x \\ s\\ t\\ c\\ v\\ z\\ f^T(:)\\ g^T(:) \end{array} \right]\end{split}$

will be used. Note $$f^T(:)$$ means the rows of $$f$$ are transposed and stacked on top of each other to form a long column vector. Table 7.2 shows the mapping between the $$\hat{x}$$ and the model variables.

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

The next step is to consider how the columns of $$A$$ is defined. Reusing the idea in Section 7.1.1 then the following pseudo code describes the setup of $$A$$.

$\begin{split}\begin{array} {ll} for & j=1:n \\ & \hat{x}_j = x_j \\ & A_{1,j} = 1.0\\ & A_{2:n+1,j} = G_{j,1:n}^T\\ & A_{n+1+j,j} = -1.0\\ & A_{2n+1+j,j} = 1.0\\ &\\ \hat{x}_{n+1} = s &\\ &\\ for & j=1:n\\ & \hat{x}_{n+1+j} = t_j\\ & A_{1+j,n+1+j} = -1.0\\ &\\ for & j=1:n\\ & \hat{x}_{2n+1+j} = c_j \\ & A_{1,2n+1+j} = m_j\\ & A_{3n+1+3(j-1)+2,2n+1+j} = 1.0\\ &\\ for & j=1:n\\ & \hat{x}_{3n+1+j} = v_j \\ & A_{3n+1+3(j-1)+1,3n+1+j} = 1.0\\ & A_{6n+1+3(j-1)+3,3n+1+j} = 1.0\\ &\\ for & j=1:n\\ & \hat{x}_{4n+1+j} = z_j \\ & A_{1+n+j,4n+1+j} = 1.0\\ & A_{1+2n+j,4n+1+j} = 1.0\\ & A_{3n+1+3(j-1)+3,4n+1+j} = 1.0\\ & A_{6n+1+3(j-1)+1,4n+1+j} = 1.0\\ &\\ for & j=1:n\\ & \hat{x}_{7n+1+3(j-1)+1} = f_{j,1} \\ & A_{3n+1+3(j-1)+1,7n+(3(j-1)+1} = -1.0\\ & \hat{x}_{7n+1+3(j-1)+2} = f_{j,2} \\ & A_{3n+1+3(j-1)+2,7n+(3(j-1)+2} = -1.0\\ & \hat{x}_{7n+1+3(j-1)+3} = f_{j,3} \\ & A_{3n+1+3(j-1)+3,7n+(3(j-1)+3} = -1.0\\ & \\ for & j=1:n\\ & \hat{x}_{10n+1+3(j-1)+1} = g_{j,1} \\ & A_{6n+1+3(j-1)+1,7n+(3(j-1)+1} = -1.0\\ & \hat{x}_{10n+1+3(j-1)+2} = g_{j,2} \\ & A_{6n+1+3(j-1)+2,7n+(3(j-1)+2} = -1.0\\ & \hat{x}_{10n+1+3(j-1)+3} = g_{j,3} \\ & A_{6n+1+3(j-1)+3,7n+(3(j-1)+3} = -1.0 \end{array}\end{split}$

The example code in Listing 7.5 demonstrates how to implement the model (14).

Listing 7.5 Code implementing model (14). Click here to download.
#include <math.h>
#include <stdio.h>

#include "mosek.h"

#define MOSEKCALL(_r,_call)  if ( (_r)==MSK_RES_OK ) (_r) = (_call)

static void MSKAPI printstr(void *handle,
const char str[])
{
printf("%s",str);
} /* printstr */

int main(int argc, const char **argv)
{
char            buf[128];
const MSKint32t n       = 3;
const double    w       = 1.0,
x0[]    = {0.0, 0.0, 0.0},
gamma   = 0.05,
mu[]    = {0.1073,  0.0737,  0.0627},
m[]     = {0.01, 0.01, 0.01},
GT[][3] = {{0.1667,  0.0232,  0.0013},
{0.0000,  0.1033, -0.0022},
{0.0000,  0.0000,  0.0338}};
double          b[3]    = {0.0, -1.0/8.0, 0.0};
double          rtemp,
expret,
stddev,
xj;
MSKenv_t        env;
MSKint32t       k,i,j,
offsetx,offsets,offsett,offsetc,
offsetv,offsetz,offsetf,offsetg;
MSKrescodee     res=MSK_RES_OK;

/* Initial setup. */
env  = NULL;
MOSEKCALL(res,MSK_makeenv(&env,NULL));

rtemp = w;
for(k=0; k<n; ++k)
rtemp += x0[k];

/* Constraints. */
sprintf(buf,"%s","budget");

for(i=0; i<n; ++i)
{
sprintf(buf,"GT[%d]",1+i);
}

for(i=0; i<n; ++i)
{
sprintf(buf,"zabs1[%d]",1+i);
}

for(i=0; i<n; ++i)
{
sprintf(buf,"zabs2[%d]",1+i);
}

for(i=0; i<n; ++i)
{
for(k=0; k<3; ++k)
{
sprintf(buf,"f[%d,%d]",1+i,1+k);
}
}

for(i=0; i<n; ++i)
{
for(k=0; k<3; ++k)
{
sprintf(buf,"g[%d,%d]",1+i,1+k);
}
}

/* Offsets of variables into the (serialized) 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. */

/* x variables. */
for(j=0; j<n; ++j)
{
for(k=0; k<n; ++k)
if( GT[k][j]!=0.0 )

sprintf(buf,"x[%d]",1+j);
}

/* s variable. */
sprintf(buf,"s");

/* t variables. */
for(j=0; j<n; ++j)
{
sprintf(buf,"t[%d]",1+j);
}

/* c variables. */
for(j=0; j<n; ++j)
{
sprintf(buf,"c[%d]",1+j);
}

/* v variables. */
for(j=0; j<n; ++j)
{
sprintf(buf,"v[%d]",1+j);
}

/* z variables. */
for(j=0; j<n; ++j)
{
sprintf(buf,"z[%d]",1+j);
}

/* f variables. */
for(j=0; j<n; ++j)
{
for(k=0; k<3; ++k)
{
sprintf(buf,"f[%d,%d]",1+j,1+k);
}
}

/* g variables. */
for(j=0; j<n; ++j)
{
for(k=0; k<3; ++k)
{
sprintf(buf,"g[%d,%d]",1+j,1+k);
}
}
if ( res==MSK_RES_OK )
{
/* sub should be n+1 long i.e. the dimmension of the cone. */

if ( sub )
{
sub[0] = offsets+0;
for(j=0; j<n; ++j)
sub[j+1] = offsett+j;

for(k=0; k<n; ++k)
{
sprintf(buf,"f[%d]",1+k);
}

for(k=0; k<n; ++k)
{
sprintf(buf,"g[%d]",1+k);
}

}
else
res = MSK_RES_ERR_SPACE;
}

#if 1
/* no log output. */
#else
#endif

#if 0
/* Dump the problem to a human readable OPF file. */
#endif

/* Display the solution summary for quick inspection of results. */
#if 1
#endif

if ( res==MSK_RES_OK )
{
expret=0.0;
stddev=0.0;

for(j=0; j<n; ++j)
{
expret += mu[j]*xj;
}

printf("\nExpected return %e for gamma %e\n",expret,stddev);
}

MSK_deleteenv(&env);

return ( 0 );
}


The example code above produces the result

Interior-point solution summary
Problem status  : PRIMAL_AND_DUAL_FEASIBLE
Solution status : OPTIMAL
Primal.  obj: 7.4390654948e-002   Viol.  con: 2e-007   var: 0e+000   cones: 2e-009
Dual.    obj: 7.4390665143e-002   Viol.  con: 2e-008   var: 2e-008   cones: 0e+000
Expected return 7,439065E-002 for gamma 5,000000E-002


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

Listing 7.6 OPF file for problem (14).
[comment]
Written by MOSEK version 7.0.0.86
Date 01-10-13
Time 08:59:30
[/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-001 'x[1]' + 7.37e-002 'x[2]' + 6.270000000000001e-002 'x[3]'
[/objective]

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

[bounds]
[b]          0 <= * [/b]
[b]               s =  5e-002 [/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 setup.