# Selecting the Best Forecasting-Implied Volatility Model Using Genetic Programming

Authors: Wafa Abdelmalek; Sana Ben Hamida; Fathi Abid
Journal: Journal of Applied Mathematics and Decision Sciences
 (2009)

---

## Abstract

The volatility is a crucial variable in option pricing and hedging strategies. The aim of this paper is to provide some initial evidence of the empirical relevance of genetic programming to volatility's forecasting. By using real data from S&P500 index options, the genetic programming's ability to forecast Black and Scholes-implied volatility is compared between time series samples and moneyness-time to maturity classes. Total and out-of-sample mean squared errors are used as forecasting's performance measures. Comparisons reveal that the time series model seems to be more accurate in forecasting-implied volatility than moneyness time to maturity models. Overall, results are strongly encouraging and suggest that the genetic programming approach works well in solving financial problems.


---

## Body

## 1. Introduction

Although the Black-Scholes (BS) formula [1] is commonly used to price derivative securities, it has some well-known deficiencies. When the BS formula is inverted to imply volatilities from reported option prices, the volatility estimates exhibit a dependency on both strike price and maturity, giving rise to what options professionals call “volatility smile” and “volatility term structure.” In many markets, prior to the October 1987 stock market crash, there appeared to be symmetry around the zero moneyness, where in-the-money (ITM) or out-of-the-money (OTM) options have higher implied volatilities than at-the-money (ATM) options. This dependency of implied volatility on the strike, for a given maturity, became known as the smile effect, although the exact structure of volatility varied across markets and even within a particular market from day to day. However, after the crash, the smile had a changed shape in many markets, particularly for traded options on stock market indexes (such as S&P500 index options), where the symmetric smile pattern (U-shape) has changed to more of a sneer. This is often referred to in the markets as the “volatility smirk,” where call option-implied volatilities are observed to decrease monotonically with strike price. This can arise when the market places a relatively greater probability on a downward price movement than an upward movement, resulting in a negatively skewed implied terminal asset distribution.The idea of the volatility smile had its genesis in the early papers documenting the systematic pricing biases of the BS option pricing model. Early tests by Black [2] found that the BS model underprices OTM stock options and overprices ITM stock options. In contrast to Black [2], Macbeth and Merville [3], studying call options, listed on the Chicago Board Options Exchange (CBOE) from December 1975 to December 1976, found evidence that the BS model underprices ITM options and overprices OTM options. Later, Rubinstein [4], studying options price data for the thirty most actively traded option classes on the CBOE between August 1976 and October 1978, found some confusing patterns. Rubinstein [4] reported a systematic mispricing pattern similar to that reported by Macbeth and Merville [3], where the BS model overprices OTM options and underprices ITM options, for a time period between August 1976 and October 1977. However, for a time period between October 1977 and October 1978, Rubinstein [4] reported a systematic mispricing pattern similar to that reported by Black [2], where the BS model underprices OTM options and overprices ITM options. Clewlow and Xu [5], studying options on stock index futures listed on the Chicago Mercantile Exchange (CME), found asymmetric patterns for the smile. Xu and Taylor [6], studying currency options traded in Philadelphia between 1984 and 1992, found empirically that the option bias was twice the size they were expecting, increasing its magnitude as maturity approaches. Duque and Paxson [7] also found the smile effect for options traded on the London International Financial Futures and Options Exchange (LIFFE) during March 1991, and conjectured that there is a possible empirical relation between time to maturity and the U-shape format of the smile. Heynen [8] also found empirical evidence for the exercise price bias when observing stock index options during 1989 on the European Option Exchange (EOE). Gemmill [9] found the same bias for options on the FTSE100 during 5-year period, although the smile showed different patterns for different days extracted from the sample. Dumas et al. [10] also found empirical smile patterns for options on the S&P500 stock index, but its shape seemed to be asymmetric and changing along time to maturity.Volatility smiles are generally thought to result from the parsimonious assumptions used to derive the BS model. In particular, the BS model assumes that security log prices follow a constant variance diffusion process. The existence of systematic bias in the model reflects a departure from the BS assumptions, with the result that the probability distribution is not lognormal. Black [2] suggested that the nonstationary behavior of volatility would lead the BS model to overprice or underprice options. Since then, a considerable body of research has attempted to explain this bias. Some authors attributed the smile to transaction costs [11]. Some other authors imputed the smile to the behavior of traders [12], others tried to explain the smile by the introduction of conditional heteroskedasticity [13], stochastic interest rate [14], jump processes [15], and stochastic volatility [16]. Literature is not unanimous in finding causes for the smile effect and the models developed in order to cover this bias have only partially solved the problem. Unfortunately, these increasingly complex models have some caveats in common: they require strong assumptions, their parameters are difficult to estimate, and they are particularly prone to in-sample overfitting. Furthermore, their stochastic volatility process itself is unobservable and must be proxied for.In contrast to parametric forecasting models, nonparametric approach has the distinct advantage of not relying on specific assumptions about the underlying asset price dynamics and is therefore robust to specification errors that might affect adversely parametric models. As mentioned in Ma et al. [17], traditional financial engineering methods based on parametric models such as the GARCH model family seem to have difficulty to improve the accuracy in volatility forecasting due to their rigid as well as linear structure. The complex and nonlinear nature of processes means that more advanced techniques are often needed. Genetic programming (GP) could effectively deal with nonlinearity and flexibility, which opens up an alternative path besides other data-based approaches.This paper intends to establish a systematic approach and a GP software tool for analysts to improve the accuracy of the forecast. This volatility's forecasting approach should be free of strong assumptions and more flexible than parametric method. Reliable volatility forecasts can greatly benefit professional option traders, market makers who need to price derivatives, and all investors with risk management concerns. Implied volatilities, generated from option markets, can be particularly useful in such contents as they are forward-looking measures of the market's expected volatility during the remaining life of an option. We attempt then to forecast the implied volatility of BS by formulating a nonlinear and nonparametric approach based on GP algorithm.It is important to notice here that we do not claim to have found the perfect implied volatility's forecasting model. Ultimately, all implied volatility's forecasting models, no matter how complex, are misspecified and the search over the universe of volatility forecasting should concentrate on finding the least misspecified model.The remainder of the paper is divided into three sections. Section2 describes the design of the GP and the methodology applied in the current research. A brief introduction to the GP is given in Section 2.1, while the data to conduct the analysis and the GP implementation are described in Section 2.2. Section 3 reports the experimental results of the GP search for volatility forecast. Selection of the best generated GP-implied volatility forecasting models is presented in Section 3.1. Implied volatility patterns are described in Section 3.2. Finally, conclusion and future directions of the research are provided in Section 4.


## 2. Research Design and Methodology

### 2.1. Overview of Genetic Programming

