Scenario generation for single-period portfolio selection problems with tail risk measures: coping with high dimensions and integer variables

In this paper we propose a problem-driven scenario generation approach to the single-period portfolio selection problem which use tail risk measures such as conditional value-at-risk. Tail risk measures are useful for quantifying potential losses in worst cases. However, for scenario-based problems these are problematic: because the value of a tail risk measure only depends on a small subset of the support of the distribution of asset returns, traditional scenario based methods, which spread scenarios evenly across the whole support of the distribution, yield very unstable solutions unless we use a very large number of scenarios. The proposed approach works by prioritizing the construction of scenarios in the areas of a probability distribution which correspond to the tail losses of feasible portfolios. The proposed approach can be applied to difficult instances of the portfolio selection problem characterized by high-dimensions, non-elliptical distributions of asset returns, and the presence of integer variables. It is also observed that the methodology works better as the feasible set of portfolios becomes more constrained. Based on this fact, a heuristic algorithm based on the sample average approximation method is proposed. This algorithm works by adding artificial constraints to the problem which are gradually tightened, allowing one to telescope onto high quality solutions.


Introduction
In the portfolio selection problem one must decide how to invest in a collection of financial instruments with uncertain returns which in some way balances portfolio return against risk. There are many ways of modeling this problem. In the robust optimization setting, the returns of the portfolio are assumed to fall within some uncertainty set, and one minimizes the worst-case of loss [BT06]. This approach is sometimes considered too conservative as it does not make effective use of available information. In stochastic programming, the user uses their knowledge and available data to explicitly model asset returns as random vectors, and then optimize some combination of expected return and risk measure. In between these two paradigms is distributionally-robust optimization [HZFF10] where the returns are modeled by a random vector whose distribution lies in some uncertainty set, and the worst-case expected loss is minimized. This approach is particularly useful when only limited or unreliable data is available to model the uncertainty, but can lead to intractable problems. Besides the optimization paradigm employed, portfolio selection problems can be further categorized into single-period problems, where only one portfolio selection is made, and multiperiod problems where the portfolio may be rebalanced several times.
The work in this paper applies to the stochastic programming single-period formulation of the portfolio selection problem. This approach is popular as it allows one to flexibly model the return of the distributions, can easily incorporate problem details such as transactions costs [LFB07], while remaining generally more tractable than other types of models. We deal in particular with a difficult variety of this problem type which involve tail risk measures and potentially integer variables.
In the typical set-up the uncertain returns are modeled by random variables, the total return of a portfolio is some linear combination of these, and riskiness is measured by a real-valued function of the total return which should in some way penalize potential large losses. This approach was first proposed by Markowitz [Mar52] who used variance as a risk measure.
Although the use of variance has remained popular [LFB07] because it leads to tractable convex programs, its use as a measure of risk is problematic for a few reasons. The foremost of these is perhaps that variance penalizes large profits as well as large losses. As a consequence, in the case where the returns of financial assets are not normally distributed, using the variance can lead to potentially bad decisions; for instance, a portfolio can be chosen in favor of one which always has higher returns (see [You98] for an example of this). This particular issue can be overcome by using a "downside" risk measure, that is one which only depends on losses greater than the mean, or some other specified threshold, for example the semi-variance [Mar59, Chapter 9], mean regret [DR99], or value-at-risk [Jor96]. More recently, much research has been given to coherent risk measures, a concept introduced in [ADEH99]. These are risk measures which have sensible properties such as subadditivity, which in particular ensures that a risk measure incentivizes diversification of a portfolio. Using a coherent risk measure in a portfolio selection problem should avoid flawed decisions, such as the one cited in the case of variance.
In this work, we are interested in portfolio selection problems involving tail risk measures. These can be thought of as risk measures which only depend on the upper tail of a distribution above some specified quantile. A canonical example of a tail risk measure is the value-at-risk (VaR) [Jor96]. The β-VaR is defined to be the β-quantile of a random variable. In portfolio selection problems this has the appealing interpretation as the amount of capital required to cover up to β × 100% of potential losses. Thus, tail risk measures, in particular those which dominate the β -VaR, are useful as they can give us some idea of the amount capital at risk in the worst (1 − β) × 100% of potential losses. Like variance, the value-at-risk is also problematic as it is not a coherent measure of risk. Specifically, it is not subadditive (see [Tas02] for example). Moreover, β -VaR leads to difficult and intractable problems when used in an optimization context. The conditional value-at-risk (CVaR), sometimes referred to as the expected shortfall, is another tail risk measure and can be roughly thought as the conditional expectation of a random variable above the β-VaR. It is both coherent [AT02], and more tractable in an optimization setting [RU00].
However, the use of risk measures, even coherent ones such as β -CVaR, is still problematic in portfolio selection problems where the asset returns are modeled with continuous probability distributions. This is because the evaluation of many risk measures for arbitrary continuously distributed returns would involve the evaluation of a multidimensional integral, and this becomes computationally infeasible when our problems involve many assets. On the other hand, the evaluation of such an integral reduces to a summation if the returns have a discrete distribution.
Scenario generation is the construction of a finite discrete distribution to be used in a stochastic optimization problem. This may involve fitting a parametric model to asset returns and then discretizing this distribution, or directly modeling them with a discrete distribution, for example via moment-matching [HKW03]. In either case, standard scenario generation methods struggle to adequately represent the uncertainty in problems using tail risk measures. This is because the value of a tail risk measure, by definition, only depends on a small subset of the support of a random variable, and typical scenario generation methods will spread their scenarios evenly across the whole support of the distribution. This means that the region on which the value of the tail risk depends, is represented by relatively few scenarios. Hence, unless there is a very large number of scenarios, the value of of tail risk measure is very unstable (see [KWVZ07] for example).
The natural remedy to this problem is to represent the regions of the distribution on which the tail risk measure depends with more scenarios. Intuition would tell us that these correspond to the "tails" of the distribution. However, for a multivariate distribution there is no canonical definition of the tails. If by tails, we simply mean the region where at least one of the components exceeds a large value, then the probability of this region quickly converges to one with the problem dimension, and thus prioritizing scenarios in this region will be of little benefit. Finding the relevant tails of the distribution is a non-trivial problem.
In our previous paper [FTW17] we addressed the problem of scenario generation for stochastic programs with an arbitrary loss function functions which use tail risk measures, and for this we defined the concept of a β-risk region. In portfolio selection, to each valid portfolio there is a distribution of losses (or returns). The β-risk region consists of all potential asset returns which lead to a loss in the β-tail for some portfolio. We have shown that under mild conditions the value of a tail risk measure in effect only depends on the distribution of returns in the risk region. Although characterizing this region in a convenient way is generally not possible, we have been able to do this for the portfolio selection problem when the asset returns are elliptically distributed. We have also proposed a sampling approach to scenario generation using these risk regions which prioritizes the generation of scenarios in the risk region. We demonstrated for simple examples that this methodology can produce scenario sets which yield better and more stable solutions than does basic sampling.
In this paper we address issues related to the application of this methodology to realistic portfolio selection problems. The first major contribution of this work is the application of the methodology to problems where the asset returns have non-elliptical distributions. In this case the distribution of returns for a portfolio will in general not have a convenient closed form, and so it is necessary to represent the asset returns with a scenario set. In order to apply the methodology, we approximate the risk region for a non-elliptical distribution with the risk region of an elliptical distribution. Moreover, we demonstrate that the methodology is effective on difficult problems which have high dimensions and integer variables.
In the paper [FTW17] it was shown that our methodology is more effective as the problem becomes more constrained. The second major contribution of this work is the proposal of an heuristic algorithm based on the stochastic average approximation (SAA) method [KSHdM01] which exploits this fact. This algorithm works by adding artificial constraints to the problem which are gradually tightened, allowing one to telescope onto high quality solutions. This algorithm is presented in a general way and could be potentially used on problems other than portfolio selection.
The paper is organized as follows. In Section 2 we formally define risk regions for general stochastic programs, and present some new results on the use of approximate risk regions. In Section 3 we define the portfolio selection problem and recall some results on the use of the risk region methodology for this problem. We also provide some new technical results related exploiting risk regions of elliptical distributions. In Section 4 we describe how risk regions are exploited for the purpose of scenario generation, and present a heuristic based on the SAA method based which uses of artificial constraints. In Section 5 we make some empirical observations on how the probability of the non-risk region, a quantity which determines the effectiveness of the methodology, varies with the type of distribution of asset returns. In Section 6 we present results for a broad range of numerical tests on the effectiveness of our sampling and reduction algorithm using distributions constructed from real asset return data. In Section 7 we demonstrate the performance of the proposed heuristic on a difficult case study problem. Finally, in Section 8 we make some concluding remarks.

