## Risk Parity Portfolio Selection using MOSEK Optimizer API¶

Portfolio selection (or portfolio management) is the art and science of making decisions about investment mix and policy, matching investments to objectives, asset allocation for individuals and institutions, and balancing risk against performance.

This tutorial demonstrates the use of exponential cones in modeling logarithmic terms using Python Fusion. It is largerly based on Bai, X., Scheinberg, K., Tutuncu, R. (2013), Least-squares approach to risk parity in portfolio selection. We assume basic familiarity with the Markowitz risk-return model.

# Risk Parity¶

Consider the problem of investing in $n$ different assets. Let $x_i$ be the percentage of money invested in assets $i$, and $\mu_i$ its expected return. Let assume we know the covariance matrix $\Sigma$ that link the assets together. We can then define two global measures of the portfolio performance, expected return:

$$\mu =\sum_{i=1}^n \mu_i x_i,$$

and risk:

$$\sigma^2 = \sum_{i=1}^n \sum_{j=1}^n \Sigma_{i,j}x_i x_j = x^T \Sigma x.$$

For different choices of $x_i$ the investor will get different combinations of $\mu$ and $\sigma^2$. The set of all possible ($\sigma^2$, $\mu$) combinations is called the $\textit{attainable set}$. The theory assumes that investors prefer to minimize risk: given the choice of two portfolios with equal returns, investors will choose the one with the least risk.

In a fully-invested portfolio, every asset $i$ has a contribution in terms of risk. Let $RC_i(x)$ be this contribution, we have that the total risk of the invested portfolio is:

$$\mathcal{R}(x) = \sum_{i=1}^n RC_i(x)$$

where

$$RC_i(x) = x_i \frac{\partial \mathcal{R}(x)}{\partial x_i}$$

Let $b=(b_1,..,b_n)$ be a vector of $\textbf{budgets}$ such as $b_i > 0$ and $\sum_{i=1}^n b_i = 1$. We define a risk parity (or risk budgeting) portfolio as the solution of the following set of equations:

\begin{aligned} &\begin{cases} RC_1(x) = b_1 \mathcal{R}(x)\\ ..\\ RC_i(x) = b_i \mathcal{R}(x)\\ ..\\ RC_n(x) = b_n \mathcal{R}(x) \end{cases} \end{aligned}

From the above equations, we can write the risk parity constraint as:

$$\frac{x_i}{b_i} \frac{\partial \mathcal{R}(x)}{\partial x_i} = \frac{x_j}{b_j} \frac{\partial \mathcal{R}(x)}{\partial x_j} \hspace{1em} \forall i,j$$

The Risk Parity portfolio requires that each asset has the same total contribution to risk, that is:

\begin{aligned} &\frac{x_i}{b_i} \frac{\partial \mathcal{R}(x)}{\partial x_i} = \frac{x_j}{b_j} \frac{\partial \mathcal{R}(x)}{\partial x_j} && \forall i,j\\ &x_i \geq 0, &&b_i > 0 \\ &\sum_{i=1}^n x_i = 1, &&\sum_{i=1}^n b_i = 1 \end{aligned}

If we set

$$\mathcal{R}(x) = \sqrt{x^T\Sigma x}$$

solving the following problem that incorporates a logarithmic barrier in the objective function is equivalent to find a Risk Parity solution:

$$$$\label{eq:11} \begin{array}{l} \min_x \frac{1}{2} x^T \Sigma x - c \sum_{i=1}^{n} b_i ln(x_i)\\ \text{s.t.}\\ x > 0 \end{array}$$$$

where $b_i > 0$, $\sum_i b_i = 1$ and $c$ is a positive constant.

### Proof¶

Since $\Sigma$ is positive semidefinite and the logarithm function is strictly concave, the objective function is $\textbf{strictly convex}$. From the first order condition, the unique solution is in corrispondence of the point where the gradient of the objective function is zero:

$$\Sigma x - c b_i x^{-1} = 0$$

Hence, at optimality we have

$$(\Sigma x)_i = \frac{c b_i}{x_i} \Rightarrow \frac{x_i(\Sigma x)_i}{b_i} = \frac{x_j(\Sigma x)_j}{b_j}, \quad \forall i,j,$$

since

$$\frac{\partial \mathcal{R}(x)}{\partial x_i} = (\Sigma x)_i.$$

# MOSEK Implementation¶

We begin by translating the problem in conic form. Assume that we have a factorization $\Sigma=F^TF$, then $x^T\Sigma x=\|Fx\|_2^2$ and we can write an equivalent model as:

$$\begin{array}{ll} \min_x & r - c b^Tt \\ \text{s.t.}& (1,r,Fx)\in \mathcal{Q}_r\\ & t_i\leq \mathrm{ln} x_i,\ i=1,\ldots,n \\ & x \geq 0 \end{array}$$

Indeed, the first constraint involving the rotated quadratic cone $\mathcal{Q}_r$ means

$$2\cdot 1\cdot r\geq \|Fx\|^2 = x^T \Sigma x$$

and moreover each inequality $t_i\leq \mathrm{ln} x_i$ can be expressed with an exponential cone as

$$(x_i, 1, t_i)\in K_{\mathrm{exp}}.$$

We proceed with implementation of this model in Fusion:

In [1]:
from mosek.fusion import *

def parityModel(n, F, c, b):
M = Model('parity')

# Define all variables appearing in the model
x = M.variable('x', n)
t = M.variable(n)
r = M.variable()

