# 7.1 Portfolio Optimization¶

This case studies is devoted to the Portfolio Optimization Problem.

The complete code can be downloaded here.

## 7.1.1 The Basic 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)

$\mathbf{E}(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}{llll} \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$$ denotes 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 the parameter $$\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.2. For a given $$G$$ we have that

$\begin{split}\begin{array}{lll} 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 Q^{n+1}.$

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

(3)$\begin{split}\begin{array}{llll} \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 formulated and solved with Fusion. 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}$

We will make use of a simple helper functions to compute the inner product of two vectors and their sum.

static double sum(std::shared_ptr<ndarray<double,1>> x)
{
double r = 0.0;
for (auto v : *x) r += v;
return r;
}

static double dot(std::shared_ptr<ndarray<double,1>> x,
std::shared_ptr<ndarray<double,1>> y)
{
double r = 0.0;
for (int i = 0; i < x->size(); ++i) r += (*x)[i] * (*y)[i];
return r;
}

static double dot(std::shared_ptr<ndarray<double,1>> x,
std::vector<double> & y)
{
double r = 0.0;
for (int i = 0; i < x->size(); ++i) r += (*x)[i] * y[i];
return r;
}


Listing 7.1 demonstrates how the basic Markowitz model (3) is implemented using Fusion.

Listing 7.1 Code implementing problem (3).
/*
Purpose:
Computes the optimal portfolio for a given risk

Input:
n: Number of assets
mu: An n dimmensional vector of expected returns
GT: A matrix with n columns so (GT')*GT  = covariance matrix
x0: Initial holdings
w: Initial cash holding
gamma: Maximum risk (=std. dev) accepted

Output:
Optimal expected return and the optimal portfolio
*/
double BasicMarkowitz
( int                                n,
std::shared_ptr<ndarray<double,1>> mu,
std::shared_ptr<ndarray<double,2>> GT,
std::shared_ptr<ndarray<double,1>> x0,
double                             w,
double                             gamma)
{
Model::t M = new Model("Basic Markowitz"); auto _M = finally([&]() { M->dispose(); });
// Redirect log output from the solver to stdout for debugging.
M->setLogHandler([](const std::string & msg) { std::cout << msg << std::flush; } );

// Defines the variables (holdings). Shortselling is not allowed.
Variable::t x = M->variable("x", n, Domain::greaterThan(0.0));

//  Maximize expected return
M->objective("obj", ObjectiveSense::Maximize, Expr::dot(mu,x));

// The amount invested  must be identical to intial wealth
M->constraint("budget", Expr::sum(x), Domain::equalsTo(w+sum(x0)));

// Imposes a bound on the risk
M->constraint("risk", Expr::vstack(gamma,Expr::mul(GT,x)), Domain::inQCone());

// Solves the model.
M->solve();

return dot(mu,x->level());
}


The source code should be self-explanatory except perhaps for

M->constraint("risk", Expr::vstack(gamma,Expr::mul(GT,x)), Domain::inQCone());


where the linear expression

$\left[ \begin{array}{c} \gamma G^T x \end{array} \right]$

is created using the Expr.vstack operator. Finally, the linear expression must lie in a quadratic cone implying

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

## 7.1.2 The Efficient Frontier¶

The portfolio computed by the Markowitz model is efficient in the sense that there no other 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 for 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

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

computes an efficient portfolio. Note that the objective maximize 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 (4) should be solved for all values $$\alpha \geq 0$$ but in practice impossible. Using the example data from Section 7.1.1, the optimal values of return and risk for several $$\alpha$$s are listed below:

Efficient frontier

alpha         return        risk
0.0000        1.0730e-01    7.2700e-01
0.0100        1.0730e-01    1.6667e-01
0.1000        1.0730e-01    1.6667e-01
0.2500        1.0321e-01    1.4974e-01
0.3000        8.0529e-02    6.8144e-02
0.3500        7.4290e-02    4.8585e-02
0.4000        7.1958e-02    4.2309e-02
0.4500        7.0638e-02    3.9185e-02
0.5000        6.9759e-02    3.7327e-02
0.7500        6.7672e-02    3.3816e-02
1.0000        6.6805e-02    3.2802e-02
1.5000        6.6001e-02    3.2130e-02
2.0000        6.5619e-02    3.1907e-02
3.0000        6.5236e-02    3.1747e-02
10.0000       6.4712e-02    3.1633e-02


