A Novel Particle Swarm Optimizer and Its Application to the Yield Curve Estimation Problem

Particle swarm optimization (PSO) has been considered as one of the main swarm intelligence algorithms for solving single-objective optimization problem. How to update the velocities and positions of a swarm of particles is key to its optimization performance. In this paper, we propose to use the moving average of the local best positions visited so far as a key information to update the velocities and positions. Further, a central learning strategy is proposed in which the center of the local best positions is computed and used to update the global best position. Combining these two strategies with the updating formulas which are inspired by the free electron model in metal conductors placed in an external electric field, we name the proposed PSO algorithm as PSO with moving average and central learning strategies (dubbed as MAPSO). We test the performance of MAPSO on the CEC 2017 test problems with 10 and 30 dimensions. Experimental results show that MAPSO significantly outperforms some well-known PSOs in general. We then apply MAPSO on a very challenging real-world problem, i.e. the yield curve estimation problem, in macroeconomics. Our experimental study on the yield curve estimation problem with Shanghai interbank offered rates shows that MAPSO can effectively solve the problem and achieve the state-of-the-art performance.


I. INTRODUCTION
Swarm intelligence is one of the main streams of evolutionary algorithm (EA), along with genetic algorithm (GA), differential evolution (DE), evolutionary strategies (ES), evolutionary programming (EP) and many others. As opposed to other EAs which are established based on the survival of the fittest principle, swarm intelligence algorithms mimic cooperative communication and information exchange among individuals within a swarm of particles (e.g. a flock of fish or a swarm of ants) [1]. It is one of the main branches of population-based stochastic optimization paradigms for solving complex optimization problem.
As a representative of swarm intelligence algorithms, particle swarm optimization (PSO) was proposed by Kennedy and Eberhart in 1995 [1]. Since then it has been widely applied The associate editor coordinating the review of this manuscript and approving it for publication was Yu-Da Lin .
to optimize complex optimization problems and train neural networks. PSO is inspired by the hunting behavior of a flock of birds, where each bird is considered as a particle. In PSO, each particle has its velocity and position. PSO mimics the bird's group search pattern through cooperative interaction among particles, and evolves the particles to search for the optimal solution of an optimization problem. Different from other EAs, such as GA, DE, ES and EP, etc., which employ crossover and/or mutation operators, PSO determines the search direction of each particle through position sharing within the swarm. Each particle updates its own position according to the positions of other particles. With such a search scheme, a global search ability can be realized [2], [3].
PSO is easy to understand and implement, and does not require a complicated tuning procedure for its control parameters in practice [4]. The main component that affects the performance of PSO lies on its updating mechanism for the velocities and positions of the particles, and the use of VOLUME 10, 2022 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ neighborhood topology for information sharing within the swarm. PSO suffers from fast convergence, lack of exploration capability, decline of solutions variety, and easily trapping in local optimum and other disadvantages. A great number of PSO variants have been developed to address these issues [2]. The development of PSO can be roughly categorized into the following research avenues. 1 One of the popular research avenues in developing efficient PSO is the introduction of the so-called 'inertia weight', which was firstly proposed by Shi and Eberhart [5]. The idea is to equip each particle with an inertia weight on its current velocity. In [6], Shi and Eberhart proposed to decrease linearly the inertia weight along with the evolution. This idea has been shown empirically to be helpful for accelerating the convergence speed of PSO. Yang et al. [7] proposed to change the inertia weight non-linearly according to the positions of the particles, which can be seen as an adaptive strategy. Bansal et al. [8] compared various inertia weight adaption strategies and presented some guidelines. Different from the original PSO in which each particle equips with one inertia weight, Taherkhani et al. [9] proposed to use different weights on different dimensions. Zhang and Li proposed an adaptive inertia weight PSO in [10], in which the speed change and momentum of each dimension of the particle are set adaptively. An adaptive inertia weight technique has been developed in [11] based on the concept of success rate usually adopted in DE [12]. Recently, researchers have developed methods to adapt the weights based on sigmoid function [13] or exponential function [14].
Secondly, the updating method for the velocities and positions of the particles has been studied extensively. In the original PSO, the updating of the velocities involves three parts, including the momentum, the cognitive and the social part. The cognitive part reflects the learning from a particle's own search experience, while the social part represents its collaboration with other particles. In the original PSO, the cognitive part is associated with the best position a particle has found until current iteration (called 'local best' or lbest), while the social part involves the best position found by all the particles (called 'global best' or gbest).
In the PSOs proposed in [15] and [16], the particle position is determined by sampling from a Gaussian distribution based on the positions of lbest and gbest. Some researchers proposed to use Levy distribution which has empirically shown better algorithmic performance than using Gaussian [17], [18]. Liang et al. [19] proposed a learningbased PSO, in which the particles to be learned are selected by setting the learning rate. This learning strategy can increase the diversity of the particles and avoid the premature convergence. Sun et al. [20] proposed to update the velocity by taking the social part as the addition of a random motion term and a drift motion term. Their development is inspired by the movement model of electrons in metal conductors. 1 Since the study on PSO is very popular, this paper does not intend to review all the literatures.
Cheng et al. [21] proposed a PSO based on human's social behavior. Tanweer et al. [22] proposed a strategy to improve the particle's velocity, in which the best particle and the other particles at the current iteration adopt different update rules. The strategy aims to speed up the exploration ability of the best particle and improve the convergence speed of other particles. Karim et al. [23] proposed to incorporate an optimal guide creation module to guide the particles escaping from local or non-optimal solution regions.
An improved PSO inspired by social learning was developed in [24] where the three best particles in the current iteration are updated using a differential mutation strategy. Li et al. [25] proposed to mix some mutation strategies to update the velocity based on the idea of dividing the total population into an elitist population and a general population. Sedighizadeh et al. proposed to add two new terms in the original PSO updating formula [26]. Xu et al. [27] proposed to use reinforcement learning to enhance the usefulness of the communication topology in PSO.
Thirdly, the multi-population strategy has been studied for improving PSO's exploration capability. In [25], the authors proposed a multi-population method with mixed mutation strategies, in which the whole population is divided into some subpopulations and each subpopulation evolves cooperatively. Peng and Chen [28] proposed a cooperative PSO with multiple strategies, in which each particle represents a fuzzy regulation and each particle in each subpopulation evolves separately to avoid trapping in local optimums.
Tanweer et al. [29] proposed a PSO with dynamic guidance and self-adjustment, in which the particles are divided into three categories according to the Euclidean distances from the particles to the global optimum. In [30], Ye et al. proposed a new multi-swarm PSO algorithm based on a dynamic learning strategy, in which the particles in each subpopulation are divided into common particles and corresponding particles. At each iteration, the common particles mainly pay attention to exploration, while the corresponding particles intend to perform exploitation. Chen et al. [31] proposed to dynamically group particles during the search procedure.
The fourth avenue is to combine PSO with evolutionary operators, such as mutation and crossover, adopted in other EAs, expecting to make use of their advantages. In [32], Angeline proposed to adopt the tournament selection method used in GA to improve the performance of PSO. Lvbjerg et al. proposed a hybrid PSO algorithm which combines the idea of particle swarm with the concept of GA [33].
To improve the diversity of particle swarms in PSO, crossover operators have been used to generate individuals with better quality. In [34], Meng et al. proposed a cross-search particle swarm optimization (CSPSO) algorithm. Unlike PSO and its variants, the particles in CSPSO are sequentially updated in each generation through the original PSO update strategy and the crossover operation. In [35], Ratnaweera et al. proposed a novel parameter automation strategy and two extended methods to update the control parameters in PSO. A time-varying acceleration coefficient is then proposed for adapting the inertia weights, which is then combined with the mutation operator in GA.
In PSO, the particles cooperate and share information with each other through a specific structure, known as neighborhood topology. The neighborhood topology structure could greatly affect the performance of PSO since it determines how the particles communicate and how the useful information spreads within the swarm. In a word, the neighborhood topology could influence the particles' behavior greatly. The original PSO adopts a star topology where each particle considers all the other particles as neighbors. The ring topology, on the other hand, takes its each particle's two immediate neighbors to form a circle [36]. The Von Neumann topology can be considered as a rectangle matrix which determines which particles are connected [37]. In dynamic topology, the neighborhood is reassigned after a certain number of iterations [38].
Kennedy and Mendes [36] tested the ''pyramid'' triangle topology based on 3D wireframe, which is a star topology with one central node and the community structure of Von Neumann. Bonyadi and Michalewicz [3] established a time adaptive topology which makes the variant of PSO locate many feasible regions in the early stage and lock the most promising regions in the later stage of the optimization process. Zeng et al. [39] proposed a dynamic neighborhood scheduling strategy to adaptively determine the neighborhood topology of particles, which is realized by monitoring the aggregation factor of particles. In [40], Zhang et al. proposed a hybrid topology, which combines the closed-loop topology with the star topology or quadrilateral topology.
Besides the aforementioned research studies, there are many other research avenues, including but not limited to, hybridizing with other evolutionary algorithms (e.g. with DE [41], whale optimization [42], quantum [43], and others), swarm initialization [44], population size controlling [45], etc. Further, PSO has been also widely studied and applied to a great number of other benchmarking problems such as constrained optimization problems [46], large-scale problems [47], and combinatorial optimization problems [48]; and real-world optimization problems such as the position of ocean-moored buoys [49], the process synthesis and design problem [50], the scale-free network problem [51], the image segmentation problem [52], and many others. Interested reader please refer to [2] for a comprehensive review.
In this paper, we propose a novel PSO focusing on the development of new updating method. The main contributions of the paper can be summarized as follows: • A new updating strategy is proposed, in which the moving average of the local best positions with correction is rigorously analyzed and applied to update the cognitive part of the velocity.
• A new central learning strategy is proposed to update the global best position; • The proposed PSO is compared with a number of PSOs on the commonly-used CEC 2017 test suite. Experimental results show that the proposed algorithm significantly outperforms the compared algorithms in general.
• The proposed algorithm is applied to solve a very challenging real-world non-convex optimization problem, namely the yield curve estimation problem. Experimental results show that the proposed algorithm significantly outperforms a variety of algorithms including some classical optimization methods and EAs. The rest of the paper is organized as follows. Section II presents the proposed strategies. Section III presents the yield curve structure estimation problem. Experimental results on the CEC 2017 test suite are presented in Section IV, while the application to the yield curve estimation problem is presented in Section V. Section VI concludes the paper.

