11.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\).

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} {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 11.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].

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

11.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 11.1 the details of the concatenation are specified.

Table 11.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 Java programming language employs 0 as index origin and that should be kept in mind when reading the example code.

Listing 11.1 Code implementing model (3). Click here to download.
package com.mosek.example;

public class case_portfolio_1 {
  static final int n = 3;

  public static void main (String[] args) {
    // Since the value infinity is never used, we define
    // 'infinity' symbolic purposes only
    double infinity = 0;
    double gamma = 0.05;
    double[]   mu = {0.1073,  0.0737,  0.0627};
    double[][] GT = {
      {0.1667,  0.0232,  0.0013},
      {0.0000,  0.1033, -0.0022},
      {0.0000,  0.0000,  0.0338}
    };
    double[] x0 = {0.0, 0.0, 0.0};
    double   w = 1.0;

    int numvar = 2 * n + 1;
    int numcon = n + 1;


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

    try ( mosek.Env env   = new mosek.Env ();
          mosek.Task task = new mosek.Task (env, 0, 0) ) {
      // Directs the log task stream to the user specified
      // method task_msg_obj.stream
      task.set_Stream(
        mosek.streamtype.log,
        new mosek.Stream()
      { public void stream(String msg) { System.out.print(msg); }});

      //Constraints.
      task.appendcons(numcon);
      for ( int i = 1; i <= n; ++i) {
        w += x0[i - 1];
        task.putconbound(i, mosek.boundkey.fx, 0., 0.);
        task.putconname(i, "GT[" + i + "]");
      }
      task.putconbound(0, mosek.boundkey.fx, w, w);
      task.putconname(0, "budget");

      //Variables.
      task.appendvars(numvar);

      int[] xindx = {offsetx + 0, offsetx + 1, offsetx + 2};
      task.putclist(xindx, mu);

      for (int i = 0; i < n; ++i) {
        for (int j = i; j < n; ++j)
          task.putaij(i + 1, offsetx + j, GT[i][j]);

        task.putaij(i + 1, offsett + i, -1.0);

        task.putvarbound(offsetx + i, mosek.boundkey.lo, 0., 0.);

        task.putvarname(offsetx + i, "x[" + (i + 1) + "]");
        task.putvarname(offsett + i, "t[" + (i + 1) + "]");
        task.putvarbound(offsett + i, mosek.boundkey.fr, 0., 0.);
      }
      task.putvarbound(offsets, mosek.boundkey.fx, gamma, gamma);
      task.putvarname(offsets, "s");

      double[] e = {1.0, 1.0, 1.0};
      task.putarow(0, xindx, e);

      //Cones.
      int[] csub = {offsets, offsett + 0, offsett + 1, offsett + 2};
      task.appendcone( mosek.conetype.quad,
                       0.0, /* For future use only, can be set to 0.0 */
                       csub);
      task.putconename(0, "stddev");

      /* A maximization problem */
      task.putobjsense(mosek.objsense.maximize);

      //task.solutionsummary(mosek.streamtype.log);

      //task.writedata("dump.opf");

      /* Solve the problem */
      try {

        task.optimize();

        double expret = 0.0, stddev = 0.0;
        double[] xx = new double[numvar];

        task.getxx(mosek.soltype.itr, xx);

        for (int j = 0; j < n; ++j)
          expret += mu[j] * xx[j + offsetx];

        System.out.printf("\nExpected return %e for gamma %e\n", expret, xx[offsets]);

      } catch (mosek.Warning mw) {
        System.out.println (" Mosek warning:");
        System.out.println (mw.toString ());
      }
    } catch ( mosek.Exception e) {
      System.out.println ("An error/warning was encountered");
      System.out.println (e.toString ());
      throw e;
    }
  }
}

The above code produce the result

Listing 11.2 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.

In the lines

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

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

      int[] xindx = {offsetx + 0, offsetx + 1, offsetx + 2};
      task.putclist(xindx, mu);

      for (int i = 0; i < n; ++i) {
        for (int j = i; j < n; ++j)
          task.putaij(i + 1, offsetx + j, GT[i][j]);

        task.putaij(i + 1, offsett + i, -1.0);

        task.putvarbound(offsetx + i, mosek.boundkey.lo, 0., 0.);

        task.putvarname(offsetx + i, "x[" + (i + 1) + "]");
        task.putvarname(offsett + i, "t[" + (i + 1) + "]");
        task.putvarbound(offsett + i, mosek.boundkey.fr, 0., 0.);
      }
      task.putvarbound(offsets, mosek.boundkey.fx, gamma, gamma);
      task.putvarname(offsets, "s");

      double[] e = {1.0, 1.0, 1.0};
      task.putarow(0, xindx, e);

