A new IPSO-SA approach for cardinality constrained portfolio optimization

Article history: Received 1 October 2010 Received in revised form 7 January 2011 Accepted 10 January 2011 Available online 14 January 2011 The problem of portfolio optimization has always been a key concern for investors. This paper addresses a realistic portfolio optimization problem with floor, ceiling, and cardinality constraints. This problem is a mixed integer quadratic programming where traditional optimization methods fail to find the optimal solution, efficiently. The present paper develops a new hybrid approach based on an improved particle swarm optimization (PSO) and a modified simulated annealing (SA) to find the cardinality constrained efficient frontier. The proposed algorithm benefits from simple and easy characteristics of PSO with an adaptation of inertia weights and constriction factor. In addition, incorporating an SA procedure into IPSO helps escaping from local optima and improves the precision of convergence. Computational results on benchmark problems with up to 225 assets signify that our proposed algorithm exceeds not only the standard PSO but also the other heuristic algorithms previously presented to solve the cardinality constrained portfolio problem. © 2011 Growing Science Ltd. All rights reserved


Introduction
Today, managers and investors confront a wide range of possible securities or assets where each has its own potential risk and rewards.The traditional portfolio optimization techniques attempt to select an appropriate portfolio of assets with a reasonable trade-off between the risk and the rate of return.The mean-variance approach introduced by Markowitz (1954) for portfolio optimization problems plays an important role in developing modern portfolio selection techniques.Markowitz model is a quadratic programming problem, which can be solved through exact methods such as active set methods, interior point techniques, etc.
Markowitz seminal work has been widely extended in different ways.Following his work, alternative models such as mean-absolute deviation (MAD) and absolute deviation were proposed for the same problem.Konno and Yamazaki (1991) proposed a linear programming model which is equivalent to Markowitz model whenever returns are multivariate and normally distributed.Ignoring the covariance matrix in the Markowitz model results in great estimation risk and could cause serious error (Simaan, 1997).A piecewise linear approximation, weighted goal programming (Lee & Chesser, 1980) and mini-max model (Young, 1998) are the other examples of these models.Some authors have added some practical constraints such as transaction costs, liquidity, buy-in threshold, cardinality, etc to Markowitz model to make it more realistic.However, adding constraints to the portfolio optimization problem make it intractable even for small instances.
A relatively new meta-heuristic approach, where its application in portfolio optimization is still limited, is particle swarm optimization (PSO).Mous et al. (2006) compared PSO and GA which was applied to portfolio selection, Xu et al. (2007) and Gao and Chu (2010) proposed improved PSO algorithms for portfolio problem with transaction costs and quality constraints.Mishra et al. (2009) studied portfolio optimization problem from multi-objective perspective and considered four wellknown multi-objective evolutionary algorithms including multi-objective PSO (MOPSO) which outperformed its counterparts for some numerical experiments.Cura (2009) presented a PSO algorithm to solve cardinality constrained portfolio.Compared with other meta-heuristics, PSO is simple, high speed, large scope and easy to be implemented by programs.However, there are still many issues in particle swarm, such as slow convergence during the latter search, poor precision and converging to local optimum.To overcome the above problems, a lot of revised PSO have emerged.
Investors generally incline to restrict the number of assets in the portfolio and purchase just a subset of them.Therefore, cardinality is a practical constraint that has to be considered in decisions.In this paper, we address a portfolio optimization problem with floor, ceiling, and cardinality constraints.This problem is a mixed integer programming with quadratic objective function where traditional optimization methods fail to find the optimal solution, efficiently.We propose a new hybrid solution approach based on an improved PSO (IPSO) and a simulated annealing (SA) procedure to optimize the cardinality constrained portfolio problem.Our proposed algorithm takes advantage of inertia weights mechanism and constriction factor approach in updating the particle's velocity to improve the PSO searching capabilities.In addition, an SA procedure is embedded into IPSO to improve the capacity of fine-tuning solution in the latter period of the search.Incorporating SA procedure can enhance the ability of IPSO for jumping out of the local optimum.Computational experiments on standard benchmark problems are carried out to assess the effectiveness of the IPSO-SA algorithm to solve cardinality constrained portfolio problem.
The rest of this paper is organized as follows.Section 2 describes the formulation of the portfolio optimization problem with cardinality constraints.In Section 3, an IPSO-SA algorithm is proposed to optimize the cardinality constrained portfolio problem.Section 4 is devoted to computational experiments and involves parameter setting for the proposed algorithm, numerical results and statistical performance analysis on benchmark data sets.Finally, Section 5 discusses the concluding remarks.

