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)