sets up the data for x variables. For instance

      int[] xindx = {offsetx + 0, offsetx + 1, offsetx + 2};
      task.putclist(xindx, mu);

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

        task.putvarname(offsetx + i, "x[" + (i + 1) + "]");
        task.putvarname(offsett + i, "t[" + (i + 1) + "]");

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

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

      //task.writedata("dump.opf");

does that and this will produce a file with the content

Listing 11.3 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.

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.

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 11.1.1, the optimal values of return and risk for several \(\alpha\)s are listed below:

Listing 11.4 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 11.5 demonstrates how to compute the efficient portfolios for several values of \(\alpha\).

Listing 11.5 Code implementing model (7). Click here to download.
package com.mosek.example;

import mosek.*;

public class case_portfolio_2 {
  static final int n = 3;

  public static void main (String[] args) {
    // Since the value infinity is never used, we define
    // 'infinity' symbolic purposes only
    double infinity = 0;
    double gamma = 0.05;
    double[]   mu = {0.1073,  0.0737,  0.0627};
    double[][] GT = {
      {0.1667,  0.0232,  0.0013},
      {0.0000,  0.1033, -0.0022},
      {0.0000,  0.0000,  0.0338}
    };
    double[] x0 = {0.0, 0.0, 0.0};
    double   w = 1.0;
    double[] 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};
    int numalphas = 12;

    int numvar = 2 * n + 1;
    int numcon = n + 1;


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

    try ( Env env  = new mosek.Env ();
          Task task = new mosek.Task (env, 0, 0) ) {

      // Directs the log task stream to the user specified
      // method task_msg_obj.stream
      task.set_Stream(
        mosek.streamtype.log,
        new mosek.Stream()
      { public void stream(String msg) { System.out.print(msg); }});


      //Constraints.
      task.appendcons(numcon);
      for ( int i = 1; i <= n; ++i) {
        w += x0[i - 1];
        task.putconbound(i, mosek.boundkey.fx, 0., 0.);
        task.putconname(i, "GT[" + i + "]");
      }
      task.putconbound(0, mosek.boundkey.fx, w, w);
      task.putconname(0, "budget");

      //Variables.
      task.appendvars(numvar);

      int[] xindx = {offsetx + 0, offsetx + 1, offsetx + 2};
      task.putclist(xindx, mu);

      for ( int i = 0; i < n; ++i) {
        for ( int j = i; j < n; ++j)
          task.putaij(i + 1, offsetx + j, GT[i][j]);

        task.putaij(i + 1, offsett + i, -1.0);

        task.putvarbound(offsetx + i, mosek.boundkey.lo, 0., 0.);

        task.putvarname(offsetx + i, "x[" + (i + 1) + "]");
        task.putvarname(offsett + i, "t[" + (i + 1) + "]");
        task.putvarbound(offsett + i, mosek.boundkey.fr, 0., 0.);
      }
      task.putvarbound(offsets, mosek.boundkey.fr, gamma, gamma);
      task.putvarname(offsets, "s");

      double[] e = {1.0, 1.0, 1.0};
      task.putarow(0, xindx, e);

      //Cones.
      int[] csub = {offsets, offsett + 0, offsett + 1, offsett + 2};
      task.appendcone( mosek.conetype.quad,
                       0.0, /* For future use only, can be set to 0.0 */
                       csub);
      task.putconename(0, "stddev");

      /* A maximization problem */
      task.putobjsense(mosek.objsense.maximize);

      //task.writedata("dump.opf");

      try {
        //Turn all log output off.
        task.putintparam(mosek.iparam.log, 0);

        System.out.printf("%-12s  %-12s  %-12s", "alpha", "exp ret", "std dev");
        for (int k = 0; k < numalphas; ++k) {
          task.putcj(offsets, -alphas[k]);

          task.optimize();

          task.solutionsummary(mosek.streamtype.log);

          double expret = 0.0, stddev = 0.0;
          double[] xx = new double[numvar];

          task.getxx(mosek.soltype.itr, xx);

          for (int j = 0; j < n; ++j)
            expret += mu[j] * xx[j + offsetx];

          System.out.printf("%-12.3e  %-12.3e  %-12.3e\n", alphas[k], expret, xx[offsets]);

        }
        System.out.println("");
      } catch (mosek.Warning mw) {
        System.out.println (" Mosek warning:");
        System.out.println (mw.toString ());
      }
    } catch ( mosek.Exception e) {
      System.out.println ("An error/warning was encountered");
      System.out.println (e.toString ());
      throw e;
    }
  }
}

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

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