Cardinality constrained portfolio optimization model
The well-known mean-variance Markowitz (1954) model which played an important role in the development of modern portfolio problems can be expressed as follows, where N is the number of distinctive assets, i x is the proportion of asset i in the portfolio, ij σ is the covariance between the returns of assets i and j, i μ denotes the mean return of asset i and λ is the risk aversion parameter.When λ is zero, the model maximizes the mean return of the portfolio, neglecting the risk.In contrast, when λ is 1, the model minimizes the risk of the portfolio regardless of the mean return.Chang et al. (2000) have extended the Markowitz model by adding two constraints: floor-ceiling and cardinality constraints.The cardinality constrained portfolio optimization model can be expressed as follows, where K is a given value which restricts the number of assets in the portfolio.The decision variable is 1, if asset i is assigned to portfolio, otherwise y i =0.Parameters ε i and δ i are respectively the floor and ceiling bounds for asset's proportion.The cardinality constrained mean-variance model is a mixed integer nonlinear program, which cannot be solved efficiently by the existing algorithms due to its NP-Hard nature (Fernandez & Gomez, 2007).In the following section, we propose a new hybrid optimization strategy to solve medium and large instances of the cardinality constrained portfolio optimization problem encountered in realworld applications.

The IPSO-SA algorithm for solving cardinality constrained portfolio optimization
Particle swarm optimization (PSO) which was originally developed by Kennedy andEberhart (1995a, 1995b) is a population based stochastic optimization technique imitating the social behaviors of bird flocking.PSO has many similarities to evolutionary computation techniques like GA; however, the typical evolution operators such as selection, crossover and mutation are not involved.PSO has been widely applied in solving combinatorial optimization problems due to its exploration capability, implementation simplicity and performance reliability (Chen et al, 2006).In PSO, individual particles move thorough the problem search space seeking an optimal solution.The new position of each particle is adjusted according to its velocity and the distance between its current position and the local best position as well as the entire swarm best known position.As the algorithm is iterated, the swarm focuses more and more on an area of the search space containing high-quality solutions (Blum & Li, 2008).The combination of global exploration and local exploitation features during the optimization process is the key success of PSO algorithm.However, the main pitfall of the standard PSO is that it may be trapped into a sub-optimal solution region.Hence, many authors have extended the standard PSO and proposed efficient variants.For example, Shi and Eberhart (1998) proposed the inertia weights approach to balance the global and local search abilities in standard PSO.Inertia weights PSO investigates the space globally in early steps and searches around the optimal solution locally in final steps.Clerc and Kennedy (2002) developed a PSO variant by using a constriction factor in velocity update equation, and improved the convergence of the algorithm.Koshino et al. (2007) incorporated the inertia weights and constriction factor approaches and rose to higher-level performance for PSO algorithm.Pram et al. (2003) presented the fitness-distance-ratio based PSO, with near neighbor interactions.This algorithm considers another particle, nbest, in velocity update, which has a higher fitness value and is nearer to the updated particle.Some authors have hybridized PSO algorithms adding features of other meta-heuristic approaches and achieved more robust and effective methods for optimizing hard problems.Examples are Tasgetiren et al. (2004), Liu et al. (2007), Shen et al. (2008), Behnamian andFatemi Ghomi (2010), andTang et al. (2010).In this paper, a new hybrid PSO algorithm is proposed for the solution of mean-variance cardinality constrained problem.The proposed algorithm is inspired by recent improved variant of PSO, called IPSO, which integrates merits of inertia weights and constriction factor approaches (IWCFA).A modified simulated annealing approach with strong local-search power is incorporated to further improvement.The proposed IPSO-SA approach makes full utilization of the exploration capabilities in IPSO and SA algorithms and offsets the weaknesses of each algorithm.

Encoding scheme
A particle is represented by a 2×N dimensions matrix in which the first row involves proportion variables x pi and the second row includes binary decision variables y pi .N reflects the number of assets in the portfolio.The initial position and the velocity populations are generated randomly containing P particles.

