Worst-case risk analysis

Covariance uncertainty

In this example we do worst-case risk analysis using CVXPY. Our setting is a single period Markowitz portfolio allocation problem. We have a fixed portfolio allocation $w \in {\bf R}^n$. The return covariance $\Sigma$ is not known, but we believe $\Sigma \in \mathcal S$. Here $\mathcal S$ is a convex set of possible covariance matrices. The risk is $w^T \Sigma w$, a linear function of $\Sigma$.

We can compute the worst (maximum) risk, over all possible covariance matrices by solving the convex optimization problem

$$ \begin{array}{ll} \mbox{maximize} & w^T\Sigma w \\ \mbox{subject to} & \Sigma \in \mathcal S, \quad \Sigma \succeq 0, \end{array} $$

with variable $\Sigma$.

If the worst-case risk is not too bad, you can worry less. If not, you'll confront your worst nightmare

Example

In the following code we solve the portfolio allocation problem

$$ \begin{array}{ll} \mbox{minimize} & w^T\Sigma_\mathrm{nom} w \\ \mbox{subject to} & {\bf 1}^Tw = 1, \quad \mu^Tw \geq 0.1, \quad \|w\|_1 \leq 2, \end{array} $$

and then compute the worst-case risk under the assumption that $\mathcal S = \left\{ \Sigma^\mathrm{nom} + \Delta \,:\, |\Delta_{ii}| =0, \; |\Delta_{ij}| \leq 0.2 \right\}$.

We might expect that $|\Delta_{ij}| = 0.2$ for all $i \neq j$. This does not happen however because of the constraint that $\Sigma^\mathrm{nom} + \Delta$ is positive semidefinite.

In [1]:
# Generate data for worst-case risk analysis.
import numpy as np
np.random.seed(2)
n = 5
mu = np.abs(np.random.randn(n, 1))/15
Sigma = np.random.uniform(-.15, .8, size=(n, n))
Sigma_nom = Sigma.T.dot(Sigma)
print("Sigma_nom =")
print(np.round(Sigma_nom, decimals=2))
Sigma_nom =
[[ 0.58  0.2   0.57 -0.02  0.43]
 [ 0.2   0.36  0.24  0.    0.38]
 [ 0.57  0.24  0.57 -0.01  0.47]
 [-0.02  0.   -0.01  0.05  0.08]
 [ 0.43  0.38  0.47  0.08  0.92]]
In [2]:
# Form and solve portfolio optimization problem.
# Here we minimize risk while requiring a 0.1 return.
import cvxpy as cp


w = cp.Variable(n)
ret = mu.T*w 
risk = cp.quad_form(w, Sigma_nom)
prob = cp.Problem(cp.Minimize(risk), 
               [cp.sum(w) == 1, 
                ret >= 0.1,
                cp.norm(w, 1) <= 2])
prob.solve()
print("w =")
print(np.round(w.value, decimals=2))
w =
[-0.01  0.13  0.18  0.88 -0.18]
In [3]:
# Form and solve worst-case risk analysis problem.
Sigma = cp.Variable((n, n), PSD=True)
Delta = cp.Variable((n, n), symmetric=True)
risk = cp.quad_form(w.value, Sigma)
prob = cp.Problem(cp.Maximize(risk), 
               [Sigma == Sigma_nom + Delta, 
                cp.diag(Delta) == 0, 
                cp.abs(Delta) <= 0.2])
prob.solve()
print("standard deviation =", cp.sqrt(cp.quad_form(w.value, Sigma_nom)).value)
print("worst-case standard deviation =", cp.sqrt(risk).value)
print("worst-case Delta =")
print(np.round(Delta.value, decimals=2))
standard deviation = 0.16889580642544394
worst-case standard deviation = 0.421993016105032
worst-case Delta =
[[ 0.    0.04 -0.2  -0.    0.2 ]
 [ 0.04  0.    0.2   0.09 -0.2 ]
 [-0.2   0.2   0.    0.12 -0.2 ]
 [-0.    0.09  0.12  0.   -0.18]
 [ 0.2  -0.2  -0.2  -0.18  0.  ]]