Sana Ben Hamida
email: sana.ben_hamida@u-paris10.fr

Wafa Abdelmalek

Fathi Abid

APPLYING DYNAMIC TRAINING-SUBSET SELECTION METHODS USING GENETIC PROGRAMMING FOR FORECASTING IMPLIED VOLATILITY

Keywords: genetic programming, implied volatility forecast, static training-subset selection, dynamic training-subset selection, mean squared errors, percentage of non-fitted observations

Volatility is a key variable in option pricing, trading, and hedging strategies. The purpose of this article is to improve the accuracy of forecasting implied volatility using an extension of genetic programming (GP) by means of dynamic training-subset selection methods. These methods manipulate the training data in order to improve the out-of-sample patterns fitting. When applied with the static subset selection method using a single training data sample, GP could generate forecasting models, which are not adapted to some out-of-sample fitness cases. In order to improve the predictive accuracy of generated GP patterns, dynamic subset selection methods are introduced to the GP algorithm allowing a regular change of the training sample during evolution. Four dynamic training-subset selection methods are proposed based on random, sequential, or adaptive subset selection. The latest approach uses an adaptive subset weight measuring the sample difficulty according to the fitness cases' errors. Using real data from S&P500 index options, these techniques are compared with the static subset selection method. Based on mean squared error total and percentage of non-fitted observations, results show that the dynamic approach improves the forecasting performance of the generated GP models, especially those obtained from the adaptive-random training-subset selection method applied to the whole set of training samples.

INTRODUCTION

Financial market volatility is a key variable in financial investment decisions and plays a central role in derivative valuation and in conducting dynamic hedging strategies. To assess the fair value of an option or to hedge market risk, an investor needs to specify his expectations regarding future volatility. Because of their forward-looking nature, option prices are especially useful for extracting such information. A number of investigations supported the idea of using implied volatility as a good predictor of future volatility [START_REF] Latané | Standard deviations of stock price ratios implied in option prices[END_REF][START_REF] Chiras | The information content of option prices and a test of market efficiency[END_REF][START_REF] Fleming | The rationality of market volatility forecasts implied by S&P 100 index option prices[END_REF][START_REF] Blair | Forecasting S& P100 volatility: The incremental information content of implied volatilities and high frequency index returns[END_REF][START_REF] Corrado | The forecast quality of CBOE implied volatility indexes[END_REF]. Assuming that an option pricing model correctly represents investors' behavior, the implied volatility can be derived from observed option prices by appropriately inverting the option pricing model. In contrast, genetic programming (GP) offers explicit formulas that can compute directly the implied volatility expressed as a function of option prices and other observable variables. This volatility's forecasting approach should be free of strong assumptions and more flexible than parametric models. GP [START_REF] Koza | Genetic Programming: On the Programming of Computers by Means of Natural Selection[END_REF]) is an evolutionary-based search technique, which is based on the principles of natural evolution. Using its basic and flexible tree-structured representation, GP is capable of solving some difficult problems without requiring the user to know or specify the form or the structure of the solution in advance.

General programming has been proven to be successful at forecasting time series (TS) volatility in different markets, such as foreign exchange and index markets. [START_REF] Neely | Using a genetic program to predict exchange rate volatility, Chapter 13[END_REF] have tested the forecasting performance of GP for USD-DEM and USD-YEN daily exchange rates against that of the generalized autoregressive conditional heteroskedasticity (GARCH) model [START_REF] Baillie | Fractionally integrated generalized autoregressive conditional heteroskedasticity[END_REF]) and a related RiskMetrics volatility forecast over different time horizons. According to various accuracy criteria, GP has produced significantly superior results. Using high frequency foreign exchange USD-CHF and USD-JPY TS, [START_REF] Zumbach | Genetic programming with syntactic restrictions applied to financial volatility forecasting[END_REF] have compared the GP forecasting accuracy to that of historical volatilities and some popular ARCH type models, notably the GARCH (1,1) model of [START_REF] Bollerslev | Generalized autoregressive conditional heteroscedasticity[END_REF], the fractionally integrated GARCH model of [START_REF] Baillie | Fractionally integrated generalized autoregressive conditional heteroskedasticity[END_REF], and the heterogeneous ARCH model of [START_REF] Müller | Volatilities of different time resolutions-analyzing the dynamics of market components[END_REF]. According to the root mean squared errors (MSEs), the generated GP volatility models have outperformed the benchmarks. Using historical returns of Nikkei 225 and S&P500 indices, [START_REF] Chen | Using genetic programming to model volatility in financial time series[END_REF] have applied a recursive genetic programming approach to estimate volatility by simultaneously detecting and adapting to structural changes. Results have shown that recursive genetic programming is a promising tool for the study of structural changes. Applying a combination of techniques such as evolutionary algorithms GA and GP, [START_REF] Ma | An engineering approach to forecast volatility of financial indices[END_REF][START_REF] Ma | Volatility forecasting using time series data mining and evolutionary computation techniques[END_REF] have proposed a systematic approach to address specifically nonlinear problems in the forecast of financial indices using intraday data of S&P100 and S&P500 indices. As a result, accuracy of forecasting has reached an average of over 75% surpassing other publicly available results on the forecast of any financial index. [START_REF] Abdelmalek | Selecting the best forecasting-implied volatility model using genetic programming[END_REF] have extended the studies mentioned earlier by forecasting the implied volatility of Black-Scholes (BS) from the S&P500 index call options instead of the integrated volatility based on historical returns. They have considered the problem of managing too large databases when training GP. As a consequence, they have proposed to split data into smaller subsets by TS and moneyness-time to maturity (MTM) classes and to train GP separately on all learning subsamples. The proposed approach is called static trainingsubset selection method. According to the total and out-of-sample MSEs, results have shown that TS models seem to be more accurate in forecasting implied volatility than MTM models. Such an approach has provided some local solutions not adaptive to the enlarged data set, especially when learning with a MTM sample. According to [START_REF] Gilli | Heuristic optimisation in financial modelling[END_REF], the relationship between in-sample fit and out-of-sample performance is not monotonous and an optimal in-sample solution might be ineffective when applied to out-of-sample data.

The present article investigates the application of a dynamic subset selection method, in which the training-subset samples change during the GP run. This allows GP to learn simultaneously on all training subsamples rather than just a single subset, which seems to have better generalization ability. This technique aims to intensify search space exploration and thus enhance the robustness of GP with large data set. Robustness is an important feature of an evolved program [START_REF] Ito | Robustness of robot programs generated by genetic programming[END_REF]. It is defined as the ability to cope with noisy or unknown situations. The use of the dynamic training-subset selection method could reduce the problem of fitting out-of-sample patterns and could improve the forecasting accuracy. The major contribution of this article is the use of the adaptive subset selection (ASS) method, which is performed in proportion to a ratio of difficulty associated to each training sample. Two other selection methods are used for purposes of comparison: the random subset selection (RSS) method where samples are selected in a random way and the sequential subset selection (SSS) method where samples are selected in a regular way. Comparative experiments are provided to show how dynamic training subset-selection methods are applied to improve the robustness of GP to generate general models relative to static training-subset selection method. Using the total MSE and percentage of non-fitted observations (NFO) as performance criteria, results show that the forecasting accuracy is improved with the ASS.