II. PARTICLE SWARM OPTIMIZER BASED ON MOVING AVERAGE AND CENTRAL LEARNING STRATEGY A. BASIC PSO
In this paper, we consider box-constrained global optimization problem: where i and u i represents the lower and upper bound of each decision variable x i , respectively.
PSO maintains a population of particles and searches the decision space by updating the positions and velocities of the particles according to some information stochastically. In PSO, at the t-th iteration, each particle is composed of a position vector x t i = x t i1 , · · · , x t id T (which can be considered as a solution to the optimization problem) and a velocity In the original PSO, let y t i be the best position found by the i-th particle (i.e. lbest), and g t be the best position (i.e. gbest) found by all the particles, until the t-th iteration. The update equations for the position and velocity are as follows: where c 1 and c 2 are the cognitive and social acceleration coefficients, r t 1 and r t 2 are two random numbers sampled within [0, 1] from the uniform distribution U(0, 1).
The first term in the right hand side of Eq. (1) is the momentum part, while the second and third part is called the cognitive and social part, respectively. The cognitive part reflects the influence of the search experience of an individual particle, while the influence of the rest particles comes from the social part. For a minimization problem, the update of y t i is as follows: It should be pointed out that PSO differs from DE mainly on the following aspects. First, in DE, each individual is not VOLUME 10, 2022 considered to be a particle. There is no velocity associated with each individual in DE. Second, the updating of each individual, usually called the mutation operator, in DE is quite different from PSO. An example DE mutation operator, named rand/1, is as follows: where F is a parameter called scaling factor, r 1 , r 2 and r 3 are random integers sampled from [1, N ]. Third, there is no crossover operator in PSO.