Risk regions for general stochastic programs
In this section we formally define the concept of risk region and present results related to these. This theory does not just apply to portfolio selection problems but more generally to stochastic programs with a tail risk measure. In Section 2.1 we recall the basic definitions and fundamental results for risk regions which appeared in our previous paper [FTW17]. In Section 2.2 we present some new results related to the approximation of risk regions.

Tail risk measures and risk regions
A risk measure is a function of a real-valued random variable representing a loss. For 0 < β ≤ 1, a β-tail risk measure can be thought of as a function of a random variable which depends only on the upper (1 − β)-tail of the distribution. The precise definition uses the generalized inverse distribution function or (lower) quantile function.
Definition 2.1 (Quantile function and β-tail risk measure). Suppose Z is a random variable with distribution function F Z . Then the generalized inverse distribution function, or quantile function is defined as follows: Now a β-tail risk measure is any function of a random variable, ρ β (Z), which depends only on the quantile function of a random variable above β.
Example 2.2 (Value at risk (VaR)). Let Z be a random variable, and 0 < β < 1. Then, the β−VaR for Z is defined to be the β-quantile of Z: . Example 2.3 (Conditional value at risk (CVaR)). Let Z be a random variable, and 0 < β < 1. Then, the β -CVaR can be thought roughly as the conditional expectation of a random variable above its βquantile. The following alternative characterization of β -CVaR [AT02] shows directly that it is a β-tail risk measure.
The observation that we exploit for this work is that very different random variables will have the same β-tail risk measure as long as their β-tails are the same.
In the optimization context we suppose that the loss depends on some decision x ∈ X ⊆ R k and the outcome of some latent random vector Y with support Y ⊆ R d , defined on a probability space (Ω, F, P), and which is independent of x. That is, we suppose our loss is determined by some function, f : X × R d → R, which we refer to as the loss function. For a given decision x ∈ X , the random variable associated with the loss is thus f (x, Y ).
To avoid repeated use of cumbersome notation we introduce the following short-hand for distribution and quantile functions: In [FTW17] we introduced the concept of a risk region for a stochastic program using a tail-risk measure. We define this now for a general stochastic program.
Definition 2.4 (Risk region). The β-risk region associated associated with the random vector Y and the feasible region X ⊆ R d is as follows: The risk region consists precisely of those outcomes of Y which have a loss in the β-tail of the loss distribution for some feasible decision. We refer to the complement of the risk region as the non-risk region and this consists of outcomes which never lead to a loss in the β-tail; it can be written as follows: The following theorem was proved in [FTW17] and states that under mild conditions the value of a tail risk measure is completely determined by the the distribution of the random vector Y in the risk region. That is, the values of the tail risk measure of any two random vectors with identical distributions in the risk region will be the same for all feasible decisions. The technical condition in (3), which we call the aggregation condition for reasons explained below, precludes certain degenerate cases. In essence, this condition ensures that there is enough mass in the set to ensure that the β-quantile does not depend on the probability distribution outside of it.
Theorem 2.5. Let R ⊇ R Y,X (β) be such that for all x ∈ X the following condition holds: IfỸ is a random vector for which the following holds: then ρ β (f (x, Y )) = ρ β f (x,Ỹ ) for all x ∈ X , for any β-tail risk measure ρ β .
With regards to scenario generation, this theorem says that any scenarios in the non-risk region can be aggregated into a single point, reducing the size of the problem, without affecting the value of the tail risk measure. This motivates the term aggregation condition for (3). The transformed random vector where all mass in a region has been concentrated into its conditional expectation plays a special role in this work. We call this the aggregated random vector.
Definition 2.6 (Aggregated Random Vector). For some set R ⊇ R Y,X the aggregated random vector is defined as follows: The conditional expectation E [ Y |Y ∈ R c ] is guaranteed to fall in the non-risk region if, for example, the loss function is convex [FTW17, Proposition 3].