The remainder of the article is organized as follows. Section 2 presents the research objectives and theoretical foundation regarding forecasting implied volatility and usefulness of subset selection methods. Section 3 illustrates the research design and methodology. Section 4 provides a description of the dynamic subset selection methods implemented. Section 5 reports and discusses the results of the comparison between static and dynamic selection methods. Finally, Section 6 summarizes and gives perspectives.

RESEARCH OBJECTIVES AND THEORETICAL FOUNDATION

This article addresses the application of GP by means of a dynamic subset selection method for the purposes of forecasting implied volatility. The objective is to improve forecasting accuracy.

In financial volatility forecasting, it is important to find accurate models fitting a maximum number of input cases from learning data. The goal is to correctly predict the volatility of new input data. By considering the assumption that there are inherently many different patterns in financial series [START_REF] Povinelli | Times series data mining: Identifying temporal patterns for characterization and prediction of time series events[END_REF], instead of using one formula (model) to explain the entire data, a better idea would be to select a set of best models that could be combined and used to forecast future implied volatility values. As a consequence, the main interest of this work lies in applying a dynamic subset selection method for training GP on the full training samples, so that the entire input data become in-learning-sample. The main contribution of this article is the use of the ASS method. This method is inspired by the dynamic subset selection proposed by [START_REF] Gathercole | Dynamic training subset selection for supervised learning in genetic programming[END_REF]. It is based on the assumption that there is a benefit in focusing the GP's abilities on difficult training samples, that is, the ones that have the highest MSE. The challenge is to make GP adaptive to all training samples and be able to generate general models. The originality of this method is to assign weights to each training-subset and update these weights through the generations. Initial weight is assigned according to the initialization method adopted (random or sequential) and increases each time an individual is not able to solve the corresponding fitness sample cases. This approach lightens the training task for GP and favors the discovery of solutions that are more robust across different data samples.

In the following subsections, background information regarding implied volatility is first introduced. Second, a discussion about the usefulness of subset selection methods is made. This is to explain the originality of the dynamic selection methods created versus the use of previous sampling methods.

Implied Volatility

There are two approaches to generate volatility forecasts. One is to extract information about the variance of future returns from their history and the second is to elicit market expectations about the future volatility from observed option prices. Options markets provide market participants and policy makers with a rich source of information for gaging market sentiment. An option contract is a derivative security that gives the holder the right to buy (call) or to sell (put) the underlying asset by a certain date for a certain price. The price in the contract is known as the exercise price or strike price. The date in the contract is known as the expiration date or maturity. American options can be exercised at any time up to the expiration date. European options can be exercised only on the expiration date itself. Options contracts can be divided into several classes according to either moneyness or term to expiration. By the term to expiration, an option contract can be short term, medium term, and long term. By the moneyness, a call option can be in-the-money (ITM) if the stock price is above the strike price (S > K), out-of-the-money (OTM) if the stock price is below the strike price (S < K), and vice versa for a put option. If the strike price is closest to the current value of the underlying stock (S D K), the option contract is said to be at-the-money (ATM). The difference between the stock price and the strike price represents the intrinsic value for a call option. It can be a positive number or zero otherwise. The total value of an option called premium is basically the sum of its intrinsic value and its time value.

In the BS framework, the option price is a function of variables, which are directly observable except for the volatility. The price of an option therefore depends on the market's opinion about the future volatility of the underlying asset upon which the option is written.

The [START_REF] Black | The pricing of options and corporate liabilities[END_REF] option pricing model assumes that volatility is constant. It was first derived for the European call option written on a non-dividend paying stock, as defined in Equation (1).

C BS D Sˆ.d 1 / Ke r£ ˆ.d 2 /
(1)

where d 1 D ln. S K /C.rC0:5¢ 2 /£ ¢ p £ ; d 2 D d 1 ¢ p £.
C BS denotes the price of a European call option, S is the market price of the underlying asset, K is the strike price of the option, r is the risk-free interest rate, £ is the time to maturity, ˆis the cumulative normal distribution function, and ¢ is the volatility.

By equating the observed market price C t of an option with the BS price C BS and implicitly solving for ¢, an implied volatility can be found.

9Š¢ BS

t .K; T/ 0; C BS S t ; K; £; ¢ BS t .K; T/ D C t .K; T/

(2)

According to the BS assumptions, this implicitly calculated volatility should be constant. However, it can be easily shown empirically that the implied volatility is not constant and changes with different option strike prices and expiry dates. For example, short-dated options will be less sensitive to implied volatility, whereas long-dated options will be more sensitive. This is based on the fact that long-dated options have more time value priced into them, whereas short-dated options have less. Besides, options with strike prices that are near the money are most sensitive to implied volatility changes, whereas options that are further ITM or OTM will be less sensitive to implied volatility changes. Implied volatility represents the expected volatility of a stock over the life of the option. As expectations change, option premiums react appropriately. Implied volatility is directly influenced by the supply and demand of the underlying options and by the market's expectation of the share price's direction. As expectations rise, or as the demand for an option increases, implied volatility will rise. Options that have high levels of implied volatility will result in high-priced option premiums. Conversely, as the market's expectations decrease, or demand for an option diminishes, implied volatility will decrease. Options containing lower levels of implied volatility will result in cheaper option prices. This is important because the rise and fall of implied volatility will determine how expensive or cheap time value is to the option.

Subset Selection Methods: A Discussion

Evolving programs is often a time-consuming task particularly in terms of fitness evaluation's effort. When using GP with the large set of training cases and a large population size, a very large number of tree evaluations must be carried out every generation. Many methods try to reduce the number of such evaluations by selecting a small subset of the training data set during fitness evaluation. These methods differ in how they choose proper subsets from the set of all fitness cases1 for evaluation. The simplest technique is to use a static subset. However, using a single learning sample might lead to a local optimum that solves only a part of the fitness cases. To reinforce learning from different parts of the search space, some solutions were proposed such as carrying out several runs using one subset for each run and selecting the appropriate model from different resulting models [START_REF] Abdelmalek | Selecting the best forecasting-implied volatility model using genetic programming[END_REF]. The historical subset selection [START_REF] Gathercole | Dynamic training subset selection for supervised learning in genetic programming[END_REF] extends the static subset selection by recording all fitness cases that are not solved by the best population member in any given generation over a small number of runs. These fitness cases become part of a static subset and are used in further GP runs.

A more flexible method is to pick a variety of subsets during the course of a training run. There are many ways to select different subsets from the training set. The goal is to pick the right subsets to allow the learning algorithm to proceed as fast and as accurately as possible. The simplest method for picking a different set for each generation is random. RSS [START_REF] Gathercole | Dynamic training subset selection for supervised learning in genetic programming[END_REF] chooses a new subset for each generation. Each learning data instance is selected independently with equal probability, which leads to varying subset sizes. Stochastic sampling [START_REF] Nordin | An on-line method to evolve behavior and to control a miniature robot in real time with genetic programming[END_REF][START_REF] Banzhaf | Genetic Programming. An Introduction[END_REF] chooses a new subset for each generation and for each individual, respectively, where all data cases having the same probability of being selected. As a result, different individuals will probably be evaluated on different data samples, which cast some doubts on the fairness of the selection step in the evolutionary algorithm (Section 3.2).

