8 Appendix

8.1 Conic optimization refresher

In this section we give a short summary about conic optimization. Read more about the topic from a mathematical and modeling point of view in our other publication, the MOSEK Modeling Cookbook [MOSEKApS21].

Conic optimization is a class of convex optimization problems, which contains and generalizes in an unified way many specific and well known convex models. These models include linear optimization (LO), quadratic optimization (QO), quadratically constrained quadratic optimization (QCQO), second-order cone optimization (SOCO), and semidefinite optimization (SDO). The general form of a conic optimization problem is

(8.1)\[\begin{split}\begin{array}{lrcl} \mbox{maximize} & \mathbf{c}^\mathsf{T}\mathbf{x} & &\\ \mbox{subject to} & \mathbf{A}\mathbf{x} + \mathbf{b} & \in & \K.\\ \end{array}\end{split}\]

where \(\K\) is a product of the following basic types of cones:

  • Linear cone:

    \[\R^n,\ \R_+^n,\ \{0\}.\]

    Linear cones model all traditionally LO problems.

  • Quadratic cone and rotated quadratic cone:

    The quadratic cone is the set

    \[\Q^n = \left\{x\in\real^n\middle\vert x_1\geq\sqrt{x_2^2+\cdots+x_n^2}\right\}.\]

    The rotated quadratic cone is the set

    \[\Q_\mathrm{r}^n = \left\{x\in\real^n\middle\vert 2x_1x_2\geq x_3^2+\cdots+x_n^2,\ x_1,x_2\geq 0\right\}.\]

    Modeling with these two cones cover the class of SOCO problems, which include all traditionally QO and QCQO problems as well. See in Sec. 8.1.2 (Traditional quadratic models) for more details.

  • Primal power cone:

    \[\POW_n^{\alpha,1-\alpha} = \left\{x\in\real^n\middle\vert x_1^\alpha x_2^{1-\alpha}\geq \sqrt{x_3^2+\cdots+x_n^2},\ x_1,x_2\geq 0\right\},\]

    or its dual cone.

  • Primal exponential cone:

    \[\EXP=\left\{x\in\real^3 \middle\vert x_1\geq x_2\exp\left(\frac{x_3}{x_2}\right),\ x_1,x_2\geq 0\right\},\]

    or its dual cone.

  • Semidefinite cone:

    \[\PSD^n=\left\{X\in\real^{n\times n}\middle\vert X\ \textrm{is symmetric positive semidefinite}\right\}.\]

    Semidefinite cones model all traditionally SDO problems.

Each of these cones allow formulating different types of convex constraints.

8.1.1 Selection of conic constraints

In the following we will list the constraints appearing in financial context in this book and show how can we convert them to conic form.

8.1.1.1 Maximum function

We can model the maximum constraint \(\mathrm{max}(x_1, x_2, \dots, x_n) \leq c\) using linear constraints by introducing an auxiliary variable \(t\):

(8.2)\[\begin{split}\begin{array}{lrcl} t & \leq & c, \\ t & \geq & x_1,\\ & \vdots &\\ t & \geq & x_n.\\ \end{array}\end{split}\]

For instance, we could write a series of constraints \(\mathrm{max}(x_i, 0) \leq c_i,\ i = 1,\dots,n\) as

(8.3)\[\mathbf{t} \leq \mathbf{c},\quad \mathbf{t} \geq \mathbf{x},\quad \mathbf{t} \geq \mathbf{0},\]

where \(\mathbf{t}\) is an \(n\)-dimensional vector.

A special case of modeling the maximum function is to model the positive part \(x^+\) and negative part \(x^-\) of a variable \(x\). We define these as \(x^+ = \mathrm{max}(x, 0)\) and \(x^- = \mathrm{max}(-x, 0)\). We can model them by the following set of constraints:

(8.4)\[\begin{split}\begin{array}{lrcl} x & = & x^+ - x^-,\\ x^+,\ x^- & \geq & 0.\\ \end{array}\end{split}\]