GP [18] is an evolutionary algorithm which extends the basic genetic algorithms [19] to process nonlinear problem structure. This optimization technique is based on the principles of natural evolution. The algorithm attempts to evolve solutions by using the principle of survival and reproduction of the fittest and genetic operators analogous to those occurring in biological species. In GP, solutions are represented as tree structures that can vary in size and shape, rather than fixed length character strings as in genetic algorithms. This means that GP can be used to perform optimization at a structural level. The advantage of the GP approach is that it allows one to be agnostic about the general form of optimal forecasting model, and to search efficiently in a nondifferentiable space of models. It also has the attractive feature that one can build into the search procedure the relevant performance criterion directly in the form of the measure of fitness. It differs from most other numerical search techniques in that it simultaneously involves a parallel search involving many points in the search space.Given the advantages of GP, it has been widely used to solve complex high-dimensional and nonconvex optimization problems (Kallel et al. [20]). In fact, the GP has been successfully applied to forecasting financial time series, including stock prices (Tsang et al. [21]), stock returns (Kaboudan [22]), and international currency exchange rates (Kaboudan [23]). In particular, the GP has proved successful at forecasting time series volatility in different markets. Zumbach et al. [24] and Neely and Weller [25] have applied GP to forecast foreign exchange rate volatility. Neely and Weller [25] have tested the forecasting performance of GP for USD-DEM and USD-YEN daily exchange rates against that of GARCH (1,1) model and a related RiskMetrics volatility forecast over different time horizons, using various accuracy criteria. While the GP rules didn't usually match the GARCH (1,1) or RiskMetrics models' MSE or R2, its performance on those measures was generally close; but the GP did consistently outperform the GARCH model on mean absolute error (MAE) and model error bias at all horizons. Overall, on some dimensions, the GP has produced significantly superior results. Using high-frequency foreign exchange USD-CHF and USD-JPY time series, Zumbach et al. [24] have compared the GP forecasting accuracy to that of historical volatilities, the GARCH (1,1), FIGARCH, and HARCH models. According to the root-mean-squared errors, the generated GP volatility models did consistently outperform the benchmarks.The GPs forecasting accuracy has been shown also in index market. Using historical returns of Nikkei 225 and S&P500 indices, Chen and Yeh [26] have applied a recursive genetic programming (RGP) approach to estimate volatility by simultaneously detecting and adapting to structural changes. Results have shown that the RGP is a promising tool for the study of structural changes. When RGP discovers structural changes, it will quickly suggest a new class of models so that overestimation of volatility due to ignorance of structural changes can be avoided. Applying a combination of theory and techniques such as wavelet transform, time series data mining, Markov chain-based discrete stochastic optimization, and evolutionary algorithms GA and GP, Ma et al. [27, 28] have proposed a systematic approach to address specifically nonlinearity problems in the forecast of financial indices using intraday data of S&P100 and S&P500 indices. As a result, accuracy of forecasting had reached an average of over 75% surpassing other publicly available results on the forecast of any financial index.These encouraging findings suggest that a GP may be a powerful tool for forecasting historical volatility based on previous returns. The present paper extends the studies mentioned earlier by forecasting the implied volatility instead of historical volatility using GP. In fact, implied volatility has been found to be a more accurate and efficient forecast of future volatility than historical volatility [29–32]. Furthermore, the implied volatility is estimated from prices of index options in a wide range of strike prices, not just ATM strikes.Our GP software is referred to as symbolic regression written in C++ language. Symbolic regression was one of the earliest applications of GP (Koza [18]), and continues to be widely studied (Cai et al. [33]; Gustafson et al. [34]). It is designed to perform genetic regression which is a function that relates a set of inputs to an output. The set data on which the program operates to determine the relationship between input parameters and output is called training set. The set of data on which the resulting formula is tested is called test set. The GP's algorithm includes the following steps (Algorithm 1).Algorithm 1:A flow chart summarizing the GP algorithm structure.
Initialize populationwhile  (termination condition not satisfied)  dobeginEvaluate the performance of each function according to the fitness criterionSelect the best functions in the population using the selection algorithmGenerate new individuals by crossover and mutationendwhileReport the best solution foundend


#### 2.1.1. Terminal and Function Sets

Two of the most important aspects of the GP algorithm are the function and terminal sets that building blocks use to construct expression trees which are solutions to the problem. Together, they define the ingredients that are available to GP to create computer programs.The terminal set corresponds to the inputs of the program. It is determined according to the domain of problems and the inputs can be constants or variables. The function set may contain standard arithmetic operations, standard programming operations, standard mathematical functions, logical functions, or domain-specific functions. The function set combined with the terminal set enables the algorithm to generate a tree-structured solution to the problem, where nodes define function set and where leaf nodes define terminal set.


#### 2.1.2. Initialization

GP starts by randomly creating an initial population of trees, which are generated by randomly picking nodes from a given terminal set and function set. A restriction on the maximum allowed depth is imposed to avoid that the generated trees to be too complex.


#### 2.1.3. Fitness Evaluation

The evolutionary process is driven by a fitness function that evaluates the performance of each individual (tree) in the population. This fitness function assigns numerical values to each individual in the population. These values are then used to select individuals for reproduction.


#### 2.1.4. Selection

After the fitness of each individual in the population is measured, GP then 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 on 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 offsprings with larger chance.


#### 2.1.5. Genetic Operators

Crossover and mutation are the two basic operators which are applied to the selected individuals in order to generate the next population. The mutation operator affects random changes in a tree by randomly altering certain nodes or subtrees, whereas the crossover operator performs swaps of randomly subtrees between two selected parents. The parameter choices for crossover and mutation are clearly critical in ensuring a successful GP application. They impact on populational diversity and the ability of GP to escape from local optima [35].After applying the genetic operators, the new individuals are evaluated through the fitness function in order to replace those of the current population.


#### 2.1.6. Termination Criterion

Once the new population has been created, the current population is replaced by the new one, and the process of fitness evaluation, selection, and the application of genetic operators are repeated until some termination criterion is satisfied. This may be the maximum number of generations or convergence level of evolution generations.


### 2.2. Data and Implementation