Example code

Listing 7.2 demonstrates how to compute the efficient portfolios for several values of $$\alpha$$ in Fusion.

Listing 7.2 Code for the computation of the efficient frontier based on problem (4).
/*
Purpose:
Computes several portfolios on the optimal portfolios by

for alpha in alphas:
maximize   expected return - alpha * standard deviation
subject to the constraints

Input:
n: Number of assets
mu: An n dimmensional vector of expected returns
GT: A matrix with n columns so (GT')*GT  = covariance matrix
x0: Initial holdings
w: Initial cash holding
alphas: List of the alphas

Output:
The efficient frontier as list of tuples (alpha,expected return,risk)
*/
void EfficientFrontier
( int n,
std::shared_ptr<ndarray<double,1>> mu,
std::shared_ptr<ndarray<double,2>> GT,
std::shared_ptr<ndarray<double,1>> x0,
double      w,
std::vector<double> & alphas,

std::vector<double> & frontier_mux,
std::vector<double> & frontier_s)
{

Model::t M = new Model("Efficient frontier");  auto M_ = finally([&]() { M->dispose(); });
M->setLogHandler([](const std::string & msg) { std::cout << msg << std::flush; } );

// Defines the variables (holdings). Shortselling is not allowed.
Variable::t x = M->variable("x", n, Domain::greaterThan(0.0)); // Portfolio variables
Variable::t s = M->variable("s", 1, Domain::unbounded()); // Risk variable

M->constraint("budget", Expr::sum(x), Domain::equalsTo(w+sum(x0)));

// Computes the risk
M->constraint("risk", Expr::vstack(s,Expr::mul(GT,x)),Domain::inQCone());

Expression::t mudotx = Expr::dot(mu,x);

for (double alpha : alphas)
{
//  Define objective as a weighted combination of return and risk
M->objective("obj", ObjectiveSense::Maximize, Expr::sub(mudotx,Expr::mul(alpha,s)));

M->solve();

frontier_mux.push_back(dot(mu,x->level()));
frontier_s.push_back((*s->level())[0]);
}
}


Note the efficient frontier could also have been computed using the code in Section 7.1.1 by varying $$\gamma$$. However, when the constraints of a Fusion model is changed the model has to be rebuild whereas a rebuild is not needed if only the objective is modified.

## 7.1.3 Improving the Computational Efficiency¶

In practice it is often important to solve the portfolio problem in a short amount of time. Therefore, in this section it is discussed 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 a positive definite diagonal matrix. Moreover, $$V$$ is a matrix with $$n$$ rows and $$p$$ columns. Such a model for the covariance matrix is called a factor model index{factor model} and usually $$p$$ is much smaller than $$n$$. In practice $$p$$ tends be a small number independent of $$n$$ say less than 100.

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

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

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

(5)$\begin{split}\begin{array}{llll} \mbox{maximize} & \mu^T x & & \\ \mbox{subject to} & e^T x + \sum_{j=1}^{n} T_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

$T_j( x_j - x_j^0)$

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

## 7.1.5 Market Impact Costs¶

If the initial wealth is fairly small and no short selling is allowed, then the holdings will be small and then the amount traded of each asset must also be small. Therefore, 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 the market impact cost for asset $$j$$ can be modelled by

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

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

$T_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

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

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

$\begin{split}\begin{array}{ll} z_j & = |x_j - x_j^0|, \\ (s_j,t_j,z_j),(z_j,1/8,s_j) & \in Q_r^3,\\ \sum_{j=1}^n T(x_j - x_j^0) & = \sum_{j=1}^n t_j. \end{array}\end{split}$

Unfortunately this set of constraints is nonconvex due to the constraint

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

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

