I'm trying to implement the Rockafellar's function described in this paper http://past.rinfinance.com/agenda/2009/yollin_slides.pdf with a risk aversion parameter for my thesis. The function to maximize should be: $$ \max (w):\, w'\mu - \lambda\Big(R_\text{VAR}+\frac{\sum ds}{(1-\beta)s}\Big) $$
The code I have implemented is:
cvarOpt= function(rmat, alpha=0.05, lambda = 10000, rmin=0, wmin=0, wmax=1, weight.sum=1) {
require(Rglpk)
n = ncol(rmat)
# number of assets
s = nrow(rmat)
# number of scenarios i.e. periods
averet = colMeans(rmat)
# creatobjective vector, constraint matrix, constraint rhs
#Amat = rbind(cbind(rbind(1, averet), matrix(data = 0, nrow = 2, ncol = s+1)), cbind(rmat, diag(s), 1))
first <- cbind(matrix(1, nrow = 2, ncol = n), rbind(1, averet), matrix(data = 0, nrow = 2, ncol = s+1))
second <- cbind(matrix(1, nrow = s, ncol = n), rmat, diag(s), 1)
colnames(first) <- 1:(1+2*n+s)
colnames(second) <- 1:(1+2*n+s)
Amat <- rbind(first, second)
objL = c(averet, rep(0, n), rep(lambda/(alpha*s), s), lambda)
bvec= c(weight.sum, rmin, rep(0, s))
# direction vector
dir.vec= c("==", ">=", rep(">=", s))
# bounds on weights
bounds = list(lower = list(ind= 1:n, val= rep(wmin, n)), upper = list(ind = 1:n, val = rep(wmax, n)))
res= Rglpk_solve_LP(obj = objL, mat = Amat, dir = dir.vec, rhs = bvec, types = rep("C", length(objL)), max = TRUE, bounds = bounds)
w = as.numeric(res$solution[1:n])
return(list(w = w, status = res$status))
}
I guess it doesn't work because with a portfolio of only two assets I get always (1, 0) as weights whatever the lambda is while with the standard CVaR formulation I get (0.96, 0.4)...
Thank you in advance.