Approximation of Risk Regions
The methodology proposed in this paper requires a characterization of a risk region which allows one to easily test membership. A convenient characterization for the exact risk region as defined in (1) is in general not possible as this set is determined by the loss function, distribution of the random vector and the problem constraints. Even for the portfolio selection which has a simple loss function, we cannot find a convenient form for arbitrary distributions of asset returns.
Recall that Theorem 2.5 applies to any set containing the risk region. Therefore, one way to circumvent the problem of finding the exact risk region would be to use a conservative risk region, that is, a set which contains the exact risk region. This approach is particularly useful for problems which have constraints which cannot be easily taken into account, such as constraints involving integer variables in the portfolio selection problem. By Definition 2.4, if X ⊆ X then R Y,X ⊆ R Y,X . Therefore, ignoring some constraints will yield a risk region which is conservative.
In the case where one cannot construct a risk region for a given loss function or distribution, it may be difficult to find a conservative risk region. Moreover, it could be the case that a conservative risk region may be too conservative to be of any use. Instead one might try to use an approximate risk region. For the portfolio selection problem we handle distributions for which we cannot conveniently characterize the exact risk region by using the risk region of a surrogate distribution which is similar to the true distribution.
Denote by R ⊂ R d an approximate risk region. When using an approximate risk region the value of the tail risk measure may be distorted for a decision x ∈ X if R does not contain all outcomes in the β-tail, that is, unless the following condition holds: If (5) holds then we say that that the approximate risk region R is valid for decision x.
We show in this section that if the approximate risk region is not valid for a particular decision, then, under mild assumptions, the values of β -VaR and β -CVaR tail risk measures are distorted downwards. In Section 3, we will exploit this observation to show that for the problems in which we are interested, if there is no distortion of the value of the tail risk measure at the optimal solution, then this solution is also optimal with respect to the true problem.
For the results in this section we employ the following notation:F x andF −1 x denote respectively the distribution and quantile functions of ψ R (Y ). We require the following conditions: Assumption (B) requires that conditional expectation of Y in the complement of approximate risk region belongs to the exact non-risk region. This means that the loss at the aggregated point will have a loss below the β-quantile of f (x, Y ) for all feasible decisions x ∈ X . Before stating and proving the key result, we require the following lemma.
Proof. Fix x ∈ X and z > f (x, E [Y |Y ∈ R]), and without loss of generality assume that The key result states that β-quantile (or β -VaR) and β -CVaR for the aggregated random vector cannot increase when using an approximate risk region under the above assumptions. The implications of this result on the portfolio selection problem are made clear in Section 3.1.
Proposition 2.8. Under assumptions (A) and (B), we have with equality if R is valid for x ∈ X (in the sense of (5)) and the aggregation condition holds. Proof.
. For the β -CVaR, recall that for a random variable Z this can be written as follows [AT02]: where 1 A denotes the indicator function of event A. Since F x is continuous we can write: On the other hand,F x could have a discontinuity atF −1 . We therefore consider two cases: In the first case,F x in continuous atF −1 x (β) by Lemma 2.7 so we can write: Note that the integrand of the first term is greater than that of the second term over the respective domain of integration. Therefore, to show that the above quantity is non-negative, it is enough to show that the domain of integration of the first term has the same probability as the second. To show this, first not that: Therefore, rearranging which gives as required.
In the second case,F x has a discontinuity atF −1 x (β), and so the β -CVaR is written as follows: Noting that, Finally, manipulation of the probabilities above yields: , as required. The fact that the inequalities hold with equality if R is valid for decision x and the aggregation condition holds follows directly from Theorem 2.5 for the special case X = {x}.

Risk Regions for Portfolio Selection
In this section we present results relating to risk regions for the portfolio selection problem. In Section 3.2 we define the problem and present general results from [FTW17] related to the risk region for this problem. The remaining two subsections deal with risk regions for elliptical distributions since these are used as approximate risk regions. Specifically, in Section 3.2 we formally define elliptical distributions and give a convenient characterization of their corresponding risk regions, and in Section 3.3 we present some new results related to testing membership to a risk region.

Problem statement and application of risk regions
We use the following basic set-up: we have a set of financial assets indexed by i = 1, . . . , d, by x i we denote how much we invest in asset i, and by Y i we denote the random future return of asset i. The portfolio return associated to a particular investment decision The loss function associated to an investment decision is thus f (x, Y ) = −x T Y , and so for a given β-tail risk measure ρ β we would like an investment with small tail risk ρ β (−x T Y ). The aim of a portfolio selection problem is to choose a decision which balances choosing a portfolio with high expected portfolio return against choosing one with small risk. This typically corresponds to solving a problem of one of the following forms: where 0 ≤ λ ≤ 1 and X ⊂ R d represents the set of feasible portfolios. This feasibility region will typically encompass a constraint which specifies the amount of capital to be invested, and may include others which, for example the exclusion of short-selling, or a limit on the amount that can be invested in certain industries.
In the case of the portfolio selection problem, the risk region is R Y,X (β) : x (β)}, that is, it is the union over all feasible portfolios, of the half spaces of points with returns above the β-quantile. We can find this region by brute force, and this is illustrated for a hypothetical discrete random vector on the left-hand side of Figure 1. Also illustrated in this figure is the set of returns where all the mass in the non-risk region has been aggregated into the conditional expectation of the random vector in the non-risk region, that is, the aggregated random vector. The figure also demonstrates, as implied by Theorem 2.5, that the β-quantile lines do not change after aggregation. Assuming E [Y |Y ∈ R c ] ∈ R c Y,X , as well as preserving the value of a tail risk measure, the aggregated random vector has the additional property of preserving the overall expected return of the original random vector. The following corollary taken from [FTW17] summarizes this result and provides sufficient conditions so that (3) holds.
and X contains at least two linearly independent elements. Then Y satisfies (3). In addition, if R c is convex thenỸ = ψ R (Y ) satisfies condition (4) and so for all x ∈ X we have: In Section 2.2 we showed that under mild conditions when using an approximate risk region, a misspecification for a particular decision will decrease the value of the β -VaR and β -CVaR. Building on this result, the following corollary gives a condition under which the optimal solution yielded by using an approximate risk region is also optimal for the true problem.
Corollary 3.2. Under assumptions (A) and (B), for the problems (P1), (P2) and (P3), suppose that replacing the solution the random vector Y with the aggregated random vector ψ R (Y ) for an approximate risk region R yields an optimal solutionx. Then, if R is valid forx and the aggregation condition holds for R thenx is also an optimal solution for the true problem.
Proof. We prove this only for P1. The proofs for the other problems are very similar. First note that x is an optimal solution to the true problem then, Hence,x is optimal with respect to the true problem.
This result guarantees a certain robustness to misspecification of the approximate risk region. Although checking whether an approximate risk region is valid for a decision could in principle be used as an optimality check, we will not use it in this way as checking directly condition (5) may be difficult. As will be seen in Section 7 we instead will rely on out-of-sample testing to verify the quality of a solution.