Data used to perform our empirical analysis are daily prices of European S&P500 index calls 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. The minimum tick for series trading below 3 is 1/16 and for all other series 1/8. 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). Option prices are set to be the average of the bid and ask prices. The risk-free interest rate is approximated by using 3-month US Treasury bill rates.The beginning sample contains 42504 daily observations of call option prices and their determinants. Four exclusion filters are applied, which yields a final sample of 6670 daily observations and this is the starting point for our empirical analysis. First, as call options with time to maturity less than 10 days may induce liquidity biases, they are excluded from the sample. 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 due to the lack of trading volume. Finally, option prices not satisfying the arbitrage restriction [14], C≥S-Ke-rτ, are not included.The full sample is sorted by time series and moneyness-time to maturity. For time series, data are divided into 10 successive samples each containing 667 daily observations. For moneyness-time to maturity, data are divided into 9 classes, classes are, respectively, OTMST (C1), OTMMT (C2), OTMLT (C3), ATMST (C4), ATMMT (C5), ATMLT (C6), ITMST (C7), ITMMT (C8), and ITMLT (C9). According to moneyness criterion: a call option is said OTM ifS/K<0.98; ATM if S/K∈[0.98,1.03[; and ITM if S/K≥1.03. According to time to maturity criterion: A call option is short-term (ST) if τ<60 days; medium term (MT) if τ∈[60,180] days; long-term (LT) if τ>180 days.The first preparatory step to run a GP is to specify the ingredients the evolutionary process can use to construct potential solutions. The terminal and function sets are described in Table1.Table 1
Terminal set and function set.
ExpressionDefinitionC/KCall price/Strike priceS/KIndex price/Strike priceτTime to maturity+Addition-Subtraction*Multiplication%Protected division:x%y=1 if y=0;   x%y=x%y otherwiselnProtected natural log:ln(x)=ln(|x|)ExpExponential function:exp(x)=exSqrtProtected square root:x=|x|NcdfNormal cumulative distribution functionΦThe search space for GP is the space of all possible parse trees that can be recursively created from the sets of terminal and function. The terminal set includes the inputs variables, notably, the call option price divided by strike priceC/K, the index price divided by strike price S/K and time to maturity τ. The function set includes unary and binary nodes. Unary nodes consist of mathematical functions, notably, cosinus function (cos), sinus function (sin), log function (ln), exponential function (exp), square root function (), and the normal cumulative distribution function (Φ). Binary nodes consist of the four basic mathematical operators, notably, addition (+), subtraction (-), multiplication (×), and division (%). The basic division operation is protected against division by zero and the log and square root functions are protected against negative arguments. The basic measure of model accuracy used is the mean squared error (MSE) between the target (yt) and forecasted (ŷt) output volatility, computed as follows:
(2.1)MSE=1N∑t=1N(yt-ŷt)2,
where yt is the implied volatility computed using the BS formula (C=SN(d1)-Ke-rτN(d2), d1=(ln(S/K)+(r+0.5σ2)τ)/στ, d2=d1-στ.),ŷt is the generated GP volatility, and N is the number of data sample.The implementation of genetic programming involves different parameters that determine its efficiency. A common approach in tuning a GP is to undertake a series of trial and error experiments before making parameter choices for the final GP runs. A two-step approach is used. First, the optimal set of genetic parameters is determined based on time series data. By varying parameters for the genetic programs when training them on the first sample data, their performance is tested on a test data set from a later date. Various genetic programs are tested. Each program is run ten times with ten different random seeds. The choice of the best genetic program is made according to the mean and median of mean squared errors (MSEs) for training and testing sets. Second, the optimal set of parameters, relative to the best genetic program selected in the first step, are used for training the genetic program separately on each of the first 9 samples using ten different seeds and testing the program's performance only on the test data from the immediately following date. This limits the effects of nonstationarities and avoid data snooping. The same genetic parameters and random seeds used for time series data are applied to moneyness-time to maturity data in an attempt to analyze the effect of data' group on the GP performance; each moneyness-time to maturity class is subdivided into training set and test set. Parameters used in this paper are listed in Table2.Table 2
Summary of GP settings and parameters.
Method of generation:Ramped half and halfPopulation size:100Offspring size:200Generations number:400Maximum depth of new individual:6Maximum depth of the tree:17Selection algorithm:Tournament selection (4)Crossover probability:60%Mutation probability:40%Branch mutation:20%Point Mutation:10%Expansion mutation:10%Replacement strategy:Comma replacementStopping criterion:Maximum number of generationsAn initial population of 100 individuals is randomly generated using random number generator algorithm with random seeds. The generative method used for initialization is the ramped half and half as detailed in Koza [18]. This method involves generating an equal number of trees using a maximum initial depth that ranges from 2 to 6. For each level of depth, 50% of the initial trees are generated via the full method and the other 50% are generated via the grow method. A maximum size of tree measured by depth is 17. This is a popular number used to limit the size of tree (Koza [18]). It is large enough to accommodate complicated formulas and works in practice. Based on the fitness criterion, the selection of the individuals for reproduction is done with the tournament selection algorithm. A group of individuals is selected from the population with a uniform random probability distribution. The fitness values of each member of this group are compared and the actual best is selected. The size of the group is given by the tournament size which is equal here to 4. The crossover operator is used to generate about 60% of the individuals in the population, while the mutation operator is used to generate about 40% of the population. 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 subtree rooted at that node with a new randomly generated subtree. Expansion mutation operator randomly selects a terminal node in the tree, and then replaces it with a new randomly-generated subtree. Branch mutation is applied with a rate of 20%; point and expansion mutations are applied with a rate of 10% each. The method of replacing parents for the next generation is comma replacement strategy, which selects the best offspring to replace the parents. Of course, it assumes that offspring size is higher than parents size. The stopping criterion is the maximum number of generations 400.


## 2.1. Overview of Genetic Programming

GP [18] is an evolutionary algorithm which extends the basic genetic algorithms [19] to process nonlinear problem structure. This optimization technique is based on the principles of natural evolution. The algorithm attempts to evolve solutions by using the principle of survival and reproduction of the fittest and genetic operators analogous to those occurring in biological species. In GP, solutions are represented as tree structures that can vary in size and shape, rather than fixed length character strings as in genetic algorithms. This means that GP can be used to perform optimization at a structural level. The advantage of the GP approach is that it allows one to be agnostic about the general form of optimal forecasting model, and to search efficiently in a nondifferentiable space of models. It also has the attractive feature that one can build into the search procedure the relevant performance criterion directly in the form of the measure of fitness. It differs from most other numerical search techniques in that it simultaneously involves a parallel search involving many points in the search space.Given the advantages of GP, it has been widely used to solve complex high-dimensional and nonconvex optimization problems (Kallel et al. [20]). In fact, the GP has been successfully applied to forecasting financial time series, including stock prices (Tsang et al. [21]), stock returns (Kaboudan [22]), and international currency exchange rates (Kaboudan [23]). In particular, the GP has proved successful at forecasting time series volatility in different markets. Zumbach et al. [24] and Neely and Weller [25] have applied GP to forecast foreign exchange rate volatility. Neely and Weller [25] have tested the forecasting performance of GP for USD-DEM and USD-YEN daily exchange rates against that of GARCH (1,1) model and a related RiskMetrics volatility forecast over different time horizons, using various accuracy criteria. While the GP rules didn't usually match the GARCH (1,1) or RiskMetrics models' MSE or R2, its performance on those measures was generally close; but the GP did consistently outperform the GARCH model on mean absolute error (MAE) and model error bias at all horizons. Overall, on some dimensions, the GP has produced significantly superior results. Using high-frequency foreign exchange USD-CHF and USD-JPY time series, Zumbach et al. [24] have compared the GP forecasting accuracy to that of historical volatilities, the GARCH (1,1), FIGARCH, and HARCH models. According to the root-mean-squared errors, the generated GP volatility models did consistently outperform the benchmarks.The GPs forecasting accuracy has been shown also in index market. Using historical returns of Nikkei 225 and S&P500 indices, Chen and Yeh [26] have applied a recursive genetic programming (RGP) approach to estimate volatility by simultaneously detecting and adapting to structural changes. Results have shown that the RGP is a promising tool for the study of structural changes. When RGP discovers structural changes, it will quickly suggest a new class of models so that overestimation of volatility due to ignorance of structural changes can be avoided. Applying a combination of theory and techniques such as wavelet transform, time series data mining, Markov chain-based discrete stochastic optimization, and evolutionary algorithms GA and GP, Ma et al. [27, 28] have proposed a systematic approach to address specifically nonlinearity problems in the forecast of financial indices using intraday data of S&P100 and S&P500 indices. As a result, accuracy of forecasting had reached an average of over 75% surpassing other publicly available results on the forecast of any financial index.These encouraging findings suggest that a GP may be a powerful tool for forecasting historical volatility based on previous returns. The present paper extends the studies mentioned earlier by forecasting the implied volatility instead of historical volatility using GP. In fact, implied volatility has been found to be a more accurate and efficient forecast of future volatility than historical volatility [29–32]. Furthermore, the implied volatility is estimated from prices of index options in a wide range of strike prices, not just ATM strikes.Our GP software is referred to as symbolic regression written in C++ language. Symbolic regression was one of the earliest applications of GP (Koza [18]), and continues to be widely studied (Cai et al. [33]; Gustafson et al. [34]). It is designed to perform genetic regression which is a function that relates a set of inputs to an output. The set data on which the program operates to determine the relationship between input parameters and output is called training set. The set of data on which the resulting formula is tested is called test set. The GP's algorithm includes the following steps (Algorithm 1).Algorithm 1:A flow chart summarizing the GP algorithm structure.
Initialize populationwhile  (termination condition not satisfied)  dobeginEvaluate the performance of each function according to the fitness criterionSelect the best functions in the population using the selection algorithmGenerate new individuals by crossover and mutationendwhileReport the best solution foundend


### 2.1.1. Terminal and Function Sets

Two of the most important aspects of the GP algorithm are the function and terminal sets that building blocks use to construct expression trees which are solutions to the problem. Together, they define the ingredients that are available to GP to create computer programs.The terminal set corresponds to the inputs of the program. It is determined according to the domain of problems and the inputs can be constants or variables. The function set may contain standard arithmetic operations, standard programming operations, standard mathematical functions, logical functions, or domain-specific functions. The function set combined with the terminal set enables the algorithm to generate a tree-structured solution to the problem, where nodes define function set and where leaf nodes define terminal set.


### 2.1.2. Initialization

GP starts by randomly creating an initial population of trees, which are generated by randomly picking nodes from a given terminal set and function set. A restriction on the maximum allowed depth is imposed to avoid that the generated trees to be too complex.


### 2.1.3. Fitness Evaluation

The evolutionary process is driven by a fitness function that evaluates the performance of each individual (tree) in the population. This fitness function assigns numerical values to each individual in the population. These values are then used to select individuals for reproduction.


### 2.1.4. Selection

After the fitness of each individual in the population is measured, GP then 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 on 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 offsprings with larger chance.


### 2.1.5. Genetic Operators

Crossover and mutation are the two basic operators which are applied to the selected individuals in order to generate the next population. The mutation operator affects random changes in a tree by randomly altering certain nodes or subtrees, whereas the crossover operator performs swaps of randomly subtrees between two selected parents. The parameter choices for crossover and mutation are clearly critical in ensuring a successful GP application. They impact on populational diversity and the ability of GP to escape from local optima [35].After applying the genetic operators, the new individuals are evaluated through the fitness function in order to replace those of the current population.


### 2.1.6. Termination Criterion

Once the new population has been created, the current population is replaced by the new one, and the process of fitness evaluation, selection, and the application of genetic operators are repeated until some termination criterion is satisfied. This may be the maximum number of generations or convergence level of evolution generations.


## 2.1.1. Terminal and Function Sets

Two of the most important aspects of the GP algorithm are the function and terminal sets that building blocks use to construct expression trees which are solutions to the problem. Together, they define the ingredients that are available to GP to create computer programs.The terminal set corresponds to the inputs of the program. It is determined according to the domain of problems and the inputs can be constants or variables. The function set may contain standard arithmetic operations, standard programming operations, standard mathematical functions, logical functions, or domain-specific functions. The function set combined with the terminal set enables the algorithm to generate a tree-structured solution to the problem, where nodes define function set and where leaf nodes define terminal set.


## 2.1.2. Initialization

GP starts by randomly creating an initial population of trees, which are generated by randomly picking nodes from a given terminal set and function set. A restriction on the maximum allowed depth is imposed to avoid that the generated trees to be too complex.


## 2.1.3. Fitness Evaluation

The evolutionary process is driven by a fitness function that evaluates the performance of each individual (tree) in the population. This fitness function assigns numerical values to each individual in the population. These values are then used to select individuals for reproduction.


## 2.1.4. Selection

After the fitness of each individual in the population is measured, GP then 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 on 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 offsprings with larger chance.


## 2.1.5. Genetic Operators

Crossover and mutation are the two basic operators which are applied to the selected individuals in order to generate the next population. The mutation operator affects random changes in a tree by randomly altering certain nodes or subtrees, whereas the crossover operator performs swaps of randomly subtrees between two selected parents. The parameter choices for crossover and mutation are clearly critical in ensuring a successful GP application. They impact on populational diversity and the ability of GP to escape from local optima [35].After applying the genetic operators, the new individuals are evaluated through the fitness function in order to replace those of the current population.


## 2.1.6. Termination Criterion

Once the new population has been created, the current population is replaced by the new one, and the process of fitness evaluation, selection, and the application of genetic operators are repeated until some termination criterion is satisfied. This may be the maximum number of generations or convergence level of evolution generations.


## 2.2. Data and Implementation

Data used to perform our empirical analysis are daily prices of European S&P500 index calls 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. The minimum tick for series trading below 3 is 1/16 and for all other series 1/8. 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). Option prices are set to be the average of the bid and ask prices. The risk-free interest rate is approximated by using 3-month US Treasury bill rates.The beginning sample contains 42504 daily observations of call option prices and their determinants. Four exclusion filters are applied, which yields a final sample of 6670 daily observations and this is the starting point for our empirical analysis. First, as call options with time to maturity less than 10 days may induce liquidity biases, they are excluded from the sample. 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 due to the lack of trading volume. Finally, option prices not satisfying the arbitrage restriction [14], C≥S-Ke-rτ, are not included.The full sample is sorted by time series and moneyness-time to maturity. For time series, data are divided into 10 successive samples each containing 667 daily observations. For moneyness-time to maturity, data are divided into 9 classes, classes are, respectively, OTMST (C1), OTMMT (C2), OTMLT (C3), ATMST (C4), ATMMT (C5), ATMLT (C6), ITMST (C7), ITMMT (C8), and ITMLT (C9). According to moneyness criterion: a call option is said OTM ifS/K<0.98; ATM if S/K∈[0.98,1.03[; and ITM if S/K≥1.03. According to time to maturity criterion: A call option is short-term (ST) if τ<60 days; medium term (MT) if τ∈[60,180] days; long-term (LT) if τ>180 days.The first preparatory step to run a GP is to specify the ingredients the evolutionary process can use to construct potential solutions. The terminal and function sets are described in Table1.Table 1
Terminal set and function set.
ExpressionDefinitionC/KCall price/Strike priceS/KIndex price/Strike priceτTime to maturity+Addition-Subtraction*Multiplication%Protected division:x%y=1 if y=0;   x%y=x%y otherwiselnProtected natural log:ln(x)=ln(|x|)ExpExponential function:exp(x)=exSqrtProtected square root:x=|x|NcdfNormal cumulative distribution functionΦThe search space for GP is the space of all possible parse trees that can be recursively created from the sets of terminal and function. The terminal set includes the inputs variables, notably, the call option price divided by strike priceC/K, the index price divided by strike price S/K and time to maturity τ. The function set includes unary and binary nodes. Unary nodes consist of mathematical functions, notably, cosinus function (cos), sinus function (sin), log function (ln), exponential function (exp), square root function (), and the normal cumulative distribution function (Φ). Binary nodes consist of the four basic mathematical operators, notably, addition (+), subtraction (-), multiplication (×), and division (%). The basic division operation is protected against division by zero and the log and square root functions are protected against negative arguments. The basic measure of model accuracy used is the mean squared error (MSE) between the target (yt) and forecasted (ŷt) output volatility, computed as follows:
(2.1)MSE=1N∑t=1N(yt-ŷt)2,
where yt is the implied volatility computed using the BS formula (C=SN(d1)-Ke-rτN(d2), d1=(ln(S/K)+(r+0.5σ2)τ)/στ, d2=d1-στ.),ŷt is the generated GP volatility, and N is the number of data sample.The implementation of genetic programming involves different parameters that determine its efficiency. A common approach in tuning a GP is to undertake a series of trial and error experiments before making parameter choices for the final GP runs. A two-step approach is used. First, the optimal set of genetic parameters is determined based on time series data. By varying parameters for the genetic programs when training them on the first sample data, their performance is tested on a test data set from a later date. Various genetic programs are tested. Each program is run ten times with ten different random seeds. The choice of the best genetic program is made according to the mean and median of mean squared errors (MSEs) for training and testing sets. Second, the optimal set of parameters, relative to the best genetic program selected in the first step, are used for training the genetic program separately on each of the first 9 samples using ten different seeds and testing the program's performance only on the test data from the immediately following date. This limits the effects of nonstationarities and avoid data snooping. The same genetic parameters and random seeds used for time series data are applied to moneyness-time to maturity data in an attempt to analyze the effect of data' group on the GP performance; each moneyness-time to maturity class is subdivided into training set and test set. Parameters used in this paper are listed in Table2.Table 2
Summary of GP settings and parameters.
Method of generation:Ramped half and halfPopulation size:100Offspring size:200Generations number:400Maximum depth of new individual:6Maximum depth of the tree:17Selection algorithm:Tournament selection (4)Crossover probability:60%Mutation probability:40%Branch mutation:20%Point Mutation:10%Expansion mutation:10%Replacement strategy:Comma replacementStopping criterion:Maximum number of generationsAn initial population of 100 individuals is randomly generated using random number generator algorithm with random seeds. The generative method used for initialization is the ramped half and half as detailed in Koza [18]. This method involves generating an equal number of trees using a maximum initial depth that ranges from 2 to 6. For each level of depth, 50% of the initial trees are generated via the full method and the other 50% are generated via the grow method. A maximum size of tree measured by depth is 17. This is a popular number used to limit the size of tree (Koza [18]). It is large enough to accommodate complicated formulas and works in practice. Based on the fitness criterion, the selection of the individuals for reproduction is done with the tournament selection algorithm. A group of individuals is selected from the population with a uniform random probability distribution. The fitness values of each member of this group are compared and the actual best is selected. The size of the group is given by the tournament size which is equal here to 4. The crossover operator is used to generate about 60% of the individuals in the population, while the mutation operator is used to generate about 40% of the population. 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 subtree rooted at that node with a new randomly generated subtree. Expansion mutation operator randomly selects a terminal node in the tree, and then replaces it with a new randomly-generated subtree. Branch mutation is applied with a rate of 20%; point and expansion mutations are applied with a rate of 10% each. The method of replacing parents for the next generation is comma replacement strategy, which selects the best offspring to replace the parents. Of course, it assumes that offspring size is higher than parents size. The stopping criterion is the maximum number of generations 400.


## 3. Findings and Results Analysis

### 3.1. Selection of the Best Genetic Programming-Implied Volatility  Forecasting Models

Selection of the best generated GP volatility model, relative to each training set, for time series and moneyness-time to maturity data, is made according to the training and test MSE. In definitive, nine generated GP volatility models are selected for time series (M1S1⋯M9S9) and similarly nine generated GP volatility models are selected for moneyness-time to maturity classification (M1C1⋯M9C9). To assess the internal and external accuracy of these generated GP volatility models, we use the MSE total computed for the enlarged sample and the MSE out-of-sample computed for external samples to the training one, as performance's measures. Figure 2 describes the evolution's pattern of the squared errors for these volatility models.Evolution of the squared errors for total sample of the best generated GP volatility models relative to time series (a) and moneyness-time to maturity (b).
(a)(b)Performance of the generated GP volatility models according to MSE total and out-of-sample for the time series samples (a) and the moneyness-time to maturity classes (b) and (c).
(a)(b)(c)It appears throughout Figure2 that for the enlarged sample, the nine generated GP volatility models relative to each data group display different patterns of forecasting errors. It is important to note that each model fits well to the data on which it is trained. Furthermore, the performance of the models is not uniform. Total errors are higher for the moneyness-time to maturity classes than for the time series samples. Time series models are adaptive not only to training samples but also to the enlarged sample. In contrast, the moneyness-time to maturity models are adaptive to training classes, but not all to the enlarged sample. The same patterns of forecasting errors are found using out-of-sample errors as performance criterion. Although we have used the same data, the same genetic parameters, and the same random seeds, the GP has generated new models adaptive to the moneyness-time to maturity classes, which are different from the time series models. This implies that the data group has an effect on the GP results. In other words, the GP fits well to the data on which is trained.Figure3 shows that, for the time series samples, the generated GP volatility model M4S4 has the smallest MSE in enlarged sample and out-of-sample. For the moneyness-time to maturity classes, the
generated GP volatility models M4C4 and M6C6 seem to be more accurate in forecasting implied volatility than the other models. They present near MSE on the enlarged sample and out-of-sample. It is interesting to notice that these models belong, respectively, to ATMST and ATMLT classes. This can be explained by the fact that ATM options are generally more actively traded than other options. Theoretically, Corrado and Miller Jr. [36] have shown that implied volatility from ATM options is the most efficient and also unbiased estimator of the exante underlying volatility when return volatility is stochastic.Volatility smile of the generated GP volatility models against BS model for short-term call options (a), medium term call options (b) and long-term call options (c) at different dates.
(a)(b)(c)Based on the MSE total and the MSE out-of-sample as performance criteria, the generated GP volatility models M4S4, M4C4, and M6C6 are selected.Table3 reports the best generated GP volatility models relative to time series samples and moneyness-time to maturity classes and shows that the time series model M4S4 seems to be more performing than moneyness-time to maturity models M4C4 and M6C6 for the enlarged and out-of-sample samples. The decoding of these models yields the following generated GP volatility forecasting formulas:
(3.1)M4S4:σGP=exp[(ln(Φ(CK))*τ-2*CK+SK)-cos(CK)],M4C4:σGP=|CK*ln(|ln(ln(ln(SK-CK)+(SK+τ))+SK)|)|,M6C6:σGP=((C/K)/(S/K))-ln(S/K)+(sin(C/K)/(S/K))Φ(τ*τ).Table 3
Selection of the best generated GP volatility models relative to time series and moneyness-time to maturity in terms of the MSE total and the MSE out-of-sample (The numbers in parentheses are standard errors of the corresponding measures of the generated GP volatility models).
ModelsMSE totalMSE out-of-sampleM4S40,001444 (0,002727)0,001539 (0,002788)M4C40,001710 (0,0046241)0,001944 (0,005108)M6C60,002357 (0,0040963)0,002458 (0,00415354)A detailed examination of the trees shows that the volatility generated by GP is function of all the inputs used, namely, the call option price divided by strike priceC/K, the index price divided by strike price S/K and time to maturity τ. The four basic mathematical operators (+,-,*,/), the square root function, and the log function appear in all genetic trees. The normal cumulative distribution function is used only in trees relative to M4S4 and M6C6 models. The exponential function, which represents the root node, and the cosinus function appear only in the tree relative to M4S4. The sinus function is used only in the tree relative to M6C6.The implied volatility is defined as the standard deviation which equates the BS model price to the observed market option price. Since there is no explicit formula available to compute directly the implied volatility, the latter can be obtained by inverting the BS model. On the contrary, the GP offers explicit formulas which can compute directly the implied volatility expressed as a function of time to maturity, index, strike, and option prices. Therefore, the generated GP volatilities take account of moneyness and time to maturity. This contradicts the BS assumption of constant volatility. Furthermore, the comparison of the GP forecasting models reveals that the volatility generated by the M6C6 model can be negative, whereas this phenomenon did not exist in the M4S4 and M4C4 models since the implied volatility is computed using the square root and the exponential functions as the root nodes. It appears that even if the variables are nulls, the implied volatility of M4S4 model is always different from zero. This means that a call option can never get a zero value.


### 3.2. Implied Volatility Patterns: Volatility's Smile, Term Structure,  and Surface

Implied volatility patterns of the generated GP volatility models are compared across both moneyness and maturity, as was done in Rubinstein [4]. The call option-implied volatilities illustrated in Figure 4 conform to the general shape hypothesized by Dumas et al. [10] for S&P500 options since the 1987 crash. ITM call options are generally trading at higher implied volatilities than OTM call options. In other words, the implied volatility, at a fixed maturity, is an increasing (decreasing) nonlinear function of moneyness (strike price). The volatility smile is consistent with the market view that the market falls more quickly than it rises. Options writers are concerned with the direction and nature of price movements. Call option writers prefer prices that creep upward and gap downward, while put option writers like the reverse. Option traders aware of a changing volatility may be able to take advantage of the knowledge. For example, if there are indications that the skew may be flattening, one strategy could involve selling at the higher volatility and reversing the position at the lower volatility. Brown [37] has examined the implied volatility smile for call and put options trading on the share price index (SPI) futures contract on the Sydney Futures Exchange (SFE), and offered a possible explanation for its shape. Implied volatilities are viewed as prices reflective of the willingness of market participants to take on and lay off the risks involved in trading volatility, and other risks not priced by the model. Figure 4 shows that the M4S4 model's implied volatility pattern smiles the least for short-term, medium-term, and long-term options. However, the M4C4, M6C6, and BS models still exhibit a significant smile for all maturities. The volatility smiles are the strongest for short-term options, indicating that short-term options are the most severely mispriced by these models. The implied volatility models yield, for each given maturity category, the same implied volatility values for some options. For instance, for long-term options, the implied volatility curves of M4C4, M6C6, and BS models intersect all at approximately the same level of moneyness. Options with moneyness of 0.98 (ATM options) are priced at about the same volatility for long-term maturity.Implied volatility of the generated GP volatility and BS models against degree of moneyness for different times to maturity at the date 27/01/2003.
(a)(b)(c)(d)Figure5 shows that the smile effect is most clear for short-term options. For long-term options, the smile effect is almost vanished. In fact, it is clear that the degree of curvature of M4S4, M4C4, M6C6, and BS models is most extreme when the options are closest to expiration and becomes more linear with longer expirations. This effect was found in [38] for DAX index options from January 1999 to February 2003. It corresponds to the empirical fact that near the expiry, the smile effect becomes stronger. The fact that the implicit volatility smile is more pronounced for short-term options contrary to long-term options is called the vanishing smile. This effect is caused by the smaller time value of short-term options which renders them more sensitive to forecasting errors. As opposed to longer term options which are strongly affected by the value of volatility, short-term options are not much affected by the value of volatility and their price cannot give much information on short-term volatility, which therefore leads to a large uncertainty on short-term volatility. A further result is the changing shape of the smile as the time to maturity increases, except the genetic model M4S4. The smiles of the M4S4 model are the flattest for all maturities, whereas the smiles of the M6C6 model are the steepest for all maturities. The implied volatilities of M4S4, M4C4, and BS models are higher for short-term maturity than for long-term maturity. The M6C6 model has an implied volatility that is higher for long-term maturity than for short-term maturity.Volatility term structure of the generated GP volatility models against BS model for OTM call options (a), ATM call options (b), and ITM call options (c).
(a)(b)(c)In contrast to the numerous studies relating models to volatility smiles, there has been relatively little research reported that relates model choice to the term structure of volatility (TSV). One reason for this is that the TSV is a far more complex phenomenon. While the smile has the same general shape for all maturities, the TSV can vary widely with the specific choice of strike price. Figure6 illustrates the volatility term structure of the generated GP volatility models against BS model for all moneyness categories. As Figure 6 shows, the implied volatility of M4S4, M4C4, and BS models is a decreasing function of time to maturity for all moneyness classes except the genetic model M6C6. Rubinstein  [4] has shown that for out-of-the money or deep out-of-the money call, the shorter the time to expiration, the higher its implied volatility. The same conclusion was found for at the money call in the second period from 24 October 1977 to 31 October 1978. For ATM call options, all the curves intersect at the same time to maturity, and in particular the curves of M4C4 and BS are very close together. The differences in ATM-implied volatilities across maturities can be attributed to the variation of implied volatility over the full range of traded strikes.Volatility surfaces of the generated GP volatility models against BS model for short-term OTM call options at 30/01/2003 (a), short-term ATM call options at 30/04/2003 (b), and short-term ITM call options at 20/03/2003 (c).
(a)(b)(c)The dynamic properties of implied volatility time series have been studied in many markets by various authors. Most of these studies either focus on the term structure of ATM-implied volatilities or study separately implied volatility smiles for different maturities [39–46]. Other studies focus on the joint dynamics of all implied volatilities quoted on the market, looking simultaneously at all available maturity and moneyness values [38, 47]. In order to investigate the generated GP volatility functions, a three-dimensional graph of implied volatility against moneyness and maturity of the option, called the volatility surface, is plotted in Figure 6. This volatility surface plays an important role in the pricing of options in practice. The implied volatility surface gives a snapshot of market prices of vanilla options: specifying the implied volatility surface is equivalent to specifying prices of all vanilla options quoted on the strike price and expiration date at a given date. While the BS model predicts a flat profile for the implied volatility surface, it is a well documented empirical fact that it exhibits both a nonflat strike and term structure [8, 10, 39, 40, 48, 49].It is, therefore, interesting to see if the genetic models reproduce realistic instantaneous profiles for the surface. Figure6 regroups all these candidate solutions on the same graph at the same date. As seen in Figure 6, the best performing volatility functions obtained have quite different surfaces. Volatility varies across the maturity and the moneyness. The implied volatility surfaces appear to display some of the features reported by Derman and Kani [50], and Corrado and Su [51], and there is some degree of negative skewness. Graphs in Figure 6 indicate that the volatility smiles are more asymmetric for the M4C4, M6C6, and BS than the M4S4 model which produces quasi-symmetric smile such as stochastic volatility models, the class of parametric models most commonly used to generate smile effects. Options are in fact capturing skewness of expectations especially for the M4C4, M6C6, and BS models. For the S&P500 index, negative skewness in the distribution appearing since the 1987 stock market crash reflects either a higher probability attached to a future crash, or higher valuation of returns in the event of a crash. In a stochastic volatility model, this skewness results from a negative correlation between shocks to price and shocks to volatility, as discussed in Hull and White [52] and Heston [16]. Jackwerth and Rubinstein [53] showed that this negative skewness has been roughly constant in S&P500 options since about one year after 1987 crash.


## 3.1. Selection of the Best Genetic Programming-Implied Volatility  Forecasting Models

Selection of the best generated GP volatility model, relative to each training set, for time series and moneyness-time to maturity data, is made according to the training and test MSE. In definitive, nine generated GP volatility models are selected for time series (M1S1⋯M9S9) and similarly nine generated GP volatility models are selected for moneyness-time to maturity classification (M1C1⋯M9C9). To assess the internal and external accuracy of these generated GP volatility models, we use the MSE total computed for the enlarged sample and the MSE out-of-sample computed for external samples to the training one, as performance's measures. Figure 2 describes the evolution's pattern of the squared errors for these volatility models.Evolution of the squared errors for total sample of the best generated GP volatility models relative to time series (a) and moneyness-time to maturity (b).
(a)(b)Performance of the generated GP volatility models according to MSE total and out-of-sample for the time series samples (a) and the moneyness-time to maturity classes (b) and (c).
(a)(b)(c)It appears throughout Figure2 that for the enlarged sample, the nine generated GP volatility models relative to each data group display different patterns of forecasting errors. It is important to note that each model fits well to the data on which it is trained. Furthermore, the performance of the models is not uniform. Total errors are higher for the moneyness-time to maturity classes than for the time series samples. Time series models are adaptive not only to training samples but also to the enlarged sample. In contrast, the moneyness-time to maturity models are adaptive to training classes, but not all to the enlarged sample. The same patterns of forecasting errors are found using out-of-sample errors as performance criterion. Although we have used the same data, the same genetic parameters, and the same random seeds, the GP has generated new models adaptive to the moneyness-time to maturity classes, which are different from the time series models. This implies that the data group has an effect on the GP results. In other words, the GP fits well to the data on which is trained.Figure3 shows that, for the time series samples, the generated GP volatility model M4S4 has the smallest MSE in enlarged sample and out-of-sample. For the moneyness-time to maturity classes, the
generated GP volatility models M4C4 and M6C6 seem to be more accurate in forecasting implied volatility than the other models. They present near MSE on the enlarged sample and out-of-sample. It is interesting to notice that these models belong, respectively, to ATMST and ATMLT classes. This can be explained by the fact that ATM options are generally more actively traded than other options. Theoretically, Corrado and Miller Jr. [36] have shown that implied volatility from ATM options is the most efficient and also unbiased estimator of the exante underlying volatility when return volatility is stochastic.Volatility smile of the generated GP volatility models against BS model for short-term call options (a), medium term call options (b) and long-term call options (c) at different dates.
(a)(b)(c)Based on the MSE total and the MSE out-of-sample as performance criteria, the generated GP volatility models M4S4, M4C4, and M6C6 are selected.Table3 reports the best generated GP volatility models relative to time series samples and moneyness-time to maturity classes and shows that the time series model M4S4 seems to be more performing than moneyness-time to maturity models M4C4 and M6C6 for the enlarged and out-of-sample samples. The decoding of these models yields the following generated GP volatility forecasting formulas:
(3.1)M4S4:σGP=exp[(ln(Φ(CK))*τ-2*CK+SK)-cos(CK)],M4C4:σGP=|CK*ln(|ln(ln(ln(SK-CK)+(SK+τ))+SK)|)|,M6C6:σGP=((C/K)/(S/K))-ln(S/K)+(sin(C/K)/(S/K))Φ(τ*τ).Table 3
Selection of the best generated GP volatility models relative to time series and moneyness-time to maturity in terms of the MSE total and the MSE out-of-sample (The numbers in parentheses are standard errors of the corresponding measures of the generated GP volatility models).
ModelsMSE totalMSE out-of-sampleM4S40,001444 (0,002727)0,001539 (0,002788)M4C40,001710 (0,0046241)0,001944 (0,005108)M6C60,002357 (0,0040963)0,002458 (0,00415354)A detailed examination of the trees shows that the volatility generated by GP is function of all the inputs used, namely, the call option price divided by strike priceC/K, the index price divided by strike price S/K and time to maturity τ. The four basic mathematical operators (+,-,*,/), the square root function, and the log function appear in all genetic trees. The normal cumulative distribution function is used only in trees relative to M4S4 and M6C6 models. The exponential function, which represents the root node, and the cosinus function appear only in the tree relative to M4S4. The sinus function is used only in the tree relative to M6C6.The implied volatility is defined as the standard deviation which equates the BS model price to the observed market option price. Since there is no explicit formula available to compute directly the implied volatility, the latter can be obtained by inverting the BS model. On the contrary, the GP offers explicit formulas which can compute directly the implied volatility expressed as a function of time to maturity, index, strike, and option prices. Therefore, the generated GP volatilities take account of moneyness and time to maturity. This contradicts the BS assumption of constant volatility. Furthermore, the comparison of the GP forecasting models reveals that the volatility generated by the M6C6 model can be negative, whereas this phenomenon did not exist in the M4S4 and M4C4 models since the implied volatility is computed using the square root and the exponential functions as the root nodes. It appears that even if the variables are nulls, the implied volatility of M4S4 model is always different from zero. This means that a call option can never get a zero value.


## 3.2. Implied Volatility Patterns: Volatility's Smile, Term Structure,  and Surface

Implied volatility patterns of the generated GP volatility models are compared across both moneyness and maturity, as was done in Rubinstein [4]. The call option-implied volatilities illustrated in Figure 4 conform to the general shape hypothesized by Dumas et al. [10] for S&P500 options since the 1987 crash. ITM call options are generally trading at higher implied volatilities than OTM call options. In other words, the implied volatility, at a fixed maturity, is an increasing (decreasing) nonlinear function of moneyness (strike price). The volatility smile is consistent with the market view that the market falls more quickly than it rises. Options writers are concerned with the direction and nature of price movements. Call option writers prefer prices that creep upward and gap downward, while put option writers like the reverse. Option traders aware of a changing volatility may be able to take advantage of the knowledge. For example, if there are indications that the skew may be flattening, one strategy could involve selling at the higher volatility and reversing the position at the lower volatility. Brown [37] has examined the implied volatility smile for call and put options trading on the share price index (SPI) futures contract on the Sydney Futures Exchange (SFE), and offered a possible explanation for its shape. Implied volatilities are viewed as prices reflective of the willingness of market participants to take on and lay off the risks involved in trading volatility, and other risks not priced by the model. Figure 4 shows that the M4S4 model's implied volatility pattern smiles the least for short-term, medium-term, and long-term options. However, the M4C4, M6C6, and BS models still exhibit a significant smile for all maturities. The volatility smiles are the strongest for short-term options, indicating that short-term options are the most severely mispriced by these models. The implied volatility models yield, for each given maturity category, the same implied volatility values for some options. For instance, for long-term options, the implied volatility curves of M4C4, M6C6, and BS models intersect all at approximately the same level of moneyness. Options with moneyness of 0.98 (ATM options) are priced at about the same volatility for long-term maturity.Implied volatility of the generated GP volatility and BS models against degree of moneyness for different times to maturity at the date 27/01/2003.
(a)(b)(c)(d)Figure5 shows that the smile effect is most clear for short-term options. For long-term options, the smile effect is almost vanished. In fact, it is clear that the degree of curvature of M4S4, M4C4, M6C6, and BS models is most extreme when the options are closest to expiration and becomes more linear with longer expirations. This effect was found in [38] for DAX index options from January 1999 to February 2003. It corresponds to the empirical fact that near the expiry, the smile effect becomes stronger. The fact that the implicit volatility smile is more pronounced for short-term options contrary to long-term options is called the vanishing smile. This effect is caused by the smaller time value of short-term options which renders them more sensitive to forecasting errors. As opposed to longer term options which are strongly affected by the value of volatility, short-term options are not much affected by the value of volatility and their price cannot give much information on short-term volatility, which therefore leads to a large uncertainty on short-term volatility. A further result is the changing shape of the smile as the time to maturity increases, except the genetic model M4S4. The smiles of the M4S4 model are the flattest for all maturities, whereas the smiles of the M6C6 model are the steepest for all maturities. The implied volatilities of M4S4, M4C4, and BS models are higher for short-term maturity than for long-term maturity. The M6C6 model has an implied volatility that is higher for long-term maturity than for short-term maturity.Volatility term structure of the generated GP volatility models against BS model for OTM call options (a), ATM call options (b), and ITM call options (c).
(a)(b)(c)In contrast to the numerous studies relating models to volatility smiles, there has been relatively little research reported that relates model choice to the term structure of volatility (TSV). One reason for this is that the TSV is a far more complex phenomenon. While the smile has the same general shape for all maturities, the TSV can vary widely with the specific choice of strike price. Figure6 illustrates the volatility term structure of the generated GP volatility models against BS model for all moneyness categories. As Figure 6 shows, the implied volatility of M4S4, M4C4, and BS models is a decreasing function of time to maturity for all moneyness classes except the genetic model M6C6. Rubinstein  [4] has shown that for out-of-the money or deep out-of-the money call, the shorter the time to expiration, the higher its implied volatility. The same conclusion was found for at the money call in the second period from 24 October 1977 to 31 October 1978. For ATM call options, all the curves intersect at the same time to maturity, and in particular the curves of M4C4 and BS are very close together. The differences in ATM-implied volatilities across maturities can be attributed to the variation of implied volatility over the full range of traded strikes.Volatility surfaces of the generated GP volatility models against BS model for short-term OTM call options at 30/01/2003 (a), short-term ATM call options at 30/04/2003 (b), and short-term ITM call options at 20/03/2003 (c).
(a)(b)(c)The dynamic properties of implied volatility time series have been studied in many markets by various authors. Most of these studies either focus on the term structure of ATM-implied volatilities or study separately implied volatility smiles for different maturities [39–46]. Other studies focus on the joint dynamics of all implied volatilities quoted on the market, looking simultaneously at all available maturity and moneyness values [38, 47]. In order to investigate the generated GP volatility functions, a three-dimensional graph of implied volatility against moneyness and maturity of the option, called the volatility surface, is plotted in Figure 6. This volatility surface plays an important role in the pricing of options in practice. The implied volatility surface gives a snapshot of market prices of vanilla options: specifying the implied volatility surface is equivalent to specifying prices of all vanilla options quoted on the strike price and expiration date at a given date. While the BS model predicts a flat profile for the implied volatility surface, it is a well documented empirical fact that it exhibits both a nonflat strike and term structure [8, 10, 39, 40, 48, 49].It is, therefore, interesting to see if the genetic models reproduce realistic instantaneous profiles for the surface. Figure6 regroups all these candidate solutions on the same graph at the same date. As seen in Figure 6, the best performing volatility functions obtained have quite different surfaces. Volatility varies across the maturity and the moneyness. The implied volatility surfaces appear to display some of the features reported by Derman and Kani [50], and Corrado and Su [51], and there is some degree of negative skewness. Graphs in Figure 6 indicate that the volatility smiles are more asymmetric for the M4C4, M6C6, and BS than the M4S4 model which produces quasi-symmetric smile such as stochastic volatility models, the class of parametric models most commonly used to generate smile effects. Options are in fact capturing skewness of expectations especially for the M4C4, M6C6, and BS models. For the S&P500 index, negative skewness in the distribution appearing since the 1987 stock market crash reflects either a higher probability attached to a future crash, or higher valuation of returns in the event of a crash. In a stochastic volatility model, this skewness results from a negative correlation between shocks to price and shocks to volatility, as discussed in Hull and White [52] and Heston [16]. Jackwerth and Rubinstein [53] showed that this negative skewness has been roughly constant in S&P500 options since about one year after 1987 crash.


## 4. Conclusion

In this study, we have used GP to forecast implied volatility from the S&P500 index options. Based on the MSE total and out-of-sample as performance criteria, the generated GP volatility models M4S4, M4C4, and M6C6 can be considered to the best selected models. According to total and out-of-sample forecasting errors, the time series model M4S4 seems to be more accurate in forecasting-implied volatility than moneyness-time to maturity models M4C4 and M6C6. In terms of volatility patterns, the M4S4 model smiles the least for all maturity options and generates a quasiflat and symmetric implied volatility surface. Our results suggest some interesting issues for further investigation. First, the 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 hedging. Finally, the GP can be extended to allow for dynamic selection techniques applied to the enlarged sample. We believe these extensions are of interest for application and will be object of our future works.



---

*Source: 179230-2009-08-31.xml*