Note that here we need additional constraints penalizing the magnitude of these two variables, so that they do not both become positive in any optimal solution. One possible way to ensure this can be used in mixed integer problems, see Sec. 8.2.1.4 (Positive and negative part).

8.1.1.2 Absolute value

We can model the absolute value constraint \(|x| \leq c\) using the maximum function by observing that \(|x| = \mathrm{max}(x, -x)\) (see Sec. 8.1.1.1 (Maximum function)):

(8.5)\[-c \leq x \leq c.\]

Another possibility is to model it using the quadratic cone:

(8.6)\[(c, x) \in \Q^{2}.\]

8.1.1.3 Manhattan norm

Let \(\mathbf{x} \in \R^n\) and observe that \(\|\mathbf{x}\|_1 = |x_1| +\cdots+ |x_n|\). Then we can model the Manhattan norm or 1-norm constraint \(\|\mathbf{x}\|_1 \leq c\) by modeling the absolute value for each coordinate:

(8.7)\[-\mathbf{z} \leq \mathbf{x} \leq \mathbf{z},\quad \sum_{i=1}^n z_i = c,\]

where \(\mathbf{z}\) is an auxiliary variable.

8.1.1.4 Euclidean norm

Let \(\mathbf{x} \in \R^n\) and observe that \(\|\mathbf{x}\|_2 = \sqrt{x_1^2+\cdots+x_n^2}\). Then we can model the Euclidean norm or 2-norm constraint \(\|\mathbf{x}\|_2 \leq c\) using the quadratic cone:

(8.8)\[(c, \mathbf{x}) \in \Q^{n+1}.\]

8.1.1.5 Squared Euclidean norm

Let \(\mathbf{x} \in \R^n\) and observe that \(\|\mathbf{x}\|_2^2 = \mathbf{x}^\mathsf{T}\mathbf{x} = x_1^2+\cdots+x_n^2\). Then we can model the squared Euclidean norm or sum-of-squares constraint \(\|\mathbf{x}\|_2^2 \leq c\) using the rotated quadratic cone:

(8.9)\[\left(c, \textstyle{\frac{1}{2}}, \mathbf{x}\right) \in \Q_\mathrm{r}^{n+2}.\]

8.1.1.6 Quadratic form

Let \(\mathbf{x} \in \R^n\) and let \(\mathbf{Q} \in \PSD^n\), i. e., a symmetric positive semidefinite matrix. Then we can model the quadratic form constraint \(\frac{1}{2}\mathbf{x}^\mathsf{T}\mathbf{Q}\mathbf{x} \leq c\) either using the quadratic cone or the rotated quadratic cone. To see this, observe that there exists a matrix \(\mathbf{G} \in \R^{n\times k}\) such that \(\mathbf{Q} = \mathbf{G}\mathbf{G}^\mathsf{T}\). Of course this decomposition can be done in many ways, so the matrix \(\mathbf{G}\) is not unique. The most interesting cases are when \(k \ll n\) or \(\mathbf{G}\) is very sparse, because these make the optimization problem much easier to solve numerically (see in Sec. 8.1.2 (Traditional quadratic models)).

