Portfolio optimization by improved NSGA-II and SPEA 2 based on different risk measures

In this study, we analyze three portfolio selection strategies for loss-averse investors: semi-variance, conditional value-at-risk, and a combination of both risk measures. Moreover, we propose a novel version of the non-dominated sorting genetic algorithm II and of the strength Pareto evolutionary algorithm 2 to tackle this optimization problem. The effectiveness of these algorithms is compared with two alternatives from the literature from five publicly available datasets. The computational results indicate that the proposed algorithms in this study outperform the others for all the examined performance metrics. Moreover, they are able to approximate the Pareto front even in cases in which all the other approaches fail.


Introduction
The portfolio selection problem can be defined as the optimal allocation of wealth among a finite number of assets that follows careful processing of all available information about both investors and markets (Meucci 2009).Markowitz's mean-variance model is by far the most popular procedure in asset allocation (Guerard 2009).Even if it is considered the cornerstone in this field, the mean-variance portfolio optimization model presents two serious drawbacks from a theoretical point of view.First, when asset returns are skewed and fat-tailed, they tend to include only a limited proportion of stochastically dominant assets in the efficient solutions, and prematurely preclude asset with negatively skewed returns.Second, risk is measured by variance, which treats both the above and below target returns equally, while investors are more concerned about the probability of investment returns falling below the target return.Consequently, risks are under-estimated and portfolios that are downside efficient are ruled out.
A solution proposed in the literature to make Markowitz's approach more effective is to replace the variance with a downside risk measure in order to model the loss aversion of investors properly.In this context, semi-variance has been studied extensively (Nawrocki 1999;Sing and Ong 2000).Two main reasons justify these efforts.First, semivariance is an approximation of the skewness for the return distribution, since it measures below-target returns.The higher is the value of semi-variance, the greater are both the degree of negative skewness and the risk of the investment.Second, semi-variance efficient portfolios closely approximate the stochastic-dominance efficient set.However, the computation of portfolio semi-variance is a difficult task owing to the endogenous nature of the portfolio co-semi-variance matrix, which depends on the weights assigned to each asset, that is, changes in the weights affect the periods in which the portfolio underperforms the target level, which, in turn, affects the evaluation of the co-semi-variance matrix itself.Some attempts have been proposed to solve this problem directly (Hogan and Warren 1972;Konno et al. 2002;Markowitz 1959;Markowitz et al. 1993).Other studies have focused on the definition of an exogenous co-semi-variance matrix that satisfactorily approaches the endogenous one (Ballestero and Pla-Santamaria 2005;Cumova and Nawrocki 2011;Estrada 2008;Nawrocki 1991).Nowadays, because of the regulatory importance of quantifying large losses in banking and insurance, another class of downside risk measures, called quantile-based measures, occupies a leading position in the risk management sector.One of the most popular examples is value-at-risk (VaR), defined as the maximum loss occurring over a given period at a given confidence level.Although VaR is apparently easy to use and intuitive, it presents several disadvantages.Primarily, it ignores losses exceeding VaR and is not sub-additive, that is, diversification of the portfolio may increase VaR (Artzner et al. 1999).From a computational perspective, VaR is difficult to use when investors want to optimize their portfolios, since it is represented by a non-linear, non-convex, and non-differentiable function with multiple local optima (Gaivoronski and Pflug 2005).Moreover, rational agents wishing to act on their decisions according to expected utility theory may be misled by the information related to the portfolio VaR (Yamai et al. 2002).To deal with these shortcomings, (Rockafellar et al. 2000) introduced the conditional value-at-risk (CVaR), which is defined as the conditional expectation of losses above the VaR (Sarykalin et al. 2008).CVaR is a coherent risk measure in the sense of (Artzner et al. 1999) and, because it is a convex function, optimization problems with CVaR as the minimization objective and/or constraints can be efficiently handled (Krokhmal et al. 2002;Larsen et al. 2002;Rockafellar and Uryasev 2002).
Increasing complexity of practical applications has led researchers to develop heuristic procedures for solving their portfolio optimization problems.These techniques require less domain information to be considered than the standard gradient-based mathematical programming methods do.Moreover, they guarantee satisfactory approximations to solutions in a fair computational time even when they deal with non-convexity, discontinuity, and integer decision variables.The approaches that have been proposed in the soft-computing literature can be categorized into the following two groups.On one hand, single objective methods optimize a weighted sum of the portfolio objectives.On the other hand, multi-objective evolutionary algorithms (MOEAs) attempt to tackle the allocation problem directly in its multi-objective form by simultaneously optimizing risk and reward.In the first case, the complete set of risk-return profiles is obtained by varying a parameter that represents the risk aversion of the investor (Chang et al. 2000;Crama and Schyns 2003;Cura 2009;Woodside-Oriakhi et al. 2011).In the second case, the complete efficient frontier is represented in a single run (Anagnostopoulos and Mamanis 2011a;Meghwani and Thakur 2017;Mishra et al. 2014).Both categories pay great attention to encoding types and constraint-handling techniques (Liagkouras and Metaxiotis 2015;Meghwani and Thakur 2017;Metaxiotis and Liagkouras 2012;Ponsich et al. 2013).Some real-life situations, which are not considered in Markowitz's model, have been analyzed recently by (Eftekharian et al. 2017) and (Meghwani and Thakur 2018).(Eftekharian et al. 2017) include as constraints some restrictions on the number of assets in the portfolio, limitations on investing in assets from a given industry, and cardinality, class, and quantity constraints.Furthermore, they developed an improved version of the NSGA II algorithm, called 2-Phase NSGA II, to solve the resulting optimization problem.Meghwani and Thakur (2018) focus on the problem of handling equality constraints, like self-financing constraints, and constraints arising from the inclusion of transaction cost models using MOEAs.Researchers have also focused on so-called swarm intelligence methods to overcome the computational difficulties of realistic portfolio designs.Unlike evolutionary algorithms that utilize the principle of natural selection, these approaches are inspired by the behavior and self-organizing interaction among agents, such as foraging of ant and bee colonies, bird flocking, and fish schooling (see Ertenlice and Kalayci (2018) for a detailed review of the subject).
Over the past few decades, machine-learning algorithms have been widely used to explore financial data and make data-driven predictions (Chao et al. 2019;Kou 2019).For instance, (Huang and Kou 2014) present a kernel entropy manifold learning algorithm to measure the relationships between two financial data points in order to describe the characteristics of a financial system by deriving the dynamic properties of the original data space.Similarly, (Huang et al. 2017) propose an information metric-based manifold learning algorithm to extract the intrinsic manifold of a dynamic financial system and to detect impending crises.(Ergu et al. 2014) focus on the use of the analytical network process in risk assessment and decision analysis of an emergent event.They propose a new consistency index to assess the key factors of risks and analyze the impacts and preferences of decision alternatives.Zhang et al. (2019) develop three soft minimum cost models based on different weighted average operators for resolving consensus decision-making problems.The usefulness of the proposed models is validated in a real-world loan consensus problem, using data from a Chinese P2P platform.Kalayci et al. (2019) review state-of-the-art methods dedicated to mean-variance portfolio optimization.
Evolutionary computations are made by iterations and, in each iteration, the weights of the portfolios are known a priori because they are generated by the algorithm.Thus, it is possible to directly determine in which periods the portfolio underperforms the target level.In this manner, the difficulty of computing the ex-post portfolio semi-variance can be fruitfully resolved by the use of an evolutionary algorithm.One of the first studies to attempt this is (Dueck and Winker 1992), who reformulate the bi-objective optimization problem into a single-objective and solve it through a local search technique, called threshold accepting.In a similar way, (Arnone et al. 1993) propose a solution to a parametric programming problem with the objective of optimizing a convex combination of mean and semi-variance.The algorithm is based on Whitley's GENITOR system (Whitley 1988) and uses a steady-state breeding strategy and elitist selection.Chang et al. (2009) apply a genetic algorithm to portfolio optimization problems in different risk measures, namely variance, semi-variance, absolute deviation, and variance with skewness.Recently, (Liagkouras and Metaxiotis 2013) address the mean-semi-variance portfolio optimization problem from a multi-objective perspective by developing an ad-hoc evolutionary algorithm.Evidence of the robustness of the algorithm is accomplished in out-of-sample testing during both bull and bear market conditions on the FTSE 100.Macedo et al. (2017) compare the non-dominated sorting genetic algorithm II (NSGA-II, Deb et al. (2002)) and the strength Pareto evolutionary Algorithm 2 (SPEA 2, Zitzler et al. (2001)) within the mean-semi-variance portfolio optimization framework.Numerical experiments indicate that NSGA-II outperforms SPEA 2 in-sample.Senhaji et al. (2016) propose to resolve the problem by combining the continuous Hopfield neural network with NSGA-II.The effectiveness of this strategy is proved using a portfolio of 20 assets.
Regarding the quantile-based minimization problems involving many financial instruments, a large number of price scenarios is necessary to estimate risk correctly.As the dimension of the problems increases, this operation can be time consuming (Lim et al. 2010) and the use of heuristics may be advisable or even necessary to detect a solution.An example is provided by the multi-purpose data-driven optimization heuristic proposed in Gilli M and Hysi (2006), which deals with different risk functions, such as VaR, CVaR, maximum loss, and Omega, as well as with practical constraints on the portfolio composition.An evolutionary computation approach is developed by (Hochreiter 2007) to solve the general scenario-based risk-return portfolio optimization problem when standard deviation, VaR, or CVaR are used to represent the risk of the investment.A variant of the NSGA-II algorithm based on differential operators is developed by (Krink and Paterlini 2011) for portfolio optimization involving real-world constraints and quantilebased risk measures.Meanwhile, a comparison of the capabilities of different MOEAs to adapt in any addition of new constraints and/or replacement of the risk function is given in (Anagnostopoulos andMamanis 2011b), and(Baixauli-Soler et al. 2012) focus on SPEA 2 performance.Recently, hybrid stock trading systems based on evolutionary metaheuristics and mean-CVaR models are proposed in Chen and Wang (2015) and (Qin et al. 2014).
The contribution of this study is twofold.On one hand, based on the promising results of (Subbu et al. 2005) and (Baixauli-Soler et al. 2010), we include in the same portfolio optimization framework the loss-averse attitude of investors as well as the capital requirements imposed by the regulator, and we investigate the relationship between semivariance and CVaR in quantifying the downside risk.On the other hand, we extend the research of (Liu et al. 2010) and investigate the effectiveness of using the uniform selection scheme, the extended intermediate crossover operator (Gen and Cheng 2000;Mühlenbein and Schlierkamp-Voosen 1993), and the Gaussian mutation (Hinterding 1995;Schwefel 1987) in the NSGA-II and SPEA 2 algorithms in order to generate the approximated Pareto fronts for the considered downside risk-based portfolio optimization problems.In addition, the performance of the proposed algorithms is compared to that of other variants of the NSGA-II and SPEA 2 algorithms that have already been shown to be highly competitive in portfolio optimization problems.The results on five publicly available datasets show that our procedures completely dominate the fronts produced by the counterparts from the literature.Furthermore, our variants of the NSGA-II and SPEA 2 algorithms can generate the entire Pareto front for large-scale problems for which the other Pareto-based approaches are unable to work properly.
The rest of the paper is organized as follows.In "The portfolio selection problem under downside risk measures" section, we introduce three preference relations based on the reward-downside risk principle and discuss the related portfolio selection problems."Pareto-based evolutionary algorithms" section presents a description of the MOEAs used for solving the resulting optimization problems."Experimental analysis" section outlines the numerical experiments and "Concluding remarks and future research" section presents concluding remarks and ideas for further research.