Constraint handling and fitness evaluation
To ensure the feasibility of a candidate solution, we employ a constraint handling approach inspired from Cura (2009).When the number of assets is fewer than K, we select the maximum c-valued asset and conclude this asset to portfolio with a probability β; or an asset is chosen randomly with (1-β) probability.On the contrary, when the number of assets is more than K, we remove the minimum cvalue asset or a random one from the portfolio.C-value provides the proportion between mean return and mean risk for a given asset, with respect to parameter λ; and can be computed as follows, We used β = 0.5 for asset reposition probability.
Besides, in order to meet constraint 5, we can easily reposition pi x with x x 1 and to guarantee the boundary constraints on x i , we provide a substitute for x i such that and * δ are the sum of d i , for Q i ∈ ( i.e.where y i is equal to 1) and d i > 0,

and *
ε are the sum of e i , where , where Q i ∈ and d i < 0 and finally γ is the sum of (−1× e i ), where Q i ∈ and e i < 0. To assess how much a solution fits to the optimization purpose, a fitness function is evaluated.The fitness value of each particle p is defined as follows

Simulated annealing procedure
Simulated annealing (SA) is a trajectory-based heuristic, which performs a stochastic neighborhood search of the solution space.SA's major advantage over classical local search methods is its ability to avoid getting trapped into local optimum while searching for a global optimum.SA exploits an analogy between the way in which a metal cools and freezes into a minimum energy crystalline structure (the annealing process) and the search for a minimum in a more general system.Kirkpatrick et al. (1983) and Cerny (1985) pioneered the use of this algorithm for combinatorial problems.
Starting from a current solution p 1 , SA generates another solution p 2 by taking a stochastic step in some neighborhood of p 1 .If the new solution improves the value of the objective function, it replaces p 1 as the new current solution.Otherwise, the new solution is accepted with a given probability.The possibility of moving to solutions with a higher cost (i.e., performing degrading moves) characterizes SA and enables it to escape from local optimum.The probabilistic acceptance of the worst solution depends on the cost difference between the two solutions and it also decreases during the search.A Boltzmann function, Eq. ( 16), inspired from thermodynamics models often defines this probability, and T n is the temperature at stage n.The temperature is kept unchanged during each stage, which consists of a constant number of iterations.After each stage, the temperature is decreased by a constant factor ) 1 , 0 ( ∈ α .Simulated annealing algorithms differ from each other with respect to the various factors such as neighborhood search, cooling (annealing) schedule and termination criterion.After evaluating the fitness value of each particle, we sort the population based on the objective function values in an ascending order to choose the global best solution of the swarm.We apply SA algorithm as a local search mechanism around this global best solution (Gbest).Thus, Gbest changes in each iteration.The other critical characteristic of SA, which determines the search quality of solution space, is the definition of neighborhood structure.Here we considered generating neighbors by modifying the weights of a subset of the assets of the current portfolio.Two neighborhood structure examined by modifying the weight of only one asset: 1. Choose randomly one asset, which is currently in portfolio.The quantity of a chosen asset is either increased or decreased by a factor ) 1 , 0 ( ∈ q and idR is also neighborhood structure (Schaerf, 2002).2. Choose one asset randomly.If it is currently in portfolio, set its quantity to zero.Otherwise, add it to portfolio and set its quantity, randomly.Applying both of the neighborhood structure to our portfolio optimization problem, the second one has better performance in reaching better solutions and escaping from local optima.The other characteristic of SA algorithm is the feasibility of the final solution.There are two different approaches about the feasibility of the produced final solution.In the first approach, throughout all iterations of the SA algorithm, feasibility must be preserved and any solution violating the constraints will not be considered.Therefore, the neighborhood of a current solution must entirely consist of feasible solutions.On the contrary, to the first approach, the second one permits the consideration of infeasible solutions but adds a penalty term to the objective function for each violated constraint.Both approaches, "all-feasible" vs. "penalty" are not equally convenient in all situations.Since penalty approach does not seem to be appropriate for cardinality constrained portfolio optimization (Crama, Y. & Schyns, 2003), we use the so called all-feasible approach.In order to satisfy the constraints at any step of the search process, the constraint handling mechanism must be performed on the produced neighbor solution at every stage of SA.The initial temperature of the proposed SA is derived from the objective function value of the initial starting solution ( 1000 / t Gbest ).2N iterations are performed at each temperature.The proposed SA procedure pseudo code is displayed in Fig. 1.

Particels' velocity and position update
Particle swarm contains two primary operators: velocity and position update.In every iteration, each particle is updated to find a better position using its own memory as well as the swarm memory of the best-found solution.The first is the best solution it has achieved so far called Pbest.The latter is the best value obtained so far by any particles in the swarm called Gbest.Considering t Pbest and , t Gbest respectively as the particle's best and the swarm's best found solutions up to iteration t, the original PSO updates the velocity and position of each particle p through the following equations, ), ( ) ( where r 1 and r 2 are two random numbers uniformly distributed in the interval [0,1].c 1 and c 2 denotes the cognitive and social scaling parameters and w t is the particle inertia weight.In this paper, each particle is moved by an improved way proposed by Koshino et al, (2007) which takes advantage of both inertia weights and constriction factor approaches.According to inertia weights approach, the iteration number changes w t , gradually.In this way, the problem space is investigated globally in the early steps, and locally near the optimal solution in the final steps.Constriction factor approach uses coefficients to the present velocity, the direction of the global best position, and the direction of the individual-best, in order to control the behavior of the swarm.Applying this improved approach, the velocity of each particle is updated as follows, )) where w max and w min are the maximum and the minimum inertia weights, respectively.t max denotes the maximum number of iterations.are the best position of the swarm on dimension x i and y i .Thus, the position of each particle on dimension x i can be updated using the Eq.23.
Since t pi y is a binary variable, a discrete variant of PSO will be utilized to update the position of the particles on dimension y i , as follows.
) 1 where , and parameter ϖ is a small number.
After moving each particle, the feasibility of the new solution must be ensured through the constraint handling process described before and the fitness of new population is evaluated.This procedure continues until the iteration counter reaches the maximum number of iterations t max .Fig. 2 depicts the entire pseudo code of the IPSO-SA algorithm proposed for the solution of cardinality constrained portfolio optimization problem.