B. THE MOVING AVERAGE OF INDIVIDUAL PARTICLE's BEST POSITIONS
In classical PSO, the velocity update formula only uses the particle's local best position and the global best position at the current iteration. In other words, each particle does not consider the positions it has visited before, which means that the particles in the present PSOs have no memory. However, we argue that such a update pattern has its limitation. Consider a situation when a particle reaches a position near the global optimum at the t-th iteration, but falls into a local optimum with a smaller objective function value at the (t +1)- x t i will be forgotten in the following iterations in the classical PSO, which means that such useful information is lost. As a result, the particle may converge to the local optimum. This situation often occurs in the search process of PSO. Moreover, since the local best position of a particle after an iteration will be updated by a better position, the original optimal position will be lost. However, the lost position may lead to a more promising area. If such a situation happens, the particle swarm will converge prematurely. As a result, the PSO's global exploration ability will be greatly affected.
In consideration of this shortcoming, we propose to maintain all the local best positions that each particle has visited during the search process. In the sequel, we use P 1:t to denote all the local best positions until the t-th iteration. Assume that the population size is N , P 1:t can then be presented as follows: where each row represents all the local best positions that a particle has visited until the t-th iteration. The global best position of all the local best positions, i.e. g t , can be derived from P 1:t as where P 1:t [i] represents the i-th column of P 1:t .
To better make use of the historical information of the particles, we propose the following formula for each particle.
That is, for each particle i, let It is seen from the equation thatỹ t i is the weighted average over time on y t i and γ ∈ [0, 1] is the weight decay rate. Let's only consider one particle's historical information to better understand the effectiveness of the above equation. Expanding Eq. (5), we havẽ From the above equation, it can be found thatỹ t i is the weighted average of all the local best positions until t. The closer to the current iteration t, the larger the weight value for the corresponding local best position.
Since the weights are exponentially decreasing along with the increasing of t and it is unwise to keep all the visited local best positions, we therefore propose to truncate the memory in case γ t < 1/e where e is the mathematical constant. Any local best positions beyond the iteration bigger than −1/log γ will be ignored. Since 1/e is the limitation of we can determine the number of iteration for truncation depending on the choice of γ . Based on the above analysis, we can regardỹ t i as its weighted average value within 1/(1 − γ ) iterations. Setting a large γ can weaken the influence of the current local best position on the entire particle swarm as much as possible. For example, let γ = 0.9,ỹ t i can be viewed as the weighted average value among the local best positions in ten consecutive iterations, i.e. y t i , y t−1 i , · · · , y t−9 i . If γ is set smaller, e.g. γ = 0.5,ỹ t i can be considered as the weighted average value of y t i and y t−1 i . That is, it is only influenced by the local best positions of two previous immediate iterations. It is seen that setting a small γ cannot compensate for the loss of valuable information. Therefore, it is suggested to set γ be close to 1.
It is worth noting that if we set the initial value of y 0 i to be 0, then y 1 i = (1 − γ )y 0 i . Since γ is nearly to 1,ỹ i 1 is not an accurate estimation for y t i . This is to say that the local best positions of the particle are completely discarded. At the beginning of the PSO iteration, it will cause the distortion of the estimated local best position. To address this problem, we propose to correct it to maintain the accuracy of the estimation by dividingỹ t i with (1 − γ t ). That is, we set To better explain the bias correction, we consider the relationship betweenỹ t i and y t i . Note Eq. (6) can be simplified as If we takeỹ t i as a random variable, its expectation can be obtained as follows: This equation indicates thatỹ t i is not an unbiased estimation of y t i . However, the approximation becomes unbiased if t increases. On the other hand, the corrected history informationỹ t i can ensure that the estimated local best position used in the current iteration is an unbiased estimate to the local best position used in Eq. (8).
Through using the corrected local best information, the referenced positions of particles found at each iteration are recorded and new positions can be generated by using the weighted average method. The method avoids particle updating by relying too much on the local best position of the current step, increases the chance that the particles jump out of the local optimum, and reduces the chance of premature convergence to a certain extent.

