Hmida Hmida
email: hhmida@gmail.com

Sana Ben Hamida

Amel Borgi
email: amel.borgi@insat.rnu.tn

Marta Rukoz
email: marta.rukoz@lamsade.dauphine.fr

A new adaptive sampling approach for Genetic Programming

Keywords: Genetic Programming, Machine Learning, Training set Sampling, adaptive sampling, sampling frequency control, KDD intrusion detection

Genetic Programming (GP) is afflicted by an excessive computation time that is more exacerbated with data intensive problems. This issue has been addressed with different approaches such as sampling techniques or distributed implementations. In this paper, we focus on dynamic sampling algorithms that mostly give to GP learner a new sample each generation. In so doing, individuals do not have enough time to extract the hidden knowledge. We propose adaptive sampling which is half-way between static and dynamic methods. It is a flexible approach applicable to any dynamic sampling. We implemented some variants based on controlling re-sampling frequency that we experimented to solve KDD intrusion detection problem with GP. The experimental study demonstrates how it preserves the power of dynamic sampling with possible improvements in learning time and quality for some sampling algorithms. This work opens many new relevant extension paths.

I. INTRODUCTION

Evolutionary Algorithms (EA) [START_REF] Yu | Introduction to Evolutionary Algorithms, ser. Decision Engineering[END_REF]- [START_REF] Ptrowski | Evolutionary Algorithms[END_REF] are meta-heuristics that complies with a wide range of problems and especially those related to optimization and machine learning. Their flexibility and expressiveness comes with two major flaws: an excessive computational cost and a problematic parameters setting. Genetic Programming (GP) [START_REF] Koza | Genetic programming -on the programming of computers by means of natural selection, ser. Complex adaptive systems[END_REF] is the most affected EA by these drawbacks.

In supervised learning field, the lack of data leads to unsatisfactory learners. This is no longer an issue with numerous data sources and high data volume that we witness in the era of Big Data. Nonetheless, this toughens up the computation problem of GP and precludes its application in data intensive problems.

Many techniques have been deployed to tackle the computation time issue that can be grouped in two main categories: hardware acceleration based techniques and software based techniques. The most affordable is software based solutions that do not require any specific hardware configuration. Sampling is the mainstream approach in this category. It relies on reducing processing time by reducing data while keeping relevant records.

In [START_REF] Hmida | Sampling methods in genetic programming learners from large datasets: A comparative study[END_REF], we exposed a review of sampling methods used with GP. Mainly, sampling methods can be classified with regard to three properties: re-sampling frequency, sampling scheme or strategy and sampling quantity. According to resampling frequency, algorithms use a unique or a renewable sample. They are called respectively static or dynamic. On the one hand, in static sampling like the Historical Subset Selection [6] and bagging/boosting [START_REF] Iba | Bagging, boosting, and bloating in genetic programming[END_REF], [START_REF] Paris | Exploring overfitting in genetic programming[END_REF], a selection of representative training set records needs to be performed. With large datasets, this poses a problem of combining downsizing and data coverage objectives. On the other hand, dynamic sampling creates samples per generation according to its selection strategy. Consequently, GP individuals do not have enough time to learn from sampled data. The population might waste some good resources for solving some difficult cases in the current training set. Otherwise, re-sampling at each GP iteration might be computationally expensive, especially when using some sophisticated sampling strategies.

We propose, in this paper, a new extension to sampling techniques in which sample renewal is controlled through a parameter that adapts the sampling to the learning process. This extension aims to preserve original sampling strategy while making an enhancement in learning robustness and/or learning time.

This paper is organized as follows: in the first section, we expose the background of this work in GP and sampling methods. Section III introduces the novel sampling approach and gives how it can extend dynamic sampling. Then, in section IV, an experimental study gives the proof of concept of adaptive sampling and traces its effect on learning process through the discussion of registered results in section V. Finally, we give some conclusions and propose further developments.

II. BACKGROUND AND RELATED WORKS

A. Genetic Programming

As any EA, GP evolves a population of individuals through iterations called generations. Each individual represents a complete mathematical expression or a small computer program. The standard GP uses a tree representation of individuals built from a function set for nodes and a terminal set for leaves. When GP is applied to a classification problem, each individual is a candidate classifier. The terminal set is composed of dataset features and some randomly generated constants and the function set contains mostly arithmetic and logic functions. As for the fitness function, it is often based on learning performance measures such as accuracy.

The main steps of GP with dynamic or active sampling are: The evaluation is the prevailing step with regard to the overall computation cost and it depends simultaneously on the sample size, the population size and each individual complexity.

B. Dynamic Sampling

In this section we present some dynamic sampling algorithms with focus on those used in the conducted experiments in section IV.

1) Random Subset Selection (RSS): Gathercole et al.

[6] proposes a simple algorithm that selects at every generation g a record i among T records in the initial dataset with a uniform probability P i (g):

∀i : 1 ≤ i ≤ T, P i (g) = S T ( 1
)
where S is the target subset size. The sampled subset has a fluctuating size around S.

2) Dynamic Subset Selection (DSS): DSS algorithm [9]- [START_REF] Gathercole | An investigation of supervised in genetic programming[END_REF] is inspired by boosting techniques and aims to bias selection to keep difficult cases (i.e. fitness cases frequently unsolved by the best solutions) and fitness cases that have not been selected for several generations. DSS computes two measures: a difficulty degree D i (g) and an age A i (g) starting with 0 at first generation and updated at every generation. The difficulty is incremented for each classification error and reset to 0 if the fitness case is solved. The age is equal to the number of generations since last selection, so it is incremented when the fitness case has not been selected and reset to 0 otherwise. Selection probabilityP i (g) in eq. ( 3) depends on each fitness weight W i (g) (eq. ( 2)).

∀i : 1 ≤ i ≤ T, W i (g) = D i (g) d + A i (g) a (2)
where d and a are given parameters denoting respectively the difficulty exponent and the age exponent.

∀i :

1 ≤ i ≤ T, P i (g) = P i (g) * S T j=1 W j (g) (3)
3) Balanced sampling: Balanced sampling [START_REF] Hunt | Sampling methods in genetic programming for classification with unbalanced data[END_REF] aims to improve classifier accuracy by correcting the original dataset imbalance within majority and minority class instances. It has some methods based on the minority class size and thus reduces the number of instances like the methods studied in this paper. Several approaches are proposed, we summarize hereafter three sampling techniques used with GP. First, Static Balanced Sampling that selects cases with uniform probability from each class without replacement until obtaining a balanced subset of the desired size. Then, Basic Under-sampling (BUSS) (resp. Basic Over-sampling) selects all minority (resp. majority) class instances and then an equal number from the majority (resp. majority) class randomly.

4) RSS-DSS and DSS-DSS: Curry et al. conceived an extension to the DSS algorithm into a 3 levels hierarchy [START_REF] Curry | Towards efficient training on large datasets for genetic programming[END_REF]. At level 0, the data set is first partitioned into blocks that are sufficiently small to reside within RAM alone. Then, at level 1, one block is chosen from this partition based on RSS or DSS. Finally, at level 2, the selected block is considered as the full dataset on which DSS was applied for several generations. Depending on the level 1 algorithm, two approaches are possible: RSS-DSS hierarchy or DSS-DSS hierarchy.

5) RSS-TBS (TBS: Topology Based Selection): Based on the same idea, Hmida et al. proposed two new variants of hierarchical sampling : the RSS-TBS and the BUSS-RSS-TBS [START_REF] Hmida | Hierarchical data topology based selection for large scale learning[END_REF]. The RSS-TBS uses the Topology Based Subset Selection at level 2 instead of RSS or DSS. The second variant BUSS-RSS-TBS extends the first variant with a Basic Under-sampling at the level 0 block creation. BUSS favors the minority class by calculating the block size according to their cardinalities. For majority class, an equal number of instances are selected randomly.

III. CONTROLLING SAMPLING FREQUENCY A. Sampling approach

Sampling frequency (f ) is a main parameter for any active sampling technique. It defines how often the training subset is changed across the learning process. When f=1, the training sample is extracted at each generation and the sampling approach is considered as a generation-wise sampling technique.

Most of the sampling techniques applied with GP belong to this category. This is the case of the techniques described in section II-B. When f is set to 1, individuals in the current population have only one generation to adapt their genetic materials to the current environment characterized by the training sample. It is very difficult even impossible for a population to solve all cases in a training set in one generation. A higher value of f corresponds to a lower number of samples to be generated. We think that the sampling frequency must be updated according to the evolution state and the difficulty of the current training set.

The fundamental idea behind adaptive sampling by controlling the sampling frequency is to add an extra parameter to sampling algorithms acting as a moderator or re-sampling regulator. While active or dynamic methods use a fixed renewal frequency equal to 1, adaptive sampling decides to generate a new sample for the subsequent generations according to a condition that must be true.