Computational experiments
In this section, computational results of the IPSO-SA algorithm to find cardinality constrained efficient frontier are discussed.The parameter setting to improve the performance of the proposed algorithm with statistical analysis is implemented.Then, the performance of fine-tuned algorithm is compared with previously proposed algorithms and its positive significant effect on results is proved through ANOVA and a multiple comparison test.To compare the results of the IPSO-SA with previously proposed algorithms, we used the benchmark instances from OR-Library (Beasley, 1996) available at http://people.brunel.ac.uk/~mastjjb/jeb/orlib/files.The data corresponds to weekly prices between March 1992 and September 1997 from different well known indices of Hang Seng in Hong Kong, DAX 100 in Germany, FTSE 100 in UK, S&P 100 in USA and Nikkei 225 in Japan.The numbers of different assets for each benchmark instance are 31, 85, 89, 98 and 225, respectively.All the results have been computed using the values K=10, i ε =0.01 and i δ =1 (i=1,…,N) for the problem formulation, and λ ∇ =0.02 for the implementation of the algorithms.With λ ∇ =0.02, the number of different λ values, denoted by ξ , is 51 and the number of particles is set to 50.4.1 Parameter setting for the proposed algorithm based on statistical analysis Each meta-heuristic has parameters which affect its solution quality and computation time.In order to improve the efficiency of an algorithm, these parameters should be tuned and set to their best values.A well-known approach to study the effect of different factors and their interactions is design of experiments (DOE) and its associated statistical tool, ANOVA.Analysis of variance (ANOVA) is a method to evaluate whether the impact of different factors on performance is significant or not.Since considering all parameters of the proposed algorithm and their combinations increases the dimension of the experiment in ANOVA, we choose to study four main parameters for tuning the algorithm.Table1 shows these parameters and their study levels.Note that cognitive and social parameters of IPSO part (c 1 and c 2 respectively) are set to their default value 2 which perform best for our algorithm.3 indicate that w max =1.5 and α =0.9 or 0.95 present the lowest error measure and subsequently the highest performance.To examine the impact of parameters on computation time of the algorithm, another ANOVA experiment is conducted with parameters in the previous study.The results of factorial ANOVA followed by SNK test suggest 2-opt as the best level for neighborhood search type to improve the efficiency of the algorithm in computation time.An analysis of the effects of these parameters on the other two error measures, variance of return error and mean return error, with the same approach determined best value for w min at 0.2 or 0.1 level.

Numerical results
Computational results of IPSO-SA algorithm to find cardinality constrained efficient frontier are presented in this section.To examine the efficiency of the proposed algorithm, the performance of the IPSO-SA has been compared with GA, TS, SA and standard PSO using three defined error measures for five benchmark data sets.