Risk regions for elliptical distributions
In order to exploit risk regions for scenario generation one has to be able to characterize these in a way which allows one to conveniently test whether or not a point belongs to it. In our previous paper, we were able to do this in the case where the asset returns have elliptical distributions. Elliptical distributions are a general class of distributions which include, among others, multivariate Normal and multivariate t-distributions. See [FKN89] for a full overview of the subject.
where ∼ means the two operands have the same distribution function. Let Y be a random vector in R d , then Y is said to be elliptical if it can be written Y = P T X + µ where P ∈ R d×d is non-singular, µ ∈ R d , and X is random vector with spherical distribution. Such an elliptical distribution will be denoted Elliptical (X, µ, P ).
This definition says that a random vector with a spherical distribution is rotation invariant, and that an elliptical distribution is an affine transformation of a spherical distribution. Elliptical distributions are convenient in the context of portfolio selection as we can write down exactly the distribution of loss of a portfolio. In particular, if Y ∼ Elliptical (X, µ, P ) and where · denotes the standard Euclidean norm and X 1 is the first component of the spherical random vector X. Therefore, the β-quantile of the loss −x T Y is as follows: For Y ∼ Elliptical (X, µ, P ), we can thus rewrite the risk region in (1) as follows: In this form it is still difficult to check whether a given pointỹ ∈ R d belongs to it. In [FTW17] we provided a more convenient characterization of the risk region for elliptical returns. This characterization makes use of the conic hull of the feasible region X ⊂ R d .
Definition 3.4 (Convex cones and conic hull). A set K ⊂ R d is a cone if for all x ∈ K and λ ≥ 0 we have λx ∈ K. A cone is convex if for all x 1 , x 2 ∈ K and λ 1 , λ 2 ≥ 0 we have λ 1 x 1 + λ 2 x 2 ∈ K. The conic hull of a set A ⊂ R d is the smallest convex cone containing A, and is denoted conic (A).
For example, suppose that our feasible region consists of portfolios with non-negative investments (i.e. no short-selling) and whose total investment is normalized to one, that is: then the conic hull of this is the positive quadrant, that is conic (X ) = R d + . The alternative characterization also makes use of projections.
Definition 3.5 (Projection). Let C ⊂ R d be a closed convex set, then for any point y ∈ R d we define its projection onto C to be the unique point We are now ready to give a characterization of the risk region. For this we use the following convenient abuse of notation: for a set A ⊂ R d and a matrix T ∈ R d×d , we write T (A) := {T y : y ∈ A}. The following result was proved in [FTW17].
Theorem 3.6. Suppose Y ∼ Elliptical (X, P, µ), X ⊆ R d is convex and let K = conic (X ). Then the risk region can be characterized exactly as follows: where K = P K is a linear transformation of the conic hull K.
When µ = 0 and P = I d the risk region can be simplified to This allows us to interpret the projection of y onto K as the portfolio which leads to the largest loss.

Testing membership to a risk region
Scenario generation algorithms which exploit risk regions rely on the ability to test membership of the risk region for randomly sampled points. The characterization of the risk region for portfolio selection problems given in (7) relies on one being able to calculate the conic hull of the set of feasible portfolios, and also the ability to project points onto a transformation of this. In Section 3.3.1 we show how one can find the conic hull of the feasible region for typical constraints of a portfolio selection problem. This conic hull is a finitely generated cone. In Section 3.3.2 we show how one can project points onto this type of cone. Finally in Section 3.3.3 we briefly discuss the computational issues for the membership tests.

Conic hull of feasible region
In portfolio problems, the feasible region is usually defined by linear constraints, that is That is, the feasible region is the intersection of a finite number of half-spaces. It is a well-known fact that any such intersection can be written as the convex hull of a finite number of points plus the conical combination of some more points (see Theorem 1.2 in [Zie08] for example). That is, there exists x 1 , . . . , x k ∈ R d and y 1 , . . . , y l ∈ R d such that The conic hull of this region is the following finitely generated cone: To express the intersection of half-spaces in the form (8), we could use Chernikova's algorithm (also known as the double description method) [Che65,LV92]. Every finitely generated cone can also be written as a polyhedral cone, that is, of the form {x ∈ R d : Dx ≥ 0}, and vice versa (see [Zie08, Chapter 1]). Chernikova's algorithm again provides a concrete method for going between these two different representations. Although these two representations are mathematically equivalent, as we shall see, they are algorithmically different. We will suppose the constraints for our portfolio selection problem have the following form: where 1 is column vector of ones and c > 0. The first of these constraints specifies the total of amount of capital to be invested, the inequalities represent other constraints such as quotas on the amount one can invest in a specific company or industry. In this case, we can describe immediately the conic hull as a polyhedral cone.
Proposition 3.7. Let X be the set defined in (9) and let Proof. Given that X is convex, to show that conic (X ) = Y, it suffices to show that We demonstrate first the forward implication. Suppose x ∈ Y \ {0}. Then, given that x > 0, we must have v := 1 T x > 0. Then, setting λ = c v , we have Since Y is a cone, we have λx ∈ Y, hence and so λx ∈ conic (X ). We now prove the backwards implication. Suppose x ∈ conic (X ) \ {0}. Then there exists λ > 0 such that λx ∈ X , that is Hence x ∈ Y as required.
x 2 x 1 Figure 2: Conic hull from simple quota constraints given x 1 + x 2 = 1 and x 1 , x 2 ≥ 0 Figure 2 shows how simple constraints in R 2 affect the conic hull of the feasible region given the total investment and positivity constraints.