C. THE CENTRAL LEARNING STRATEGY
The motivation of the proposed central learning strategy comes from the observation of the swarm's search behavior.
To demonstrate the trajectory change of the local best position y t i , we show the optimization process of the classical PSO on the two-dimensional Rastrigrin function in Fig. 1. From Fig. 1(a) to (d), the local best position of each particle is shown at the 0-th, 10th, 40th, and 80th iteration, respectively. In the figure, the black circle represents the local best particle, the blue diamond stands for the weighted center of these particles, and the red asterisk represents the global best position. We can find that the distribution of the local best position changes constantly as the search goes. However, the center of the local best (i.e. the blue diamond) is not far from the global optimum (the red asterisk), which indicates that the center of the local best positions reflects some sort of global information.
Further, we show an example search process when optimizing the 30-dimensional Rastrigrin function in Fig. 2, in which the distances between the global best position and the centre of the local best positions to the global optimum solution are plotted against the iteration. From the figure, we can see that during the search process, the center of the local best positions is closer to the global optimum than that of the global best position, but they gradually converge to the same position as the search iterates. Therefore, it may be useful to take the center of the local best positions as an essential information for the velocity updating.
Based on the above observation, we propose the following central updating strategy. That is, a particle is set to explicitly visit the center of the local best position at each iteration. Particularly, at the t-th iteration, we set the N -th particle as the central particle: Different from other particles which are updated according to the proposed strategy, the central particle will be updated independently. It will be used only to update the global best position, but not to be used for the updating of other particles' velocities.

D. PARTICLE SWARM OPTIMIZER BASED ON THE MOVING AVERAGE AND CENTRAL LEARNING STRATEGY 1) THE FREE ELECTRON MODEL
Besides aforementioned strategies, our algorithm is built upon RDPSO [20] which is inspired by the free electron model for conductors that are placed on external electric VOLUME 10, 2022 field [53]. According to this model, the movement of an electron is the addition of the thermal and the drift motions. The thermal or random motion always exists, while the drift motion makes the electron particle move towards the opposite direction of the external electric field. The electron's movement will lead it to the position of the minimum potential energy in the conductor.
Since there may exist many local minimum potential energies, the drift motion itself may drive the electron to local minimum, which is analogical to the situation that an optimizer converges to a local minimum. On the other hand, the random motion can help the electron escape a local minimum potential energy. Therefore, the movement of the electron can be considered to be a process of minimizing its potential energy. Particularly, the thermal motion guides the global search of the particle, while the drift motion realizes the local search.

