Efficient CVaR Scenario Optimization for a Large number of Scenarios
Great- Thanks again for all of your help
On Feb 16, 2012, at 12:55 PM, alexios ghalanos <alexios at 4dscape.com> wrote:
Not sure but this sounds to me like the Lower Partial Moment problem (but i could be mistaken)...if it is, you should read up on the literature on this measure as it has been extensively covered (for the problem of the partial moment to the power 1 there is an LP representation, for the power 2 a QP representation, and for the rest or general LPM problem you can definitely represent it as a proper NLP). Regards, Alexios On 16/02/2012 15:55, Robert Harlow wrote:
Alexios,
I have one more question on this topic. With CVaR, one generally
knows the percentile, and then has to find the quantile that that
percentile corresponds to (VaR). Suppose that my objective function is
slightly different. In the following case, I have a fixed hurdle
return, and I would like to minimize the average return that is less
than that hurdle return multiplied by the likelihood of having a return
that is lower than the hurdle. Basically, this is just the expected
value of the return from -Inf to the hurdle (rather than -Inf to
positive Inf for the total expected value). This seems like it should
be an easier problem to solve than the CVaR one, but for some reason I
find it difficult. Is there a better way to write my objective function
(cvarStar) such that it nicely follows the process you already laid out
in the CVaR case?
Thanks again,
-Bob
library(mvtnorm)
mu <- c(.04, .06, .1)
sig <- c(.04, .09, .2)
cor.mat <- rbind(c( 1, .6, .2),
c(.6, 1, .1),
c(.2, .1, 1)
)
v.cov <- sig%o%sig*cor.mat
set.seed(100)
n <- 1000
samps <- rmvnorm(n, mean = mu, sigma = v.cov)
cvarStar <- function(w, samps, hurd){
z <- samps %*%w
tail.v <- z[z<=hurd]
if(length(tail.v) == 0) tail.v <- 0
-(length(tail.v)/length(z))*mean(tail.v)
}
cvarStar(rep(1/3,3), samps, -.05)
On Wed, Feb 15, 2012 at 6:01 PM, Robert Harlow <rharlow86 at gmail.com
<mailto:rharlow86 at gmail.com>> wrote:
Wow, this is fantastic Alexios. Thank you very much.
-Bob
On Wed, Feb 15, 2012 at 5:11 PM, alexios <alexios at 4dscape.com
<mailto:alexios at 4dscape.com>> wrote:
Hi Bob,
There is a way to work with nonlinear solvers to solve this type
of problem (without resorting to expensive global optimization
approaches), once you work out the derivatives and decide on how
to deal with the min/max type functions so that the problem is
'made convex'.
A complete example follows:
# We work with a smooth approximation to the max function:
epsilon = 1e-8
.max.smooth = function(x)
{
(sqrt(x*x + epsilon ) + x) / 2
}
# the objective function (CVaR)
func.cvar = function(w, Data, alpha, mu, rmin)
{
# we reverse the sign since we are minimizing
-w[1] + mean( .max.smooth(w[1] - Data %*% w[-1]) )/alpha
}
# the gradient of the objective function
grad.cvar = function(w, Data, alpha, mu, rmin)
{
epsilon2 = epsilon^2
m = length(w)
N = dim(Data)[1]
g = vector(mode = "numeric", length = m)
port = Data %*% w[-1]
# first the derivative wrt to VaR
g[1] = ( sum( ( port - w[1] )/( 2*N * sqrt( epsilon2 + (
port - w[1] )^2 ) ) ) - 0.5 )/alpha + 1
for(i in 2:m){
g[i] = ( sum(Data[,i-1]/(2*N)) - sum(
(Data[,i-1] * ( port - w[1]) )/( ( 2*N ) * (epsilon2 + ( port -
w[1] )^2 )^0.5 ) ) )/alpha
}
# reverse the sign since we are minimizing
return(-1 * g)
}
# the budget constraint:
func.eq = function(w, Data, alpha, mu, rmin){
sum(w[-1]) - 1
}
# gradient of equality
grad.eq = function(w, Data, alpha, mu, rmin){
m = length(w)-1
return( matrix(c(0, rep(1, m)), ncol = m+1, nrow = 1) )
}
# Your inequality (weighted forecast >= rmin)
# but the solver we will use, sets: g(x) <= 0
# so we set: (rmin - weighted forecast) <= 0
func.ineq = function(w, Data, alpha, mu, rmin){
wM = sum( mu * w[-1] )
return( rmin - wM )
}
# gradient inequality
grad.ineq = function(w, Data, alpha, mu, rmin){
m = length(w)
Mu = matrix( c(0, mu), ncol = m )
return(-Mu)
}
## Main Program Example (replicating your data)
library(mvtnorm)
mu <- c(.04, .06, .1)
sig <- c(.04, .09, .2)
cor.mat <- rbind(c( 1, .6, .2),
c(.6, 1, .1),
c(.2, .1, 1)
)
v.cov <- sig%o%sig*cor.mat
set.seed(100)
n <- 10000
samps <- rmvnorm(n, mean = mu, sigma = v.cov)
rmin = .08
alpha = 0.05
wmin = 0
wmax = 1
# I am going to use the nloptr solver (get it from CRAN)
library(nloptr)
.slsqp.ctrl = function(control = list())
{
ctrl = list()
ctrl$algorithm = "NLOPT_LD_SLSQP"
if( is.null(control$ftol_rel) ) ctrl$ftol_rel = 1e-12
else ctrl$ftol_rel = control$ftol_rel
if( is.null(control$xtol_rel) ) ctrl$xtol_rel = 1e-10
else ctrl$xtol_rel = control$xtol_rel
if( is.null(control$maxeval) ) ctrl$maxeval = 5000 else
ctrl$maxeval = as.integer( control$maxeval )
if( is.null(control$maxtime) ) ctrl$maxtime = 10000 else
ctrl$maxtime = control$maxtime
if( is.null(control$print_level) ) ctrl$print_level =0
else ctrl$print_level = as.integer( control$print_level )
if( is.null(control$local_opts$__algorithm) )
ctrl$local_opts$algorithm = "NLOPT_LD_MMA"
if( is.null(control$local_opts$__ftol_rel) )
ctrl$local_opts$ftol_rel = 1e-12 else ctrl$local_opts$ftol_rel
= control$local_opts$ftol_rel
if( is.null(control$local_opts$__xtol_rel) )
ctrl$local_opts$xtol_rel = 1e-10 else ctrl$local_opts$xtol_rel
= control$local_opts$xtol_rel
if( is.null(control$local_opts$__print_level) )
ctrl$local_opts$print_level = 0 else
ctrl$local_opts$print_level = as.integer(
control$local_opts$print_level )
return(ctrl)
}
# set control parameters and subsolvers (use SQP)
ctrl = .slsqp.ctrl(control = list())
# set starting values
q1 = quantile(samps %*% rep(1/3,3), alpha)
# we are estimating [VaR AND the weights]
sol = nloptr(
x0 = c(q1,rep(1/3,3)),
eval_f = func.cvar,
eval_grad_f = grad.cvar,
lb = c(-1, rep( wmin, 3) ),
ub = c( 0, rep( wmax, 3) ),
opts = ctrl,
eval_g_ineq = func.ineq,
eval_jac_g_ineq = grad.ineq,
eval_g_eq = func.eq,
eval_jac_g_eq = grad.eq,
Data = samps, alpha = 0.05, mu = colMeans(samps), rmin = 0.08)
# [VaR w1 w2 w3]
# nloptr[-0.1125099 1.104523e-14 0.480181 0.519819]
# glpk[0.1125283 0.000000 0.480181 0.519819]
# end example
HTH.
Regards,
Alexios
On 15/02/2012 21:05, Robert Harlow wrote:
Hi,
I am running into a memory issue when I try to run a
mean-CVaR
optimization on a large number of scenarios. I am running R
32-bit on
windows (see sessionInfo() output below). I have attached
code to display
the problem below (using a sample from a random
multi-variate normal
distribution as an example, though that is not what I am
actually doing.)
I am using the methodology proposed by Uryasev (thanks to
Guy Yollin for
the cvarOpt function) to formulate the problem as an LP.
When I have a
small number of asset classes, I can arrive at good (though
sometimes non
"optimal" for a given set of scenarios) solutions using
different
optimizers (DEoptim, BBoptim, Rsolnp) where I don't run into
the scenario
size constraint. However, when I increase the number of
asset classes, the
non-LP solvers do a worse and worse job of finding a stable
"global"
optimum, with the exception of DEoptim which will arrive at
a good
solution, but will take a very long time for, say, 25 asset
classes and
10,000 scenarios. The LP, on the other hand, seems much
more sensitive to
number of scenarios than to the number of asset classes.
Thanks in advance,
-Bob
## Code
library(Rglpk)
library(mvtnorm)
mu<- c(.04, .06, .1)
sig<- c(.04, .09, .2)
cor.mat<- rbind(c( 1, .6, .2),
c(.6, 1, .1),
c(.2, .1, 1)
)
v.cov<- sig%o%sig*cor.mat
set.seed(100)
n<- 1000 ##works
n<- 10000##doesn't work due to memory error
samps<- rmvnorm(n, mean = mu, sigma = v.cov)
cvarOpt<- function(rmat, alpha = .05, rmin = .05, wmin = 0,
wmax = 1,
weight.sum =1){
nAss = ncol(rmat) # number of assets
s = nrow(rmat) # number of scenarios i.e. periods
averet = colMeans(rmat)
##create objective vector, constraint matrix,
constraint rhs
Amat =
rbind(cbind(rbind(1,averet),__matrix(data=0,nrow=2,ncol=s+1)__),cbind(rmat,diag(s),1))
objL = c(rep(0,nAss), rep(-1/(alpha*s), s), -1)
bvec = c(weight.sum,rmin,rep(0,s))
##direction vector
dir.vec = c("==",">=",rep(">=",s))
##bounds on weights
bounds = list(lower = list(ind = 1:nAss, val =
rep(wmin,nAss)),upper =
list(ind = 1:nAss, val = rep(wmax,nAss)))
res = Rglpk_solve_LP(obj=objL, mat=Amat, dir=dir.vec,
rhs=bvec,
types=rep("C",length(objL)), max=T, bounds=bounds)
w = as.numeric(res$solution[1:__nAss])
return(list(w=w,status=res$__status))
}
ans.cvar<- cvarOpt(rmat = samps, rmin = .08)
## session info
R version 2.13.2 (2011-09-30)
Platform: i386-pc-mingw32/i386 (32-bit)
locale:
[1] LC_COLLATE=English_United States.1252
[2] LC_CTYPE=English_United States.1252
[3] LC_MONETARY=English_United States.1252
[4] LC_NUMERIC=C
[5] LC_TIME=English_United States.1252
attached base packages:
[1] stats graphics grDevices utils datasets
methods base
other attached packages:
[1] mvtnorm_0.9-9991 Rglpk_0.3-5 slam_0.1-22
loaded via a namespace (and not attached):
[1] tools_2.13.2
[[alternative HTML version deleted]]
_________________________________________________
R-SIG-Finance at r-project.org
<mailto:R-SIG-Finance at r-project.org> mailing list
https://stat.ethz.ch/mailman/__listinfo/r-sig-finance
<https://stat.ethz.ch/mailman/listinfo/r-sig-finance>
-- Subscriber-posting only. If you want to post, subscribe
first.
-- Also note that this is not the r-help list where general
R questions should go.