Projection onto a finitely generated cone
First, suppose that we can represent the conic hull of the feasible region X ⊂ R d as a finitely generated cone with k generators, that is K = {Ay : y ≥ 0} where A ∈ R k×d . By definition, the projection of a point x 0 ∈ R d can be found by solving the following quadratic program: In particular, if y * is the optimal solution then p K (x 0 ) = Ay * . By formulating the KKT conditions [BV04, Chapter 5] of this problem, it can be seen that this problem is equivalent to solving the following linear complementarity problem (LCP): Find y, z ∈ R d such that If (y, z) is a solution to the above problem, then the required projection is p K (x 0 ) = Ay. LCPs can be solved by more specialized algorithms than standard quadratic programs such as Lemke's algorithm [CPS92]. Now, suppose instead we have a polyhedral characterization of the conic hull, that is a cone of the form: The projection of a point x 0 ∈ R d onto the polyhedral cone in (11) is the solution of the following quadratic program: Although the former problem in (10) can often be solved more efficiently using specialized algorithms, we will in practice use both approaches. For conic hulls with a small number of extremal rays, for example K = R d + we will use use the former method. As we add more constraints to the problem, we have found from experience that the number of extremal rays can exponentially increase, which for the former approach leads to cumbersomely large LCP problems. In this case we will use the polyhedral representation for projection.

Computational issues
Given that testing membership of the risk region (for elliptically distributed returns) involves solving a small LCP or quadratic program, using this methodology could potentially become computationally expensive, especially if used to construct large scenario sets for high dimensional problems. However, this issue can be mitigated in a few ways. Firstly, the membership test for a point can be conducted independently from that of another point, which means that membership tests for a large number of points is naturally parallelizable. Secondly, for the case where the case K ⊆ R d + the loss function y → −x T y is monotonic. Therefore, if we have set of points y 1 , . . . , y k which are in the risk region, we know that a pointỹ is also in the risk region if it is dominated by any of these points. A similar shortcut exists for testing ifỹ is in the non-risk region. Finally, the membership test could also be made more efficient, by directly testing the condition p K (y) ≤ α without calculating the full projection p K (y). For example, the quadratic program used to calculate the projection could be solved only to an accuracy sufficient to test this condition. This could be easily implemented through a callback function in the quadratic program solver.

Scenario generation
In this Section we show how risk regions can be exploited for the purposes of scenario generation. In Section 4.1 we present two specific methods which work essentially by prioritizing the construction of scenarios in the risk region. In Section 4.2 we propose a new heuristic algorithm based on the SAA method [KSHdM01]. This heuristic boosts the performance of the proposed sampling algorithm through the addition of artificial constraints to the problem.

Aggregation sampling and reduction
In [FTW17] we proposed two methods to exploit risk regions. The first of these allows the user to specify the final number of scenarios in advance. The algorithm, which is called aggregation sampling, samples scenarios, aggregating all samples in the non-risk region and keeping all in the risk region, until we have the required number of risk scenarios, that is the required number of scenarios in the risk region. This is described in Algorithm 1.
Let q = P Y ∈ R c Y,X be the probability of the non-risk region, and n be the number of risk scenarios required. Define N (n) to be the effective sample size from aggregation sampling, that is, the number of draws until the algorithm terminates 1 . The quantity N (n) is a random variable: where N B(N, q) denotes a negative binomial random variable. Recall that a negative binomial random variable N B(n, q) is the number of failures in a sequence of Bernoulli trials with probability of success q until n successes have occurred. The expected effective sample size of aggregation sampling is thus as follows: The expected effective sample size can be thought of as the sample size required for basic sampling to produce the same number of scenarios in the risk region. The difference between the desired number of risk scenarios, and expected effective sample size is proportional to the ratio q 1−q . In particular, as the probability of the non-risk region approaches one, this gain tends to infinity.
The converse to aggregation sampling is sampling a set of a given size n and then aggregating all scenarios in the risk region of the underlying distribution. We call this aggregation reduction. This can be viewed as a sequence of n Bernoulli trials, where success and failure are defined in the same way as input : R ⊂ R d approximate risk region, N R number of required risk scenarios output: Algorithm 1: Aggregation sampling described above. The number of scenarios in the reduced sample, R(n) is as follows: where B(n, q) denotes a binomial random variable. The expected reduction in scenarios in aggregation reduction is thus nq − 1.
The reason why aggregation sampling and aggregation reduction work is that, for large samples, they are equivalent to sampling from the aggregated random vector. Suppose that Y 1 , Y 2 , . . . is a sequence of independently identically distributed (i.i.d.) random vectors with the same distribution as Y , then ψ R (Y 1 ), ψ R (Y 2 ), . . . is a sequence of i.i.d. random vectors with the same distribution as the aggregated random vector ψ R (Y ). Denote byρ n,β (x) the value of the tail-risk measure for the decision x ∈ X for the sample ψ R (Y 1 ), . . . , ψ R (Y n ), and byρ n,β the analogous function the scenario set constructed by aggregation sampling. The following result, adapted from [FTW17] for the portfolio selection problem, gives precise conditions under which aggregation sampling is asymptotically valid.
Theorem 4.1. Suppose the following conditions hold: (i) For each x ∈ X , F x is strictly increasing and continuous in some neighborhood of Then, with probability 1, for n large enoughρ n,β ≡ρ N (n),β .
See [FTW17, Section 4] for a full proof of the consistency of these algorithms.