2) THE PROPOSED ALGORITHM
The convergence analysis [20] proved that the i-th particle converges to its local attractor p t i = (p t i1 , p t i2 , · · · , p t id ) which can be written as follows: T is a random number vector and each element is uniformly sampled in [0, 1], and means element-wise multiplication.
Moving towards the local attractor p t i reflects the local search ability of PSO. It is seen that only the information at the current iteration is taken into consideration in Eq. (12). From the discussion in Section II-B, we know that using previous information, i.e. the moving average of the local best positions of each particle, can help the particles escape from local optimum. Therefore, we propose to update the local attractor for the first N − 1 particles as follows: for i = 1 to N − 1, define whereỹ t i is defined as in Eq. (5) and Eq. (8). The drift motion of a particle can be considered as the movement towards the local attractor, which is analogical to the combination of the cognitive and social recognition parts. In the proposed algorithm, to analogy the electron's movement, we use the random motion to replace the monument part in the classical PSO. By combing all the motion terms together, including the random term and the drift term, the velocity of a particle i at the j-th dimension in the proposed algorithm can be written as follows: where vr t ij and vd t ij represent the random and drift term, respectively. In the sequel, we use to represent the vector form of the movement of particle i Further, we assume that vr t ij follows the Maxwell velocity distribution law [54], which means vr t ij follows a zero-mean normal distribution as follows: where σ t ij is the standard deviation of the distribution. vr t ij can be sampled as follows: where ψ t ij is a random number sampled from the standard Gaussian distribution. To make σ t ij adapt to the optimization process of the particle swarm, we propose to adjust it by using previously collected information.
where each component is defined as: The standard deviation σ t ij can then be updated as follows: where α is the thermal coefficient. Then, vr t+1 ij is updated as: It is expected that the drift velocity component vd t ij realizes the local search of particles. As mentioned above, this can be achieved by moving towards the local attractors. To realize this, we express vd t ij as follows: In summary, the updated formula for each particle i, 1 ≤ i ≤ N − 1 at the j-th dimension is written as follows: while the N -th particle is updated as shown in Eq. (11). By combining the aforementioned updating strategies, the proposed algorithm, named PSO with moving average and central learning strategies (MAPSO), can be summarized in Alg. 1. In the algorithm, first a swarm of particles are randomly initialized as X 0 (line 1). We then set the local best y 0 i , 1 ≤ i ≤ N − 1 to be the initialized particle (line 2) and compute the central position (line 4). The algorithm iterates from line 6 to 19 until the termination criterion has been met. At each iteration, the particles are updated from line 7 to 11. With the updated particles, the local best positions and the central position are computed respectively in line 12 and 13. The new particle swarm is then formed in line 15. The new corrected local best position is updated in Compute p t i and C t j by Eq. (12) and (18), respectively, for 1 ≤ i ≤ N − 1; 4: Compute the N -th particle x 0 N by Eq. (11); 5: Set t ← 0; 6: while t < MaxIter do 7: for i = 1 → N − 1 do 8: Compute v t i according to Eq. (15); 9: 10: Evaluate f (x t+1 i ); Compute p t i and C t j by Eq. (12) and (18), respectively, for 1 ≤ i ≤ N − 1; 18: t ← t + 1; 19: end while 20: Return x * = arg max f (X t ).
line 16 while the local attractor and central position is updated in line 17. After termination, the best solution found so far is returned (line 20).

III. THE YIELD CURVE ESTIMATION PROBLEM
In this paper, we intend to apply the proposed algorithm to a real-world problem, namely the yield curve estimation problem. In this section, we briefly introduce the problem's background and its optimization modeling.

A. BACKGROUND
The three main components of bank liabilities include customer deposits, interbank liabilities, and borrowings from the central bank. The primary source of bank liabilities is customer deposits, while interbank liabilities include the certificates of deposit and financial debts, both of which are important supplementary sources of liabilities for commercial banks, and they also are the source of financing with the highest degree of interbank marketization. As a certification for time deposit issued by depository financial institutions in the national interbank market, their investment and transaction entities are fund management companies, fund products, and members of the national interbank lending market. Within the issuance and filing quota in the current year, deposit-taking financial institutions may determine the terms and amounts of every interbank certificate, but the total amount issued cannot be less than a certain amount, e.g. RMB 50 million in China. In most cases, the time limit of the certificate is within a year.
Interbank certificates of deposit can be circulated in the interbank market before expiration, and the market quotation is determined by a group of banks. In China, eighteen commercial banks with higher credit ratings form a Shibor (Shanghai Interbank Offered Rate) quotation bank group and they determine the official quotations of different varieties of Shibor. Shibor has a total of eight varieties, including overnight, one week, two weeks, one month, three months, six months, nine months, and one year. Taking the official Shibor offer on May 14, 2021, as an example, its quotation form is shown in Fig. 3:

B. THE OPTIMIZATION MODEL FOR THE YIELD CURVE ESTIMATION PROBLEM
Macroeconomists and financial market practitioners are very much interested in modeling the term structure of government bonds and interbank offered rates. It is critical for bond and derivatives pricing, risk management, and revealing market expectations, which are critical for monetary policy decision-making.
A yield curve is a graphical depiction of the term structure of interest rates which is a function of time. It is a remarkably meaningful financial tool used in investment and policy decisions. Yield curve estimation is basically a nonlinear parameter optimization problem. Market expectations and conditions are reflected in the term structure of interest rates. Therefore, its estimate should be as accurate as possible to ensure that the information it describes reflects the true market conditions and prospects.
Various efficient structural models have been proposed to estimate the term-rate structure, such as the Nelson-Siegel (NS) model [55], the SV model [56] and so on. Among them, the NS model is one of the most popular term structure models in academia, practitioners, and central banks. The solution accuracy of these nonlinear models directly affects the level of returns in practical applications.
In this paper, we adopt the NS model to transform the yield curve estimation problem into an optimization problem.