More efficient criteria can be used to guide the selection of the new subset based on fitness-case topology or the performance of the current GP population. Fitness-case topology-based sampling [START_REF] Lasarczyk | Dynamic subset selection based on a fitness case topology[END_REF]) relies on creating a dynamic weight for each couple of fitness case that is updated by the evolutionary system according to the number of GP solutions able to solve both of the fitness cases of each couple. The learning subset is constructed by a random selection from the set of couples with small weights. Dynamic subset selection [START_REF] Gathercole | An investigation of supervised learning in genetic programming[END_REF]Ross 1994, 1997;[START_REF] Gathercole | An investigation of supervised learning in genetic programming[END_REF]) makes use of the difficulty of each training case, that is, how often it is misclassified, and its age, that is, how many generations since it was last selected. This has worked well on some large classification problems, using less computer resources to produce better results than standard GP.

However, the resampling procedure during the course of a training run increases the complexity of the GP system and requires the entire data set to be stored in the main memory. With GP systems, using the full learning data might be impossible when the input sample does not fit within the main memory and could cause serious problems in the realization of predictors. In this case, data reduction through the partitioning of the data set into smaller subsets seems to be a good approach [START_REF] Abdelmalek | Selecting the best forecasting-implied volatility model using genetic programming[END_REF].

In this article, we suggest to construct subset data with a fixed size from the full database using data division schemes and to apply dynamic training-subset selection methods to select subsets already built up for learning process. With dynamic subset selection methods, a new subset is picked in each g generations (g is the number of generations to change sample). These methods can be used to prevent a biasing influence of subset selection in evolution. These methods differ from the existing approaches as they don't extract a fixed number of fitness cases from the training set but select a subset from a set of subsamples data already built up using the data division scheme, which avoids increasing the complexity of the GP engine. In this article, we proposed four dynamic training-subset selection methods: RSS, SSS, adaptive-sequential subset selection (ASSS) method, and adaptiverandom subset selection (ARSS) method. The RSS and SSS allow the GP to learn on all training samples sequentially (SSS or randomly RSS). However, with these methods, there is no certainty that GP will focus on the samples that are difficult to learn. Then, the ASSS and the ARSS, which are variants of ASS, are introduced to focus the GP's attention onto the difficult samples, that is, having the greatest MSE and then to improve the learning algorithm. An adaptive weight is associated to each subset (adaptive subset weight (ASW)) and updated each generation according to the average fitness of all the cases in the corresponding sample.

RESEARCH DESIGN AND METHODOLOGY

To achieve the research goals presented earlier, we followed three steps. The initial stage is devoted to the data preparation, to the implementation of GP2 and to the subset selection procedures. The second step is devoted to the learning process using static and dynamic subset selection methods. The last step is dedicated to results comparison and the selection of the best forecasting models. The following research steps are summarized in the following scheme.

Data Preparation

The data used to perform our empirical analysis are daily prices of European S&P500 index calls options, from the Chicago Board of Options Exchange (CBOE), for the sample period running from January 2, 2003 to August 29, 2003. The database includes the time of the quote, the expiration date, the exercise price, and the option price. Similar information for the underlying S&P 500 index is also available on a daily basis. The daily bid and ask quotes for the call options are obtained from the CBOE. Following a standard practice, we use the average of an option's bid and ask price as a stand-in for the market value of the option. Strike price intervals are 5 points and 25 for far months. The expiration months are three near-term months followed by three additional months from the March quarterly cycle (March, June, September, and December). The risk-free interest rate is approximated by using 3-month U.S. Treasury bill rates.

Two preparation procedures have been applied to the data before use, preprocessing, and division. The data preprocessing serves the purpose of "smoothing" the raw data and removing what is not essential before the machine learning algorithm is applied. It is widely accepted that preprocessing is usually beneficial and has positive effects on the learning process [START_REF] Chen | Genetic Programming and Financial Trading How Much about "What We Know[END_REF]. The original training set contains 42,504 daily observations of call option prices and their determinants. To reduce the likelihood of errors, data screening procedures are used (Harvey andWhaley 1991, 1992). Then, four exclusion filters are applied to construct the final sample. First, call options with time to maturity less than 10 days are excluded from the sample. This can be explained by the fact that implied volatilities of short-term options are very sensitive to small errors in the call price and may convey liquidity-related biases. Second, call options with low quotes are eliminated to mitigate the impact of price discreteness on option valuation. Third, deep-ITM and deep-OTM option prices are also excluded because of the lack of trading volume. Finally, option prices not satisfying the arbitrage restriction [START_REF] Merton | Theory of rational option pricing[END_REF], C S Ke r£ , are not included. The final sample contains 6670 daily option quotes, with ATM, ITM, and OTM options, respectively, taking up 37%, 34%, and 29% of the total sample.

For the data subsampling procedure, two schemes were used. For the first division scheme, the full sample is sorted by TS and for the second division scheme, by MTM. For TS, data are divided chronologically into successive samples (S 1 ; S 2 : : : S 10 ), each containing 667 daily observations. These samples will be used simultaneously for training and test steps. For MTM, data are divided into nine classes with respect to moneyness and time to maturity. 3 Each class C i is divided into a training set C L i and a test set C T i , which produce, respectively, nine training and nine test MTM subclasses. Figure 1 illustrates the two division schemes.

The Design of Genetic Programming

There are several GP techniques that might be used to deal with financial forecasting. People are familiar with regression analysis. One problem with regression analysis is that the results of the analysis depend very much on the skills and inventiveness of the experimenter. Furthermore, in many application areas, there is a considerable tradition of using only linear or quadratic models even when the data might be better fitted by a more complex model. Symbolic regression attempts to go beyond this. It consists in finding a function that fits some target values without making any assumptions about the structure of that function. Symbolic regression was one of the earliest applications of GP [START_REF] Koza | Genetic Programming: On the Programming of Computers by Means of Natural Selection[END_REF]) and has continued to be widely studied [START_REF] Keijzer | Scaled symbolic regression[END_REF][START_REF] Gustafson | On improving genetic programming for symbolic regression[END_REF][START_REF] Cai | Heat transfer correlations by symbolic regression[END_REF][START_REF] Lew | Identification of response surface models using genetic programming[END_REF]. The steps necessary to implement the GP's symbolic regression are summarized in Algorithm 1. and function sets The standard GP tree is a simple structure made by terminal (or leaf) nodes and nonterminal (or function) nodes with branches. Terminal and function sets, which are described in Table 1, define the ingredients that GP can use to create function models and to construct potential solutions.

The terminal set includes input variables, mainly, the call option price divided by strike price C K , the index price divided by strike price S K , and time to maturity £. The predictive target output is the implied volatility ¢ BS t computed using the BS formula. The function set includes basic mathematical operators and BS components. The mathematical operators we use are the basic arithmetic operators together with the cosine functions and the sine functions. The BS components involve the log function (ln), the exponential function (exp), the square root function p , and the normal cumulative distribution function .ˆ/, which may be useful for implied volatility models.

Initialization General programming starts by randomly creating an initial population of trees, which are generated by randomly picking nodes from a given terminal set and function set (Table 1). The initialization scheme used in this article is the ramped half-and-half method [START_REF] Koza | Genetic Programming: On the Programming of Computers by Means of Natural Selection[END_REF], which is a combination of the full and grow initialization methods. This method involves generating an equal number of trees using a maximum initial depth that ranges from 2 to 6, as specified in Table 2. For each level of depth, half of the initial trees are generated via the full method, and the other half is generated via the grow method.