Ghost constraints
We noted above that the performance of our methodology improves as the probability of the non-risk region increases. In particular, the expected effective sample size in aggregation sampling increases as the probability of the non-risk risk region increases. Now, by its definition (2) the non-risk region grows as the problem becomes more constrained. This suggests that it may be helpful to add constraints to our problem which shrink the set of feasible portfolios, but which are not themselves active, in the sense that their presence does not affect the set of optimal solutions. We will refer to a constraint added to a problem to boost the performance of our methodology, loosely, as a ghost constraint.
Finding non-active constraints to add to our problem is non-trivial as it relies on some knowledge of the optimal solution set. Moreover, even verifying whether or not a particular constraint is active is difficult in general for stochastic programs. For a deterministic objective function which is convex and for which all constraints are convex (and the optimal solution is unique) a constraint {x : g(x) ≤ 0} is active if and only if it is binding at the optimal solution x * , that is g(x * ) = 0. For a stochastic program, we are typically solving a scenario-based approximation and so a constraint which is not binding with respect to the scenario-based approximation may be binding with respect to the true problem and vice versa. A rigorous test of whether a ghost constraint is active in the sense above is beyond the scope of this paper. We simply promote the idea here that ghost constraints may be a useful way of finding better solutions.
We resort to heuristic rules to choose ghost constraints. For example, one could constrain our set of feasible portfolios to some neighborhood of a good quality solution. This suggests an iterative procedure whereby one samples scenario sets using aggregation sampling, solves the resulting problems, adjusts the problem constraints and then resamples. In Algorithm 2 we propose such a heuristic procedure based on the sample average approximation (SAA) method of [KSHdM01]. We call this procedure the SAA method with ghost constraints. Like with the original SAA method, this algorithm is presented in a very general form, as the update rules, such as how the bounds are adjusted in each iteration, can be implemented in many different ways. In Section 7 we test this algorithm on a realistic and difficult problem.
Add constraints l ≤ x ≤ u to problem; Construct risk region R for problem; for i = 1, . . . , M do Generate a sample of size N using aggregation sampling with risk region R; Solve corresponding problem with objective value ν m N and optimal solutionx m N ; Estimate optimality gap of solution and the variance of the optimality gap; end Increase size of N ; Adjust bounds l and u; until Optimality gap and variance are sufficiently small for some m; Choose the best solutionx among all candidate solutions x m N produced, using a screening and selection procedure.
Algorithm 2: Sample average approximation method with ghost constraints

Probability of the non-risk region
The benefit of aggregation sampling and reduction depends on the probability of the non-risk region.
As was observed in [FTW17] the probability of the non-risk region tends to decrease as the problem dimension increases, but increases as we tighten our problem constraints, and as we increase β, the level of the tail risk measure. In this section we make some empirical observations on how this probability varies with heaviness of the tails, and the correlations of the distribution. The first observation is that in the presence of positivity constraints, the probability of the non-risk region increases as the the correlation between random variables increases. This can be seen in Figure  3 which plots the probability of the risk region as a function of correlation for some two-dimensional distributions. An intuitive explanation for this type of behavior is that in the case of positive correlations there is much more overlap in the risk regions of the individual portfolios.
The extent to which probabilities vary with correlation seems to be much greater in higher dimensions. In Figure 4 we have plotted for Normal returns and a range of dimensions, the probabilities of the nonrisk region for a particular type of correlation matrix: Λ (ρ) ∈ R d×d where Λ(ρ) ij = ρ for i = j and ρ > 0. In the case of ρ = 0, the probability decays very quickly to zero as the dimension increases, whereas as when ρ is close to one, the probability of the non-risk region approaches β for all dimensions.
Our next observation is that the probability of the non-risk region seems to increase as the tails of the distributions become heavier. In Figure 5 are plotted the probabilities of risk regions for some spherical distributions and a range of dimensions. Note that multivariate t-distributions have heavier tails than the multivariate Normal distribution, but the tails get lighter as the degrees of freedom parameter increases. This phenomenon can also be observed in Figure 3.
The observations made in this section suggest that that the application of our methodology will be particularly effective when applied to real stock data tend to be positively correlated and have heavy tails.

Numerical tests
In this Section we test the numerical performance of our methodology for realistic distributions. There are three parts to these tests: the calculation of the probability of the non-risk region for a range of distributions and constraints, the performance of aggregation sampling, and the performance of aggregation reduction. To allow us to measure the quality of the solutions yielded by the proposed methods, the majority of the tests are performed for elliptical distributions on problems without integer variables. These problems can be solved exactly using other methods. In Section 6.1 we describe our experimental set-up, in particular we justify the distributions constructed for these experiments. The remaining three sections detail the individual experiments and their results.

Experimental set-up
For robustness we will use several randomly constructed distributions for each family of distribution and each dimension we are testing. We construct these by fitting parametric distributions to real data. We use real data rather than arbitrarily generated problem parameters for two reasons. Firstly, generating parameters which correspond to well-defined distributions can be problematic. For example, for the moment matching algorithm, there may not exist a distribution which has a given set of target moments (see [KME00] and [JR51] for instance). Secondly, as was observed in Section 5, the probability of the non-risk region can vary widely, and so it is most meaningful to test the performance of our methodology for distributions which are realistic for portfolio selection problems.
We construct our distributions from monthly return data from between January 2007 and February 2015 for 90 companies in the FTSE 100 index. For each dimension in the test, we randomly sample five sets of companies of that size, and for each of these sets fit Normal, t distributions and Skew-t distributions to the associated return data. Figure 6 shows for two stocks the contours of the fitted density functions overlaying the historical return data. For the t distributions we fix the degrees of freedom parameter to 4.0. This is so that we can more easily compare the effect of heavier tails on the results of our tests. We allow the corresponding parameter for the skew-t distributions to be fitted from the data.
These three distributions are fitted to the data through maximum likelihood estimation, weighing each observation equally; our aim here is not to construct distributions which accurately capture the uncertainty of future returns, but to simply construct distributions which are realistic for this type of problem. We also use scenario sets constructed using the moment-matching algorithm. For each random set of companies, we calculate all the required marginal moments and correlations from their historical returns, and use these as input to the moment-matching algorithm. To allow us to compare results, the same constructed distributions are used across the three sets of numerical tests.
Throughout this section we use the β -CVaR as our tail risk measure. This is not only because the β -CVaR leads to tractable scenario-based optimization problems, but also for elliptically distributed returns we can evaluate the β -CVaR exactly which provides us with a means to evaluate the true performance of the solutions yielded by the approximate scenario-based problems. In addition, to ensure that the non-risk region does not have negligible probability, we will assume that we always have positivity constraints on our investments (i.e. no short selling).
The first class of non-elliptical distributions we use in this paper are known as multivariate Skew-t Skew t Figure 6: Contour plots of distributions fitted to financial return data for two assets distributions [AC03]. This class of distributions generalizes the elliptical multivariate t-distributions through the inclusion of an extra set of parameters which regulate the skewness. In this case we approximate the risk region with the risk region of the corresponding t-distribution.
The second class we use are discrete distributions constructed using the moment matching algorithm of [HKW03]. These distributions have been applied previously to financial problems [KWVZ07]. This algorithm constructs scenario sets with a specified correlation matrix and whose marginals have specified first four moments. This algorithm works by first taking a sample from a multivariate Normal distribution, and then iteratively applying transformations to this until the difference between its marginal moments and correlation matrix are sufficiently close to their target values. Since the algorithm is initialized with a sample from an elliptical distribution, the final distribution is near elliptical and we approximate the risk region for these distributions with the risk region of a multivariate normal distribution with the same mean and covariance structure.