1) THE NELSON-SIEGEL MODEL
The NS model was firstly proposed by Nelson and Siegel in 1987 [55]. It involves four parameters with economic meaning, and can describe the changes and states of the interest rate curve profoundly. In the following, we briefly introduce the NS model. First of all, the forward instantaneous interest rate function can be derived as follows: where θ = (β 0 , β 1 , β 2 , τ 1 ) is the model parameter. Since the spot rate and the forward rate have the following relationship: and r(0) = R(0), r(∞) = R(∞). The spot rate can then be calculated as: and lim t→0 In the estimated parameters, the slope factor of the parameter β 0 is 1. It does not decay over time, i.e., the effect of β 0 on all term rates is the same and its movement will cause the yield curve to move up and down horizontally. The slope factor of the parameter β 1 is 1−e − t τ 1 t τ 1 and τ 1 is varied from 1 to 10. Thus, This shows that the initial value of the slope factor of β 1 is 1 and it decays to 0 at a rather fast rate along the increase of time, which indicates that β 1 has a greater impact on short-term interest rate. It represents the slope factor and affects the slope of the interest rate curve. On the other hand, the slope factor of the parameter β 2 is 1−e − t which is a function that increases with time and then gradually decays to 0. This parameter has a weak influence on the short and long ends of the interest rate curve structure and has a greater impact on the mid-end interest rate. This character makes β 2 be very suitable for representing the slope of the interest rate curve. τ 1 is the decay degree of the slope factors of β 1 and β 2 . The smaller τ 1 , the faster of the decay rate of the slope factor and vice versa.

2) OPTIMIZATION MODELING
By optimizing the yield to maturity (YTM) in the past period, the four parameters of the NS model are updated so that the yield curve is expected to be close to the actual term interest rate structure of the secondary market. According to the updated NS model of the short-term arbitrage, the specific pricing process can be divided into the following three steps: • Applying the NS model to fit the interest rate term structure and obtain the initial values of the parameters of the NS model according to the quotation reference of Shibor's official secondary market circulation; • Using the transaction data of interbank deposit receipts in the secondary market on a certain day to build the objective function and optimize the parameters of the NS model; • The required YTM can be updated using an NS model with the estimated parameters to perform short-term arbitrage. The parameters of the NS model are obtained by fitting it to the official quotation of the interbank offered rate, while the quotation of the actual transaction certificate of deposit is not used. In other words, we assume that there is a certain difference between the official quotation and the actual transaction quotation. Our goal is to correct the four parameters of the NS model through the transaction data of the actual certificate of deposit. The objective function can then be written as follows: where ω i is the weight of a single interbank deposit receipt traded in the total market volume, pr * i is the actual price and pr i is the theoretical price, its discount function is calculated as follows, Substituting the above equations, we obtain the final optimization problem, which is shown as follows:

IV. EXPERIMENTAL RESULTS
In this section, we first compare the proposed algorithm with some well-known PSOs. Then the effectiveness of the proposed strategies is studied.
In the experiments, the single-objective functions in the 2017 IEEE Congress on Evolutionary Computing (CEC2017) are used as benchmarks. The function set contains a total of 30 test problems with different characteristics. In our experiment, F2 is excluded since the function is not computationally stable.