Fig. 1 depicts this approach. GP based learner interacts with sampling process by providing adequate data about the learning process in order to perform the underlying selection strategy. For example, DSS needs to know misclassified cases to update the difficulty value. Then, the sampling algorithm delivers a new sample generated according to the updated data. With adaptive sampling, a predicate controls the re-sampling decision. We assume that any data required to evaluate this predicate is available within the data dispensed by the GP engine. Various predicates can be used for adaptive sampling that based on other parameters such as:

• generations number,

• population mean fitness,

• mean fitness improvement rate,

• best fitness improvement rate,

• etc. The straightforward approach is to define a threshold per measure and the predicate is then a comparison of the current value to the corresponding threshold.

For example, if we define a threshold of 0.002 to best fitness improvement rate then GP will continue to use the same sample if the best fitness of the current generation is better than the previous generation with 0.2%. Otherwise a new sample must be created. In a more complex approach, threshold can be auto-adapted to learning process.

We propose below some new techniques to adapt the sampling frequency. After a brief recall of the classic method considering the sampling frequency as an EA input parameter with a fixed value usually set to 1, we introduce the deterministic sampling frequency which the value evolves according to predefined schema. We then present the adaptive sampling frequency based on the current state of the population. It uses either the evolution of the mean fitness or the number of resolved cases to decide whether to create a new sample or to carry on learning with the previous one.

B. Fixed Sampling Frequency

The sampling frequency f is set before starting the GP run like any GP parameter. This value remains unchanged till the last generation. This can be represented by the following algorithm: Input: f {sampling frequency} 1: for all generation g < g max do 2:

if g mod f = 0 then An increasing frequency starts with samples having short lifetime. By the end of the run, samples are used for a larger number of generations. This can be achieved through the following steps:

1: for all generation g < g max do

2: f = (C × g) α {C, α ∈ R} 3: if g mod f = 0 then 4: re-sample 5:
end if 6: end for

The opposite process (i.e. decreasing frequency) uses the same steps but updates frequency with a decreasing function such as: f = (C × (g max -g)) α .

D. Adaptive sampling frequency

In this technique, sampling frequency value is adjusted to be more suitable to the current learning state. Therefore, it can increase or decrease by a varying amount. We assume, that less performing learners need more time to improve their performance and symmetrically efficient learners on a particular sample need to see different data from a new sample.

Adaptive sampling can rely on various learning performance indicators. It retrieves these indicators from the GP engine. Hereafter, two examples of adapting techniques.

The first uses a threshold for the population mean fitness to detect if the population is making improvements or not. In the latter case, a new sample is generated since the old one is fully exploited. Input: r{mean fitness variation rate: (mean fitness(g)-mean fitness(g-1))/mean fitness(g-1)} Input: t{threshold} 1: for all generation g < g max do 2:

if r > t then 3:

re-sample 4:

end if 5: end for The second example is based on measuring the mean number of individuals that have resolved each record in the sample. When this value reaches a designated value then new records are selected in a new sample.

In the following sections, we give details about the settings used for the conducted experiments and implementation of adaptive sampling over some dynamic sampling algorithms discussed in section II-B. Then we expose the experimental results and discuss them to analyze the effect of sampling frequency and adaptive sampling on GP performance in resolving KDD intrusion detection problem.

IV. EXPERIMENTAL SETTINGS

A. Cartesian Genetic Programming

Cartesian Genetic Programming [START_REF] Miller | Cartesian genetic programming[END_REF] is a GP variant where individuals represent graph-like programs. It is called "Cartesian" because it uses a two-dimensional grid of computational nodes implementing directed acyclic graphs. Each graph node encodes a function from the function set.

CGP shows several advantages over other GP approaches. Unlike trees, there are more than one path between any pair of nodes. This enables the reuse of intermediate results. A genotype can also have multiple outputs which make CGP able to solve many types of problems and classification problems in particular [START_REF] Harding | Implementing cartesian genetic programming classifiers on graphics processing units using gpu.net[END_REF]. Otherwise, CGP has the great advantage of counteracting the bloating effect (genotype growth), frequent phenomena with other GP representations. CGP is easy to implement, and it is highly competitive compared to other GP methods.

1) CGP settings: The design of CGP parameters used in this work is summarized in Table I. In this work, the parameter tuning is not fully explored.

B. Dataset

The KDD-99 intrusion detection problem consists in classifying connections into normal or attack classes. It uses a large dataset called 10% KDD-99 dataset [START_REF]Kdd cup[END_REF]. The dataset is already divided into training and test sets which are presented in Table III. Each fitness case is described by 41 features. The original data is preprocessed with the following steps:

• Transforming discrete nominal attributes to numeric values,

• Scaling data using MinMax scaler:

X sc = X -X min X max -X min ,
• Binarization of attack classes: the problem is converted to a binary classification problem with a 'Normal'class and 'Attack'class. The original four attack types (Dos, Probe, R2L and U2R) are fused in a single class.

C. Performance Metrics

We recorded, for each run, its accuracy (eq. ( 4)) and False Positive Rate (FPR) (eq. ( 3)) to measure the learning performance on both training and test sets. We also recorded the learning time measuring the computational cost.

Accuracy =

T rue P ositives + T rue N egatives T otal patterns

F P R = F alse P ositives F alse P ositives + T rue N egatives .

(5)

D. Framework

Software framework: Among several evolutionary computation frameworks, Sean Luke's ECJ [START_REF] Luke | Ecj homepage[END_REF] was used in this work to implement and test the CGP. It's an open source framework written in Java and benefits of many contribution packages as the one used here for implementing Cartesian GP developed by David Oranchak [START_REF] Cgp | Cartesian gp website[END_REF]. This framework provides a very flexible API using parameter files well documented in the ECJ owner's manual.

Hardware framework: Experiments are performed on an Intel i7 -4810M Q (2.8GHZ) workstation with 8GB RAM running under Windows 8.164 -bit Operating System.

E. Sampling settings

In the first set of experiments, we tested six values for the sampling frequency: 1, 10, 20, 30, 40 and 50 on four sampling methods BRSS, BUSS, DSS and RSS. BRSS is Balanced RSS. It is an RSS variant where the random sample is balanced according to a given ratio between problem classes.

In the second part, we implemented four different techniques to control sampling frequency. Two deterministic techniques and two adaptive ones as follows:

• Deterministic+: deterministic controlling with the increasing function f = (2 × g) 0.5 , • Deterministic-: deterministic controlling with the decreasing function f = (2 × (g max -g)) 0.5 , • Average Fitness: adaptive controlling based on population average fitness with a threshold of 0.001, • Min Resolved: adaptive controlling based on the average population proportion of the population representing the individuals that resolved sample records. We use a minimum threshold of 0.5. The underlying active sampling algorithms have their own parameters described in Table IV.

V. RESULTS AND DISCUSSION

This study is organized in 2 types of experiments. The first aims to study the impact of sampling frequency variation on learning time and learning performance indicators.

For each value of re-sampling frequency and controlling technique, 21 runs of each sampling algorithm are conducted. We report the mean learning time of each configuration and accuracy and FPR of the best individual.

A. Effect of sampling frequency Fig. 2 illustrates learning time variation on re-sampling frequency for the studied algorithms. Fig. 3 and Fig. 4 show the effect of sampling frequency on two learning quality measures accuracy and FPR. The shape of the curves in Fig. 2 reveals two distinct behaviors when the sampling frequency increases. The first concerns the BUSS, BRSS and RSS algorithms that have recorded an insignificant decreasing variation of the average learning time. The second behavior is that of DSS with a more remarkable decrease.

Indeed, for a frequency equal to 1, 200 sampling operations are made and for a frequency of 50 this number goes down to 4. Time saving depends on the time needed to perform sample creation with respect to the time spent for a whole generation. This is why the decrease in time is not very important since population evaluation is the predominant step in the learning time for GP.

DSS algorithm differs from other algorithms by updating certain parameters (age and difficulty). But this does not affect the learning time since it is performed at each generation regardless of the re-sampling frequency. In DSS, the selection of a learning case requires the calculation of a probability based on the age and difficulty values of the whole dataset and subsequently a greater time than the other techniques. This is at origin of the difference in learning time saving. From Fig. 3, the variation of the rate of accuracy, contrary to expected results, did not make any improvement for the 2 sets. Although, there is an increase of Accuracy at frequency 50, for DSS and BRSS for the test set. However, this remains irregular and can not be generalized.

The most noticeable shift is that of the FPR (Fig. 4). However, no empirical correlation with the variation of the re-sampling frequency for BRSS, DSS and RSS can be made. Only BUSS realizes a decrease in FPR value decrease when re-sampling frequency increases for training and test sets. dynamic algorithms (Dynamic).

B. Adaptive sampling