Probability of non-risk region with quota constraints
We first estimate the probability of the non-risk region for a range of distributions, dimensions and constraints. We calculate these probabilities only for the Normal and t distributions as skew-t distributions and moment matching scenario sets use surrogate risk regions based on these distributions. The main purpose of this is to provide intuition about under what circumstances the methodology is effective: there is little to be gained from aggregating scenarios in a non-risk region of negligible probability.
For each distribution we sample 2000 scenarios and calculate the proportion of points in the nonrisk region for different levels of β and constraints. In particular, for each dimension we calculate for β = 0.95, and β = 0.99, and for a range of quotas. The feasible region corresponding to quota 0 < q < 1 Quotas are quite a natural constraint to use in the portfolio selection problem as they ensure that a portfolio is not overexposed to one asset. The quotas may also be viewed as ghost constraints to be used in cases where the probability of the non-risk region with only positivity constraints is too small.  Figure 7 for each each dimension we have tested we have plotted the results of one trial. The full results can be found in Appendix A. The first important observation from these is that the proportion of scenarios in the non-risk region, as compared to the uncorrelated case in Figure 4, is surprisingly high; even for β = 0.95 and dimension 30, this proportion is non-negligible. As expected, the proportion of scenarios in the non-risk region increases as we tighten our quota. However, for higher dimensions the quotas need to be a lot tighter to make a significant difference. The plots also provide further evidence that the t distribution has non-risk regions with higher probabilities than the lighter-tailed Normal distribution. In Figure 7, the non-risk region for the t-distributions has probability around 0.05 to 0.1 higher for dimensions 5 and 10, and around 0.1 to 0.2 higher for dimensions 20 and 30.

Aggregation sampling
In this section we compare the quality of solutions yielded by sampling and aggregation sampling by observing the optimality gaps of the solutions that these scenario generation methods yield. For this, we use the following version of the portfolio selection problem.
where µ is the mean of the input distribution (rather than scenario set), τ is the target return and u is a vector of asset quotas. The primary reason for using this formulation over those in (P2) and (P3) is that given a distribution of asset returns it is easy to select an appropriate expected target return τ . For simplicity, in our tests we set τ = 1 n n i=1 µ i which ensures that the constraint is feasible but not trivially satisfied. Notice that in the above formulation we use the deterministic constraint, x T µ ≥ τ rather than E x T Y ≥ τ . This is because the latter constraint depends on the scenario set. Therefore, the solution from a scenario-based approximation might be infeasible with respect to the original problem, which makes measuring solution quality problematic.
In this experiment, we test the performance of the aggregation sampling algorithm for three families of distributions: the Normal distribution, the t-distribution and the skew-t distribution.
For each distribution and problem dimension we run five trials using our constructed distributions (as described in Section 6.1). Each trial consists of generating 50 scenario sets via sampling and aggregation sampling, solving the corresponding scenario-based problem for each of these sets, and calculating the optimality gap for each solution which is yielded. For each scenario generation method we then calculate the mean and standard deviation (S.D.) of the optimality gap. For the skew-t distributions, although we are able to evaluate the objective function value for any candidate solution, to find the true optimal solution value (or one close to it), we resort to solving the problem for a very large sampled set of size 200000.
The full results for this experiment can be found in Appendix B. In Figure 8 we have plotted for one trial the raw results for dimensions 10 and 30. We observe that there is a consistent improvement in solution quality in using aggregation sampling over basic sampling. In addition the solution values are much more stable. The improvement in solution quality and stability is particularly big for the t-distributions. This is because the probability of the non-risk region is greater for heavier-tailed distributions as observed in Section 5. Aggregation sampling even lead to consistently better solutions for the skew-t distributions where we are approximating the risk region with a risk region for a t-distribution.

Aggregation reduction
The aim of these tests is to quantify the error induced through the use of aggregation reduction. In particular, we calculate the error induced in the optimal solution value. For a given scenario set, we aggregate the non-risk scenarios, solve the problem with respect to this reduced set, and calculate the optimality gap of this solution with respect to the original scenario set.
For these tests we use the same problem as in Section 6.3 and run tests for Normal, t and moment matching distributions. As explained in Section 4, we use the risk region of a Normal distribution to approximate the risk region for moment-matched scenario sets. For each family of distributions and problem dimension we again run five trials for different instances of the distribution. In each trial for different initial scenario set sizes, n = 100, 200, 500, and two different levels of tail risk measure β = 0.95, 0.99, we calculate the reduction error for 30 different scenario sets and report the mean error.
The full results can be found in Appendix C. These show that the reduction error is generally very small, in fact for almost all problems using β = 0.95, there is no error induced. For β = 0.99, and the smallest scenario set size n = 100, there is a small amount of error (< 0.01) for the Normal distributions, slightly larger errors for the heavier-tailed t distribution (< 0.02), and the largest errors (0.1-0.5) occur for reduced moment-matching scenario sets whose risk regions have been approximated with that of a Normal distribution. However, as the scenario set size is increased, all errors are reduced, and for the largest scenario set size n = 500, there is no error induced for almost all problems.
Comparing the reduction errors with the corresponding non-risk region probabilities in Appendix A, we see that the larger errors generally occur for the higher dimensional distributions whose non-risk region has a larger probability. This is to be expected as the larger the non-risk region the more scenarios that are aggregated. In Table 33 in Appendix C we have also included the proportions of reduced scenarios for moment matching scenario sets for which we approximated the risk region with that of a Normal distribution. The proportions of reduced scenarios in this case are generally slightly higher than that of the corresponding Normal distributions. This might suggest that the surrogates for the risk region are slightly too small, but this could equally be explained by the fact that moment matching scenario sets generally have heavier tails than the corresponding Normal distribution, which, as we observed in Section 5, also leads to non-risk regions of higher probabilities. In either case, the larger errors which are induced by reducing small moment-matching scenario sets could be explained by these increased probabilities.