The portfolio selection problem under downside risk measures
Consider a financial market consisting of n risky assets, indexed from 1 to n, and modeled by a probability space ( , F, P).Let agents act their investment decisions over a one-period horizon with respect to the following set of feasible, or admissible, portfolios defined according to the budget constraint and no-short selling: where x i represents the proportion of capital to be allocated to asset i, with i = 1, . . ., n.
The random rate of return of asset i at the end of the investment period, denoted as R i , is assumed to have a continuous probability density function (pdf ).Hence, the rate of return of a portfolio x = (x 1 , . . . where E(•) denotes the expectation and E(R i ) is the expected rate of return of asset i.Given a level z for the rate of return, the cumulative distribution function (cdf ) of R(x) is defined as It is also assumed that F R(x) is continuous and strictly increasing with respect to z ∈ R 1 .The portfolio loss distribution is defined as the negative of the portfolio return distribution, that is, L(x) = −R(x).
To identify the portfolio in X that guarantees the "best" rate of return, a model for preferences under uncertainty needs to be defined.We adopt the so-called reward-risk approach that relates the portfolio selection problem to a multi-objective optimization problem in two steps.First, a set of objectives that the investor perceives as beneficial is identified, and second, a set of objectives he or she considers damaging is identified in relation to R(•).Then, a preference relation is defined based on these criteria as follows.
Definition 1 Let f 1 , . . . ,f p : X → R and g 1 , . . ., g q : X → R be the reward-type and the risk-type objective functions, respectively.Then for all x, y ∈ X , we say that R(x) dominates (is preferred to) R(y) if and only if f i (x) ≥ f i (y) for all i = 1, . . ., p and g j (x) ≤ g j (y) for all j = 1, . . ., q with at least one strict inequality.Alternatively, we can say that portfolio x dominates portfolio y.
According to this definition, a portfolio x * ∈ X is (Pareto) optimal if and only if it is non-dominated with respect to X , that is, another x ∈ X that dominates x * does not exist.Thus, an efficient portfolio in the reward-risk model is a Pareto optimal solution of the following multi-objective problem: minimize F(x) = −f 1 (x), . . ., −f p (x), g 1 (x), . . ., g q (x) subject to x ∈ X . (3) The solutions to Problem (3) form the so-called efficient set, or Pareto optimal set, of which image in the reward-risk space is called efficient frontier, or Pareto front.Under certain smoothness assumptions, it can be induced from the Karush-Kuhn-Tucker conditions that, when the p + q objectives are continuous, the efficient frontier defines a piecewise continuous (p + q − 1)-dimensional manifold in the decision space (Li and Zhang 2009).Therefore, the efficient frontier of a continuous bi-objective portfolio optimization problem is a piecewise continuous 1-D curve in R 2 , while, in the case of a problem with three objectives, it is a piecewise continuous 2-D surface in R 3 .We now specialize this general reward-risk framework by considering investors that focus on the mean rate of return (2) as the reward criterion and employ different measures to assess the downside risk of their investments.

The mean-semi-variance model
One of the most well-known downside risk measures of a portfolio x is its semi-variance, already introduced in the seminal paper of (Markowitz 1959) and defined as where [u] − = max{0, −u} and b represents a given benchmark, for example, the portfolio expected rate of return, the Roy's safety first criterion, or a Treasury rate of return.In our analysis, we set the target level b equal to 0 in order to estimate the variance of portfolio losses.As proved in Fishburn (1977), this downside risk measure is consistent with utility theory An investor that operates her or his decisions according to the mean-semi-variance model is interested in maximizing the expected rate of return of her or his portfolio and, simultaneously, minimizing its semi-variance: The mean-CVaR model The CVaR of a portfolio x at the confidence level α, with α ∈ (0, 1), represents the average of losses in the 100(1 − α)% worst scenarios (Acerbi and Tasche 2002).According to our assumptions about the cdf of R(x) and the definition of the portfolio loss L(x) as the opposite of R(x), this downside risk measure is defined as where is the value-at-risk at the confidence level α of the portfolio x.
An agent whose preferences are dictated by the mean-CVaR model would select her or his portfolio weights according to the following bi-objective optimization problem: The mean-semi-variance-CVaR model Using only one criterion to model risk in a portfolio selection model may provide a restricted picture of the assessment process (Steuer et al. 2005).It is advisable to include multiple risk measures in order to cover the different facets of an investor's preferences.
In this study, we analyze the combined effects of CVaR and semi-variance.The related tri-objective optimization problem is where V − and CVaR α are the same as the definitions in ( 4) and ( 6), respectively.In this manner, the preferences are modeled as follows.
Definition 2 Portfolio x is preferred to portfolio y if and only if This model may produce improved solutions when a mean-CVaR efficient portfolio has an excessively large semi-variance or when a mean-semi-variance efficient portfolio has an excessively large CVaR.

Scenario-based framework for portfolio optimization
A practical way to handle the random variables that represent assets and portfolio rates of return is to treat them as discrete.Therefore, let = {ω 1 , . . ., ω S }, P(ω s ) = p s , with s = 1, . . ., S and S s=1 p s = 1, F be a given σ -field and assume the following table of scenarios is available where S represents the number of involved scenarios, r s = (r 1s , . . ., r ns ) T is the n-vector of rates of return for the s-th scenario, and p s is the associated probability of occurrence, with s = 1, . . ., S. Assume that p s = 1/S for all s.Thus, the rate of return of x ∈ X may assume r is x i .Analogously, assuming that the distribution of the portfolio rate of return R(x) is such that P R(x) = r p s (x) = 1/S with s = 1, . . ., S, since the expected rate of return for the i-th security can be computed as the mean rate of return over the S scenarios (i.e., E(R Now it is possible to reformulate the previously introduced downside risk-based portfolio allocation models in terms of the scenarios conveyed in (9).Based on the findings of (Cumova and Nawrocki 2011), we estimate V − (R(x)) by means of the so-called co-semi-variance matrix, henceforth denoted as C − = C − ij , which is defined element-wise as for all i, j = 1, . . ., n.It can be noted that C − is an asymmetric matrix, since in general, assets i and j do not have the same below-target returns during the same period, that is, The semi-variance of portfolio x can then be calculated as In this manner, Problem (5) translates into the maximization of a linear objective function, given by Eq. ( 10), and the simultaneous minimization of a quadratic objective function, given by Eq. ( 13), over the unit standard n-dimensional simplex.Regarding Problem ( 7), the CVaR expression ( 6) can be evaluated as an arithmetic mean over scenarios as follows.Let l p s (x) = −r p s (x) be the portfolio loss in the s-th scenario, after sorting losses in ascending order, that is, (x), the VaR α can then be defined as l p ( αS ) (x) and CVaR α can be estimated as where I(•) represents the indicator function.In this case, the investor maximizes the expected rate of return as given in Eq. ( 10) and minimizes the risk calculated according to Eq. ( 14).Finally, Problem (8) simultaneously exploits Eqs. ( 10), (13), and ( 14) in the optimization process.

Pareto-based evolutionary algorithms
The capabilities of MOEAs to generate reasonably good approximations of the Pareto front in a single run and within a limited computational time have already been shown in the literature on the mean-variance portfolio selection problem.
In this study, we focus on two variants of popular algorithms, NSGA-II and SPEA 2, and investigate their effectiveness in solving the portfolio optimization problems formulated in the previous section.These algorithms belong to the family of Pareto-based MOEAs and include a two-level ranking scheme to guide the search toward the true Pareto front (Emmerich and Deutz 2018).The first ranking is provided by the Pareto dominance relation, while the second concerns the diversity of the solutions and applies to the individuals that share the same position in the first ranking.However, the methods by which NSGA-II and SPEA 2 approximate the true Pareto front differ and the corresponding procedures are described in the following overview 2 .

NSGA-II
The NSGA-II procedure is described in Algorithm 1.The crowding distance mechanism is employed to preserve the diversity of solutions.It evaluates the volume of the hyper-rectangle defined by two nearest neighbors (Zitzler and Thiele 1999).New solutions, called offsprings, are generated using a selection mechanism and a set of variation operators.Based on the values provided by the ranking scheme, the best individuals from the combination of the current population P k and the offspring pool Q k are detected and those with lower rank and higher crowding distance are saved in the next population P k+1 .In case some candidate solutions are of the same rank and not all of them enter the next population, the less crowded individuals from that particular rank are selected to fit the next population, thereby ensuring elitism.
Algorithm 1: NSGA-II (Deb et al. 2002) As reported in Algorithm 2, SPEA 2 first initializes a population of candidate solutions P k , then stores the best solutions in an explicit archive A k , separate from the population.To emphasize non-dominated individuals, SPEA 2 uses a combination of the dominance count and the dominance rank methods.Each individual is assigned a raw fitness value depending on both the number of individuals it dominates and the number of individuals by which it is dominated.The density information is expressed as a function of the k-th smallest Euclidean distance in the objective space to the k-th nearest neighbor.The non-dominated individuals from the union of the archive and the current population are then updated.In particular, if the number of non-dominated individuals is less than the pre-established archive size, some dominated individual from the current pool form part of the archive.Otherwise, some individuals are removed from the archive using a truncation operator.This procedure recursively removes individuals based on the nearest neighbor Euclidean distance.If there is more than one candidate solution with the same minimum distance, then the decision is made by considering the second nearest neighbor, and so forth.The mating pool used to generate the next population P t+1 is filled by the individuals of the updated archive selected on the basis of a given selection mechanism.The offsprings are then generated by a set of variation operators as in the algorithm NSGA-II.

Solution approaches
NSGA-II and SPEA 2 cannot be directly implemented in their standard format for the solution of real-world portfolio selection problems and special treatments need to be adopted in relation to a number of issues, like the encoding type, genetic operators, and constraint-handling procedures.

Solution representation and initialization
We represent a portfolio by a real vector of size n, denoted by x = (x 1 , . . ., x n ) T ∈ R n , where x i is the proportion of budget invested in asset i, with i = 1, . . ., n.Each candidate solution in the initial population is generated by following Rubistein's (1982) procedure, which consists of the following steps: It is evident that, in this manner, all the individuals in the initial population are feasible.

Reproduction process
Two variants of the reproduction process are compared in this study, both employing a selection procedure, a crossover, and a mutation operator.
The first configuration is based on that proposed by (Anagnostopoulos and Mamanis 2011a) 3 .
1. Binary tournament selection is used for selecting the parents.In this process, two individuals are chosen randomly from the population and compete against each other.The individual with the highest fitness wins and is included as one parent for the next steps of the reproduction.2. The uniform crossover operator is then applied to the population of parents to produce the offspring population.Two parents generate a single child and its value for each array is selected with equal probability from one or another parent.3. The set of children is finally subject to the Gaussian mutation.A percentage P mut of these individuals is selected.Then, each member x of this subset has a probability μ m that a gene x i mutates according to the following rule: where σ m is the mutation step size, and rand 0 , rand 1 represent two independent random numbers from a standard normal distribution, with i = 1, . . ., n.
In the numerical comparisons, the algorithms using this set up are denoted by NSGA-IIb and SPEA 2b.The alternative we propose is as follows.
1. Uniform selection is used to generate two sub-populations of parents.The first of these sets contains the fraction P cross of the original population that is involved into the recombination step.The second set contains the fraction P mut of the original population that is subject to mutation.These sets cannot be separated, in that the same individual can enter the crossover stage and can mutate.2. The intermediate crossover is applied to the first sub-population of parents to generate a first set of offsprings.Differently from the uniform crossover, this operator generates two children for each pair of parents x 1 , x 2 , as follows where csf i , called the crossover scaling factor, is a random number in the interval is a parameter to be tuned, and x 1 , x 2 are the associated children, with i = 1, . . ., n.The second population of parents is modified by the Gaussian mutation as in the first configuration.3. The sub-populations of children are finally gathered together to form the offspring population.
The variants of the considered MOEAs involving this reproduction design are denoted by NSGA-IIa and SPEA 2a.

Constraint-handling procedure
There are several techniques proposed to handle constraints in MOEAs (see Mezura-Montes and Coello (2011) and references therein).Since the feasible set in our portfolio optimization problems is represented by the unitary simplex (1), the best procedure is to use the following two-step repair mechanism (Liagkouras and Metaxiotis 2015).
(a) Each candidate solution x ∈ R n is first clamped by projecting it onto [ 0, 1] n : with i = 1, . . ., n.In this manner, x = ( x 1 , . . ., x n ) T satisfies the lower bound constraints.(b) The projected vector x is then normalized through the transformation After this step, the individual x = x 1 , . . ., x n T also verifies the budget constraint.
This procedure makes all the individuals involved in the search feasible.The corresponding objective function values can now be computed and the MOEAs described above can be applied to identify an approximation of the Pareto set.

Computational complexity
Let us assume that the size of the population and the size of the external archive are equal to N, the number of objectives is m, the number of generations is G, and the dimension of the vector of decision variables is n.Based on the analysis presented by (Deb et al. 2002) for one generation, and taking into account the different definition of the parent population for the two configurations, the NSGA-IIb algorithm requires O(m(2N) 2 ) computations while the NSGA-IIa algorithm has a running time bounded by O m (1 + 2P cross + P mut ) 2 N 2 .In a similar manner, the time complexity of SPEA 2b is O (2N) 3 and that of SPEA 2a is O (1

Experimental analysis
In this section, we assess the effectiveness of the two variants of NSGA-II and SPEA 2 in solving the proposed instances of the portfolio optimization problem.

Description of the datasets
The experiments are based on five public datasets provided by (Bruni et al. 2016).They include the weekly linear returns of the constituents of the following stock markets: Dow Jones Industrial Average, Fama-French Market Industry, NASDAQ 100, S&P 500, and NASDAQ Composite.As shown in Table 2, the test problems range from 28 assets for the smaller case study to 1203 assets for the bigger one.From Table 3, it can be pointed out that, in all the datasets, an asset with a moderately or highly negative skewed return distribution (Skew < −0.5) presents a value of the standard semi-deviation (Stsd), that is, the squared root of the semi-variance, which is higher than the corresponding value of the standard deviation (Std).Conversely, an asset with a moderately or highly positive skewed return distribution (Skew > 0.5) has a Stsd value less than the value of the Std.Assets with a skewness between − 0.5 and 0.5 present approximately symmetric return distributions, resulting in very similar Stsd and Std values.
Consequently, loss-averse investors can use standard semideviation (or, equivalently, semi-variance) to capture the downside risk conveyed in skewness properly.In this sense, the portfolio selection models designed in the previous section are more appealing than the mean variance model for acting investment choices.

Performance metrics and statistical testing
It is non-trivial to evaluate the quality of the solution sets of MOEAs for two main reasons.First, the presence of multiple conflicting goals makes the definition of a "better algorithm" vague.Second, the stochastic nature of these optimizers suggests that comparisons based on the approximation sets from a single run of each algorithm are not correct.
Regarding the first problem, we can identify three major performance criteria in multiobjective optimization: the capacity of a given algorithm to generate an appropriate number of non-dominated solutions, the convergence of the solution set to the true Pareto front, and the diversity of the solutions in the objective space.Accordingly, optimal solution sets with a large number of non-dominated solutions, approaching the true Pareto front and even scattering are generally desirable (Jiang et al. 2014).
The metrics proposed in the literature usually synthesize the information content of the solution set by considering only one or two criteria at a time.Therefore, it is advisable to employ a combination of different performance metrics in order to provide a complete analysis of the experimental results.In this sense, we consider the following four metrics.(i) Schott's spacing metric (S ) (Schott 1995) The spacing metric measures how evenly the solutions are distributed in the approximate efficient front A and it is expressed as where d i represents the minimum value of the sum of the absolute difference in objective function values between the i -th solution and any other solution in the obtained non-dominated set and d is the mean value of these distance measures.
When the solutions are equidistantly spaced, the corresponding distance measure is small.Thus, an algorithm finding a set of non-dominated solutions with smaller spacing is better.
(ii) Generalized spread metric ( * ) (Zhou et al. 2006) The generalized spread metric is a generalization of the well-known metric of (Deb et al. 2002) and simultaneously gauges the distribution and spread of an optimal solution set A for high dimensional multi-objective optimization problems.It takes the form: * (A, P) = (iii) Inverted generation distance (IGD ) (Zhang and Li 2007) The IGD index has the following formulation where d i is the Euclidean distance (in the objective space) between the i -th member of the true Pareto front P and the closest solution in A. If |P| is large enough to represent the Pareto front very well, IGD(A, P) could measure both the diversity and convergence of A.
A low value of IGD(A, P) means A is very close to the Pareto front and does not miss any part of it.Thus, an algorithm with a lower value of IGD is better.
(iv) Hypervolume (HV ) (Zitzler and Thiele 1999) Similar to the IGD index, the hypervolume indicator evaluates both diversity and convergence of an approximation set A. It is defined as the size of the portion of objective space that is dominated by at least one point of A relative to a reference set R, which is formed by points worse than (or equal to) every point in A in every objective.Formally, we define it as where v i is the hypercube constructed with the reference set R and the solution a i ∈ A as the diagonal corners.Large values of HV indicate the approximate solutions are closer to the true Pareto front and, at the same time, scattered more evenly in the objective space.Thus, an algorithm with a large value of the HV metric is desirable.
The second question posed at the beginning of this section can be resolved by comparing the algorithms through a sample of approximation sets from multiple runs.In this manner, we obtain a sample of values for each performance metric and for each portfolio optimization problem.A rigorous comparison between the algorithms can then be performed based on non-parametric inference testing (Coello et al. 2007).
Since the true Pareto front is a (continuous) manifold in the decision space, we approximate it by a surrogate front, given by the non-dominated points from the union of the solution sets of all the algorithms for all the simulations.This set is used to calculate the * and IGD metrics.Furthermore, before conducting the performance analysis, we normalize each objective function value F i (x), i = 1, . . ., p + q, in (3) according to the following transformation where F min i and F max i are the minimum and maximum values of the i-th objective function for the corresponding surrogate Pareto front.In this manner, the objectives contribute almost equally to each performance metric and the results from the different algorithms can be fairly compared.

Parameter settings
Let us assume that the distribution of historical returns acts as a good proxy of the returns faced over the next holding period.Then, the historical simulation method is used to compute the financial scenarios r is , with i = 1, . . ., n and s = 1, . . ., S, as well as the co-semi-variance matrix C − .
All the variants of NSGA-II and SPEA 2 use identical population and archive sizes of 250, and 400 iterations for all the experiments as in (Anagnostopoulos and Mamanis 2011a).For the remaining parameters, the algorithms NSGA-IIb and SPEA 2b use the values proposed by the same authors.Instead, for NSGA-IIa and SPEA 2a, the parameters are tuned over the DowJones, FF49Industries, and NASDAQ100 datasets for each of the three portfolio optimization problems, for a total of nine test problems.For each pair of parents the intermediate crossover produces two offsprings while the uniform crossover provides a single child, and thus, to make fair comparisons among algorithms, we select the percentage of crossover P cross from the set {0.35, 0.45} and the parameter d, related to the crossover scaling factor, from {0, 0.5, 1}.The mutation percentage P mut varies in {0.3, 0.5}, the mutation probability μ m can assume a value in the set {0.1, 0.3}, while the mutation step size σ m is checked in {0.1, 0.15, 0.2}.In this manner, there are 72 possible combinations of these parameters.All the experiments are conducted 20 times for each algorithm and the hypervolume values corresponding to the final populations are recorded.Finally, since we search for a parameter set up that is widely applicable in finding suitable approximation fronts for all the test problems, the best group of parameter settings is selected based on the average ranking of the Friedman test for the values of the HV index over the simulations (Derrac et al. 2011).Tables 4 and 5 report the configurations with the highest average rankings for the proposed MOEAs.
Overall, the results show that the best choices are higher values for P cross and d and, at the same time, lower values for μ m and σ m .However, SPEA 2a needs a larger sample of mutants to achieve the best rankings with respect to NSGA-IIa.The configurations in bold are the best parameter settings and are used in the experiments.

Computational results and discussion
We perform a multi-problem analysis in which the four MOEAs are tested on 15 optimization problems (the three instances of the portfolio selection problem over the five datasets introduced in Table 2).To check the robustness of the results, 20 simulations for each algorithm and for each test problem are used.The algorithms are implemented in MATLAB R2018b and the experiments are carried out on a 2.2 GHz Intel Core i7 laptop with 4 GB RAM.
Table 6 reports the average number of nondominated solutions provided by the algorithms in the final approximated optimal sets, divided by problem and dataset involved.It can be noticed that including the proposed configuration in the reproduction process increases the search capabilities of both NSGA-II and SPEA 2 with respect to the standard configuration.In particular, NSGA-IIa and SPEA 2a can identify sets of optimal solutions with a size almost equal to the pre-defined threshold of 250 individuals in all the test problems, while NSGA-IIb and SPEA 2b produce much smaller optimal sets as the complexity of the problems increases.
Tables 7, 8 and 9 focus on the distribution of the values for the performance metrics, averaged over the 20 simulations for each portfolio optimization problem.Overall, the superiority of NSGA-IIa and SPEA 2a is evident with respect to all the criteria used in the comparisons.Regarding the S metric, for small-to medium-sized optimization problems, NSGA-IIa produces approximate efficient fronts with the smallest spacing.However, as the dimension of the problems increases, the standard versions of the algorithms become competitive.In terms of the generalized spread metric, NSGA-IIa confirms its good performance while SPEA 2a guarantees a better uniformly distributed set of non-dominated solutions for large-scale problems.Similar analysis can be undertaken for the inverted generation distance and for the hypervolume indicator.In these cases, NSGA-IIa and SPEA 2a are the best choices, and SPEA 2a is preferred for large-scale problems.
We check whether the differences in the performance are significant by using the R package called "scmamp" (Calvo and Santafé 2016).First, we apply the Friedman aligned omnibus test to detect if at least one of the algorithms performs differently to the others.The results reported in Table 10 show that the differences are significant for three out of   Thus, the analysis can proceed to characterize the differences found for the * , IGD, and HV indicators using a post-hoc procedure.As suggested in (Derrac et al. 2011), we use the Friedman aligned post-hoc test with the correction of (Bergmann and Hommel 1988).Table 11 lists the five hypotheses of equality among the four algorithms and the corresponding adjusted p-values.With a 1% level of significance, we find an improvement in the performance of NSGA-II and SPEA 2 when they involve the second selection and reproduction scheme.
The contrast between the medians of the samples of the results is finally computed considering all the pairwise comparisons (Garcia et al. 2010).This test points out the quantitative difference between two algorithms over multiple test problems, allowing us to estimate how much an algorithm outperforms another one.Table 12 lists the comparisons of the proposed algorithms, divided by the performance metric used.Notably, an estimator with a negative value for * and IGD, and a positive value for HV indicates that the algorithm in that row outperforms the algorithm in the corresponding column.The estimators highlight that SPEA 2a is the best performing algorithm, although the difference with NSGA-IIa is negligible.For illustrative purposes, Figs. 1 and 2 show a comparison of the Pareto fronts obtained by the four algorithms analyzed for the DowJones, NASDAQ100, and NASDAQComposite datasets corresponding to the simulations that attain the highest values of the hypervolume metric.The proposed versions of the algorithms NSGA-II and SPEA 2 can identify frontiers with more evenly scattered and spread out points than the standard alternatives in all the test problems.The differences among the algorithms are already marked for problems with 83 decision variables.As the size of the investible universe increases, NSGA-IIb and SPEA 2b no longer cover the shape of the Pareto front and tend to produce highly inefficient solutions, with higher risks and lower expected returns.Moreover, for the three-objective portfolio optimization problem, the charts on the left of Fig. 2 highlight the relationship between Stsd and CVaR, which is almost linear for higher values of risk.
Figure 3 displays the mean processing time of the four algorithms for each portfolio model and for each dataset.The procedures implementing the proposed configuration demand more computational resources owing to the use of larger offspring populations.The difference becomes more evident as the number of decision variables increases.Moreover, the total run time is greatly influenced by the time spent evaluating the CVaR objective function.The time spent on the management of the archive for the SPEA 2 algorithms is another cost to take into account.
In summary, balancing the results from the statistical analysis and those concerning the computational complexity, we suggest that the NSGA-II algorithm with the proposed configuration should be used to solve the portfolio optimization problems involving downside risk measures.In fact, it produces approximated sets similar to the SPEA 2a algorithm but takes less time to do so.

Concluding remarks and future research
In this study, we described three instances of the portfolio selection problem designed to handle the downside risk of an investment properly.A flexible multi-objective reward-torisk framework was presented in which expected returns, semi-variance, and CVaR of a   portfolio can be optimized simultaneously.These problems were tackled using two nondominated sorting algorithms, namely, NSGA-II and SPEA 2, which have already showed competitive performance for the mean variance problem.In particular, we proposed a novel combination of operators for the selection and reproduction phases to be included in both algorithms.A comparative analysis was undertaken with respect to a second variant of the same algorithms, involving another configuration design.We used five publicly available datasets ranging from small-to large-sized portfolio optimization problems.The capabilities of the procedures were assessed in terms of four performance metrics.Finally, a set of statistical tests checked the robustness of these findings.Overall, the numerical experiments showed that the proposed algorithms outperformed the others with respect to all the criteria.Even if the algorithms with the novel variation configuration demanded the use of more computational time as the dimension of the problems increased, they nonetheless yielded reasonable results for the cases in which the other algorithms failed to capture the shape of the Pareto front properly.These research findings can be put in practice to improve the risk management infrastructure of an investment company.The inclusion of several risk measures in the portfolio optimization process can increase the capabilities of the system to describe the risk, providing more attractive investment opportunities.
Future research work on the topic includes the analysis of out-of-sample effectiveness for this type of strategy, which is expected to be improved by the incorporation of other simulation techniques to estimate semi-variance and CVaR.We are also interested in exploring mechanisms that adaptively exploit several selection schemes and reproduction operators to accelerate convergence and allow the search to be stopped automatically when a suitable level of quality for the approximated set has been attained.

Endnotes
1 Under these assumptions, the distribution function F R(x) is such that no jumps and no flat parts occur, implying that the equation F R(x) (z) = α has a unique solution for any α ∈ (0, 1), say z * = F −1 R(x) (α), where F −1 R(x) denotes the inverse of F R(x) . 2 The terms solution and individual will be used interchangeably, since individuals in the population represent solutions to the problem that is being optimized.

do 4 7
Calculate fitness for P k and assign rank based on Pareto dominance 5 Perform selection on P k to fill mating pool 6 Apply recombination and mutation operators to obtain the offspring population Q k Select the best N non-dominated solutions from P t ∪ Q k by the two-step procedure to form P k+1 8 Set k = k + 1 9 end 10 return P k SPEA 2

Algorithm 2 :
SPEA 2(Zitzler et al. 2001) 1 Set k = 0 2 Initialize P 0 and set A 0 = ∅ 3 while k < G do 4 Calculate fitness for P k and A k 5 Copy all non-dominated solutions from P t ∪ A k to A k+1 6 Use the truncation operator to remove elements from A t+1 , if its capacity has been extended 7 Use dominated individuals from A k to fill A k+1 , if the capacity of A k+1 has not been exceeded 8 Perform selection on A k+1 to fill mating pool 9 Apply recombination and mutation operators to obtain P k+1 10 Set k = k + 1 11 end 12 return A k where d i denotes the Euclidean distance between neighboring solutions with the mean value d.The term d P m is the distance between the extreme (bounding) solutions of A and of the true Pareto front P corresponding to the m-th objective function.An algorithm finding a smaller value of * generates a better uniformly distributed set of non-dominated solutions.

Fig. 1 Fig. 2
Fig. 1 Approximated Pareto fronts for the Mean-semi-variance (on the left) and Mean-CVaR (on the right) portfolio optimization problems corresponding to the simulations with the highest HV values obtained by NSGA-IIb, NSGA-IIa, SPEA 2b and SPEA 2a for the datasets DowJones (charts at the top), NASDAQ100 (charts in the middle) and NASDAQComposite (charts below)

Fig. 3
Fig. 3 Mean total CPU times (in seconds) for the Mean-semi-variance (charts at the top), Mean-CVaR (charts in the middle) and Mean-semi-variance-CVaR (charts below) portfolio optimization problems for 100,000 function evaluations of the algorithms NSGA-IIb, NSGA-IIa, SPEA 2b and SPEA 2a for the considered datasets over 20 simulations

Table 1
Run time complexity for the algorithms

Table 3
Classification of the assets belonging to each dataset in terms of skewness (Skew), standard semideviation (Stsd) and standard deviation(Std)

Table 5
Average rankings achieved by the Friedman test for the SPEA 2a algorithm with different parameter settings

Table 6
Average number of nondominated solutions in the approximated Pareto front (Mean) with the corresponding standard deviation (Std) for each algorithm

Table 7
Mean, standard deviation (Std), median, minimum (Min) and maximum (Max) values of the performance metrics for the compared algorithms for the Mean-SV portfolio optimization problem

Table 8
Mean, standard deviation (Std), median, minimum (Min) and maximum (Max) values of the performance metrics for the compared algorithms for the Mean-CVaR portfolio optimization problem a Numbers in bold represent best results among the algorithms with respect to the corresponding performance metric

Table 9
Mean, standard deviation (Std), median, minimum (Min) and maximum (Max) values of the performance metrics for the compared algorithms for the Mean-semi-variance-CVaR portfolio optimization problem a Numbers in bold represent best results among the algorithms with respect to the corresponding performance metric

Table 10
Statistics and related p-values for the Friedman aligned omnibus test for each of the four performance metrics

Table 11
Adjusted p-values for the Friedman aligned post-hoc test with Bergmann and Hommel's correction for multiple comparisons among the four algorithms

Table 12
Contrast estimation results for the performance metric