(7)$z_j \geq |x_j - x_j^0|.$

which is equivalent to

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

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

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

cannot hold for an optimal solution.

Now given that the optimal solution has the property that (9) holds then the market impact costs 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. It is assumed this is not the case and hence the models (6) and (7) are equivalent.

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

The revised budget constraint

$e^T x = w + e^T x^0 - m^T t$

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

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

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

Example code

Listing 7.3 demonstrates how to compute an optimal portfolio when market impact cost are included using Fusion.

Listing 7.3 Implementation of model (10).
/*
Description:
Extends the basic Markowitz model with a market cost term.

Input:
n: Number of assets
mu: An n dimmensional vector of expected returns
GT: A matrix with n columns so (GT')*GT  = covariance matrix'
x0: Initial holdings
w: Initial cash holding
gamma: Maximum risk (=std. dev) accepted
m: It is assumed that  market impact cost for the j'th asset is
m_j|x_j-x0_j|^3/2

Output:
Optimal expected return and the optimal portfolio

*/
void MarkowitzWithMarketImpact
( int n,
std::shared_ptr<ndarray<double,1>> mu,
std::shared_ptr<ndarray<double,2>> GT,
std::shared_ptr<ndarray<double,1>> x0,
double      w,
double      gamma,
std::shared_ptr<ndarray<double,1>> m,
std::vector<double> & xsol,
std::vector<double> & tsol)
{
Model::t M = new Model("Markowitz portfolio with market impact");  auto M_ = finally([&]() { M->dispose(); });
M->setLogHandler([](const std::string & msg) { std::cout << msg << std::flush; } );

// Defines the variables. No shortselling is allowed.
Variable::t x = M->variable("x", n, Domain::greaterThan(0.0));

Variable::t t = M->variable("t", n, Domain::unbounded());
Variable::t z = M->variable("z", n, Domain::unbounded());
Variable::t v = M->variable("v", n, Domain::unbounded());

//  Maximize expected return
M->objective("obj", ObjectiveSense::Maximize, Expr::dot(mu,x));

// Invested amount + slippage cost = initial wealth

// Imposes a bound on the risk
M->constraint("risk", Expr::vstack( gamma,Expr::mul(GT,x)),
Domain::inQCone());

// z >= |x-x0|
M->constraint("sell", Expr::sub(z,Expr::sub(x0,x)),Domain::greaterThan(0.0));

// t >= z^1.5, z >= 0.0. Needs two rotated quadratic cones to model this term
M->constraint("ta", Expr::hstack(v,t,z),Domain::inRotatedQCone());
M->constraint("tb", Expr::hstack(z,Expr::constTerm(n,1.0/8.0),v),
Domain::inRotatedQCone());

M->solve();

xsol.resize(n);
tsol.resize(n);
auto xlvl = x->level();
auto tlvl = t->level();

std::copy(xlvl->flat_begin(), xlvl->flat_end(), xsol.begin());
std::copy(tlvl->flat_begin(), tlvl->flat_end(), tsol.begin());
}


The major new feature compared to the previous examples are

M->constraint("ta", Expr::hstack(v,t,z),Domain::inRotatedQCone());


and

M->constraint("tb", Expr::hstack(z,Expr::constTerm(n,1.0/8.0),v),
Domain::inRotatedQCone());


In the first line the variables $$v$$, $$t$$ and $$z$$ are stacked horizontally which corresponds to creating a list of linear expressions where the $$j$$‘th element has the form

$\begin{split}\left[ \begin{array}{c} v_j\\ t_j\\ z_j \end{array} \right]\end{split}$

and finally each linear expression are constrained to be in rotated quadratic cone i.e.

$2 v_j t_j \geq z_j^2 \mbox{ and } v_j,t_j \geq 0.$

Similarly the second line is equivalent to the constraint

$\begin{split}\left[ \begin{array}{c} z_j\\ 1/8\\ v_j \end{array} \right] \in Q_r^3\end{split}$

or equivalently

$2 z_j \frac{1}{8} \geq v_j^2 \mbox{ and } z_j \geq 0.$