Most common ways to compute the matrix \(\mathbf{G}\) are the following:

  • Cholesky decomposition: \(\mathbf{Q} = \mathbf{C}\mathbf{C}^\mathsf{T}\), where \(\mathbf{C}\) is a lower triangular matrix with nonnegative entries in the diagonal. From this decomposition we have \(\mathbf{G} = \mathbf{C} \in \R^{n\times n}\).

  • Eigenvalue decomposition: \(\mathbf{Q} = \mathbf{V}\mathbf{D}\mathbf{V}^\mathsf{T}\), where the diagonal matrix \(\mathbf{D}\) contains the (nonnegative) eigenvalues of \(\mathbf{Q}\) and the unitary matrix \(\mathbf{V}\) contains the corresponding eigenvectors in its columns. From this decomposition we have \(\mathbf{G} = \mathbf{V}\mathbf{D}^{1/2} \in \R^{n\times n}\).

  • Matrix square root: \(\mathbf{Q} = \mathbf{Q}^{1/2}\mathbf{Q}^{1/2}\), where \(\mathbf{Q}^{1/2}\) is the symmetric positive semidefinite “square root” matrix of \(\mathbf{Q}\). From this decomposition we have \(\mathbf{G} = \mathbf{Q}^{1/2} \in \R^{n\times n}\).

  • Factor model: If \(\mathbf{Q}\) is a covariance matrix of some data, then we can approximate the data series with the combination of \(k \ll n\) common factors. Then we have the decomposition \(\mathbf{Q} = \EBeta\mathbf{Q}_F\EBeta^\mathsf{T} + \mathbf{D}\), where \(\mathbf{Q}_F \in \R^{k\times k}\) is the covariance of the factors, \(\EBeta \in \R^{n\times k}\) is the exposure of the data series to each factor, and \(\mathbf{D}\) is diagonal. From this, by computing the Cholesky decomposition \(\mathbf{Q}_F = \mathbf{F}\mathbf{F}^\mathsf{T}\) we have \(\mathbf{G} = [\EBeta\mathbf{F}, \mathbf{D}^{1/2}] \in \R^{n\times (n+k)}\). The advantage of factor models is that \(\mathbf{G}\) is very sparse and the factors have a direct financial interpretation (see Sec. 5 (Factor models) for details).

After obtaining the matrix \(\mathbf{G}\), we can write the quadratic form constraint as a sum-of-squares \(\frac{1}{2}\mathbf{x}^\mathsf{T}\mathbf{G}\mathbf{G}^\mathsf{T}\mathbf{x} \leq c\), which is a squared Euclidean norm constraint \(\frac{1}{2}\|\mathbf{G}^\mathsf{T}\mathbf{x}\|_2^2 \leq c\).

We can choose to model this using the rotated quadratic cone as

(8.10)\[\left(c, 1, \mathbf{G}^\mathsf{T}\mathbf{x}\right) \in \Q_\mathrm{r}^{k+2},\]

or we can choose to model its square root using the quadratic cone as

(8.11)\[\left(\sqrt{c}, \mathbf{G}^\mathsf{T}\mathbf{x}\right) \in \Q^{k+1}.\]

Whether to use the quadratic cone or the rotated quadratic cone for modeling can be decided based on which is more natural. Typically the quadratic cone is used to model 2-norm constraints, while the rotated quadratic cone is more natural for modeling of quadratic functions. There can be exceptions, however, see for example in Sec. 2.3 (Conic formulation).

8.1.1.7 Power

Let \(x \in \R\) and \(\alpha > 1\). Then we can model the power constraint \(c \geq |x|^\alpha\) or equivalently \(c^{1/\alpha} \geq |x|\) using the power cone:

(8.12)\[(c, 1, x) \in \POW_3^{1/\alpha,(\alpha-1)/\alpha}.\]

8.1.2 Traditional quadratic models

Optimization problems involving quadratic functions often appear in practice. In this section we show how to convert the traditionally QO or QCQO problems into conic optimization form, because most of the time the solution of these conic models is computationally more efficient.

8.1.2.1 Quadratic optimization

The standard form of a quadratic optimization (QO) problem is the following:

(8.13)\[\begin{split}\begin{array}{lrcl} \mbox{minimize} & \frac{1}{2}\mathbf{x}^\mathsf{T}\mathbf{Q}\mathbf{x} + \mathbf{c}^\mathsf{T}\mathbf{x} & &\\ \mbox{subject to} & \mathbf{A}\mathbf{x} & = & \mathbf{b},\\ & \mathbf{x} & \geq & \mathbf{0}.\\ \end{array}\end{split}\]

The matrix \(\mathbf{Q} \in \R^{n\times n}\) must be symmetric positive semidefinite, otherwise the objective function would not be convex.

Assuming the factorization \(\mathbf{Q} = \mathbf{G}\mathbf{G}^\mathsf{T}\) with \(\mathbf{G}\in \R^{n\times k}\), we can reformulate the problem (8.13) as a conic problem by applying the method described in Sec. 8.1.1.6 (Quadratic form):