For the mean learning time, results from Fig. 5 confirm the previous remarks about this behaviour in section V-A. It is only learning time of DSS which is clearly affected by introducing any controlling technique whether it is deterministic or adaptive. This effect is more or less strong according to reducing samples number used per GP run.

As regards learning performance on the test set, results do not rise a general response to controlling sampling frequency by means of the 4 studied techniques.

For BUSS algorithm, accuracy records an infinitesimal improvement with controlled frequency techniques. In addition to that, FPR made a similar change in the same direction (Fig. 6).

For each sampling method, there is at least one frequency control technique that allows it to improve, largely or slightly, its learning performance. For example, for the RSS, the 'Deterministic+' and 'Min Resolved' techniques have allowed the accuracy to move from values around 80% to values greater than 90%. Similarly, the two deterministic methods (deterministic + and deterministic-) and the Adaptive method 'Min Resolved' have improved significantly the accuracy values for DSS sampling, whose values increased by 10% to 12%. However, a deterioration of the FPR has been recorded.

For the BUSS, only the adaptive controlling technique 'Average fitness' improved the accuracy, but the other techniques do not impact its performance comparing to the generationwise frequency. In fact, the control of the sampling frequency, if it does not allow an improvement of the results, does not generate a deterioration of the performance, except for some cases for the FPR measure.

VI. CONCLUSION

The idea of adaptive sampling has its origins in the assumption that a population needs enough time (generations) to adapt to learning data. We have formalized this principle in the form of a re-sampling predicate that is based on a condition to decide whether to modify the current sample or to keep it.

Experiments are limited to testing adaptive sampling by controlling sampling frequency with simple predicates. The results showed a slight effect on learning time. This effect is in the direction of a decrease but with different degrees depending on the sampling method.

With adaptive sampling, computational cost can be improved according to the underlying dynamic sampling algorithm only if the selection process is time consuming as it is for DSS.

The used predicates did not deteriorate the learning accuracy. However, it did not have a proven positive and generalizable effect on the quality of learning. Thus, they need to be refined.

Many new research paths emerges from this study that are worthy of further investigation. A first path is the exploration of other predicates that take into account the characteristics of the training dataset and the underlying problem to find more relevant predicates for GP classifier improvements. A second one is to extend the scope of adaptive sampling to other sample properties. For instance, an adaptive sampling can downsize or upsize the sample instead of generating a new one. We may also combine several sampling strategies and algorithms in a single method. Then, according to the learning evolution, a sample is generated using the suitable strategy in an interleaved way: we use a different algorithm at each time we need to create a new sample.

Fig. 1 .

 1 Fig. 1. Adaptive vs Dynamic sampling

  end for C. Deterministic sampling frequency When a deterministic control of the sampling frequency is applied, the frequency f gets different values during the GP run. It is guided by a function that delivers the same values at each run. Thus, frequency can increase, decrease or have a complex curve.

Fig. 2 .

 2 Fig. 2. Variation of the mean learning time with the re-sampling frequency.

Fig. 3 .Fig. 4 .

 34 Fig. 3. Variation of the best individual accuracy with the re-sampling frequency.

Figures 5 Fig. 5 .

 55 Figures 5 and 6 report the experimental results of 4 sampling algorithms extended with 4 different sampling frequency techniques (section IV-E): deterministic based on an increasing or decreasing function (Deterministic+ and Deterministic-)

Fig. 6 .

 6 Fig.6. Variation of the accuracy and FPR measures according to the re-sampling frequency controlling strategy.

TABLE I CGP

 I PARAMETERS.

	Parameter	Value
	Population size	256
	Sub-populations number	1
	Generations number	200
	CGP nodes	300
	Inputs	49
	Outputs	1 (2 classes)
	Tournament size	4
	Crossover probability	0.9
	Mutation probability	0.04
	Fitness	Minimize classification error
	2) Terminal and function sets: The terminal set includes
	41 features of the benchmark KDD-99 dataset. The function
	set includes basic arithmetic, comparison and logical operators
	reaching 17 functions (table II).