Analysis of results
In this section, we statistically compare the performance of our proposed algorithm with GA, TS, SA and standard PSO using analysis of variance.To test the equality of the performance for different meta-heuristic algorithms, a randomized complete block experiment is designed.In this experiment, the effect of benchmark problems is blocked to reduce experimental error and to provide more reliable results (Coffin & Saltzman, 2000).
The ANOVA for randomized complete block design is performed with five meta-heuristic algorithms as treatments, five benchmark problems as blocks and one observation per cell.The results of ANOVA comparative study of the algorithms are presented in Table 5.In this analysis, the We used SNK test for multiple comparisons of mean responses to determine the most efficient algorithms with respect to error measures.SNK grouping revealed that our proposed IPSO-SA algorithm (group A) performs better than GA, TS, SA and standard PSO (group B).

Concluding remarks
This paper studied a realistic portfolio optimization problem with floor, ceiling, and cardinality constraints, which plays an important role in financial engineering.A new IPSO-SA algorithm based on an improved particle swarm optimization and a modified simulated annealing was proposed to find the cardinality constrained efficient frontier.The proposed algorithm benefits from simple and easy characteristics of PSO with an adaptation of inertia weights and constriction factor approaches.
In addition, incorporating SA procedure into IPSO helps escape from local optimum and improve the precision of the convergence.
Computational experiments on five benchmark data sets for up to 225 assets were executed to assess the effectiveness of the IPSO-SA algorithm to solve the cardinality constrained portfolio problem.A parameter setting method with statistical analysis was employed to support the algorithm.Then, the performance of fine-tuned algorithm was compared with previously proposed algorithms and its positive significant effect on results was proved through ANOVA and a multiple comparison test.
The results indicate that our proposed algorithm can find high quality solutions with near-zero errors in reasonable amount of time.

Fig. 2 .
Fig. 2. IPSO-SA algorithm for the solution of cardinality constrained portfolio optimization problem In order to test the effectiveness of the proposed algorithm, we compute the heuristic efficient frontier and compare it with standard efficient frontiers.Three error measures are used for the evaluation of the heuristic efficient frontier relative to the standard efficient frontier: mean Euclidean distance, variance of return error and mean return error.Error measures can be quantified as follows: Let ( s i s i r v , ) (i=1, …, 2000) be the variance and the mean return of the point in the standard efficient frontier, and let ( h j h j r v , ) (j=1,…, ξ ) be the variance and the mean return of the point in the heuristic efficient frontier.Let (

Fig. 3 .
Fig. 3. IPSO-SA efficient frontier vs. the standard efficient frontier far for particle p on dimensions x i and y i , respectively.
t pi vx reflects the velocity of particle p on dimension x i , and t pi vy indicates the velocity of particle p on dimension y i at iteration t.

Table 1
Experimental design parameters general factorial experiment with four factors, each having three levels and two observations per combination is designed.The ANOVA and multiple comparison tests were performed with SAS software version 9.1 for Hang Seng benchmark data set.Similar results are expected for other benchmark instances.Table2shows the results of ANOVA test for parameter tuning based on mean Euclidean distance error values.These results indicate the significant effects of w max and α factor on the performance of the proposed model.In order to determine the statistically similar and different means, the Student-Newman-Keuls (SNK) test is used to group them.The results of SNK test for significant parameters in Table A

Table 3
Results from SNK test for Hang Seng data set Table4shows the comparative results, and Fig.3illustrates the comparison of efficient frontiers.It is worth to mention that we implement IPSO-SA algorithm in MATLAB 7.6.0using a core 2 Dou, 2.1 GHz computer with 2 GB of memory.The execution time of IPSO-SA was 183 seconds for Hang Seng, 563 seconds for DAX 100, 767 seconds for FTSE 100, 848 seconds for S&P 100 and about 72 minutes for Nikkei.

Table 4
Experimental results of five heuristics on portfolio selection

Mean return Variance of return Standard Efficient Frontier IPSO-SA Efficient Frontier Mean return Variance of return Standard Efficient Frontier IPSO-SA Efficient Frontier Variance of return Mean return Standard Efficient Frontier IPSO-SA Efficient Frontier Variance of return Mean return Standard Efficient Frontier IPSO-SA Efficient Frontier Mean return Variance of return Standard Efficient Frontier IPSO-SA Efficient Frontier performance
of the algorithms are measured with mean Euclidean distance error.According to ANOVA test, there is a significant difference among the performance of different algorithms.