Efficient CVaR Scenario Optimization for a Large number of Scenarios

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]]