(8.14)\[\begin{split}\begin{array}{lrcl} \mbox{minimize} & t + \mathbf{c}^\mathsf{T}\mathbf{x} & &\\ \mbox{subject to} & \mathbf{A}\mathbf{x} & = & \mathbf{b},\\ & \mathbf{x} & \geq & \mathbf{0},\\ & (t, 1, \mathbf{G}^\mathsf{T}\mathbf{x}) & \in & \Q_\mathrm{r}^{k+2}. \end{array}\end{split}\]

8.1.2.2 Quadratically constrained quadratic optimization

Consider the quadratically constrained quadratic optimization (QCQO) problem of the form

(8.15)\[\begin{split}\begin{array}{lrcl} \mbox{minimize} & \frac{1}{2}\mathbf{x}^\mathsf{T}\mathbf{Q}_0\mathbf{x} + \mathbf{c}_0^\mathsf{T}\mathbf{x} + a_0 & &\\ \mbox{subject to} & \frac{1}{2}\mathbf{x}^\mathsf{T}\mathbf{Q}_i\mathbf{x} + \mathbf{c}_i^\mathsf{T}\mathbf{x} + a_i & \leq & 0,\quad i=1,\dots,m.\\ \end{array}\end{split}\]

The matrices \(\mathbf{Q}_i \in \R^{n\times n}\), \(i=0,\dots,m\) must all be symmetric positive semidefinite, otherwise the optimization problem would not be convex.

Assuming the factorization \(\mathbf{Q}_i = \mathbf{G}_i\mathbf{G}_i^\mathsf{T}\) with \(\mathbf{G}_i\in \R^{n\times k_i}\), we can reformulate the problem (8.15) as a conic problem by applying the method described in Sec. 8.1.1.6 (Quadratic form) for both the objective and the constraints:

(8.16)\[\begin{split}\begin{array}{lrclr} \mbox{minimize} & t_0 + \mathbf{c}_0^\mathsf{T}\mathbf{x} + a_0 & & &\\ \mbox{subject to} & t_i + \mathbf{c}_i^\mathsf{T}\mathbf{x} + a_i & \leq & 0, & i=1,\dots,m,\\ & (t_i, 1, \mathbf{G}_i^\mathsf{T}\mathbf{x}) & \in & \Q_\mathrm{r}^{k_i+2}, & i=0,\dots,m.\\ \end{array}\end{split}\]

8.1.2.3 Practical benefits of conic models

The key step in the model conversion is to transform quadratic terms \(\mathbf{x}^\mathsf{T}\mathbf{Q}\mathbf{x}\) using the factorization \(\mathbf{Q} = \mathbf{G}\mathbf{G}^\mathsf{T}\), where \(\mathbf{G}\in \R^{n\times k}\). Assuming \(k \ll n\), it results in the following benefits:

  • The storage requirement \(nk\) of \(\mathbf{G}\) can be much lower than the storage requirement \(n^2/2\) of \(\mathbf{Q}\).

  • The amount of work to evaluate \(\mathbf{x}^\mathsf{T}\mathbf{Q}\mathbf{x}\) is proportional to \(n^2\) whereas the work to evaluate \(\|\mathbf{G}^\mathsf{T}\mathbf{x}\|_2^2\) is proportional to \(nk\) only.

  • No need to numerically validate positive semidefiniteness of the matrix \(\mathbf{Q}\). This is otherwise difficult owing to the presence of rounding errors that can make \(\mathbf{Q}\) indefinite.

  • Duality theory is much simpler for conic quadratic optimization.

In summary, the conic equivalents are not only as easy to solve as the original QP or QCQP problems, but in most cases also need less space and solution time.

8.2 Mixed-integer models

Mixed integer optimization (MIO) is an extension of convex optimization, which introduces variables taking values only in the set of integers or in some subset of integers. This allows for solving a much wider range of optimization problems in addition to the ones convex optimization already covers. However, these problems are no longer convex. Handling of integer variables make the optimization NP-hard and needs specific algorithms, thus the solution of MIO problems is much slower. In many practical cases we cannot expect to find the optimal solution in reasonable time, and only near-optimal solutions will be available.