Fitness function

The evolutionary process is driven by a fitness function that evaluates the performance of each individual (tree) in the population. The fitness criterion used for the selection of the best individuals is the MSE between the target output volatility ¢ BS t and the generated GP volatility .O ¢ t /, computed as follows:

MSE D 1 N N X t D1 ¢ BS t O ¢ t 2 (3)
where, N is the number of fitness cases in the learning sample. At the end of evolution, each individual is evaluated according to the MSE computed with Equation (3) using a test sample that must be different from the learning sample.

Selection

Based on fitness measure, GP probabilistically selects the fitter individuals from the population to act as the parents of the next generation. Selection determines which individuals of the population will have all or some of their genetic material passed to the next generation. In general, parents displaying a higher level of performance are more likely to be selected with the hope that they can produce better offspring with larger chance. The most commonly used method for selecting individuals in GP is tournament selection. In tournament selection, a number of individuals, called the tournament size, are selected randomly from the population and they compete with each other. The best is to be selected. As specified in Table 2, the tournament size used for experiments is equal to 4.

Genetic operators

Crossover and mutation are the two basic operators that are applied to the selected individuals in order to generate new individuals for the next generation. They are needed to explore the search space.

Crossover operator.

The most commonly used form of crossover is sub-tree crossover. Given two parents, subtree crossover randomly selects a crossover point (a node) in each parent tree. Then, it creates the offspring by replacing the sub-tree rooted at the crossover point in a copy of the first parent with a copy of the sub-tree rooted at the crossover point in the second parent. All non-terminal nodes (except the root node) have the same probability to be selected. The generated offspring should not surpass the fixed size. As indicated in Table 2, the crossover operator is used to generate about 60% of the individuals in the population.

The maximum tree size (measured by depth) allowed after the crossover is 17. This is a popular number used to limit the size of the tree [START_REF] Koza | Genetic Programming: On the Programming of Computers by Means of Natural Selection[END_REF]. It is large enough to accommodate complicated formulas and works in practice.

Mutation operator.

The basic role of the mutation operator in the evolutionary process is to ensure diversity in the population. It affects small random changes in a tree by randomly altering nodes or subtrees to create a new offspring and continue the search process. Many mutation operators are used in GP. The most commonly used form of mutation in GP is sub-tree (or branch) mutation. It replaces a randomly selected sub-tree with another randomly created sub-tree [START_REF] Koza | Genetic Programming: On the Programming of Computers by Means of Natural Selection[END_REF]. Another common form of mutation is point mutation (or node replacement), which randomly changes a node in the individual and replaces it with another node with the same arity [START_REF] Mckay | Using a tree structural genetic algorithm to perform symbolic regression[END_REF]. Expansion mutation randomly selects a terminal node in the tree and then replaces it with a new randomly generated sub-tree. As indicated in Table 2, Branch mutation is applied with a rate of 20%; point and expansion mutations are applied with a rate of 10%, respectively.

The parameter choices for crossover and mutation are clearly critical in ensuring a successful GP application. They impact on population diversity and the ability of GP to escape from local optima [START_REF] Yin | Genetic programming for dynamic environments[END_REF]).

Replacement

Once the new population has been created, the current population (parents) is replaced by the new one (offspring). We use a comma replacement method to replace parents for the next generation [START_REF] Schwefel | Numerical Optimization of Computer Models[END_REF]. This method selects the best offspring to replace the parents. If is the population size and is the number of new individuals (which can be larger than ), the population is constructed using the best out of the new individuals.

Termination criterion

The termination criterion we use is the maximum number of generations to be run. We take 400 and 1,000 for static and dynamic training-subset selection, respectively. In the dynamic training-subset selection approach, the maximum number of generations is increased to allow GP to train on the maximum of samples simultaneously. Typically, the single best-so-far individual is then obtained and designed as the result of the run.

Parameters of the GP design used in this work are summarized in Table 2. The optimal set of genetic parameters is determined based on a series of trial and error experiments.

DESCRIPTION OF DYNAMIC TRAINING-SUBSET SELECTION METHODS

To implement the dynamic training-subset selection methods, two supporting decisiondesigns are necessary to be fixed. First, the frequency of the replacement of the current training subset; such a frequency is designed by a fixed iteration number (g) added as a parameter to the GP system. The value of this parameter must be chosen such as the GP system has enough time to adapt the genetic material in the population to make it able to solve the current subset. However, g does not have a very high value in order to keep the population diversity necessary for the next learning step.

Second, the procedure of replacement and selection of the new training sample is determined. To design the subset replacement protocol, two approaches are available. In the first approach, all subsamples are treated equally. Subsamples are then selected with a uniform probability (RSS) or in a regular way by taking the samples in turn (SSS). In the second approach, subset selection is performed proportionally to a predefined ratio of ASS method. It involves an ASW updated at each g generations, according to the fitness of the best current model computed according to all observations in each subset. The lower the sample performance is, the higher the selection pressure is. Samples having a great weight are selected more frequently in the learning process. The goal of this approach is to focus on the samples having the fitness cases the most difficult to learn. The GP is then guided to adapt its models to the more difficult samples. Only the individuals with the desirable characteristics that are well adapted to the environmental change will survive. Thereby, the problem of overfitting that can be encountered with the static subset selection can be prevented. We have two variants of ASS: ASSS and ARSS methods, differing on the initialization procedure of the sample weights. The operating principle of the dynamic selection is as follows.

Let S be the set of training samples S i.iD1:::k/ , where k is the total number of samples. A selection probability P (S i / is assigned to each sample S i , which is changed in each g generation (g is the number of generations to change sample) according to the dynamic training-subset selection method used. Once a new training sample is selected, the best individuals are used as population for the next training samples. This procedure is repeated until the maximum number of generations is reached. This permits GP to adapt its generating process to changing data in response to feedback from the fitness function.

To decide the training-subset to select with the ASS method, an ASW is computed for each subset according to the MSE values obtained along the last g generations. The ASW and the four proposed methods, RSS, SSS, ASSS, and ARSSS, are described in the following.

Adaptive Subset Weight

For ASS method, the selection probability depends on the subset weights computed proportionally to the sample's average fitness. After g generations, the weight of the learning sample S i is updated as follows:

W.S i / D g P tD1 M P jD1 f.X j / M g (4)
where M is the population size, g is the number of generations to change sample, and f .X j / is the MSE of the individual X j , where X j 2 P t .P t is the current population).

If several individuals in the population have difficulty to solve some fitness cases in a sample S i , then this sample will have a high subset weight and a great probability to be selected for the next learning step.

The update of the subset weights does not increase the complexity of the GP program and does not need additional computational cost as some dynamic subset selection methods yet proposed in the literature [START_REF] Gathercole | Small populations over many generations can beat large populations over few generations in genetic programming[END_REF][START_REF] Lasarczyk | Dynamic subset selection based on a fitness case topology[END_REF]). Indeed, the individuals' fitness (used in Equation ( 4)) is computed by the GP system and the subsets selection probabilities are computed only by each g generation.

Random Training-Subset Selection Method:

This method randomly selects the training samples with replacement. All the samples from S have the same probability to be selected: P (S i / D 1=k; 1 Ä i Ä k. Figure 2 illustrates an example of the best fitness (MSE) curve along evolution using the RSS method. As selection of training samples is random, the performance of the current population changes with the training sample used for evolving the genetic program.

Sequential Training-Subset Selection Method

This method selects all the training samples in turn. All the learning subsets are used during the evolution in an iterative way. If, at generation g 1, the current training sample is S i; then at generation g: P (S j / D 1, with j D iC1 if i < k or j D 1 if i D k. Figure 3 illustrates an example of the best fitness (MSE) curve along evolution using the SSS method. It shows that all the learning subsets are used during the evolution in an iterative way.

Adaptive Training-Subset Selection Method

Instead of selecting a training-subset data in a random or sequential way, one can use an adaptive approach to dynamically select "difficult" training-subset data having high fitness errors (MSE). This approach is inspired by the dynamic subset selection method of [START_REF] Gathercole | Dynamic training subset selection for supervised learning in genetic programming[END_REF], which is based on the idea of dynamically selecting instances, not training samples, which are difficult and/or have not been selected for several generations. ASS simplifies this method by selecting training samples, not instances, containing unsolved fitness cases.

Selection is made according to the ASW. After g generations, training samples are reordered, so that most difficult training samples, those having higher ASW, will be moved to the beginning of the ordered training list, and the easiest training samples, those having smaller ASW, will be moved to the end of the ordered training list.

Adaptive-sequential training-subset selection method

The initial weights are initialized with a constant and the selection of samples is done in an iterative way: W .S i / D C; 1 Ä I Ä k. Then, for the k first generations, selection of training samples is made in the order using the SSS method (Section 4.1). Later, from the generation kC1, samples are selected for the next step according to the adaptive approach based on the reordering procedure (Equation ( 4)). Figure 4 illustrates an example of the best fitness (MSE) curve along evolution using the ASSS method.

Adaptive-random training-subset selection method

The ARSS method uses the same procedure as the ASSS method, except that the initial weights are generated randomly at the start of running rather than initialized with a constant:

For t D 0; W.S i / D Q P i ; Q P i 2 OE0; 1; 1 Ä i Ä k.
Then, for the few first generations, samples are selected using the RSS method (Section 4.2) after the selection of samples is made using the adaptive approach based on the reordering procedure (Equation ( 4)). Figure 5 illustrates an example of the best fitness (MSE) curve along evolution using ARSS method.

FINDINGS AND RESULTS ANALYSIS

Experiments

The experiments were performed in two major phases: static subset selection experiments (phase 1) and dynamic subset selection experiments (phase 2). 3 1 7 3 4 6 1 3 7 1 1 4 2 3 2 2 9 2 5 Phase 1. First, the genetic program is trained separately on each of the first nine TS subsamples .S 1 ; : : : ; S 9 / using 10 different seeds and is tested on the subset data from the immediately following date .S 2 ; : : : ; S 10 /. Second, using the same genetic parameters and the random seeds applied on TS data, GP is trained separately on each of the first nine MTM subclasses C L 1 ; : : : ; C L 9 and is tested on the second nine MTM subclasses C T 1 ; : : : ; C T 9 . Actually, each subsample is independently evolvable by GP, and the best individual generated from each subsample is selected.

Phase 2. First, the genetic program is trained on the first nine TS subsamples simultaneously .S 1 ; : : : ; S 9 / using 10 different seeds and it is tested only on the tenth subsample data (S 10 ). Second, GP is trained on the first nine MTM subclasses simultaneously C L 1 ; : : : ; C L 9 , and it is tested on the second nine MTM subclasses regrouped in one test sample data C T 1 C C T 2 : : : C C T 9 . Third, GP is trained on both the nine TS subsamples and the nine MTM subclasses simultaneously S 1 ; : : : ; S 9 I C L 1 ; : : : ; C L 9 , and it is tested on one test sample data composed of the TS and MTM test data S 10 C C T 1 C C T 2 : : : C C T 9 . Table 3 summarizes the training and test data samples used for static and dynamic trainingsubset selection methods, respectively.

For each case, the best individual (tree function) is selected according to the MSE's fitness measure. Selected models are then analyzed and compared with each other according to different measures as described in the following.

GP Model Selection

Results analysis focuses on the comparison between GP solutions given by static and dynamic training-subset selection methods in terms of their ability to forecast implied volatility. GP models subject to the comparative study are selected as follows.

First, selection of the best generated GP volatility model, relative to each training set, for TS, MTM, and both TS and MTM classifications, is made according to the training and test MSE. For static training-subset selection method, nine generated GP volatility models (M 1 S 1 : : : M 9 S 9 / are selected for TS samples and similarly, nine generated GP volatility models (M 1 C 1 : : : M 9 C 9 / are selected for MTM classes.

Second, for dynamic training-subset selection methods (RSS, SSS, ASSS, and ARSS), four generated GP volatility models are selected for TS classification (MSR, MSS, MSAS, and MSAR). Similarly, four generated GP volatility models are selected for MTM classi- fication (MCR, MCS, MCAS, and MCAR) and four generated GP volatility models are selected for global classification using both TS and MTM classes (MGR, MGS, MGAS, and MGAR). Table 4 summarizes the 30 volatility models selected for the comparative study.

Results Analysis

To assess the accuracy of the training-subset selection methods applied as well as the generated GP volatility models selected, two measures are used. First, the MSE total is computed using the same formula as the basic MSE (Equation 3) but according to the enlarged data sample. The best model is the model that provides the smallest forecasting error. MSE total aims to measure the generalization ability of the GP generated models. Second, the number of NFO for a given data sample is used to describe how well a model fits the sample observations; otherwise, it measures the model's ability to solve all fitness cases in the corresponding subset. An observation is supposed to be well fitted if the corresponding error (The absolute difference between target and forecasted output volatility) is less than 0.1 (value determined according to the experimental results). NFO aims to compare the final solutions in terms of adaptation level to all input fitness cases.

Figure 6 describes the performance of the best generated GP volatility models, using static and dynamic training-subset selection methods, according to the MSE total for all data samples.

Figure 6 (a) presents the MSE total relative to the 18 generated GP volatility models, using static training-subset selection method, selected for TS samples and for MTM classes. Some extreme MSE values for MTM data are not shown in this figure. As shown, the performance of the static models is not uniform. Total errors are higher for the MTM classes than for the TS samples. Indeed, the MSE exceed 1 with some fitness cases of MTM classes, and it doesn't reach 0.006 for all TS sample cases. Thus, it seems that TS models are more general than MTM models. The difference in accuracy between GP applied on TS samples and GP applied on MTM classes is very striking. With MTM classes, GP was unable to find satisfactory models with high forecasting ability, which might be caused by insufficient search intensity. Furthermore, Figure 6 (a) shows that the generated GP models M4S4 and M4C4 have the smallest MSE in enlarged sample, for the TS and for the MTM classes, respectively. They seem to be more accurate in forecasting implied volatility than the other models.

Figure 6(b) illustrates the MSE relative to the 12 generated GP volatility models, using dynamic training-subset selection methods (RSS, SSS, ASSS, and ARSS), selected for TS samples, for MTM classes and for global classification using both TS and MTM classes. It appears throughout Figure 6(b) that the generated GP volatility models, relative to each dynamic subset selection method, are performing on the enlarged sample and present forecasting errors, which are small and much close. The MSE relative to these models don't reach 0.003 for all sampling scheme data, except the MCS model generated using the SSS method for MTM classes. Figure 6(b) shows that, for the TS samples, the MSS model generated using the SSS method has the smallest MSE on the enlarged sample. For the MTM classes, the MCAR model generated using ARSS method outperforms the other models generated using the other methods. For both TS and MTM data, the MGAR model generated using ARSS method presents the highest accuracy in enlarged sample. Overall, the best forecasting's performance is achieved by the ARSS method. This can be explained by the fact that this method permits to generate more general models adaptive to all sample data. In order to demonstrate more explicitly the improvements accomplished by the dynamic subset selection, the percentage of NFO in the whole learning sample is computed for the best generated GP models using static and dynamic training-subset selection methods. Results are illustrated in Figure 7 for TS samples, MTM classes, and both TS and MTM samples.

Figure 7(a) describes the performance of the 18 generated GP volatility models relative to TS samples and MTM classes, using static training-subset selection method, according to the percentage of NFO.

It appears throughout Figure 7(a) that the TS models are best-fit patterns compared with MTM models. In fact, the percentages of NFO given by TS models are markedly lower than those achieved by MTM models for the enlarged sample. The NFO percentages given by TS models do not exceed 12%, but they reach 62% with MTM models. Indeed, training is more homogeneous with TS samples than with MTM classes. This leads us to affirm that the static subset selection depends on the sampling scheme. Otherwise, it's important to note that each model (from TS models or MTM models) fits well to the data on which is trained. Actually, the generated GP volatility models M 4 S 4 and M 4 C 4 present the lower percentages of NFO for TS and MTM samples, respectively.

Results show that the gap between the models' performances is more remarkable for MTM classes. This can be explained by the fact that the majority of models fit well their training classes but not other out-of samples. GP is most generally not efficient when the training data pattern is different from the one relative to out-of-sample data [START_REF] Chen | Genetic Programming and Financial Trading How Much about "What We Know[END_REF]. robust models with high forecasting accuracy. For MTM training, the NFO percentage is relatively high and varies from 3.36% (MCR) to 10.25% (MCS). It could be observed that the dynamic methods (SSS, ASSS, and ARSS) have the highest performance with global classification data (both TS and MTM) than with TS samples or MTM classes. In fact, the GP volatility models generated using these methods have NFO percentages lower for global classification (MGS, MGAS, and MGAR) than for TS samples (MSS, MSAS, and MSAR) and MTM classes (MCS, MCAS, and MCAR), respectively. The RSS method gives high training quality for the MTM classes. This can be attributed to the randomness of the learning sample order, which allows a better adaptation of the population to the environmental changes. Although the SSS method provides the best accuracy with global classification (NFO D 0.63%) and less to TS samples (NFO D 1.2%), it was unable to provide high performance for MTM classes (NFO D 10.25%). This can be explained by the fact that the sample selection scheme with the SSS method is steady and unchanged along evolution.

Comparison between static (Figure 7(a)) and dynamic training-subset selection methods (Figure 7(b)) reveals that the percentages of NFO are in most cases lower for the models generated using dynamic training-subset selection methods than for the models generated using static training-subset selection method. The NFO percentage is reduced for most samples, in particular, the MTM classes when using dynamic training-subset selection methods.

Overall, according to the MSE total and the NFO percentage, the generated GP models using dynamic training-subset selection methods exhibit a very high accuracy relative to that using static training-subset selection method. This observation is confirmed by the measures in Table 5, which illustrates the average of MSE total and NFO percentage for all the models obtained with static and dynamic subset selection for each set of learning data samples (TS series, MTM classes, and both TS series and MTM classes). The dynamic subset selection was able to achieve the desired goal and improve the GP research process in order to fit better the whole learning data. It presents the smallest averages of MSE total and NFO percentage. Otherwise, we can note that the diversity of the input samples in the case of dynamic selection applied on both TS and MTM samples makes GP more robust in supervised learning and so as the generated forecasting models outperform all other models.

The last step of the present work is the selection of the most accurate models from the 30 GP generated models as listed in Table 4.

Based on the MSE total and the percentage of NFO as performance criteria, the generated GP volatility models M 4 S 4 and M 4 C 4 are selected for static training-subset selection method. Similarly, the generated GP volatility models MSS, MCAR, and MGAR are According to MSE total, the MCAR and MGAR models generated using the ARSS method seem to outperform the MSS model generated using the SSS method. Although the latter has the lowest percentage of NFO, it seems to be less performing than the MCAR and MGAR models. This can be explained by two points. First, the TS model MSS presents the highest MSE relative to the other models. Second, the total of error values computed for the non-fitted observations is higher for MSS than the other models even it presents the lowest percentage of NFO.

Comparison between models reveals that the best models generated, respectively, by static .M 4 S 4 / and dynamic selection methods (MCAR and MGAR) present small and very close total MSE values. Although the generated GP volatility models M 4 S 4 and MCAR have a total MSE smaller than the MGAR model, the latest seems to be more accurate in forecasting implied volatility than the other models. This can be explained by the fact that, on one hand, the difference between forecasting errors is small, and on the other hand, the MGAR model is more general than MCAR and M 4 S 4 models because it is adaptive to all TS and MTM classes simultaneously. According to the percentage of NFO, the MGAR model presents a percentage of NFO relatively higher than the M 4 S 4 model, trained only on TS data, and relatively smaller than the MCAR model, trained only on MTM classes.

The decoding of these models yields the following GP volatility forecasting formulas:

M4S4: ¢ GP D exp " ln  ˆÂ C K Ãà r £ 2 C K C S K ! cos  C K à # MCAR: ¢ GP D £ q C K ln. S K C£/ S K 2 S K C £ MGAR: ¢ GP D v u u t C K S K 6 C S K 5 £

CONCLUSION

This article presents a GP-based technique to generate implied volatility forecasting models from S&P500 index options. We have demonstrated that the accuracy of the generated models depends on the training sample, especially when learning with MTM classes. To reduce the gap between in-sample fit and out-of-sample performance, we introduced the dynamic training that aims to enlarge the training set to the whole input data. Four techniques of dynamic subset selection are proposed: RSS, where samples are selected in random way; SSS, where samples are selected in a regular way; and ASSS and ARSS, which use a training instance weight to enhance learning on the "difficult" fitness cases. These techniques are applied on the TS samples and on MTM classes and compared with the static training-subset selection method using a single sample for the learning process.

Experiments indicate that using the dynamic training with GP yields better results than applying the static training, especially when learning on TS and MTM samples simultaneously. Otherwise, based on the MSE total and the percentage of NFO as performance criteria, three generated GP volatility models are selected: M4S4 generated using the static training-subset selection method, MCAR generated using the ARSS method applied on MTM classes, and MGAR generated using the ARSS method applied on times series and MTM classes regrouped. However, the MGAR seems to be more accurate in forecasting implied volatility than MCAR and M4S4 models. This means that the MGAR model is more general than MCAR and M4S4 models because it is adaptable to all TS and MTM classes simultaneously.

Our results revealed some interesting issues for further investigation. First, the dynamic training GP can be used to forecast implied volatility of other models than BS model, notably stochastic volatility models and models with jump. Second, this work can be reexamined using data from individual stock options, American-style index options, options on futures, currency, and commodity options. Third, the performance of the generated GP volatility models can be measured in terms of trading and hedging. Finally, the GP approach can be applied to extract risk-neutral densities, which provide valuable information about market expectations. We believe that these extensions are of interest for application and will be the object of our future works.

FIGURE 1 .

 1 FIGURE 1. Data division schemes.

FitnessFIGURE 2 .

 2 FIGURE 2. Example of fitness curve of the best individuals generated by GP using RSS method for TS samples.

FIGURE 3 .

 3 FIGURE 3. Example of fitness curve of the best individuals generated by GP using SSS method for MTM classes.

FIGURE 4 .

 4 FIGURE 4. Example of curve fitness of the best individuals generated by ASSS method for TS samples.

FIGURE 5 .

 5 FIGURE 5. Example of curve fitness of the best individuals generated by ARSS method for MTM classes.

FIGURE 6 .

 6 FIGURE 6. Performance of the generated GP volatility models, using static and dynamic training-subset selection methods, according to the MSE total for TS samples (S 1 : : : S 9 ), MTM classes C L 1 : : : C L 9 and both TS and MTM samples S 1 : : : S 9 ; C L 1 : : : C L 9 .

FIGURE 7 .

 7 FIGURE 7. Performance of the generated GP volatility models, using static and dynamic training-subset selection methods, according to the percentage of non-fitted observations (NFO) for TS samples (S 1 : : : S 9 ), MTM classes C L 1 : : : C L 9 and both TS and MTM data S 1 : : : S 9 ; C L 1 : : : C L 9 .

TABLE 1 .

 1 Terminal Set and Function Set.

	Expression	Definition
	Terminal set	C/K Call price/strike price S/K Index price/strike price £ Time to maturity
		C	Addition
			Subtraction
			Multiplication
		%	Protected division: x%y D 1 if y D 0I x%y D x%y otherwise
	Function set	cos Sin ln	Cosinus function Sinus function Protected natural log: ln.x/ D ln.jxj/
		Exp Exponential function: exp.x/ D e x Sqrt Protected square root: p x D p jxj
		Ncdf Normal cumulative distribution function	Terminal

TABLE 2 .

 2 Summary of GP Parameters.

	Population size:	100
	Offspring size:	200
	Generations' number for static method:	400
	Generations' number for dynamic method:	1000
	Generations' number to change sample	20-100
	Maximum depth of new individuals:	6
	Maximum depth of the tree:	17
	Tournament size:	4
	Crossover probability:	60%
	Mutation probability:	40%
	Branch mutation:	20%
	Point mutation:	10%
	Expansion mutation:	10%

TABLE 3 .

 3 Definition of Training and Test Data Samples for Static and Dynamic Training-Subset Selection Methods.

	Subset selection	Learning data sample	Test data sample
	Static subset	1. S i 2 TS samples (S 1 ; : : : ; S 9 /	The successive TS sample S j ; j D i C 1
	Selection	(1 subset for a run)
		2. C L i 2 MTM training samples C L 1 ; : : : ; C L 9 (1 subset for a run)	The corresponding MTM test sample C T i
	Dynamic subset 1. TS samples S 1 ; : : : ; S 9	The last subset in TS samples set (S 10 /
	selection	(9 subsets for a run)
	(RSS/SSS/	2. MTM training samples	The nine MTM test samples
	ASSS/ARSS)	C L 1 ; : : : ; C L 9 (9 subsets for a run)	C T 1 C C T 2 : : : C C T 9
		3. TS samples C MTM samples The last TS sample with the nine MTM test
		S 1 ; : : : ; S 9 ; C L 1 ; : : : ; C L 9 (18 subsets for a run)	samples S 10 C C T 1 C C T 2 : : : C C T 9

TABLE 4 .

 4 Definition of GP Generated Models Given by Static and Dynamic Training-Subset Selection Methods.

	Subset selection	Learning data		GP volatility models
	Static Subset TS samples (S Selection MTM classes C L 1 ; : : : ; C L 9	M 1 C 1 : : : M 9 C 9 (MTM models)
	Dynamic Subset		RSS	SSS	ASSS	ARSS
	Selection	TS series (S 1 ; : : : ; S 9 /	MSR MSS MSAS MSAR
		MTM classes C L 1 ; : : : ; C L 9	MCR MCS MCAS MCAR
		TS series C MTM classes	MGR MGS MGAS MGAR
		S 1 ; : : : ; S 9 ; C L 1 ; : : : ; C L 9

1 ; : : : ; S 9 / M 1 S 1 : : : M 9 S 9 (TS models)

TABLE 5 .

 5 Average of MSE Total and NFO Percentage for Static and Dynamic Subset Selection Applied to the Different Learning Data Sets. The numbers in parentheses are the standard deviation corresponding to MSE values of all observations in each sample set.

				Methods
		Static subset selection	Dynamic subset selection
		TS	MTM	TS	MTM TS series C MTM
	Measures	series	classes	series	classes	classes
	Average of MSE total	0.002599	0.416320	0.002372 0.003600	0.002033
		(0.064383) (56.067)	(0.003894) (0.126530)	(0.003506)
	Average of NFO percentage 4.27%	20.29%	3.41%	5.67%	1.77%

TABLE 6 .

 6 Selection of the Best Generated GP Volatility Models, Using Static and Dynamic Training-Subset Selection Methods, in Terms of MSE Total and Percentage of NFO. Bold rows denote the selected models. selected for dynamic training-subset selection method. Table 6 reports the performance of these selected models. Results show that according to MSE total and NFO percentage, the TS model M 4 S 4 seems to be more performing than the MTM model M 4 C 4 for static training-subset selection method. It presents the lowest MSE and percentage of NFO.

	Models	MSE total	Percentage of NFO (%)
	M4S4	0.001444 (0.002727)	0.67
	M4C4	0.001710 (0.004624)	1.66
	MSS	0.002076 (0.004044)	1.20
	MCAR 0.001424 (0.003527)	3.40
	MGAR 0.001599 (0.003590)	1.45

A fitness case is an input/output pair, which measures how well an evolved individual predicts the output(s) from the input(s).

GP system is built around the evolving object library, which is an ANSI-CCC evolutionary computation framework (EO library).

According to moneyness criterion: A call option is said OTM if S=K < 0:98; ATM if S=K 2 OE0:98; 1:03OE; and ITM if S=K 1:03. According to time to maturity criterion: A call option is short term if £ < 60 days; medium term if £ 2 OE60; 180 days; and long term if £ > 180 days.

APPENDIX: SYMBOLS AND ABBREVIATIONS

Genetic programming symbols

Financial symbols GP: genetic programming C: call option price GA: genetic algorithms K: strike price NFO: non-fitted observations S: the index price MSE: mean squared error £: time to maturity ITM: in-the-money Dynamic subset selection symbols OTM: out-of-the-money ATM: at-the-money SSS: sequential subset selection method ST: short term RSS: random subset selection method MT: medium term ASS: adaptive subset selection method LT: long term ASSS: adaptive-sequential subset selection method TS: times series ARSS: adaptive-random subset selection method MTM: moneyness-time to maturity


Adnane Saoud
email: adnane.saoud@l2s.centralesupelec.fr

Antoine Girard
email: antoine.girard@l2s.centralesupelec.fr

Laurent Fribourg
email: fribourg@lsv.ens-cachan.fr

Contract based Design of Symbolic Controllers for Vehicle Platooning

Keywords: Symbolic control, Contract based design, Vehicle platooning

come

INTRODUCTION

Vehicle platoons are groups of autonomous vehicles traveling closely. Platooning makes it possible to reduce traffic congestion while increasing safety and fuel efficiency. In this paper, we present an approach to vehicle platooning based on the use of symbolic control techniques and assume-guarantee contracts. Due to scalability issues, vehicle platooning can only be solved using compositional techniques. Existing work [START_REF] Borri | Decentralized symbolic control of interconnected systems with application to vehicle platooning[END_REF], applies to discrete-time models and specifications and assumes that sampling of time in all components is synchronous. In comparison, our approach makes it possible to take into account the system behavior between sampling instants and to consider asynchronous controllers.

In Section 2, we describe the model of a platoon of vehicles and formulate the control objective. In Section 3, we briefly present the general framework introduced in [START_REF] Saoud | On the composition of discrete and continuous-time assume-guarantee contracts for invariance[END_REF], for compositional reasoning and assign an assume-guarantee contract to each vehicle. In Section 4, we explain how symbolic control techniques can be used to synthesize controllers enforcing the specified contracts. Finally, Section 5 presents numerical simulations of vehicle platooning on circular roads to show the effectiveness of the proposed approach even in the presence of asynchronous sampling and of heterogeneous vehicles in the platoon. (Due to space constraints, the simulations of the vehicle platooning on the straight road are omitted).

PROBLEM FORMULATION

In this work, we deal with two platoon configurations corresponding to straight and circular roads (see Figure 1). On a straight road, the platoon of m vehicles is made of a leader and m -1 followers. The velocity v 1 of the leader is considered as a disturbance and assumed to belong to the set [0, v max ] with v max > 0. The dynamics of follower i ∈ I = {2, . . . , m} is given by the nonlinear system:

d i = v i -v i-1 M v i = F i -f 0 -f 1 v i -f 2 v 2 i (1)
where d i ≤ 0 represents the relative distance between vehicle i and the preceding vehicle i -1, v i its velocity and F i its control input. On a circular road, there is no leader and the dynamics of each vehicle i ∈ I = {1, . . . , m} is given by (1) with the convention that v 0 = v m .

Our control objective is to synthesize controllers, giving values of input F i , for all vehicles of a platoon, whatever the number of the vehicles, (except for the leader on a straight road) such that the velocity of each vehicle remains between 0 and v max , and the

∀i ∈ I, ∀t ∈ R + 0 , v i (t) ∈ [0, v max ] and d i (t) ∈ (-∞, d max ]. (2)
The computational complexity of existing symbolic approaches grows exponentially with the dimension and thus with the number of vehicles in the platoon. A way to tackle this computational complexity is the use of continous-time assume-guarantee contracts.

ASSUME-GUARANTEE CONTRACTS

A continuous-time system is a tuple Σ = (W , X , Y , T) where W ⊆ R m , X ⊆ R n and Y ⊆ R p , are the sets of inputs, states, and outputs; and T is a set of continuous-time trajectories (w, x, y) :

R + 0 → W × X × Y . A continuous-time assume-guarantee contract for Σ is a tuple C = (A W , G X , G Y )
where A W ⊆ W is a set of assumptions on the external input and G X ⊆ X , G Y ⊆ Y are sets of guarantees on the state and output, respectively. Intuitively, an assume-guarantee contract is strongly satisfied when if the system's input belongs to A W up to a time instant t, then the system's state belongs to G X at least until t, and the system's output belongs to G Y until t + δ with δ > 0, see [START_REF] Saoud | On the composition of discrete and continuous-time assume-guarantee contracts for invariance[END_REF].

The dynamics of the platoon on straight or circular roads can be seen as a cascade composition of m - 1). Therefore, in view of Theorems 1 and 2 in [START_REF] Saoud | On the composition of discrete and continuous-time assume-guarantee contracts for invariance[END_REF] , the control objective can be achieved, on both straight and circular roads, by assigning assume-guarantee contracts C = (A W , G X , G Y ) to systems Σ i where

1 systems Σ 2 || c . . . || c Σ m or feedback-cascade composition of m systems (Σ 1 || c . . . || c Σ m ) f , respectively (see Figure
A W = G Y = [0, v max ] and G X = (-∞, d max ] × [0, v max ]. (3)
Using the previous results, if one can synthesize a controller for each vehicle in the platoon to enforce the strong satisfaction of the contract (3), then the control objective is achieved. In the following, we provide a solution to the control objective using symbolic control techniques.

SYMBOLIC CONTROLLER SYNTHESIS

Given a dynamical system described by differential equations, a symbolic model or discrete abstraction is a finite-state dynamical system related to the original system by some formal behavioral relationship (e.g. simulation, bisimulation and their alternating/ approximate versions), which makes it possible to refine a symbolic controller, designed for the abstraction, into a concrete one that can be used on the original system.

For the vehicle platooning problem, we compute a symbolic abstraction for each vehicle under assumptions of the contract, similarly to [START_REF] Meyer | Compositional abstraction and safety synthesis using overlapping symbolic models[END_REF] by discretizing state and input spaces and by computing over-approximations of the reachable sets. Then, we show that the dynamics of the vehicle is related to its abstraction at sampling instants by an alternating simulation relation [START_REF] Tabuada | Verification and control of hybrid systems: a symbolic approach[END_REF]. Using a fixed-point algorithm, we can synthesize symbolic safety controllers, which after refinement enforce the strong satisfaction of the contract for the original system. Let us remark that even if the symbolic controller acts at discrete sampling instants; the strong satisfaction of the contracts can be certified by construction.

NUMERICAL RESULTS

We consider a platoon on a circular road consisting of three vehicles with different sampling periods and different abstraction parameters. Figure 2 shows the simulation results. One can see that the specification is enforced (d max = -10m and v max = 15m/s). Interestingly, the vehicles do not distribute uniformly on the road. On this simulation, one can see that the vehicles with larger sampling period need to keep a larger distance to the front vehicle, which can be explained by the fact, that they need more time to react.

CONCLUSION

In this paper, we have presented a compositional approach to the control of a platoon of vehicles, based on a notion of continuoustime assume-guarantee contracts. This approach makes it possible to deal with heterogeneous vehicles and asynchronous sampling.

In future work, we will deal with additional quantitative objectives (comfort maximization, fuel consumption minimization, etc.).

Figure 1 :

 1 Figure 1: A platoon of 4 vehicles on a straight road (top) or on a circular road (bottom)

Figure 2 :

 2 Figure 2: Simulation results for vehicle platooning on a circular road with asynchronous sampling: inter-vehicle distance (top), velocities (bottom).

* This work has been supported by Labex DigiCosme (project ANR-11-LABEX-0045-DIGICOSME) operated by ANR as part of the program "Investissement d'Avenir" Idex Paris Saclay (ANR-11-IDEX-0003-02).