Minimising Conditional Value-at-Risk (CVaR)
This note was motivated by a question on Stack Exchange. We will look at minimising conditional Value-at-Risk (aka Expected Tail Loss). We start with the implementation suggested by Rockafellar/Uryasev:
@ARTICLE{,
author = {R. Tyrrell Rockafellar and Stanislav Uryasev},
title = {Optimization of Conditional Value-at-Risk},
journal = {Journal of Risk},
year = 2000,
volume = 2,
number = 3,
pages = {21--41},
doi = {10.21314/JOR.2000.038}
}
Their description of the model as a linear programme (LP) is given after Equation 17 in their paper. We start with a scenario set, i.e. a sample of possible return realisations. Note that in the paper, Rockafellar/Uryasev work with losses, i.e. negative returns. A 'bad' quantile is thus a high one, say the 90th. We store the scenarios in a matrix \(R\). It has \(n_{\mathrm{A}}\) columns (one for each asset) and \(n_{\mathrm{S}}\) rows (one row for each scenario).
The variables in the model are the actual portfolio weights \(x\), plus auxiliary variables \(u\) (one for each scenario), plus the VaR level \(\alpha\). The solver may choose the VaR level and the weights together in such a way that the CVaR is minimised. The objective function only has auxiliary variables and the VaR level in it; the actual portfolio weights enter model model only through the constraints.
The weights in the objective function look as follows: \[ \left[ {\begin{array}{ccccccc} \underbrace{\begin{matrix}\phantom{000}1\phantom{000} \end{matrix}}_{\alpha}& \underbrace{\begin{matrix}0 & \cdots & 0\phantom{0}\end{matrix}}_{x} & \underbrace{\begin{matrix}\frac{1}{(1-\beta) n_{\mathrm{S}}} & \cdots & \frac{1}{(1-\beta) n_{\mathrm{S}}}\end{matrix}}_{u}\\ \end{array} } \right] \]
The constraints matrix looks as follows: \[ {\begin{array}{ccccccccc} 0 & 1 & \cdots & 1 & 0 & 0 & \cdots & 0 & = &1 \\ 1 & r_{1,1} & \cdots & r_{1, n_\mathrm{A}} & 1 & 0 & \cdots & 0 & \geq & 0 \\ 1 & r_{2,1} & \cdots & r_{2, n_\mathrm{A}} & 0 & 1 & \cdots & 0 & \geq & 0 \\ \vdots \\ \underbrace{\phantom{00}1\phantom{00}}_{\alpha} & \underbrace{\phantom{00}r_{n_{\mathrm{S}},1}\phantom{00}}_{x_1} & \cdots & \underbrace{r_{n_\mathrm{S}, n_\mathrm{A}}}_{x_{n_\mathrm{A}}} & \underbrace{0}_{u_1} & \underbrace{0}_{u_2} & \cdots & \underbrace{1}_{u_{n_\mathrm{S}}} & \geq & 0 \\ \end{array} } %\alpha & x_1 & \cdots & x_{n_\mathrm{A}} & u_1 & \cdots & u_{n_\mathrm{S}} \\ \]
Note that the first line is the budget constraint. Other than that line, the matrix consists of a column of ones (for the VaR variable), the scenario matrix \(R\), and an identity matrix of dimension \(n_\mathrm{S}\). There are non-negativity constraints for all \(x\) and \(u\), but they are not explicitly shown. Many solvers automatically enforce them.
We can try an implementation in R. We load the packages required for the examples.
library("Rglpk") library("NMOF") ## https://github.com/enricoschumann/NMOF library("neighbours") ## https://github.com/enricoschumann/neighbours
We start with a small dataset, so that we can look at the objective function and the constraints matrix: only 3 assets and 10 scenarios.
options(width=200) set.seed(2476) ns <- 10 ## number of scenarios na <- 3 ## number of assets R <- randomReturns(na, ns, sd = 0.01) ## an array of size ns x na b <- 0.8 ## beta in the original paper
The objective function gives zero weights to the x.
f.obj <- c(alpha = 1,
x = rep(0, na),
u = 1/rep((1 - b)*ns, ns))
f.obj
alpha x1 x2 x3 u1 u2 u3 u4 u5 u6 u7 u8 u9 u10 1.0 0.0 0.0 0.0 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5
The constraints matrix gains one column for each scenario.
options(digits=3) C <- cbind(1, R, diag(nrow(R))) C <- rbind(c(alpha = 0, x = rep(1, na), u = rep(0, nrow(R))), C) C
alpha x1 x2 x3 u1 u2 u3 u4 u5 u6 u7 u8 u9 u10
[1,] 0 1.000000 1.000000 1.00000 0 0 0 0 0 0 0 0 0 0
[2,] 1 0.000183 0.029174 0.00293 1 0 0 0 0 0 0 0 0 0
[3,] 1 -0.001776 -0.001673 0.00225 0 1 0 0 0 0 0 0 0 0
[4,] 1 -0.009948 -0.007892 0.01129 0 0 1 0 0 0 0 0 0 0
[5,] 1 0.008299 -0.005601 -0.00144 0 0 0 1 0 0 0 0 0 0
[6,] 1 0.005766 0.000521 -0.00940 0 0 0 0 1 0 0 0 0 0
[7,] 1 -0.017110 0.016782 -0.00122 0 0 0 0 0 1 0 0 0 0
[8,] 1 -0.008334 0.017317 0.00498 0 0 0 0 0 0 1 0 0 0
[9,] 1 -0.004077 -0.009600 0.01568 0 0 0 0 0 0 0 1 0 0
[10,] 1 -0.000532 -0.000201 0.00267 0 0 0 0 0 0 0 0 1 0
[11,] 1 -0.005090 0.002318 0.00368 0 0 0 0 0 0 0 0 0 1
Let us run it on a larger model.
ns <- 5000 na <- 20 R <- randomReturns(na, ns, sd = 0.01, rho = 0.5) b <- 0.75 f.obj <- c(alpha = 1, x = rep(0, na), u = 1/rep(( 1 - b)*ns, ns)) C <- cbind(1, R, diag(nrow(R))) C <- rbind(c(alpha = 0, x = rep(1, na), u = rep(0, nrow(R))), C) const.dir <- c("==", rep(">=", nrow(C) - 1)) const.rhs <- c(1, rep(0, nrow(C) - 1)) sol.lp <- Rglpk_solve_LP(f.obj, C, const.dir, rhs = const.rhs, control = list(verbose = TRUE, presolve = TRUE)) lp.weights <- sol.lp$solution[2:(1+na)]
GLPK Simplex Optimizer, v4.65
5001 rows, 5021 columns, 110020 non-zeros
Preprocessing...
4736 rows, 4756 columns, 104190 non-zeros
Scaling...
A: min|aij| = 7.231e-08 max|aij| = 1.000e+00 ratio = 1.383e+07
GM: min|aij| = 3.425e-03 max|aij| = 2.919e+02 ratio = 8.523e+04
EQ: min|aij| = 1.173e-05 max|aij| = 1.000e+00 ratio = 8.523e+04
Constructing initial basis...
Size of triangular part is 4736
0: obj = 0.000000000e+00 inf = 1.705e+02 (2507)
2775: obj = 1.894482283e-02 inf = 0.000e+00 (0) 27
* 4295: obj = 9.268423618e-03 inf = 0.000e+00 (0) 15
OPTIMAL LP SOLUTION FOUND
Now let us solve the same model with a heuristic. We can use one of the simplest methods, a (stochastic) Local Search. And I will also leave the implementation here simple (it could be improved to gain speed).
Local Search is straightforward: we start with a portfolio, e.g. one with equal weights for all assets. Then we evaluate this portfolio – see how good it is. So we write an objective function that, given the data, maps a portfolio to a number – the CVaR.
CVaR <- function(w, R, b) { Rw <- R %*% w ## compute portfolio loss under scenarios mean(Rw[Rw >= quantile(Rw, b)]) }
And now we run the heuristic: it will loop through the space of possible portfolios and accept portfolios that are better than the current one, but reject those portfolios that are worse. For this loop, we need a second ingredient (the first was the objective function): the neighbourhood function. The neighborhood function takes a portfolio as input and returns a slightly-modified copy of this portfolio.
nb <- neighbourfun(0, 1, length = na, type = "numeric", stepsize = 0.05)
We run the Local Search.
sol.ls <- LSopt(CVaR,
list(x0 = rep(1/na, na),
neighbour = nb,
nI = 1000),
R = R, b = b)
We can compare the objective function values. Note that the objective-function definitions are not exactly equivalent, since the LP may choose the VaR and the weights together, whereas the Local Search only varies the weights and then, in the objective function, imposes a way to compute the VaR.
CVaR(sol.ls$xbest, R, b) CVaR(lp.weights, R, b)
[1] 0.00946 [1] 0.00955
As to instability: We can run the method several times (we always should; see http://enricoschumann.net/R/remarks.htm). Let us run it 20 times.
options(width=200) sols.ls <- restartOpt(LSopt, n = 20, OF = CVaR, list(x0 = rep(1/na, na), neighbour = nb, nI = 1000, printBar = FALSE, printDetail = FALSE), R = R, b = b) summary(sapply(sols.ls, `[[`, "OFvalue"))
Min. 1st Qu. Median Mean 3rd Qu. Max. 0.00946 0.00946 0.00946 0.00946 0.00946 0.00946
A nice thing about the heuristic is that we can easily change it. Suppose you did not want CVaR any more; but now we preferred to minimise a partial moment, say. Then all you would have to do is write another objective function.
PM <- function(w, R, exp = 2, ...) { Rw <- R %*% w ## compute portfolio loss under scenarios pm(Rw, xp = exp, lower = FALSE) ## we work with losses } sol.ls <- LSopt(PM, list(x0 = rep(1/na, na), neighbour = nb, nI = 1000, printBar = FALSE), R = R, exp = 2)
Local Search. Initial solution: 2.78e-05 Finished. Best solution overall: 2.76e-05