# Objective r-cb^Tt
M.objective(ObjectiveSense.Minimize, Expr.sub(r, Expr.mul(c, Expr.dot(b,t))))

# Bound on risk - construct the vector (1,r,Fx)
M.constraint(Expr.vstack(1, r, Expr.mul(F,x)), Domain.inRotatedQCone())

#Logarithmic bounds - all together in matrix form
M.constraint(Expr.hstack(x, Expr.constTerm(n,1.0), t), Domain.inPExpCone())

# That's all
return M


Now, we create some sample data.

In [2]:
import numpy as np

n = 20
k = 8

F = np.random.sample([k,n])
Sigma = np.dot(F.transpose(), F)   # Not used

b = np.random.sample(n)
b = b/np.sum(b)

c = 1.0

PM = parityModel(n, F, c, b)
import sys
PM.setLogHandler(sys.stdout)
PM.solve()

Problem
Name                   : parity
Objective sense        : min
Type                   : CONIC (conic optimization problem)
Constraints            : 70
Cones                  : 21
Scalar variables       : 111
Matrix variables       : 0
Integer variables      : 0

Optimizer started.
Presolve started.
Linear dependency checker started.
Linear dependency checker terminated.
Eliminator started.
Freed constraints in eliminator : 21
Eliminator terminated.
Eliminator started.
Freed constraints in eliminator : 0
Eliminator terminated.
Eliminator - tries                  : 2                 time                   : 0.00
Lin. dep.  - tries                  : 1                 time                   : 0.00
Lin. dep.  - number                 : 0
Presolve terminated. Time: 0.01
Problem
Name                   : parity
Objective sense        : min
Type                   : CONIC (conic optimization problem)
Constraints            : 70
Cones                  : 21
Scalar variables       : 111
Matrix variables       : 0
Integer variables      : 0

Optimizer  - solved problem         : the primal
Optimizer  - Constraints            : 8
Optimizer  - Cones                  : 21
Optimizer  - Scalar variables       : 70                conic                  : 70
Optimizer  - Semi-definite variables: 0                 scalarized             : 0
Factor     - setup time             : 0.00              dense det. time        : 0.00
Factor     - ML order time          : 0.00              GP order time          : 0.00
Factor     - nonzeros before factor : 36                after factor           : 36
Factor     - dense dim.             : 0                 flops                  : 3.11e+03
ITE PFEAS    DFEAS    GFEAS    PRSTATUS   POBJ              DOBJ              MU       TIME
0   1.5e+01  1.3e+00  2.5e+00  0.00e+00   1.534945180e+00   0.000000000e+00   1.0e+00  0.02
1   9.5e-01  8.0e-02  1.6e-01  6.74e-02   3.799364477e+00   4.071561966e+00   6.2e-02  0.03
2   1.6e-01  1.3e-02  2.6e-02  8.65e-01   3.913169811e+00   3.949278388e+00   1.0e-02  0.03
3   6.1e-02  5.1e-03  1.0e-02  1.08e+00   3.709041715e+00   3.717025112e+00   4.0e-03  0.03
4   1.1e-02  9.2e-04  1.8e-03  1.03e+00   3.639500611e+00   3.639166289e+00   7.1e-04  0.03
5   1.0e-03  8.5e-05  1.7e-04  1.00e+00   3.638115939e+00   3.638070284e+00   6.6e-05  0.03
6   7.0e-05  5.9e-06  1.2e-05  1.00e+00   3.638231209e+00   3.638228036e+00   4.5e-06  0.03
7   7.4e-06  6.2e-07  1.2e-06  1.00e+00   3.638236463e+00   3.638236125e+00   4.8e-07  0.03
8   9.8e-07  8.2e-08  1.6e-07  1.00e+00   3.638236859e+00   3.638236815e+00   6.3e-08  0.03
9   1.3e-07  1.1e-08  2.2e-08  1.00e+00   3.638236836e+00   3.638236830e+00   8.5e-09  0.03
10  2.8e-08  9.4e-09  8.9e-10  1.00e+00   3.638236825e+00   3.638236824e+00   3.5e-10  0.04
11  4.6e-09  1.6e-09  1.5e-10  1.00e+00   3.638236816e+00   3.638236816e+00   5.9e-11  0.04
Optimizer terminated. Time: 0.04

Interior-point solution summary
Problem status  : PRIMAL_AND_DUAL_FEASIBLE
Solution status : OPTIMAL
Primal.  obj: 3.6382368157e+00    nrm: 6e+00    Viol.  con: 2e-09    var: 0e+00    cones: 8e-10
Dual.    obj: 3.6382368138e+00    nrm: 2e+00    Viol.  con: 0e+00    var: 3e-09    cones: 0e+00


We can access, normalize and plot the optimal solution

In [3]:
xx = PM.getVariable('x').level()
xx = xx/sum(xx)

%matplotlib inline
import matplotlib.pyplot as plt

plt.bar(np.arange(n), xx)
plt.show()


We can also verify that the total contribution of each assed towards total risk is proportional to the budgeting vector $b$, as required:

In [4]:
RC = xx * (np.dot(Sigma,xx))
plt.bar(np.arange(n), RC/b)
plt.show()


This work is licensed under a Creative Commons Attribution 4.0 International License. The MOSEK logo and name are trademarks of Mosek ApS. The code is provided as-is. Compatibility with future release of MOSEK or the Fusion API are not guaranteed. For more information contact our support.