Linear Programming Formulation for CVaR Portfolio Optimization

# CVAR alternatives for optimization


# 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).

## Answer by Enrico Schumann (score 7, accepted)

https://quant.stackexchange.com/a/53111

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. This objective function looks at the left tail of the returns distribution. Minimising CVaR means shifting this left tail to the right, so we need a minus in front of the `mean`. (For loss scenarios, i.e. negative returns, we could have written `mean(Rw[Rw >= quantile(Rw, b)])`.)

```
CVaR <- function(w, R, b) {
      Rw <- R %*% w   ## portfolio return under scenarios
   ## mean(Rw[Rw >= quantile(Rw, b)])      ## ==> for losses
     -mean(Rw[Rw <= quantile(Rw, 1 - b)])  ## ==> for returns
}
```

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 the same, 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. Also, the are several valid ways to compute quantiles (see `?quantile`).

```
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 https://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   ## portfolio return under scenarios
      pm(Rw, xp = exp)
}
sol.ls <- LSopt(PM,
                list(x0 = rep(1/na, na),
                     neighbour = nb,
                     nI = 1000),
                R = R,
                exp = 2)
```