Case study
In this section we demonstrate the application of the presented methodology on a difficult problem which may occur in practice. The problem is characterized by a high-dimension, a heavy-tailed non-elliptical distribution of asset returns, and the presence of integer variables. We compare the performance of the SAA method with ghost constraints algorithm proposed in Section 4.2 against the standard SAA method using basic sampling and aggregation sampling without ghost constraints.

Problem construction
The following problem is used: This problem is similar to that used in Section 6.3 except that we now use binary decision variables to limit the number of assets in which one can invest. The extra constraints involving integer variables may change the conic hull of feasible portfolios, however, the method presented in Section 3.3.1 for calculating conic hulls of feasible regions cannot handle these. We therefore ignore these constraints when constructing a risk region. This is acceptable as the resulting conic hull will contain the true conic hull. The problem is solved for β = 0.99, d = 50 assets, and maximum number of assets l = 10. The random vector Y of asset returns is constructed by fitting Skew-t distributions to historical monthly return data for companies from the FTSE100 stock index. The risk region used for this distribution is constructed as described at the end of Section 6.1.

Details of the SAA method
Estimation of optimality gap In each iteration of the SAA method the optimality gap for each of the found solutions is calculated using the procedure described in [KSHdM01]. Let N be the scenario set size and M be the number of replications used in each iteration. For m = 1, . . . , M denote by g m N and ν m N respectively the objective function and optimal solution value corresponding to the m-th constructed scenario set. In this case, g m N (x) = β -CVaR(−x T Y m ) where Y m denotes the discrete random vector corresponding to the m-th constructed scenario set.
Estimators for the optimal solution value and objective function are given by: Now, an α-level confidence interval for the optimality gap estimator by the solution x is given bȳ whereS 2 M is the standard deviation of g m N (x)−ν m N over m = 1, . . . , M , and Φ is the cumulative distribution function for a standard Normal distribution. Note that other procedures for estimating the optimality gap exist which only require the solution of one or two problems [BM06], [SB13].
Sample sizes, replications and bounds For this experiment, the initial sample size used to solve this problem is N = N 0 = 200 and at each iteration this is increased by N 0 /2 = 100. The number of replications used in each iteration is fixed at M = 10. These update rules have been chosen for simplicity; more sophisticated rules for updating the sample sizes and number of replications can be found in, for example, [RS13].
For the ghost constraints, the following heuristic rule is used. First denote byx m N for m = 1, . . . , M the solutions found for each replication. At the end of each iteration the bounds l ≤ x ≤ u are updated as follows: wherex M N andσ M N are the element-wise mean and standard deviation of the solutionsx m N over m = 1, . . . , M , and the parameter 0 < α < 1 controls how aggressively the ghost constraints are tightened. In this experiment we use α = 0.99.

Solution validation
Given the potential dangers in approximating the risk region, and misspecifying ghost constraints, it is important to verify the quality of a solution by the calculation of its corresponding out-of-sample value [KW07]. That is, after we calculate the β-CVaR for all candidate solutions with respect to a large independently sampled scenario set for all solutions yielded by the final iteration of the SAA methods. For this experiment a sample size of 100000 is used for validation.

Results
The results to this experiment are shown in Figure 9. In Figure 9a is shown the best optimality gap found at the end of each iteration of the SAA method, in Figure 9b are shown box plots with the out-of-sample values of the final solutions yielded by each method, and to aid our interpretation of the results. In Figure 9c we have plotted the evolution of the probability of the non-risk region used in the SAA method with ghost constraints.
In terms of the optimality gap and the final out-of-sample values, both aggregation sampling methods signifiantly outperform basic sampling. For the smallest sample size N = 200, the best optimality gap of aggregation sampling with and without ghost constraints is similar, since at this point no ghost constraints have been added. As the sample size increases, the ghost constraints become tighter and the probability of the non-risk region increases. This leads to a much smaller best optimality gap for aggregation sampling with ghost constraints. Since we cannot verify that the ghost constraints added are valid, we must view this gap with some caution. However, the out-of-sample validation after the final iteration reveals that aggregation sampling with ghost constraints is indeed producing higher quality solutions than aggregation sampling without ghost constraints.

Conclusions
In the paper [FTW17] we proposed a general approach to scenario generation using risk regions for stochastic programs with tail risk measure. As proof-of-concept we demonstrated how this applied for portfolio selection problems for elliptically distributed returns. In this work, we have presented how this methodology may be used for more realistic portfolio selection problems, in particular those which are high-dimensional, have non-elliptical assets returns and integer decision variables.
The main issue in applying the methodology to more realistic problems was its extension to nonelliptical distributions. In order to do this, we proposed the use of approximate risk regions, and derived results which indicated our approach would be robust against small misspecifications of risk region. Although this paper was focused on portfolio selection, the results were derived for general stochastic programs, which means that using approximate risk regions could work for other problems with tail risk measures.
We tested the performance of our methodology for solving realistic problems where the return distributions were fitted from real financial return data. Aggregation sampling generally outperformed basic sampling in terms of solution quality and stability. We also showed that aggregation reduction induces almost no error in the solution for reasonably sized scenario sets. These results not only held for elliptical distributions, but also non-elliptical distributions for which we used approximate risk regions.
The effectiveness of using risk regions for scenario generation depends upon the probability of the risk region: the greater the probability of the non-risk region, the more scenarios that can be aggregated. It follows directly from the definition of risk regions that this probability decreases as the problem becomes more constrained. Based on this observation, a heuristic based on the SAA method was proposed, which adds artificial constraints, called ghost constraints, to the problem. As the algorithm progresses, the ghost constraints are tightened which in effect allows one to focus in on high quality solutions. This algorithm was demonstrated on a difficult case study problem, and was shown to significantly outperform basic sampling and aggregation sampling without ghost constraints. This algorithm was presented in a non-problem specific way and could potentially be applied to other stochastic programs with tail risk measure.

A Reduction proportion tables
The following tables list the estimated probabilities of the non-risk region for a variety of distributions constructed from real data. See Section 6.2 for details. Each table corresponds to a family of distributions at a given dimension, and each row gives the proportions for a given set of companies. In addition, the distributions corresponding the i-th row of each table of dimension d have been fitted using the same set of companies.         Table 5: Proportion of reduced scenarios for t 4.0 distributed returns and d = 5