A general mixed integer optimization problem has the following form:

(8.17)\[\begin{split}\begin{array}{lrcl} \mbox{maximize} & \mathbf{c}^\mathsf{T}\mathbf{x} & &\\ \mbox{subject to} & \mathbf{A}\mathbf{x} + \mathbf{b} & \in & \K,\\ & x_i & \in & \mathbb{Z},\quad i \in \mathcal{I}. \end{array}\end{split}\]

where \(\K\) is a cone and \(\mathcal{I} \subseteq \{1,\dots,n\}\) contains the indices of integer variables. We can model any finite range for the integer variable \(x_i\) by simply adding the extra constraint \(a_i \leq x_i \leq b_i\).

8.2.1 Selection of integer constraints

In the following we will list the integer constraints appearing in financial context in this book and show how to model them.

8.2.1.1 Switch

In some practical cases we might wish to impose conditions on parts of our optimization model. For example, allowing nonzero value for a variable only in presence of a condition. We can model such situations using binary variables (or indicator variables), which can only take \(0\) or \(1\) values. The following set of constraints only allow \(x\) to be positive when the indicator variable \(y\) is equal to \(1\):

(8.18)\[\begin{split}\begin{array}{rcl} x & \leq & My, \\ y & \in & \{0, 1\}.\\ \end{array}\end{split}\]

The number \(M\) here is not related to the optimization problem, but it is necessary to form such a switchable constraint. If \(y=0\), then the upper limit of \(x\) is \(0\). If \(y=1\), then the upper limit of \(x\) is \(M\). This modeling technique is called big-M. The choice of \(M\) can affect the solution performance, but a nice feature of it is that the problem cannot accidentally become infeasible.

If we have a vector variable \(\mathbf{x}\) then the switch will look like:

(8.19)\[\begin{split}\begin{array}{rcl} \mathbf{x} & \leq & \mathbf{M}\circ\mathbf{y},\\ \mathbf{y} & \in & \{0, 1\}^n,\\ \end{array}\end{split}\]

where \(\circ\) denotes the elementwise product. We accounted for a possibly different big-M value for each coordinate.

8.2.1.2 Semi-continuous variable

A slight extension of Sec. 8.2.1.1 (Switch) is to model semi-continuous variables, i. e. when \(x \in \{0\} \cup [a, b]\), where \(0 < a \leq b\). We can model this by

(8.20)\[ay \leq x \leq by,\quad y \in \{0, 1\}.\]

8.2.1.3 Cardinality

We might need to limit the number of nonzeros in a vector \(\mathbf{x}\) to some number \(K<n\). We can do this with the help of an indicator vector \(\mathbf{y}\) of length \(n\), which indicates \(|\mathbf{x}| \neq \mathbf{0}\) (see Sec. 8.2.1.1 (Switch)). First we add a big-M bound \(\mathbf{M}\circ\mathbf{y}\) to the absolute value, and model it based on Sec. 8.1.1.2 (Absolute value). Then we limit the cardinality of \(\mathbf{x}\) by limiting the cardinality of \(\mathbf{y}\):

(8.21)\[\begin{split}\begin{array}{rcl} \mathbf{x} & \leq & \mathbf{M}\circ\mathbf{y},\\ \mathbf{x} & \geq & \mathbf{M}\circ\mathbf{y},\\ \mathbf{1}^\mathsf{T}\mathbf{y} & \leq & K,\\ \mathbf{y} & \in & \{0, 1\}^n.\\ \end{array}\end{split}\]

8.2.1.4 Positive and negative part

We introduced the positive part \(x^+\) and negative part \(x^-\) of a variable \(x\) in Sec. 8.1.1.1 (Maximum function). We noted that we need a way to ensure only \(x^+\) or \(x^-\) will be positive in the optimal solution and not both. One such way is to use binary variables:

(8.22)\[\begin{split}\begin{array}{lrcl} x & = & x^+ - x^-,\\ x^+ & \leq & My,\\ x^- & \leq & M(1-y),\\ x^+,\ x^- & \geq & 0,\\ y & \in & \{0, 1\}.\\ \end{array}\end{split}\]

Here the binary variable \(y\) allows the positive (negative) part to become nonzero exactly when the negative (positive) part is zero.

Sometimes we need to handle separately the case when both the positive and the negative part is fixed to be zero. Then we need to introduce two binary variables \(y^+\) and \(y^-\), and include the constraint \(y^+ + y^- \leq 1\) to prevent both variables from being \(1\):

(8.23)\[\begin{split}\begin{array}{lrcl} x & = & x^+ - x^-,\\ x^+ & \leq & My^+,\\ x^- & \leq & My^-,\\ x^+,\ x^- & \geq & 0,\\ y^+ + y^- & \leq & 1,\\ y^+,\ y^- & \in & \{0, 1\}.\\ \end{array}\end{split}\]

8.3 Quadratic cone and risk-free security

In Sec. 6 (Transaction costs) we consider portfolio optimization problems with risk-free security. In this case there is a difference in the computation of the efficient frontier, if we model using the quadratic cone instead of the rotated quadratic cone. Namely, the optimization problem will not have solutions that are a mix of the risk-free security and risky securities. Instead, it will only have either the 100% risk-free solution, or a portfolio of only risky securities.

Consider the following simple model:

(8.24)\[\begin{split}\begin{array}{lrcll} \mbox{maximize} & \EMean^\mathsf{T}\mathbf{x} + r^\mathrm{f}x^\mathrm{f} - \tilde{\delta} \sqrt{\mathbf{x}^\mathsf{T}\ECov\mathbf{x}} & & \\ \mbox{subject to} & \mathbf{1}^\mathsf{T}\mathbf{x} + x^\mathrm{f} & = & 1,\\ & \mathbf{x}, x^\mathrm{f} & \geq & 0, \end{array}\end{split}\]

where \(r^\mathrm{f}\) is the risk-free rate and \(x^\mathrm{f}\) is the allocation to the risk-free asset.

We can transform this problem using \(x^\mathrm{f} = 1 - \mathbf{1}^\mathsf{T}\mathbf{x}\), such that we are only left with the variable \(\mathbf{x}\):

(8.25)\[\begin{split}\begin{array}{lrcll} \mbox{maximize} & (\EMean-r^\mathrm{f}\mathbf{1})^\mathsf{T}\mathbf{x} + r^\mathrm{f} - \tilde{\delta} \sqrt{\mathbf{x}^\mathsf{T}\ECov\mathbf{x}} & & \\ \mbox{subject to} & \mathbf{1}^\mathsf{T}\mathbf{x} & \leq & 1,\\ & \mathbf{x} & \geq & 0. \end{array}\end{split}\]

The feasible region in this form is a probability simplex. The solution \(\mathbf{x} = \mathbf{0}\) means that everything is allocated to the risk-free security, and the hyperplane \(\mathbf{1}^\mathsf{T}\mathbf{x} = 1\) has all the feasible solutions purely involving risky assets.

Let us denote the objective function value by \(\mathrm{obj}(\mathbf{x})\). Then the directional derivative of the objective along the direction \(\mathbf{u} > \mathbf{0}\), \(||\mathbf{u}|| = 1\) will be

\[\partial_\mathbf{u}\mathrm{obj}(\mathbf{x}) = (\EMean-r^\mathrm{f}\mathbf{1})^\mathsf{T}\mathbf{u} - \tilde{\delta}\frac{\mathbf{x}^\mathsf{T}\ECov\mathbf{u}}{\sqrt{\mathbf{x}^\mathsf{T}\ECov\mathbf{x}}}.\]

This does not depend on the norm of \(\mathbf{x}\), meaning that \(\partial_\mathbf{u}\mathrm{obj}(c\mathbf{u})\) will be constant in \(c\):