TABLE II TERMINAL

 II AND FUNCTION SETS FOR GP.

	Function (node) set
	Arithmetic operators:	+, -, * , %
	Comparison operators:	<, >, <=, >=, =
	Logic operators:	AND, OR, NOT, NOR, NAND
	Other:	NEGATE, IF (IF THEN ELSE),
		IFLEZE(IF <=0 THEN ELSE)
		Terminal set
	KDD-99 Features	41
	Random Constants	8 in [-2, 2[

TABLE IV COMMON

 IV SAMPLING PARAMETERS.

	Method	Parameter	Value
	All (except BUSS)	Target Size	5000
	BRSS	Balancing method	Full dataset distribution
	BUSS	Target size	416
	DSS	Difficulty exponent Age exponent	1 3.5


Tristan Cazenave
email: cazenave@lamsade.dauphine.frsbenhami@u-paris10.fr

Sana Ben Hamida

Forecasting Financial Volatility Using Nested Monte Carlo Expression Discovery

We are interested in discovering expressions for financial prediction using Nested Monte Carlo Search and Genetic Programming. Both methods are applied to learn from financial time series to generate non linear functions for market volatility prediction. The input data, that is a series of daily prices of European S&P500 index, is filtered and sampled in order to improve the training process. Using some assessment metrics, the best generated models given by both approaches for each training sub sample, are evaluated and compared. Results show that Nested Monte Carlo is able to generate better forecasting models than Genetic Programming for the majority of learning samples.

I. INTRODUCTION

Derivative markets provide market participants and policymakers working in the financial markets with a rich source of information for gauging market sentiment. Option prices are especially useful for extracting such information. Due to their forward-looking nature, they efficiently encapsulate market perceptions. The price of an option therefore depends on the market opinion about the future volatility of the underlying asset upon which the option is written. The volatility that makes the theoretical option price equal to the observed option price is called the implied volatility. It may contain useful information about the market participant expectations regarding the future volatility of the underlying asset. It has been conjectured that the implied volatility is a good predictor of the future volatility of the underlying asset [START_REF] Corrado | The forecast quality of cboe implied volatility indexes[END_REF]. This effectively makes volatility forecasting a big challenge for many financial institutions around the word.

Traditional parametric methods have limited success in estimating and forecasting the implied volatility as they are dependent on restrictive assumptions and difficult to estimate. Several machine learning techniques have been recently used to overcome these difficulties. Genetic programming (GP) has been often applied to forecast financial time series and recently, it was successfully applied for the prediction of the implied volatility [START_REF] Abdelmalek | Selecting the best forecasting-implied volatility model using genetic programming[END_REF], [START_REF] Hamida | Applying dynamic training-subset selection methods using genetic programming for forecasting implied volatility[END_REF]. This paper investigates the application of a new metaheuristic for the implied volatility forecasting, that is the Nested Monte Carlo Search, and compares it to Genetic Programming.

Monte Carlo Tree Search (MCTS) is a general search algorithm that made a breakthrough in computer Go [START_REF] Coulom | Efficient selectivity and backup operators in Monte-Carlo tree search[END_REF]. The standard MCTS algorithm for multi-player games is Upper Confidence bounds applied to Trees (UCT) [START_REF] Kocsis | Bandit based Monte-Carlo planning[END_REF]. A good MCTS algorithm for single player games is Nested Monte Carlo Search (NMCS) [START_REF] Cazenave | Nested Monte-Carlo Search[END_REF].

The principle of Monte Carlo Search algorithms is to play a lot of random games in order to evaluate a state. A random game is also called a playout. When using UCT the average of the playouts that start with a move is used to select the next move of the simulation. The key principle of UCT is to balance exploration and exploitation. Exploitation favors the moves with the greatest means while exploration favors the moves that have not been tried much.

The principle of NMCS is different from the principle of UCT. It still uses a lot of playouts in order to find good sequences of moves but instead of using the average of the playouts, it memorizes the best sequence of moves found so far and follows it. It uses nested levels of search and tries every possible moves once. Memorizing the best sequence of moves enables it to intensify the search, while trying all possible moves once enables it to diversify the search.

NMCS has been applied to puzzles such as SameGame and Morpion Solitaire and to optimization problems. It has found world records for difficult problems such as the Snake-In-The-Box problem [START_REF] Kinny | A new approach to the snake-in-the-box problem[END_REF] or the Weak Schur problem [START_REF] Eliahou | Investigating monte-carlo methods on the weak schur problem[END_REF]. In order to find these records the random playout policy has been replaced with a more informed policy that still includes some randomness.

NMCS has also been applied to transportation problems such as the Bus Regulation problem [START_REF] Cazenave | Monte-Carlo bus regulation[END_REF] or the Traveling Salesman with Time Windows problem [START_REF] Rimmel | Optimization of the nested monte-carlo algorithm on the traveling salesman problem with time windows[END_REF]. It has also been used for software engineering so as to generate structured test data with specific properties [START_REF] Poulding | Generating structured test data with specific properties using nested monte-carlo search[END_REF].

A related algorithm is Nested Rollout Policy Adaptation (NRPA) that learns a playout policy online using nested levels of learning [START_REF] Rosin | Nested rollout policy adaptation for Monte Carlo tree search[END_REF]. NRPA has beaten human world records that were standings for more than thirty years at Morpion Solitaire. It has also been successfully applied to optimization problems such as the Traveling Salesman with Time Windows problem [START_REF] Edelkamp | Algorithm and knowledge engineering for the TSPTW problem[END_REF] and the Multiple Sequence Alignment problem [START_REF] Edelkamp | Monte-carlo tree search for the multiple sequence alignment problem[END_REF].

The remainder of the paper is divided into six sections. Section 2 gives a brief introduction to the symbolic regression and its application to the prediction of the implied volatility in finance. Section 3 describes the design of the Nested Monte Carlo search. Section 4 reminds the basics of the Genetic 978-1-4799-7560-0/15/$31 c 2015 IEEE Regression. Section 5 and 6 report the experimental settings and the results.

II. SYMBOLIC REGRESSION FOR VOLATILITY

FORECASTING

There are two approaches to generate volatility forecasts. One is to extract information about the variance of future returns from their history, the second is to elicit market expectations about the future volatility from observed option prices. The Black&Scholes(BS) [START_REF] Black | The pricing of options and corporate liabilities[END_REF] is the most commonly used model in the second category. It uses a backward induction technique to derive the implied volatility from a set of observable data for an option contract. 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.

Whilst it is the most commonly used in financial markets, the BS pricing model has some well-known deficiencies. Indeed, it assumes that volatility is 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. Black [START_REF] Black | The pricing of options and corporate liabilities[END_REF] suggested that the non-stationary behavior of volatility could lead the BS model to over-price or under-price options... To overcome the BS deficiencies, several parametric models have been introduced in financial engineering based on the realized volatility. They can be classified into three categories: models based on Past standard Deviations (such as the Random walk model, the Simple Regression method, ...), ARCH Class Conditional Volatility models (such as the well known GARCH model family [START_REF] Engle | Autoregressive conditional heteroscedasticity with estimates of the variance of united kingdom inflation[END_REF], [START_REF] Bollerslev | Generalized autoregressive conditional heteroskedasticity[END_REF]) and the stochastic models (such as Monte Carlo models).

All these models are parametric models that need some assumptions. However, as mentioned in [START_REF] Ma | An engineering approach to forecast volatility of financial indices[END_REF], the volatility forecasting approaches should be free of strong assumptions and more flexible than parametric methods.

Recently, to provide a new issue for implied volatility forecasting, several researchers have studied the application of supervised machine learning techniques. The aim is to offer explicit formulas which could compute directly the implied volatility expressed as a function of option prices and other observable variables.

The main objective of the supervised machine learning is to draw a mapping between input and output values in a training data set. The hope is that the resulting pattern (model) from the mapping process be able to correctly predict unknown examples of the same type.

In the case of learning regression models, the task is to discover a target function f (X) that map a vector of realvalued inputs X to real valued target variable Y . For financial volatility forecasting, the input set X is a financial time series data (X = X t1 , X t2 , ..., X tN ), and Y is the BS implied volatility computed from the observed data. Each input X ti is a vector of variables for one observed option at the date ti.

The variables could be the strike price, the maturity, the call or put price, ...

III. MONTE CARLO EXPRESSION DISCOVERY

Monte Carlo sampling can be used to generate expressions. Expressions are seen as stacks of atoms with a maximum size. At each step of a sample (also called a playout), a random atom is chosen and added to the stack. In order to stop the playout a maximum number of atoms are allowed in the stack. When the number of atoms of the stack plus the number of open branches is equal to the maximum, only leaves are added to the stack so that when the stack is complete it contains less than the maximum number of atoms.

Nested levels of Monte Carlo search can be used in order to improve on random sampling [START_REF] Cazenave | Nested Monte-Carlo Search[END_REF]. The principle of NMCS is to use different levels of search. At level zero it does a playout. At greater levels it performs informed playouts. At each step of a playout, it generates the possible moves and for each move it does a level -1 search. It keeps the move that is associated to the best level -1 search and plays it. It continues to choose moves like this until the end of the informed playout. The NMCS algorithm is given in algorithm 1. For a playout of level n it chooses the move that gives the best score when followed by a playout of level n -1. A search at level n calls a search at level n -1 to direct its search.

NMCS was applied with success to multiple expression discovery problems [START_REF] Cazenave | Nested monte-carlo expression discovery[END_REF], [START_REF]Monte-carlo expression discovery[END_REF]:

• The prime generating polynomial problem where the goal is to discover a polynomial that generates as many different primes in a row as possible. Cartesian Genetic Programming has also been applied to the problem of finding polynomials that generate prime numbers [START_REF] Walker | Predicting prime numbers using cartesian genetic programming[END_REF].

• The finite algebra problem where the goal is to find a discriminator term for an algebra. The discriminator term is a formula that gives special results for the possible values of the variables of the algebra. Genetic Programming has also been applied to finding a discriminator term for finite algebras [START_REF] Spector | Genetic programming for finite algebras[END_REF].

• The parity problem which is a classical problem for GP [START_REF] Koza | Genetic Programming: On the Programming of Computers by Natural Selection[END_REF]. The goal is to find a formula that tells if the number of bits of the input is odd or even. The problem can be easy or difficult depending on the atoms used to discover expressions and the size of the input.

• The N-prisoners puzzle. N prisoners wear a hat and each prisoner can see the other's hats but not his hat. A hat is associated to either a zero or a one. All prisoners are asked to find the color of their hat. They can pass and give no answer. If all prisoners that have not passed have found the right number for their hats, all the prisoners are free, else they stay in prison. A basic strategy is to have one prisoner randomly choosing its answer. It wins 50 percent of the time. Much more elaborate strategies exist. Genetic Programming has also been applied to the N-prisoners puzzle [START_REF] Burke | A puzzle to challenge genetic programming[END_REF].

• The MAX problem is to find an expression that evaluates as great as possible with a limited expression size. The authorized atoms are 0.5, + and *. In [START_REF] Langdon | An analysis of the max problem in genetic programming[END_REF] the limit was on the depth of the corresponding tree.

For NMCS the limit was on the number of atoms of the generated expression, not on the depth of the tree.

• The target number problem is derived from a popular television game where the goal is to find a an expression that evaluates as close as possible to a target number given the four basic operations as atoms as well as a given set of predefined numbers.

A good property of NMCS for the generation of expressions is that it uses a maximum expression size. It avoids a problem of Genetic Programming called bloat which consists in generating bigger and bigger expressions that are very specific to the learning data and that hinder generalization of the learned expression. Another good property of NMCS is that it there is no parameter to tune except the maximum expression size. Even without tuning it has outperformed optimized Genetic Programming systems for the problems previously cited. It is still possible to tune NMCS for a problem by improving the playout policy. For most problems that were addressed with NMCS, big improvements came from using rules during the playouts so as to avoid bad moves. It could be possible to do this for Nested Monte Carlo expression discovery, however that would be a problem specific enhancement. In future work we will consider domain specific playout policies.

For a problem of height h, the NMCS algorithm of level l has a complexity of O(|branchingF actor| n h l+1 ) [START_REF] Méhat | Combining UCT and Nested Monte Carlo Search for single-player general game playing[END_REF].

IV. GENETIC REGRESSION Genetic Programming, that is an offshoot of Genetic Algorithms, is essentially applied to develop efficient computer programs. For the implied volatility forecasting, we used a variant of Genetic Programming called Genetic Regression that implement the Symbolic Regression. This is a numerical optimization tool to select a model which best fits a time series. With Genetic Regression, the desired program is a function that relates a set of inputs to one output. In our case, the output is the implied volatility.

Genetic Regression was one of the earliest applications of GP published by John Koza on 1992 [START_REF] Koza | Genetic Programming: On the Programming of Computers by Natural Selection[END_REF], and has continued to be widely studied. The GP system takes as input a set of explanatory or lagged dependent variables (defined as terminals in GP), and a set of mathematical operators (defined as non-terminals in GP) and then randomly selects and combines independent variables and operators in search for a model specification that would satisfy some user-defined fitness function. A GP solution (individual) is a symbolic regression program able to generate predicted values of the dependent variable. The individual possessing the lowest predictive error is characterized as fittest and is kept in memory.

Several works have applied the Genetic Regression for financial times series mining. For example, for volatility forecasting, Chen and Yeh [START_REF] Chen | Using genetic programming to model volatility in financial time series: the cases of nikkei 225 and s&p 500[END_REF] have applied a recursive genetic programming (RGP) approach to estimate volatility historical returns of Nikkei 225 and S&P500 index. Zumbach and al. [START_REF] Zumbach | Genetic programming with syntactic restrictions applied to financial volatility forecasting[END_REF] and Neely and Weller [START_REF] Neely | Predicting exchange rate volatility: Genetic programming vs. garch and risk metrics[END_REF] have applied GP to forecast foreign exchange rate volatility. Recently, Abdelmalek et al. [START_REF] Abdelmalek | Selecting the best forecasting-implied volatility model using genetic programming[END_REF] and Ben Hamida et al. [START_REF] Hamida | Applying dynamic training-subset selection methods using genetic programming for forecasting implied volatility[END_REF] have extended the studies mentioned earlier by forecasting the implied volatility from the S&P500 index call options instead of historical volatility. For this work, we use the same Genetic Regression implementation as in [START_REF] Abdelmalek | Selecting the best forecasting-implied volatility model using genetic programming[END_REF], [START_REF] Hamida | Applying dynamic training-subset selection methods using genetic programming for forecasting implied volatility[END_REF].

The steps necessary to implement the Genetic Regression are summarized in algorithm 2.

Algorithm 2 Abstract Genetic Regression Algorithm

Randomly create an initial population of models while termination condition not satisfied do Evaluate the performance of each function according to the fitness criterion Select the best functions in the population using the selection algorithm Generate new individuals by crossover and mutation end while Report the best solution found

The starting GP population is generated using a ranked half and half method [START_REF] Koza | Genetic Programming: On the Programming of Computers by Natural Selection[END_REF]. This method involves generating an equal number of trees using a pre-established maximum depth. Based on the fitness criterion, the selection of the individuals for reproduction is done with the tournament selection algorithm.

The main genetic operator involved in GP are:

• Crossover: It create an offspring function by recombining randomly chosen parts from two selected parent functions.

• Mutation: Different mutation operators are used. Point mutation operator consists of replacing a single node in a tree with another randomly-generated node of the same arity. Branch mutation operator randomly selects an internal node in the tree, and then it replaces the sub-tree rooted at that node with a new randomlygenerated sub-tree. Expansion mutation operator randomly selects a terminal node in the tree, and then replaces it with a new randomly-generated sub-tree. Each mutation is applied with a rate given in the GP parameters (table II).

The method of replacing parents for the next generation is comma replacement strategy, which selects the best offspring to replace the parents. It assumes that offspring size is higher than parents' size. The stopping criterion is the maximum number of generations (table II).

To run a Genetic Regression system, four ingredients need to be specified: the terminal set, the function set, the fitness measure and the GP parameters. Definition and details of these ingredient are given in the section V. When using the symbolic regression technique, developers and users of learning meta-heuristics believe that a small sample size is usually sufficient, especially with GP. Thus, the full input sample is sorted by time series and divided chronologically into nine successive sub-samples (S 1 , S 2 , • • • , S 9 ) each containing 667 daily observations. These samples will be used simultaneously for training and test steps. Note that the training set is the data sub-sample on which the program operates to determine the relationship between input parameters and output, and the test set is the data sub-sample on which the resulting formula is tested.

If a sub-sample S i is used as a training set, then the subsample of the the immediately following date S i+1 (S 1 if i = 9) is used as test set.

B. Settings

Both GP and NMCS used the same terminal and function sets to generate and evolve the tree based models. The terminal set includes input variables: the call option price divided by strike price C/K, the index price divided by strike price S/K and time to maturity τ (table I). The function set includes basic mathematical operators, such as arithmetic operators, and Black&Scholes components such as the normal cumulative distribution function (N cf d) [START_REF] Black | The pricing of options and corporate liabilities[END_REF], which may be useful for implied volatility models (table I).

1) NMCS parameter:

There is only one NMCS parameter which is the maximum number of nodes in a tree representing an expression. We used a maximum of 40 nodes.

2) GP parameters: The implementation of genetic programming involves different parameters that determine its efficiency. A common approach in tuning GP is to undertake a series of trial to make parameter choices for the final GP runs. We used this approach to determine some parameters such as the crossover and mutation probabilities.