A. COMPARISON WITH WELL-KNOWN PSOs
The proposed algorithm, named PSO with moving-average and central learning strategies (MAPSO), is compared with several known PSO algorithms. The compared algorithms include • FIPS [36] which uses all the information in a particle's neighborhood to update itself based on different topologies; • DNLPSO [38] which proposes a dynamic neighborhood learning strategy to improve the swarm's diversity; • SRPSO [22] which uses different strategies to update the best particle and the rest particles' velocities; • RDPSO [20] which is inspired by electron's motions in conductor.
• EPSO [57] which is a recently proposed PSO with six strategies to address the shortcomings of the classical PSO. The parameters of these compared algorithms are listed in Table 1, which are the same as those set in the original references. To compare these algorithms fairly, the same experimental parameters are chosen. Specifically, the population size N = 40, MaxFES = 10, 000 × D, where D is the dimension of the test problem, the search interval is [−100, 100]. The termination is that the number of fitness evaluations is greater than MaxFES. The difference between the optimal value found by the algorithm and the actual optimal value at the end of the iteration is called the error value. The error is considered to be 0 when the error value is less than 10 −8 .
All the compared algorithms are run 50 times independently on each function and the particles' velocities are randomly initialized for each run. The means and standard deviation values of the error values are computed and shown in the following tables, where the smallest mean values are typeset in bold. The Wilcoxon's rank-sum test is used to reflect the statistically-significant differences between the compared algorithms with MAPSO at a significance level of 5%. In this paper, we use symbols '+', '−' and '≈' to represent that MAPSO performs worse than, better than, and similar to the compared algorithms as suggested by the rank sum hypothesis test. Tables 2 and 3 show the experimental results of the compared algorithms on the CEC2017 test problems with 10D and 30D, respectively. From these tables, we can find that for almost all the tested functions, the mean values of MAPSO are smaller than that of the compared algorithms, which   indicates that MAPSO achieves a better solution than the compared algorithms in general.
In Table 2, for the 10D problems, it is seen that MAPSO clearly outperforms the compared algorithms on F1. It achieves a 100% success rate on F3, F6 and F9. Further, it can obtain better performance than the compared algorithm on solving simple multimodal functions. In addition, it outperforms the compared algorithms on F12-F19, F22, F24-F26 and F28 which are difficult hybrid and composite functions. The hypothesis test suggests that MAPSO achieves significant improvement over the compared algorithm on at least half of the test functions.
In comparison with EPSO, we see from Table 2 and 3 that EPSO can perform well on the unimodal and multimodal functions (F1-F10), but for the hybrid and composite functions, its performance is clearly worse than MAPSO except few functions (F17, F18 and F21 in 10D, and F24 in 30D). Note that these functions are usually considered to be more challenging than F1-F10. This clearly indicates that MAPSO improves the exploration capability of PSO but not decreasing its exploitation capability.
In comparison with RDPSO, we see that MAPSO outperforms RDPSO significantly in general. The indicates that the proposed strategies can indeed be beneficial to the performance of MAPSO since MAPSO is built upon RDPSO. Since it is claimed in [20] that the exploration ability of RDPSO is limited, the better performance of MAPSO implies that the proposed strategies are able to improve its exploration ability which justifies that the proposed strategies are effective. We will carry out more detailed studies in the next subsection to justify the effectiveness of the proposed strategies.
To better demonstrate the performances of the compared algorithms, we apply the average performance score (APS) [58] to give an overall performance on the test function set. If there are totally n compared algorithms and M functions in the set, and algorithm A i (i ∈ {1, · · · , n}) outperforms algorithm A j (j ∈ {1, · · · , n} \ {i}) on the m-th function f m (m ∈ {1, · · · , M }) with a statistical significance, then set δ i,j = 0. Otherwise, δ i,j = 1. When comparing to all other algorithms, the performance score for A i on the test function f m equals the sum of δ i,j . That is, The APS value of A i on the entire function set is APS( It is seen that a smaller value of APS for A i means that the algorithm A i has a better overall performance.
The overall ranking on the test suite in terms of APS is represented in Table 4. It is seen from the table that MAPSO obtains the highest APS values both on 10D and 30D test functions. Thus we may conclude that MAPSO performs the best among the compared algorithms.
To make the comparison more clearer, we show the convergence properties of the compared algorithms in Fig. 4 and Fig. 5, respectively, on the CEC 2017 test functions with 10 and 30 dimensions. From the figure, it is seen that MAPSO displays an obvious advantage over other algorithms in both early and final iterations in general. As shown in Fig. 4 for F5, F8, F10, F16 and F30, MAPSO can gradually escape from a local optimum once trapped. This confirms its exploration capability. In Fig. 5, such an exploration capability can be easily seen on functions F5, F8, F10, and F21. This again confirms the improvement of the exploration capability of MAPSO against the classical PSO and RDPSO.

B. COMPONENT ANALYSIS
In this section, we investigate the effectiveness of the proposed strategies in terms of exploration and exploitation. To carry out the study, we compare the following algorithms. The first algorithm is the proposed algorithm without the central learning strategy, which is named as MAPSO/wcCL. The second algorithm is the proposed algorithm without the moving average strategy, which is named as MAPSO/wcMA. These algorithms are compared with MAPSO on the CEC 2017 test suite with 10 dimensions.
Tables 5 summarizes the obtained results. In the table, the means and standard deviations of the obtained results averaged over 50 independent runs are listed. We also include the rank sum hypothesis test results obtained by comparing with MAPSO and RDPSO. The reason to compare with RDPSO is that MAPSO is built upon RDPSO. Comparing with RDPSO enables us to see which strategy is more effective for the improvement of RDPSO. To save space, we omit the obtained results by MAPSO and RDPSO in the table since these results are available in Table 2.
The same as previous experimental settings, symbols '+', '−' and '≈' indicate that MAPSO performs worse than, better than and similar to the compared algorithms according to the Wilcoxon's rank sum hypothesis test at a significance level of 5%. Furthermore, we use symbols ' †', ' ‡' and '≈' to indicate the performance comparison with RDPSO in terms of the rank sum hypothesis test.
From Table 5, it is seen that according to the rank sum hypothesis test, without the central learning strategy,  /wcMA on the 29 test problems with 10 dimensions, where the mean values and  the standard deviations are shown. Symbols '+', '−' and '≈' (resp. ' †', ' ‡' and '≈') indicate that MAPSO (resp. RDPSO) performs worse than, better than and similar to the compared algorithms according to the the Wilcoxon's rank sum hypothesis test at a significance level of 5%. FIGURE 6. The corrected yield curve structure obtained by using the obtained parameters to correct the NS model.
MAPSO/wcCL achieves only three significantly better results than MAPSO, but 23 significantly better results than RDPSO. On the other hand, without the moving average strategy, MAPSO/wcMA achieves five significantly better results than MAPSO, but 25 better results than RDPSO. This clearly indicates that each proposed strategy can improve the performance of RDPSO. We notice that the disadvantage of RDPSO lies on its exploration capability since the random term itself may be not effectively enough to help the search escape from local minimum. In consideration of the better performance of MAPSO/wcCL and MAPSO/wcMA, we may conclude that the proposed strategies can indeed improve the exploration capability.
The results in Table 5 also indicate that each strategy alone cannot beat the usage of the combined strategies. Particularly, it is seen that for the unimodal and multimodal functions (F1-F10), MAPSO/wcMA and MAPSO/wcCL perform similarly. However, for the hybrid functions (F11-F20), MAPSO/wcMA achieves one more better better than MAPSO/wcMA in comparison with MAPSO. For the composition functions (F21-F30), however, MAPSO/wcCL achieves one more significantly better, and five similar results than MAPSO/wcMA. Since in MAPSO/wcCL, only the moving average strategy is used, we may conclude that this strategy is more beneficial to improve the performance of PSO than the central learning strategy.

V. APPLICATION ON THE YIELD CURVE ESTIMATION PROBLEM
In this section, we apply MAPSO on the yield curve estimation problem taking Shanghai interbank offered rate (Shibor) as an example, and compare it with a variety of algorithms, including GA [59], DE [60], SRPSO [22] and some classical optimization algorithms. The classical algorithms include the sequential least squares programming (SL-SQP) [61], L-BFGS-B [62] and Nelder-Mead [63]. Table 6 shows the obtained results including the means and standard deviations obtained on the problem with 50 different weight settings. The compared algorithms were allowed to run for 10 seconds. The Wilcoxon rank sum test was carried out at a significance level of 5%.
The hypothesis test suggests that MAPSO perform significantly better than all the compared algorithms. It is seen from the table that MAPSO obtains not only smaller mean values, but rather small variances. This implies that MAPSO is very stable in solving the yield curve estimation problem. It performs even better than the classical algorithms, such as SL-SQP, BFGS, and L-BFGS-B. Notice that these classical optimization algorithms are often considered to be more effective than population-based algorithms. Using the obtained optimal parameters by MAPSO to correct the NS model, we can obtain the new yield curve structure as shown in Fig. 6. In the figure, the red dots represent the official Shibor price, while the blue dots are the actual transaction data in the previous day. The green line represents the interest rate time structure corrected by using the actual transaction data. From the figure, we see that the overnight, one week and fortnight Shibor prices (red dots) usually are lower than the actual transaction interest rates (blue dots). However, the term structure of interest rates after our correction fit the actual transaction better for the secondary market. The corrected term structure provides very useful value reference for recent short-term arbitrage operation, which is exactly the motivation for us to find the global optimum of Problem 29.

VI. CONCLUSION
In this paper, we proposed a novel PSO algorithm, named MAPSO, based on the moving average strategy and central learning strategy, inspired by the movement of electron in conductor in an external electric field. Firstly, the moving average of the local best positions was proposed to take effect on the random movement of particles aiming to improve its exploration capability, while the centre position of all the local best positions was used to update the global best position. In the experiments, MAPSO was compared with several known PSOs on the widely-used CEC 2017 test suite with 10 and 30 dimensions. The experimental results revealed that MAPSO performs significantly better than the compared algorithms in general. Finally, MAPSO was applied to a challenging problem in macro-economy, the yield curve estimation problem. The comparison with a variety of algorithms showed that MAPSO achieved the state-of-the-art performance.
JIAQI ZHANG received the B.Sc. degree in human resource management from Shaanxi Normal University, Xi'an, China, in 2015, and the M.Sc. degree in management from the University of Bristol, U.K., in 2016. She is currently working with the Xi'an University of Finance and Economics, Xi'an. Her research interests include evolutionary algorithms and its application in finance and economics.
BO SHI received the Ph.D. degree in economics from Nanjing University, Nanjing, China. He is a Full Professor of economics with the Department of Mathematical Economics and Statistics, Northwest University, Xi'an, China. He has authored over 100 publications on mathematical economics and statistics. His research interests include macro-economics, energy economics, and developing political economics. VOLUME 10, 2022