11.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 11.2 shows the mapping between the \(\hat{x}\) and the model variables.

Table 11.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 11.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 11.6 demonstrates how to implement the model (14).

Listing 11.6 Code implementing model (14). Click here to download.
package com.mosek.example;

public class case_portfolio_3 {
  static final int n = 3;

  public static void main (String[] args) {
    // Since the value infinity is never used, we define
    // 'infinity' symbolic purposes only
    double infinity = 0;
    double gamma = 0.05;
    double[]   mu = {0.1073,  0.0737,  0.0627};
    double[][] GT = {
      {0.1667,  0.0232,  0.0013},
      {0.0000,  0.1033, -0.0022},
      {0.0000,  0.0000,  0.0338}
    };
    double[] x0 = {0.0, 0.0, 0.0};
    double   w = 1.0;
    double[] m = {0.01, 0.01, 0.01};

    int offsetx = 0;
    int offsets = offsetx + n;
    int offsett = offsets + 1;
    int offsetc = offsett + n;
    int offsetv = offsetc + n;
    int offsetz = offsetv + n;
    int offsetf = offsetz + n;
    int offsetg = offsetf + 3 * n;

    int numvar = offsetg + 3 * n;

    int offset_con_budget = 0;
    int offset_con_gx_t = offset_con_budget + 1;
    int offset_con_abs1 = offset_con_gx_t + n;
    int offset_con_abs2 = offset_con_abs1 + n;
    int offset_con_f =  offset_con_abs2 + n;
    int offset_con_g =  offset_con_f + 3 * n;

    int numcon = 1 + 3 * n + 2 * 3 * n;

    try ( mosek.Env env  = new mosek.Env ();
          mosek.Task task = new mosek.Task (env, 0, 0) ) {

      // Directs the log task stream to the user specified
      // method task_msg_obj.stream
      task.set_Stream(
        mosek.streamtype.log,
        new mosek.Stream()
      { public void stream(String msg) { System.out.print(msg); }});

      //Set up constraint bounds, names and variable coefficients
      task.appendcons(numcon);
      for ( int i = 0; i < n; ++i) {
        w += x0[i];
        task.putconbound(offset_con_gx_t + i, mosek.boundkey.fx, 0., 0.);
        task.putconname(offset_con_gx_t + i, "GT[" + (i + 1) + "]");

        task.putconbound(offset_con_abs1 + i, mosek.boundkey.lo, -x0[i], infinity);
        task.putconname(offset_con_abs1 + i, "zabs1[" + (i + 1) + "]");

        task.putconbound(offset_con_abs2 + i, mosek.boundkey.lo, x0[i], infinity);
        task.putconname(offset_con_abs2 + i, "zabs2[" + (i + 1) + "]");

        for (int j = 0; j < 3; ++j) {
          task.putconbound(offset_con_f + 3 * i + j, mosek.boundkey.fx, 0., 0.);
          task.putconname(offset_con_f + 3 * i + j, "f[" + (i + 1) + "," + (j + 1) + "]");

          task.putconbound(offset_con_g + 3 * i + j, mosek.boundkey.fx, 0., 0.);
          task.putconname(offset_con_g + 3 * i + j, "g[" + (i + 1) + "," + (j + 1) + "]");
        }

        task.putconbound(offset_con_g + 3 * i + 1, mosek.boundkey.fx, -1. / 8, -1. / 8.);
      }
      // e x = w + e x0
      task.putconbound(offset_con_budget, mosek.boundkey.fx, w, w);
      task.putconname(offset_con_budget, "budget");

      //Variables.
      task.appendvars(numvar);

      //the objective function coefficients
      int[] xindx = {offsetx + 0, offsetx + 1, offsetx + 2};
      task.putclist(xindx, mu);

      double[] one_m_one = {1.0, -1.0};
      double[] one_one = {1.0, 1.0};

      //set up variable bounds and names
      for ( int i = 0; i < n; ++i) {
        task.putvarbound(offsetx + i, mosek.boundkey.lo, 0., infinity);
        task.putvarbound(offsett + i, mosek.boundkey.fr, infinity, infinity);
        task.putvarbound(offsetc + i, mosek.boundkey.fr, infinity, infinity);
        task.putvarbound(offsetz + i, mosek.boundkey.fr, infinity, infinity);
        task.putvarbound(offsetv + i, mosek.boundkey.fr, infinity, infinity);
        for (int j = 0; j < 3; ++j) {
          task.putvarbound(offsetf + j + i * 3, mosek.boundkey.fr, infinity, infinity);
          task.putvarbound(offsetg + j + i * 3, mosek.boundkey.fr, infinity, infinity);
        }
        task.putvarname(offsetx + i, "x[" + (i + 1) + "]");
        task.putvarname(offsett + i, "t[" + (i + 1) + "]");
        task.putvarname(offsetc + i, "c[" + (i + 1) + "]");
        task.putvarname(offsetz + i, "z[" + (i + 1) + "]");
        task.putvarname(offsetv + i, "v[" + (i + 1) + "]");
        for (int j = 0; j < 3; ++j) {
          task.putvarname(offsetf + j + i * 3, "f[" + (i + 1) + "," + (j + 1) + "]");
          task.putvarname(offsetg + j + i * 3, "g[" + (i + 1) + "," + (j + 1) + "]");
        }

        for ( int j = i; j < n; ++j)
          task.putaij(offset_con_gx_t + i, j, GT[i][j]);

        task.putaij(offset_con_gx_t + i, offsett + i, -1.0);

        task.putaij(offset_con_budget, offsetx + i, 1.0);
        task.putaij(offset_con_budget, offsetc + i, m[i]);

        // z_j - x_j >= -x0_j
        int[] indx1 = {offsetz + i, offsetx + i};
        task.putarow(offset_con_abs1 + i, indx1, one_m_one);
        // z_j + x_j >= +x0_j
        int[] indx2 = {offsetz + i, offsetx + i};
        task.putarow(offset_con_abs2 + i, indx2, one_one);

        int[] indxf1 = { offsetv + i, offsetf + i * 3};
        task.putarow(offset_con_f + 3 * i, indxf1, one_m_one);
        int[] indxf2 = {offsetc + i, offsetf + i * 3 + 1};
        task.putarow(offset_con_f + 1 + 3 * i, indxf2, one_m_one);
        int[] indxf3 = {offsetz + i, offsetf + i * 3 + 2};
        task.putarow(offset_con_f + 2 + 3 * i, indxf3, one_m_one);

        int[] indxg1 = {offsetz + i, offsetg + i * 3};
        task.putarow(offset_con_g + 3 * i, indxg1, one_m_one);

        task.putaij(offset_con_g + 3 * i + 1, offsetg + i * 3 + 1, -1.);

        int[] indxg3 = {offsetv + i, offsetg + i * 3 + 2};
        task.putarow(offset_con_g + 3 * i + 2, indxg3, one_m_one);
      }
      task.putvarbound(offsets, mosek.boundkey.fx, gamma, gamma);
      task.putvarname(offsets, "s");

      //Cones.
      int conecount = 0;

      int[] csub = {offsets, offsett + 0, offsett + 1, offsett + 2};
      task.appendcone(mosek.conetype.quad, 0.0,  csub);
      task.putconename(conecount, "stddev");
      ++conecount;

      for (int j = 0; j < n; ++j, ++conecount) {
        int[] coneindx = {offsetf + j * 3 , offsetf + j * 3 + 1,  offsetf + j * 3 + 2};
        task.appendcone(mosek.conetype.rquad, 0.0, coneindx);
        task.putconename(conecount, "f[" + (j + 1) + "]");
      }

      for (int j = 0; j < n; ++j, ++conecount) {
        int[] coneindx = {offsetg + j * 3 , offsetg + j * 3 + 1,  offsetg + j * 3 + 2};
        task.appendcone(mosek.conetype.rquad, 0.0, coneindx);
        task.putconename(conecount, "g[" + (j + 1) + "]");
      }
      /* A maximization problem */
      task.putobjsense(mosek.objsense.maximize);

      /* Solve the problem */
      try {
        //Turn all log output off.
        //task.putintparam(mosek.iparam.log,0);

        //task.writedata("dump.opf");

        task.optimize();

        task.solutionsummary(mosek.streamtype.log);

        double expret = 0.0, stddev = 0.0;
        double[] xx = new double[numvar];

        task.getxx(mosek.soltype.itr, xx);

        for (int j = 0; j < n; ++j)
          expret += mu[j] * xx[j + offsetx];

        System.out.printf("Expected return %e for gamma %e\n\n", expret, xx[offsets]);

      } catch (mosek.Warning mw) {
        System.out.println (" Mosek warning:");
        System.out.println (mw.toString ());
      }
    } catch ( mosek.Exception e) {
      System.out.println ("An error/warning was encountered");
      System.out.println (e.toString ());
      throw e;
    }
  }
}

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