The final design of GP parameters used in this work is summarized in Table II

C. Evaluation metrics

Each formula given by NMCS or GP is evaluated to test whether it can accurately forecast the output value (implied volatility) for all entries in the training set. To assign a fitness measure to a given solution, we compute the Mean Squared Error (MSE) between the estimated volatility ( ŷi ) and the target volatility (y i ):

M SE = 1 N 1 N (y i -ŷi ) 2
where N is the number of entries in the training data sample.

Beside the (MSE), we used another measure to evaluate the training accuracy and the generalization ability of the obtained models. This measure is the percentage of Poor Fitted Observations (PFO). It compute the proportion of input entries where the gap between the estimated and the target values is greater than a given level (l).

When using mining techniques in the field of automatic learning, two difficulties are often encountered: the over-fitting and the under-fitting. Over-fitting happens when the data is learned by rote and when the learned expression does not generalize past the learned data. Under fitting is encountered when the regression model does not fit correctly a part of the learning data, and the forecasting error of the corresponding cases is relatively high. The predictive ability of the model is then not satisfactory, mainly for finance prediction. So it is necessary to introduce an other metric able to detect these case more effectively than the mean square error. The PFO is then introduced in this purpose and is computed as follows.

P F O = ( N i=1 z i )/N ; z i = { 1 if f e( ŷi ) > l, 0 otherwise (1)
f e( ŷi ) = | ŷi -y i y i |; (2)
where N is the size of the input data and l is the maximum forecasting error level allowing for the prediction of one observation X i .

The PFO value is used to estimate the percentage of variables in the input data that a given formula was unable to forecast without exceeding a fixed maximum error margin that could be accepted in the prediction. It is related to the closeness of the forecasted values against the target ones. The lower the values obtained by this metric, the better is the forecasting technique.

VI. RESULTS AND DISCUSSION

For the comparative study, a set of eighteen generated models is constructed.

S1 S2 • • • S9 NMCS model M N1 M N2 • • • M N9 GP model M G1 M G2 • • • M G9
To compare both forecasting and generalization ability of each model, the MSE according to the whole input data ("MSE Total") is computed for each best solution given by NMCS and GP using each sub-sample in the learning set. Results are shown in figure 3. It can be observed from this figure that the forecasting performance of all models is quite similar, except for the three first models. Indeed, training on the three first subsamples S 1 , S 2 and S 3 don't give solutions with high accuracy, and this is for both methods.

Otherwise, a comparison of the "MSE Total" of the two methods shows that NMCS is little better than GP, except for patterns obtained with the training set S8. Performance of the generated NMC and GP volatility models according to "MSE total" Most of the models given by the two methods have a similar (and good) performance according to the MSE metric. It is therefore interesting to see if they present the same performance according to the PFO metric. Figure 4 and 5 illustrates these results for two maximum error levels:l = 0.5 and l = 0.3. Performance of the generated NMC and GP volatility models according to the PFO percentage for l = 0.5. It appears also, through the figures 4 and 5, that the MSE metric is not sufficient to compare in detail the generated models. Some solutions might have very low MSE, but they present a poor performance when a high precision for the prediction step is required, which is the case generally in finance. For example, the model MN1 presents high performance according to the MSE Total and it seems better than the model MN7. However, its PFO percentage for both levels l = 0.5 and l = 0.3 is much higher than that of MN7. In this case, for more accurate prediction, MN7 is preferred to MN1.

For a complete comparative study between NMCS and GP, the PFO percentage for all generated models are computed according to different maximum error levels (l) varying from 0.3 to 1. These values are illustrated in Figure 6. It is clear that the percentages of PFO according to the enlarged data sample are lower for the NMCS models than GP models for all values of l in the e interval test. We can already make a first conclusion that the learning capacity of NMCS from financial time series is better than that of GP. For a more detailed comparative study, only four models are selected: the two NMCS best models and the two GP best models according to the MSE Total and the PFO percentage. The selected models are given in The performance of the four best models is compared against three criteria: 1) PFE percentage computed according to three different maximum error levels (0.3, 0.4 and 0.5); 2) Shape of the forecasting error of all observations (options) in the input data; 3) Number of PFO in each sub-samples for the three levels 0.3, 0.4 and 0.5.

The figure 7 shows PFO percentage for the first comparison criterion. Since the selected models are the best given by the two methods, their prediction accuracy are quite similar, specially MN4 and MG4. For the second comparison criterion, we have represented in figure 8 the pattern of the forecasting errors of the retained models. Each point represents the squared error in implied volatility units for one option (observation) in the data set. More monotonous the squared error pattern is, more accurate the corresponding model is. It is therefore clear that the dispersion of the squared errors given by the models MN7 and MG7 presents more inhomogeneity than that of the models MN4 and MG4. They have some difficulties in forecasting options which observed date is far from the training options' dates. This observation could be confirmed by the figures 9, 10 and 11. The detailed comparison of the best models showed that the model MN4 given by NMCS is the more accurate for the volatility forecasting. The MG4 could also be considered as an accurate model, but with a little less performance than the MN4.