\[\partial_\mathbf{u}\mathrm{obj}(c\mathbf{u}) = (\EMean-r^\mathrm{f}\mathbf{1})^\mathsf{T}\mathbf{u} - \tilde{\delta} \sqrt{\mathbf{u}^\mathsf{T}\ECov\mathbf{u}}.\]

Thus the objective is linear along the 1-dimensional slice between \(\mathbf{x} = \mathbf{0}\) and any \(\mathbf{x} > \mathbf{0}\), meaning that the optimal solution is either \(\mathbf{x} = \mathbf{0}\) or some \(\mathbf{x} > \mathbf{0}\) satisfying \(\mathbf{1}^\mathsf{T}\mathbf{x} = 1\).

Furthermore, along every 1-dimensional slice represented by a direction \(\mathbf{u}\), there is a \(\tilde{\delta}\) threshold

\[\tilde{\delta}_\mathbf{u} = \frac{(\EMean - r^\mathrm{f}\mathbf{1})^\mathsf{T}\mathbf{u}}{\sqrt{\mathbf{u}^\mathsf{T}\mathbf{Su}}}.\]

This is the Sharpe ratio of all portfolios of the form \(c\mathbf{u}\). If \(\tilde{\delta} > \tilde{\delta}_\mathbf{u}\) then \(\partial_\mathbf{u}\mathrm{obj}(\mathbf{u})\) turns negative.

In general, the larger we set \(\tilde{\delta}\), the fewer directions \(\mathbf{u}\) will exist for which \(\tilde{\delta} < \tilde{\delta}_\mathbf{u}\), i. e., \(\partial_\mathbf{u}\mathrm{obj}(\mathbf{u}) > 0\), and the optimal solution is a combination of risky assets. The largest \(\tilde{\delta}\) for which we still have such a direction is \(\tilde{\delta}^* = \mathrm{max}_\mathbf{u}\ \tilde{\delta}_\mathbf{u}\), the maximal Sharpe ratio. The corresponding optimal portfolio is the one having the smallest risk while consisting purely of risky assets. For \(\tilde{\delta} > \tilde{\delta}^*\) we have \(\partial_\mathbf{u}\mathrm{obj}(\mathbf{u}) < 0\) in all directions, meaning that the optimal solution will always be \(\mathbf{x} = \mathbf{0}\), the 100% risk-free portfolio.

8.4 Monte Carlo Optimization Selection (MCOS)

We can use the following procedure to compare the accuracy of different portfolio optimization methods.

  1. Inputs of the procedure are the estimated mean return \(\EMean\) and estimated covariance \(\ECov\). Treat this pair of inputs as true values.

  2. Compute the optimal asset allocation \(\mathbf{x}\) for the original input pair \((\EMean, \ECov)\).

  3. Repeat \(K_\mathrm{sim}\) times:

    1. Do parametric bootstrap resampling, i. e., draw a new \(N\times T\) return data sample \(\mathbf{R}_k\) and derive a simulated pair \((\EMean_k, \ECov_k)\).

    2. De-noise the covariance matrix \(\ECov_k\) using the method in Sec. 4.2.2 (Covariance de-noising).

    3. Compute the optimal portfolio \(\mathbf{x}_k\) using the optimization method(s) that we wish to analyze.

  4. To estimate the error for each security in the optimal portfolio, compute the standard deviation of the difference vectors \(\mathbf{x}_k-\mathbf{x}\). Then by taking the mean we get a single number estimating the error of the optimization method.

Regarding the final step, we can also measure the average decay in performance by any of the following:

  • the mean difference in expected outcomes: \(\frac{1}{K_\mathrm{sim}}\sum_k (\mathbf{x}_k-\mathbf{x})^\mathsf{T}\EMean\)

  • the mean difference in variance: \(\frac{1}{K_\mathrm{sim}}\sum_k (\mathbf{x}_k-\mathbf{x})^\mathsf{T}\ECov(\mathbf{x}_k-\mathbf{w})\)

  • the mean difference in Sharpe ratio or other metric computed from the above statistics.