## 7.1.6 Transaction Costs¶

Now assume there is a cost associated with trading asset $$j$$ and the cost is given by

$\begin{split}T_j(\Delta x_j) = \left\lbrace \begin{array}{lll} 0, & \Delta x_j=0, & \\ f_j + g_j |\Delta x_j|,& \mbox{otherwise}. & \end{array}\right.\end{split}$

Delta x_j is the change in the holding of asset $$j$$ i.e.

$\Delta x_j = x_j - x_j^0.$

Hence, whenever asset $$j$$ is traded a fixed cost of $$f_j$$ has to be paid and a variable cost of $$g_j$$ per unit traded. Given the assumptions about transaction costs in this section then problem (5) may be formulated as

(11)$\begin{split}\begin{array}{lllll} \mbox{maximize} & \mu^T x & &\\ \mbox{subject to} & e^T x + \sum_{j=1}^n (f_j y_j + g_j z_j) & = & w + e^T x^0, &\\ & [\gamma;G^T x] & \in & Q^{n+1}, & \\ & z_j & \geq & x_j - x_j^0, & j=1,\ldots,n,\\ & z_j & \geq & x_j^0 - x_j, & j=1,\ldots,n,\\ & z_j & \leq & U_j y_j, & j=1,\ldots,n,\\ & y_j & \in & \{0,1\}, & j=1,\ldots,n, \\ & x & \geq & 0. & \end{array}\end{split}$

First observe that

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

and hence $$z_j$$ is bounded below by $$|\Delta x_j|$$. $$U_j$$ is some a prior chosen upper bound on the amount of trading in asset $$j$$ and therefore if $$z_j>0$$ then $$y_j = 1$$ has to be the case. This implies that the transaction costs for the asset $$j$$ is given by autonomous

$f_j y_j + g_j z_j.$

Example code

The following example code demonstrates how to compute an optimal portfolio when transaction costs are included.

Listing 7.4 Code solve problem (11). Click here to download.
/*
Description:
Extends the basic Markowitz model with a market cost term.

Input:
n: Number of assets
mu: An n dimmensional vector of expected returns
GT: A matrix with n columns so (GT')*GT  = covariance matrix
x0: Initial holdings
w: Initial cash holding
gamma: Maximum risk (=std. dev) accepted
f: If asset j is traded then a fixed cost f_j must be paid
g: If asset j is traded then a cost g_j must be paid for each unit traded

Output:
Optimal expected return and the optimal portfolio

*/
std::shared_ptr<ndarray<double,1>> MarkowitzWithTransactionsCost
( int n,
std::shared_ptr<ndarray<double,1>> mu,
std::shared_ptr<ndarray<double,2>> GT,
std::shared_ptr<ndarray<double,1>> x0,
double                             w,
double                             gamma,
std::shared_ptr<ndarray<double,1>> f,
std::shared_ptr<ndarray<double,1>> g)
{
// Upper bound on the traded amount
std::shared_ptr<ndarray<double,1>> u(new ndarray<double,1>(shape_t<1>(n),w+sum(x0)));

Model::t M = new Model("Markowitz portfolio with transaction costs");  auto M_ = finally([&]() { M->dispose(); });
//M->setLogHandler(new java.io.PrintWriter(System.out));

// Defines the variables. No shortselling is allowed.
Variable::t x = M->variable("x", n, Domain::greaterThan(0.0));

Variable::t z = M->variable("z", n, Domain::unbounded());
// Binary varables
Variable::t y = M->variable("y", n, Domain::binary());

//  Maximize expected return
M->objective("obj", ObjectiveSense::Maximize, Expr::dot(mu,x));

// Invest amount + transactions costs = initial wealth
Domain::equalsTo(w+sum(x0)));

// Imposes a bound on the risk
M->constraint("risk", Expr::vstack( gamma,Expr::mul(GT,x)),
Domain::inQCone());

// z >= |x-x0|
M->constraint("sell", Expr::sub(z,Expr::sub(x0,x)),Domain::greaterThan(0.0));
// Alternatively, formulate the two constraints as