The formula of the best model MG4 is:

M G4 = exp (ln(φ( C K )) × τ -2 × C K + S K ) -cos( C K ) (3)
A detailed examination of the formula for MG4 shows that the implied volatility is function of all the inputs used, namely the option price divided by strike price ( C K for calls and C K for puts), the index price divided by strike price ( C K )and time to maturity (τ ). The implied volatility generated for calls and puts cannot be negative since they are computed using the square root and the normal cumulative distribution functions as the root nodes.

A. Discussion

Both GP and NMCS were able to give accurate predictions of the implied volatility. However, the "MSE Total" plot suggest that NMCS seems to produce slightly less error than GP. Indeed, the total MSE is higher for GP models than for NMCS models, excepted for M G 8 . According to the total MSE, the model built using NMCS method outperforms the model built using GP.

However, to select a model for prediction, the analysis should be related to efficiency, generalization, stability, and prediction ability of all models. Thus, models are compared also according to the PFO values. This comparison shows that the PFO scores are better for NMCS than for GP.

Another interesting aspect between GP and NMCS regression is related to the trained model generalization capabilities. As we have stated before, the purpose of a forecasting model is to capture important features of real phenomena in order to perform a good prediction about its future behavior. In this sense, models are compared according to the second evaluation metric, that is the PFO. According to this criteria, NMCS could generalize with higher efficacy than GP should we use this measure during the learning phase.

Otherwise, regardless the good quality of the models obtained by the two methods, the use of NMCS for the implied volatility forecasting was simpler and faster than the use of Genetic Programming. Indeed, the latter requires a phase of parameter adjustment. This step is not needed with NMCS. One must note that a choice of unsuitable values for GP parameters could have a negative effect on the quality of the results. This risk does not exist with NMCS.

VII. CONCLUSION

In this paper Nested Monte Carlo Search is applied to the volatility forecasting problem and compared to Genetic Programming.

Results presented show that NMCS excels in the volatility forecasting and outperforms Genetic Programming for this task. It could generate efficient models with high prediction ability. However, model evaluation can also use measures other than the Mean Squared Error. Using the Poor Fitted Observations measure on expressions learned with the Mean Squared Error fitness function, NMCS gives better results than GP.

The ability of NMCS to avoid over-fitting by limiting the size of a learnt expression is a key factor for improving its generalization ability over GP. Moreover it is a simple and parameter free method that gives good results.

In future work we will apply NMCS to other problems in finance and economy and improve it so as to cope with large amounts of data. We will also investigate the use of playout policies for expression discovery in order to improve even more the effectiveness of NMCS. Furthermore, the method will be compared, besides the GP, with the classic Black&Scholes method.

Fig. 1 .

 1 Fig. 1. At each step of the principal playout shown here with a bold line, a NMCS of level n performs a NMCS of level n -1 (shown with wavy lines) for each available move and selects the best one. At level 0, a simple pseudo-random playout is used.

  Data used to perform our empirical analysis are daily prices of European S&P500 index call options, from the CBOE for the sample period January 02, 2003 to August 29, 2003. S&P500 index options are among the most actively traded financial derivatives in the world.

Figure 2 Fig. 2 .

 22 Figure 2 illustrates the evolution of the mean values of the S&P500 index implied volatility computed with the Black&Scholes formula [15].

  It contains the NMCS and GP best solutions for each training sub-sample (S 1 , • • • , S 9 ), selected according to the test M SE. Obtained models are noted M N 1 , • • • M N 9 for NMCS and M G 1 , • • • , M G 9 for GP as described in the following table.

Subsample

Fig. 3 .

 3 Fig. 3.Performance of the generated NMC and GP volatility models according to "MSE total"

Fig. 4 .

 4 Fig. 4.Performance of the generated NMC and GP volatility models according to the PFO percentage for l = 0.5.

Fig. 5 .

 5 Fig. 5.Performance of the generated NMC and GP volatility models according to the PFO percentage for l = 0.3.The figures 3, 4 and 5 show clearly that the NMCS outperform the GP for the volatility forecasting problem. Indeed,

Fig. 6 .

 6 Fig. 6. Evolution of the total PFO percentage for l = varying from 0.3 to 1.

Fig. 7 .

 7 Fig. 7. Variation of the PFO percentage of the four best generated NMC and GP volatility models for a maximum error level equal to 0.3, 0.4 and 0.5.

Fig. 8 .

 8 Fig. 8. Pattern of the forecasting errors of the models M N 4 , M N 7 , M G 4 and M G 7 . Each point represents the squared error in implied volatility units for one option (observation) in the data set.

Fig. 9 .

 9 Fig. 9. Variation of the PFO values of the four best generated NMC and GP volatility models computed for each sub-sample with l = 0.5.

Fig. 10 .

 10 Fig. 10. Variation of the PFO values of the four best generated NMC and GP volatility models computed for each sub-sample with l = 0.4.

Fig. 11 .

 11 Fig. 11. Variation of the PFO values of the four best generated NMC and GP volatility models computed for each sub-sample with l = 0.3.

TABLE I .

 I TERMINAL AND F FUNCTION FOR NMCS AND GP TREES.

  .

	Population size	100
	Offspring size	200
	Maximum number of generations	500
	Maximum depth of new individual	6
	Maximum depth of the tree	17
	Tournament size	4
	Crossover probability	0.6
	Mutation :
	Branch mutation	0.2
	Point mutation	0.1
	Expansion mutation	0.1

TABLE II .

 II GP PARAMETERS.

  table III with their corresponding "MSE Total".

	Model	M N4	M N7	M G4	M G7
	MSE Total	0.000905	0.001702	0.001444	0.0019

TABLE III .

 III "MSE TOTAL" OF THE FOUR BEST MODELS GIVEN BY BOTH METHODS NMCS AND GP.


Marcel Bourgade
email: mbourgade@hotmail.com

Contribution for a revision of the taxonomic list of Spherodactyls (Ti-Mabouya) of Martinique by

Keywords: Eastern Caribbean, Martinique, dispersion center, Herpetofauna, Squamates, taxonomic list, Spherodactyls, Dwarf gecko, Sphaérodactylus, S. vincenti

Contribution for a revision of the taxonomic list of Spherodactyls (Ti-Mabouya) of Martinique

INTRODUCTION

Geological history and geosciences are key to understanding the general evolution of biodiversity. The herpetofauna of the islands of the eastern Caribbean and particularly the small amphibians and squamates living in ecological niches on soils requiring litter, are then among the first occupants of land surfaces emerged. And so the geological evolution -tectonic and volcanic -of Martinique and the Caribbean archipelago is decisive for the process of dispersions and therefore of speciations of this herpetofauna originating in the area of the current Orinoco basin. With advances in genetics, research in molecular systematics makes species of this micro herpetofauna in particular indicators allowing revelation and precision of dating of geological phenomena and events. A phylogeny of a herpetofauna which is therefore revealing of a geological history of the Eastern Caribbean where Martinique stands out as exceptional and central on the scale of this Grenadian basin where it is located; playing a main role as a refuge island from -45 My, then a center of dispersion towards -12 My for an endemic herpetofauna older than previously thought. A proto-Martinique island refuge in this almost unique Grenadian basin (with proto-Barbados) maintained not submerged from -25 My, which was quite naturally center of dispersion, first towards a proto-Saint Lucia in geological phase re-emergence around -12 My.

Research in herpetology and in particular the first results of work in molecular systematics [START_REF] Hedges | Molecular phylogeny and biogeography of West Indian frogs of the genus Leptodactylus (Anura, Leptodactylidae)[END_REF][START_REF] Surget-Groba | Journal of Biogéography n°40 -School of Biological sciences , College of Natural Sciences, Bangor University -UK -Supporting information A likelihood framework analysis of an island radiation: comparative phylogeography of the Lesser Antillean gecko Sphaerodactylus vincenti[END_REF][START_REF] Hedges | Tree of Life Reveals Clock-Like Speciation and Diversification[END_REF] relating to this Martinican herpetofauna are shaking up the geosciences and obliging a revision of the historical scenario of the basically tectonic geological formation and evolution of the Eastern Caribbean with emerged land (forming an isthmus since the Eocene, well over -45 My), then almost all submerged around -25 My. This isthmus at the top of '' a cordillera (volcanic arc), located on the southeastern edge of the so-called Caribbean plate, at the junction of the three lithospheric films (or plates), Caribbean, North America and South America, formed one of the two bridges, with that of the Aves ride, between the Orinoco basin and the Greater Antilles. (See Fig. 3)

