Efficient CVaR Scenario Optimization for a Large number of Scenarios

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

Efficient CVaR Scenario Optimization for a Large number of Scenarios

Thread (6 messages)