The value-at-risk (VaR) of a random variable $X$ at risk level $\alpha$ is defined as $$\textrm{VaR}_\alpha(X) = \min_\eta\left(\mathbb{P}[X\geq\eta]\leq 1-\alpha\right)$$ and the conditional value-at-risk (CVaR) is the conditional expected value $$\textrm{CVaR}_\alpha(X)=\mathbb{E}(X~|~X\geq\textrm{VaR}_\alpha(X)).$$ In the standard setting where the random variable $X$ represents a loss (on an investment) and we take $\alpha=0.9$, CVaR can be interpreted as the expected loss assuming we end up in the $10\%$ worst scenarios. The concept was studied in the paper Optimization of Conditional Value-at-Risk, Rockafellar and Uryasev, Journal of Risk, 2002, where we also find an equivalent formula $$\textrm{CVaR}_\alpha(X)=\min_\eta\left(\eta+(1-\alpha)^{-1}\mathbb{E}([X-\eta]^+)\right)$$ where $[x]^+=\max(x,0)$.
Some authors consider generalizations $$\rho^f_\alpha(X)=\min_\eta\left(\eta+(1-\alpha)^{-1}f^{-1}\left(\mathbb{E}(f(X-\eta))\right)\right),$$ where $f$ is any increasing, convex function on $\mathbb{R}_+$ which is identically zero on $\mathbb{R}_-$. Examples include:
In this notebook we implement the models described in Section 5 of (Vinel, Krokhmal 2017) with various risk measures. We consider the problem of investing in $n$ stocks given $m$ historical return vectors $r_1,\ldots,r_m\in\mathbb{R}^n$ taken with probabilities $p_1,\ldots,p_m$. We consider the following conditions:
We used a preprocessed file with closing prices from NYSE for 1675 assets over 2050 trading days. Each historical scenario is defined by computing returns of all assets over some period of 10 trading days They are returned in a martix $r$ with $n$ rows and $m$ columns.
%matplotlib inline
import matplotlib
import matplotlib.pyplot as plt
from mosek.fusion import *
import numpy as np
import sys, os, csv
allPrices = np.array(list(csv.reader(open("allPrices.txt"), lineterminator='\n', quoting=csv.QUOTE_NONNUMERIC)))
n,m = allPrices.shape
p1 = allPrices[:,0:m-10]
p2 = allPrices[:,10:m]
r = (p2-p1)/p1
print(r.shape)
(1675, 2040)
We first implement the linear CVaR model with $f(x)=[x]^+$. We can then write the problem as
$$ \begin{array}{ll} \mbox{minimize} & \eta + (1-\alpha)^{-1} p^Tw\\ \mbox{subject to} & 1^Tx\leq 1, \\ & x^Trp\geq \bar{r}, \\ & w_j\geq \max\{0, -r_j^Tx-\eta \}, \quad j=1,\ldots,m, \end{array} $$with variables $x\in\mathbb{R}^n, w\in\mathbb{R}^m, \eta\in\mathbb{R}$. Below is the Fusion implementation.
def modelCVaR(alpha, rbar, p, r, m, n):
M = Model("CVaR")
# Variables
X = M.variable('X', n, Domain.greaterThan(0.0))
W = M.variable('W', m, Domain.greaterThan(0.0))
eta = M.variable('eta')
# The bounds w_j + r_j^Tx + eta \geq 0
M.constraint(Expr.add([W, Expr.mul(np.transpose(r), X), Var.vrepeat(eta, m)]), Domain.greaterThan(0.0))
# Minimal risk
M.constraint(Expr.dot(X, np.matmul(r,p)), Domain.greaterThan(rbar))
# Budget
M.constraint(Expr.sum(X), Domain.lessThan(1.0))
# Set the objective
M.objective(ObjectiveSense.Minimize, Expr.add(eta, Expr.mul(1/(1-alpha), Expr.dot(p, W))))
return M
Next we have $f(x)=\exp([x]^+)-1$. Since $\sum_j p_j=1$, the model can equivalently be formulated as
$$ \begin{array}{ll} \mbox{minimize} & \eta + (1-\alpha)^{-1} t\\ \mbox{subject to} & 1^Tx\leq 1, \\ & x^Trp\geq \bar{r}, \\ & w_j\geq \max\{0, -r_j^Tx-\eta \}, \quad j=1,\ldots,m, \\ & t\geq \log\left(\sum_j p_j\exp(w_j)\right). \end{array} $$The last bound (log-sum-exp constraint) can be expressed using an auxiliary variable $\xi\in\mathbb{R}^m$:
$$ \begin{array}{l} p^T\xi \leq 1, \\ \xi_j\geq \exp(w_j-t), \end{array} $$where the last inequality is conic-representable as $(\xi_j,1,w_j-t)\in K_\mathrm{exp}$.
def modelLogExpCR(alpha, rbar, p, r, m, n):
M = Model("LogExpCR")
# Variables
X = M.variable('X', n, Domain.greaterThan(0.0))
W = M.variable('W', m, Domain.greaterThan(0.0))
eta = M.variable('eta')
xi = M.variable('xi', m)
t = M.variable('t')
# The log-sum-exp constraint
M.constraint(Expr.dot(p, xi), Domain.lessThan(1.0))
# Stack allcones together - every row has the form (xi_j, 1, w_j-t)
# The matrix notation means that every row belongs to the cone
M.constraint(Expr.hstack(xi, Expr.constTerm(m, 1.0), Expr.sub(W, Var.vrepeat(t, m))), Domain.inPExpCone())
# Continue with the linear constraints
M.constraint(Expr.add([W, Expr.mul(np.transpose(r), X), Var.vrepeat(eta, m)]), Domain.greaterThan(0.0))
M.constraint(Expr.dot(X, np.matmul(r,p)), Domain.greaterThan(rbar))
M.constraint(Expr.sum(X), Domain.lessThan(1.0))
# Set the objective
M.objective(ObjectiveSense.Minimize, Expr.add(eta, Expr.mul(1/(1-alpha), t)))
return M
The model using $f(x)=([x]^+)^q$ is the same as LogExpCR, except that the last constraint becomes
$$t\geq\left(\sum_j p_jw_j^q\right)^{1/q}.$$This $q$-norm cone constraint is equivalent to
$$ \begin{array}{l} t = 1^T\xi,\\ t^{1-1/q}\xi_j^{1/q}\geq |p_j^{1/q}w_j| \end{array} $$where the last row is precisely a conic constraint involving the $3$-dimenional power cone with parameters $(1-1/q,1/q)$.
def modelHMCR(alpha, rbar, p, r, m, n, q):
M = Model("HMCR")
# Variables
X = M.variable('X', n, Domain.greaterThan(0.0))
W = M.variable('W', m, Domain.greaterThan(0.0))
eta = M.variable('eta')
xi = M.variable('xi', m)
t = M.variable('t')
# The power cones representing the q-norm constraint
M.constraint(Expr.sub(t, Expr.sum(xi)), Domain.equalsTo(0.0))
M.constraint(Expr.hstack(Var.vrepeat(t, m), xi, Expr.mulElm(np.power(p, 1.0/q), W)), Domain.inPPowerCone(1.0-1.0/q))
# Continue with the linear constraints
M.constraint(Expr.add([W, Expr.mul(np.transpose(r), X), Var.vrepeat(eta, m)]), Domain.greaterThan(0.0))
M.constraint(Expr.dot(X, np.matmul(r,p)), Domain.greaterThan(rbar))
M.constraint(Expr.sum(X), Domain.lessThan(1.0))
# Set the objective
M.objective(ObjectiveSense.Minimize, Expr.add(eta, Expr.mul(1/(1-alpha), t)))
return M
We solve a few problems with $n=100$ stocks and $m=2000$ sample historical scenarios.
n, m = 100,2000
alpha = 0.9
rbar = 0.005
p = [1.0/m]*m # Uniform probability
MCVaR = modelCVaR(alpha, rbar, p, r[0:n,0:m], m, n)
MLogExpCR = modelLogExpCR(alpha, rbar, p, r[0:n,0:m], m, n)
MHMCR = modelHMCR(alpha, rbar, p, r[0:n,0:m], m, n, q=5.9)
for M in [MCVaR, MLogExpCR, MHMCR]:
M.setSolverParam("numThreads", 4)
M.solve()
print('{0: >8}: status={1} time={2:.2f}s'.format(M.getName(), M.getProblemStatus(), M.getSolverDoubleInfo("optimizerTime")))
plt.plot(M.getVariable('X').level(), 'o')
plt.show()
CVaR: status=ProblemStatus.PrimalAndDualFeasible time=0.20s
LogExpCR: status=ProblemStatus.PrimalAndDualFeasible time=0.69s
HMCR: status=ProblemStatus.PrimalAndDualFeasible time=0.64s
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.