The history of biological endemism on the surface of emerged land therefore does not begin there only in the Miocene at "-24 My" as it is generally accepted. It is in such a context that the Spherodactyls of Martinique (Ti-Mabouya, in Martinique language), are emblematic of this evolution of biodiversity in a logic of cause and effect with the geological history of this Martinique, obviously intrinsic to that of the eastern Caribbean.

The base of Martinique consists of a fragment of the Aves wrinkle located east of the Caribbean tectonic film which evolved 120 My ago in the Pacific Ocean. This Aves wrinkle, with its two consecutive folds of the west to east displacement of the Caribbean lithospheric film encrusted between the two most important American lithospheric films, is part of a set called the Caribbean Mesozoic Arc . [START_REF] Stephan | Paleogeodynamic maps of the Caribbean[END_REF].

Martinique is the only island in the eastern Caribbean to present three phases of volcanic activity juxtaposed with a regular displacement of these activities from east to west since 25 My (Fig. 4 and 5, Germa, Aurélie 2008) . Located remarkably on a position in the far east of this Caribbean arc, "where a North-South tangent passes", precisely at the latitude of the splitting of American lithospheric films. (Fig. 1: Bernard Mercier de Lépinay 2010) 1) In Martinique, the products of the ancient arc (gray) outcrop on the peninsulas of La Caravelle (east) and Sainte Anne (south-east). 2 -3) Then the volcanic front moved slightly towards the west. An underwater chain was built between 17 and 8 My (Vauclin-Pitault chain, purple) then emerged to build the peninsula of Trois Ilet until 7 My (Volcanism of the Southwest, green). 4) The monogenic volcanoes of Trois-Ilets were active from 2.4 to 0.34 My (orange). 5) Le Morne Jacob (dark blue, see below) is a shield volcano which emerged around 5.5 My and was active until 1.5 My. 6) The Carbet complex began to build on the western flank of Le Morne Jacob about 1 My ago. A flank collapse affected its entire western flank and was immediately followed by the eruption of the Pitons du Carbet in the depression. 7 -8) Mont Conil (light blue) emerged about 550 ka ago. Since then, volcanic activity has never really stopped in the northern block of the island. Around 126 ka, a flank collapse destroyed about 25 km3 of the building, then the volcanic activity of Mount Pelée (pink) proper became centralized in the depression formed. Germa Aurélie (2008).

II.

MATÉRIEL ET MÉTHODE

Cette contribution pour une révision de la liste taxinomique des Sphérodactyles de la Martinique, plaidant pour élever au rang d'espèce le Ti Mabouya de la péninsule de Sainte-Anne du sud-est de la Martinique, le Sphaérodactylus psammius, se fonde sur des publications datant du début du XXème siècle, mais aussi de publications récentes traitant de biogéographie, de géologie et de systématique : [START_REF] Barbour | A contribution to the zoogeography of the West Indies, with especial reference to Amphibians and Reptiles[END_REF][START_REF] Barbour | Memoirs of the Museum of Comparative Zoology at Harvard College[END_REF][START_REF] Barbour | A list of Antillean reptiles and amphibians[END_REF][START_REF] Barbour | A second list of Antillean reptiles and amphibians[END_REF][START_REF] Kl | Third list of Antillean reptiles and amphibians[END_REF][START_REF] King | Systematics of Lesser Antillean lizards of the genus Sphaerodactylus[END_REF][START_REF] Schwartz | A review of Sphaerodactylus vincenti in the southern Windward Islands[END_REF][START_REF] Schwartz | A new subspiecies of the gecko Sphaérodactylus microlepis[END_REF], Westerkamp 1975[START_REF] Schwartz | A check-list of West Indian Amphibians and Reptiles[END_REF][START_REF] Andreieff | Histoire géologique de la Martinique[END_REF][START_REF] Stephan | Paleogeodynamic maps of the Caribbean[END_REF][START_REF] Hedges | The origin of West Indian amphibians and reptiles. Contributions to West Indian herpetology: a tribute to Albert Schwartz[END_REF][START_REF] Malhotra | Reptiles & Amphibians of the Eastern Caribbean[END_REF][START_REF] Breuil | Atlas des amphibiens et reptiles de la Martinique : illustration et cartographie provisoires (IVème partie) -Histoire géologique et différenciation géologique des populations de sphérodactyles de la Martinique (Sphaerodactylus vincenti) -Laboratoire des Reptiles et Amphibiens Muséum national d'Histoire naturelle -France -Thorpe[END_REF], 2009, 2011, 2013[START_REF] Steinberg | Catalogue of American Amphibians and Reptiles. Sphaerodactylus vincenti[END_REF][START_REF] Germa | Evolution volcano-tectonique de l'île de la Martinique (arc insulaire des Petites Antilles): nouvelles contraintes géochronologiques et géomorphologiques[END_REF] The taxonomic list of Spherodactyls of Martinique here established includes the three species identified in historical times; it being agreed today by IUCN (Anonymous 2012) that this limit of these so-called historical times must not be earlier than the year 1500 of the Judeo-Christian calendar. The three taxa declined in the list of Spherodactyls proposed are endemic to Martinique; two of them S. vincenti and S. festus should be considered as subendemic to Martinique [START_REF] Lescure | Biog2ographie et insularité : l'exemple des Petites Antilles[END_REF]; since they are taxa common to the islands located south of Martinique, Saint-Vincent, Béquia (Grenadines), Saint Lucia, for S. vincenti, and north of Martinique, Dominica for S. festus.

III.

TAXINOMIC IV.

COMMENTS

The historical scenario of evolution of the genus Sphaérodactylus in Martinique presents similarities to that declined for the genus Anolis of Martinique. Moreover this historical scenario would be globally valid with chronological differences, for several genera of the micro herpetofauna the most formerly dispersed from the basin of the Orinoco towards the geographical position of the current banks of islands of the eastern Caribbean from the Eocene, for more than 40 My.

The common ancestors of the Ti Mabouya, dwarf geckos, Squamates of the genus Sphaerodactylus from Martinique, were able to disperse almost terrestrially, over the Caribbean isthmus (see Fig. 2 and Fig. 3) formed more than -40 My [START_REF] Stephan | Paleogeodynamic maps of the Caribbean[END_REF] 2013 ). This global scenario is presumed to become more precise as the work in molecular biology proceeds, and to be the same for Anurans, Amphibians of the genus Eleutherodactylus, as well as concerning the snake reptiles of the genus Tetracheilostoma.

Sphaérodactylus vincenti

Photo by A. Fong

The species was first observed and described on the island V.

CONCLUSION

In the current state of knowledge, and since 2013, only two taxa are recognized for the Martinique Spherodactyl species present in the eastern Caribbean (Dominica, Saint Lucia, Saint Vincent, Béquia): Sphaérodactylus festus and Sphaérodactylus vincenti.

The geographical position of the current Sainte-Anne peninsula in the south-east of Martinique corresponding with certainty to the zone of refuge then of center of dispersion of the Sphaérodactylus which dispersed in the Eastern Caribbean from Dominica to the Grenadines; it is clearly relevant to raise to the rank of new species S. psammius, the recognized line of Sphaérodactylus from south-eastern Martinique, of the Sainte-Anne peninsula, under the current taxon of subspecies S. v. psammius.

Fig. 1 :

 1 Fig. 1 :

