# CVAR alternatives for optimization

Are there some alternatives to the CVaR measure for portfolio optimization, which are easier to implement for ex. with a linear program? They can be just approximations of CVaR or measures representing similar concepts.

Edit: I would prefer some exact (analytical) methods over (meta-)heuristics, as they seem to be quite sensitive to changes (for my case). If not, other approaches to make heuristics more reliable/robust are also welcome.

Edit2: If there are no real approximations to CVaR, or if a more detailed explanation (of all the steps) of optimizing CVaR can be provided, that is also on topic, as that would allow us to use standard and analytical methods (such as LP).

• IIRC, CVaR can be implemented as a linear programme. Look for 'Optimization of conditional value-at-risk' by Rockafellar/Uryasev. But in general, you should differentiate between the model you want to solve (in your case, the specific objective function); and methods to solve the model. There might be several methods, and some are simpler/faster/whatever than others. – Enrico Schumann Apr 2 '20 at 6:39
• @EnricoSchumann Yes, CVaR can be implemented as a linear programme, however I have not been able to do so because the underlying mathematical model is too complex for me, hence why I was looking for alternatives or approximations. I am looking to find the global optimum of some objective function, hence why I suggested linear programming, as I suspect it will probably involve this. But if not, other solutions are also welcome. – user2974951 Apr 2 '20 at 7:01
• In general, the sensitivity of the solutions is not a result of using heuristics; it is result of the properties of financial data. In papers.ssrn.com/sol3/papers.cfm?abstract_id=2698114 we explicitly looked into this; see Figures 3 and 4. The randomness that comes from a (properly working) heuristic is negligible compared with the randomness caused by changes in the data setup, such as using different window lengths or adding a few observations to the sample. – Enrico Schumann Apr 2 '20 at 7:52
• If any optimization expert has a lot of free time nowadays, a clear "how to" writeup about implementation of Rockafellar/Uryasev CvaR optimization by LP would be very useful, maybe with some pseudocode or step by step description of the algo ;). As presented now it is not vet clear, and there have been other Q's about it in this forum. – noob2 Apr 2 '20 at 12:15
• @noob2: I have drafted some notes on implementing the CVaR model: enricoschumann.net/notes/minimising-conditional-var.html . But I suppose that would be off-topic to the question. But If the OP considers it on-topic, I could post this here as well. – Enrico Schumann Apr 6 '20 at 20:28

Following the comments and the edits to the question, I'll try to show how conditional Value-at-Risk (aka Expected Tail Loss) can be minimised for a portfolio. 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. The chosen quantile is called $$\beta$$ in the paper. 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:

$$\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}$$

That is, a vector of a 1, followed by $$n_{\mathrm{A}}$$ zeros, and then $$n_{\mathrm{S}}$$ times a constant.

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

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


(Disclosure: I am the maintainer of packages NMOF and 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.

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    ...    u9   u10
##   1.0  0.0   0.0   0.0   0.5   0.5   0.5    ...   0.5   0.5


The constraints matrix gains one column for each scenario.

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))
## GLPK Simplex Optimizer, v4.65
## 5001 rows, 5021 columns, 110020 non-zeros
## ...
## OPTIMAL LP SOLUTION FOUND


We store the resulting weights in a variable lp.weights.

lp.weights <- sol.lp$solution[2:(1+na)]  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, type = "numeric", stepsize = 0.05)  To give an idea, suppose we had a portfolio consisting of only three assets, and each asset has a weight of one-third. Then we could compute neighbours as follows: nb(rep(1/3, 3)) ## 0.3122 0.3544 0.3333 nb(rep(1/3, 3)) ## [1] 0.3272 0.3333 0.3394  Now that we have an objective and a neighbourhood function, we can 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


So both implementations give very similar results.

As for instability: We can run the heuristic several times (we always should; see http://enricoschumann.net/R/remarks.htm). Let us run it 20 times.

sols.ls <- restartOpt(LSopt,
n = 20,
OF = CVaR,
list(x0 = rep(1/na, na),
neighbour = nb,
nI = 1000),
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


So we get very similar results in all runs.

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 we 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),
R = R,
exp = 2)

• This was by far the simplest and easiest to understand example on the topic that I've found. Very nice. – user2974951 Apr 8 '20 at 8:32
• An update: The development version of the NMOF package now has a function minCVaR that implements the LP approach. – Enrico Schumann May 24 '20 at 7:57