Fig. 2 :

 2 Fig. 2 : Extracts from cartographic works by Stephan, J.F. et al. (1990).

  , Hedges & Kumar 2009, Daltry 2009, Henderson & Powel 2009, Powell & Henderson 2012, Pichot 2012, Barat 2013, Surget-Groba & Thorpe 2013, Sepulchre et al. 2014, Sacco 2015, Allard Saint-Albin 2015, Hedges, Marin, Suleski, Paymer & Kumar 2015, Mauvois 2015, Saffache et al. 2018, Bourgade & Fong 2018, Dewynter 2018, Bourgade 2019, Dewynter, De Massary, Bochaton, Bour, Ineich, Vidal, & Lescure 2019.

  It is the particular configuration of the proto Martinique (South-East) geologically, that is to say the existence in the basement, of large limestone formations out of water, compared to that of the other islands of the Eastern Caribbean (including Barbados), to have found itself as the only possible refuge for fauna and flora, with banks of islands tectonically remained emerged between -24 My and -12 My. It is this situation which explains this scenario of which we have scientific certainty since 2008, of a Martinique center of dispersion on the scale of the Eastern Caribbean, for the Anolis(Johansson, Surget-Groba & Thorpe 2008), for the Sphaérodactylus (Surget-Groba & Thorpe

  Photo by M. BourgadeThe species was first observed and described in Martinique byBarbour in 1915. It has been reduced to S. v. festus from the 1960s to the rank of subspecies subservient to S. vincenti, since it was mistakenly believed these Spherodactyla present in Martinique, originating in Saint-Vincent.Many alleged subspecies have been described in Martinique particularly, in Dominica, in Saint Lucia, Béquia, by[START_REF] Schwartz | A review of Sphaerodactylus vincenti in the southern Windward Islands[END_REF][START_REF] Lazell | A new Sphaerodactylus (Sauria, Gekkonidae) from Bequia, Grenada Bank, Lesser Antilles[END_REF] and[START_REF] Breuil | Composition et endémisme de l'herpétofaune martiniquaise : histoire géologique et différenciation intra-insulaire du sphérodactyle[END_REF], under the taxon of S. vincenti ssp. ; In reality molecular biology tells us that the lineages of Spherodactyls present in the north of Martinique (and according to our own observations in the north of the peninsula of Southeast Martinique, north of the Montagne du Vauclin), are precisely originating from this southeast of Martinique, from where they dispersed from 12 My ago (Surget-Groba & Thorpe 2013). The Sphaérodactylus present in Dominica are descendants of a line of Spherodactyls whose specimens and / or propagules have undoubtedly drifted on rafts of plants towards Dominica (thanks to climatic events and phenomena or geological upheavals), from from the North Caribbean coast of Martinique, at the latitude of the Pitons du Carbet. « S. v. psammius » => Sphaérodactylus psammius

LIST OF MARTINIQUE SPHERODACTYLS ABORIGINAL SPECIES SAUROPSIDA Huxley

	SQUAMATA Oppel, 1811 …………………………………………………………………… SQUAMATES
	«SAURIA» Brongniart, 1800 ……………………………………………………………………. SAURIENS
	SPHAERODACTYLIDAE Underwood, 1954 …………………………………… SPHERODACTYLS
	Sphaerodactylus Wagler, 1830 …………………………………………………………. Spherodactyl
	• Sphaerodactylus vincenti Boulenger, 1891 ………………………… The Spherodactyl of the péninsula
	…………….. from the southwest of Martinique scattered towards Saint-Vincent, Béquia, St-Lucia
	• Sphaerodactylus festus Barbour, 1915 ……………………. The Spherodactyl with ocelli (eyespots)
	……………………………….. dispersed towards the North of Martinique, then towards the Dominica

, 1864 …………………………………………………. SAUROPSIDES • Sphaerodactylus psammius Schwartz, 1964 …………………………… The Spherodactyl of South East ………………………………………………………………………. of Martinique ; of the peninsula of Sainte-Anne • S. p. adamas Schwartz, 1964 ……… The Spheroctactyl present on the Diamond Rock

  ; surviving on the island benches of the Martinique proto of the extinction linked to the submersion of this Caribbean isthmus which took place from -40 My to -25 My, isolated, they have therefore genetically diverged since that time from their close relatives both from the South, ancestors of Sphaerodactylus moléi Boettger 1894 (Trinidad and South America), as well as those from the North, singularly ancestors of Sphaerodactylus nicholsi Hass 1991 (Puerto Rico) and generally of the Greater Antilles. The first speciations of this Ti Mabouya on Martinique proto began -30 My ago[START_REF] Hedges | Tree of Life Reveals Clock-Like Speciation and Diversification[END_REF]; then other speciation processes and therefore divergences started from -12 My (Surget-Groba & Thorpe 2013); with dispersion towards new land, proto Saint Lucia, on Gros Ilet, and Maria Islands,

generating S. microlepis Reinhardt & Lütken, 1862, and S. festus towards Caravelle (North Martinique).

Acknowledgments -I warmly thank Ansèl Fong, and David Maceira, biologists from Santiago de Cuba, who since 2013 have encouraged me to persevere in my approach to herpetological and geological research. My thanks also go to Roger Mauvois for this basic information on geology, as well as to Garcin Malsa for his readings and reviews of my publications and manuscripts.


ON THE EXPECTED RESULT FOR THE SECOND MOMENT OF TWISTED L-FUNCTIONS

The fourth moment of Dirichlet L-functions is M 4 (q) := χ∈X * (q)

L(χ, 1/2) 4 where X * (q) stands for the set of primitive Dirichlet characters of modulus q whereas the second moment of twisted L-functions is

M 2, f (q) := χ∈X * (q) L( f × χ, 1/2) 2
where f is a fixed holomorphic primitive cusp form of level D f 1 coprime with q, nebentypus ψ f of modulus D f satisfying ψ f (-1) = (-1) k f and integer weight k f 1 (see appendix [RiRo] for the automorphic background). Note that card X * (q) := ϕ * (q) = q p||q 1 -2

p p 2 |q 1 - 1 p 2
according to [START_REF] Iwaniec | Analytic number theory[END_REF]Equation (3.7) Page 46]. Finding an asymptotic formula for M 4 (q) as q → +∞ with a power saving in the error term should be philosophically as difficult as the corresponding question for M 2, f (q) since

M 4 (q) = M 2, ∂E (1/2,.) ∂s (1/2) (q)
where E (z, s) is the real-analytic Eisenstein series on the modular curve X 0 (1) associated to the cusp ∞. Note that M 4 (q) is itself the q-analog of the fourth moment of the Riemann zeta function in the t -aspect (see [In] and [START_REF] Heath-Brown | The fourth power moment of the Riemann zeta function[END_REF]) but this analogy will not be developed here. D. Heath-Brown ( [HB]) for M 4 (q) and T. Stefanicki ( [St]) for M 2, f (q) proved the following analogous results.

Theorem. If q goes to infinity then

M 4 (q) = 1 2π 2 p|q (1 -p -1 ) 3 (1 + p -1 ) log 4 (q)ϕ * (q) + O 2 ω(q) q log 3 q
where ω(q) stands for the number of prime divisors of q and

M 2, f (q) = 6(4π) k f πΓ(k) || f || 2 p|q (1 -p -1 ) 2 (1 -p -2 ) 1 - λ f (p 2 ) -2 p + 1 p 2 log (q)ϕ * (q)
+ O 2 ω(q) q log δ q where || f || stands for Petersson's norm of f and 0 < δ < 1 is any real number satisfying

∀ε > 0, p x λ f (p) (δ + ε) x log x .
Remark 0.1. It is possible to choose δ = ( 2 + 3 3)(5 2) -1 according to [Ra].

Date: Version of July 21, 2009.

1 Remark 0.2. The formula for M 4 (q) is an asymptotic formula for almost all integers but it turns out that if ω(q) > log -1 2 log 2 q then the error term is not smaller than the main term. Similarly, if ω(q) > (1-δ) log -1 q log 2 q then the error term in M 2, f (q) is not smaller than the main term such that the formula for M 2, f (q) is an asymptotic formula only for a set of integers of zero density according to [HR].

These results have been improved by K. Soundararajan ( [So]) for M 4 (q) and in [GKR] for M 2, f (q) as follows.

Theorem. If q goes to infinity then

M 4 (q) = 1 2π 2 p|q (1 -p -1 ) 3 (1 + p -1 ) log 4 (q)ϕ * (q) 1 + O log -1/2 2 (q)
and if ω(q) exp 300 -1 log 2 q log -1 Young ( [Yo]) got an asymptotic formula with a power saving in the error term for M 4 (q) in the prime modulus case.

3 q then M 2, f (q) = 6(4π) k f πΓ(k) || f || 2 p|q (1 -p -1 ) 2 (1 -p -2 ) 1 - λ f (p 2 ) -2 p + 1 p 2 log (q)ϕ * (q) × 1 + O log -1 2 (q) . M.
Theorem. If q goes to infinity among the prime numbers then M 4 (q) = P (log q)ϕ * (q) + O ϕ * (q)q -1 80 + θ 40 where P is an explicit polynomial of degree 4 and θ = 7/64. Remark 0.3. Note that -1/80 + θ/40 < 0. The same result should hold without the restriction q prime but it remains an open question so far. It should be very difficult to extend this result to any integer following the proof given in [Yo] since many technical difficulties would occur if q is composite.

Remark 0.4. This particular value of θ is the best currently known approximation towards Ramanujan-Petersson-Selberg's conjecture in GL(2) according to [Ki] and [KH].

It turns out that the analogous result for M 2, f (q) is still an open question and the purpose of this note is to identify the underlying analytic issue, which occurs in the second question. Let us state the expected result.

Expected Result. There exists some absolute constant α > 0 such that if q goes to infinity then M 2, f (q) = P (log q)ϕ * (q) + O ϕ * (q)q -α where P is an explicit polynomial of degree 1.

Proving this requires some new input in order to solve unbalanced shifted convolution problems. Such new input would have many other interesting applications.

Acknowledgements. The author would like to thank Karim Belabas, Hendrik W. Lenstra and Don B. Zagier for inviting him in Mathematisches Forschungsinstitut Oberwolfach on the occasion of the workshop "Explicit Methods in Number

Theory".