Introduction

In the last decades, many stochastic volatility models (SVMs) for option pricing have arisen as an alternative (or generalization) to the standard Black-Scholes-Merton (BSM) model (^{Black and Scholes, 1973}; and ^{Merton, 1973}). Among the main limitations of the BSM model is the assumption of constant volatility. Stochastic volatility models undertake this problem under the assumption that prices and their volatility are stochastic processes affected by several risk factors. In this regard, SVMs handle more consistent prices. Unlike the BSM model, the SVMs considers stylized facts, such as implied volatility smile and skewness, as observed in derivative markets. The latter handle the leverage effect, i.e., volatility rises after price declines. That is, the asset price and its volatility are negatively correlated. The best known SVMs are those from ^{Hull and White (1987)}, ^{Stein and Stein (1991)}, ^{Heston (1993)}, and ^{Schöbel and Zhu (1999)}. In particular, ^{Hull and White’s (1987)} approach expand the asset price with a Taylor series function up to terms of third order, and ^{Heston’s (1993)} model provides a closed formula that requires numeric methods.¹

A main issue of option pricing models is their capability of generating volatility structures that are consistent with the expiration date and the strike price in order to capture the behavior of the underlying asset return distribution. There are three important features to be considered once a model is proposed: tail risk, option pricing, and implied volatility. Tail risk is empirically defined as the probability that the underlying asset return will move more than three standard deviations from its mean. The BSM model underestimates tail risk since the forecasted probability for extreme variations is low. This is because in the BSM model the tails of the asset return distribution are driven by a normal distribution. Thus, in this context the models of ^{Hull and White (1987)}, ^{Stein and Stein (1991)}, and ^{Heston (1993)} exhibit higher performance than the BSM model allowing estimating the probabilities for unexpected falls in the asset price (^{Gulisashvili and Stein, 2010}).

Volatility is an unobservable variable that measures the relative changes in an asset price. Often the volatility is also seen as the instantaneous standard deviation of the asset returns. It important to point out that through periods of low volatility, asset prices do not change very much; however, in periods of high volatility, significant variations are shown for the asset prices. In the BSMmodel, volatility is constant and positive until the expiration date. In contrast, a SVM, volatility is described by a non-negative stochastic process. For example, in ^{Hull and White’s (1987)} model, volatility dynamics follows a geometric Brownian motion. In ^{Stein and Stein’s (1991)} model, volatility follows an Ornstein-Uhlenbeck process or the absolute value. While in ^{Heston’s (1993)} model, volatility is driven by the square root of the ^{Cox-Ingersoll-Ross’ (1985)} process.

The term “implied volatility” appeared for the first time in the work of ^{Latané and Rendleman (1976)}, called “implied standard deviation.” These authors studied the standard deviations of returns on assets implied in observed call options prices when an investor values options prices under the standard BSM model. When implied volatility is obtained with options price data, and the strike price and expiration date are plot, the implied volatility surface is not flat, and changes along with the strike price and expiration date. This is what is known as implied volatility smile.

In relation to the estimation of the parameters of stochastic volatility, there are three main statistical techniques: maximum likelihood (ML), the generalized method of moments (GMM), and the market-based (MB) method. The latter minimizes the spread or market volatilities obtained from the theoretical model and the data market, e.g., using a loss function (^{Bakshi, Cao, and Chen, 1997}; and ^{Bams, Lehnert, and Wolff, 2009}). In this context, innovations in the computational science field offer considerable advantages to the financial industry, e.g., evolutionary computation, which has developed optimization methods based on evolutionary algorithms, representing an alternative to traditional optimization methods (either deterministic or stochastic) whenever the latter fail or the computational cost is high. Genetic algorithms are computational algorithms based on Darwinian evolution, which use operators such as reproduction, mutation, and selection to evolve a set of candidates or possible solutions to a given optimization problem. ^{Gilli and Schumann (2011)} provided a method using differential evolution to calibrate the Heston model and proposed only one loss function with artificially generated data. A notable finding is the utilization of alternative methods to decrease the risk model and the selection criteria of parameters. On the other hand, ^{Cui, del Baño and Germano (2017)} propose an efficient algorithm to calibrate the Heston model with nonlinear least squares using the ^{Levenberg (1944)} method by using the alternative representation of the characteristic function and avoiding discontinuities.

^{Soler and Kizys (2017)} review some literature related to the use of evolutionary algorithms to solve classic and emerging problems in the field of finance. A non-exhaustive list of examples includes portfolio optimization, index tracking, enhanced indexation, credit risk, stock investments, financial project scheduling, option pricing, feature selection, bankruptcy and financial distress prediction, and credit risk assessment, they conclude that the use evolutionary algorithms in finance is becoming increasingly more popular among researchers from different communities in order to make efficient decisions. In ^{Liu et al (2019)} a method to calibrate stochastic volatility models using artificial neural networks is proposed that includes three phases: training, prediction and calibration. Some relevant results are the efficiency and precision of the calculations when compared with other methods such as differential evolution that avoids getting stuck in a local minimum. In ^{Palmer (2019)} it is addressed a number of problems in computational finance, characterized by their high demand for computation time and resources. One of them is the calibration of the Heston model with variants of evolutionary algorithms such as particle swarm optimization (PSO) and differential evolution (DE). State-of-the-art DE algorithms were found to provide excellent calibration performance, and the prior use of rudimentary DE in the literature underestimates the use of these methods.

This paper is aimed at calibrating the parameters of Heston’s SVM with several loss functions by using two different methods: the first one minimizes a nonlinear function of several variables with constraints, and the second one employs the differential evolution algorithm. Both methods are applied to a set of implied volatilities on the Mexican Stock Exchange Index (IPC Spanish acronym for Indice de Precios y Cotizaciones) with four maturities and twenty-eight strike prices. The selection criterion for the parameters is minimizing the value of the mean square error of the implied volatility. It important to point out that the mean squared loss function weighted by the square of the vega reaches the minimum value by the first method with less error and time. Moreover, the empirical results show that for both methods the adjustment of implied volatilities is better for options with long-term maturities than those for short-term maturities, and the value of rho is negative, which means that the slope of the implied volatility curve is negative and is closely related to the asymmetry of the underlying asset returns distribution. The main conclusion is that a change in the slope of implied volatility reflects the supply and demand effects of options. Thus, high premium will move in the direction that market participants expect the underlying to move.

The structure of the paper is as follows: the next section provides a short review of SVMs; section 3 describes briefly the Heston model; section 4 describes the calibration of parameters using several loss functions; through section 5, the differential evolution algorithm to calibrate Heston’s model parameters is defined; in section 6, the Heston model parameters are calibrated by using loss functions according to the two methods of the previous sections, both methods are implemented for a set of implied volatilities of the IPC. Finally, in section 7, conclusions are given.

A Short Review of Stochastic volatility models

According to ^{Garman’s (1976)} theoretical proposal, an option price with stochastic volatility satisfy a second-order partial differential equation (PDE) with two fundamental variables: strike price and volatility. Volatility is not a trading asset, so arbitrage cannot eliminate volatility risk. Therefore, the market price is included in the PDE. Usually, the risk premium is zero or a constant volatility proportion. Final and boundary conditions for the option contracts state the premium. In these option pricing models; underlying volatility is explicit and its parameters are obtained by statistical methods. The PDE solution from which the option price is obtained engender additional problems, since most of the time closed formulas are not available, and the least it can be asked for is an analytical solution as a prerequisite for the empirical work.

^{Wiggins’ (1987)} work solves a PDE with a finite difference method for the logarithmic transformation of a state variable. ^{Hull and White (1987)} worked with a simple power series approximation technique based on the conditional underlying distribution of the average value of the stochastic variance. ^{Stein and Stein (1991)} and ^{Heston (1993)} proposed an analytic approach using the inverse Fourier transform. Of these three methods, the PDE numerical solution is the best known, but also the most computationally expensive. ^{Hull and White’s (1987)} and ^{Stein and Stein’s (1991)} approaches depend on the average stochastic variance distribution and assume instant correlation of state variables equal to zero, and ^{Hull and White’s (1987)} approach is more practical; however, ^{Heston (1993)} claimed that the ^{Stein and Stein (1991)} model works with the inverse Fourier transform even if the correlation is different from zero. On the other hand, ^{Schöbel and Zhu (1999)} developed a model with two correlated Brownian motions and estimated the characteristic function of the logarithmic price through a system of stochastic differential equations. Also, ^{Zhu (2010)} proposed a double root stochastic volatility model akin to ^{Longstaff´s (1989)} model. ^{Date and Islyaev (2015)} proposed an SVM based on a deterministic time structure modified by an escalated random variable. In this case, a closed approximation is obtained for a European-type option premium using higher order Greeks from its volatility.

Revisiting Heston’s SVM

An important feature in the Heston model is that the characteristic functions of the risk neutral probabilities are the solution to a second-order PDE (^{Venegas-Martínez, 2008}). Using these probabilities, a formula is obtained to value a European call option under the assumption that volatility is stochastic. In this case, there is a correlation between the volatility and the underlying price, while the strike price is stated by the put-call parity. There are some complications with Heston’s model because it includes complex variables. Therefore, the option price is obtained by calculating the probability that the call option expires in the money. This probability is obtained by calculating the inverse of the characteristic function for the logarithm of the underlying price. The formalization of the Heston model is as follows. The underlying price Χt is driven by the following stochastic differential equation:

where α denotes the annual average return and B1t is a Brownian motion or Wiener process (normal independent increments with zero mean and variance linearly depending of time). It is important noticing that Heston model uses the variance vt instead of the volatility vt Indeed, the variance follows the ^{Uhlenbeck and Ornstein (1930)} process for volatility, denoted by gt=vt. Thus, dgt=-βgtdt+δdB2t, where B2t is a Wiener process correlated to B1t, i.e., CovdB1t,dB1t=ρ. Using Itô´s lemma with vt=gt2, the variance dynamics is obtained, and in this case, according to ^{Cox, Ingersoll, and Ross (1985)}, the mean reversion process for the variance is:

where κ=2β,θ=δ2/2β, and σ=2δ. If and κ,θ, and σ satisfy the conditions 2κθ>σ2 and v0>0. The above stochastic process is always positive, and the distribution of (1) is a non-central chi-squared. Moreover, the variance always has a mean reversion structure, where κ is the speed of adjustment toward the mean, θ>0, which is the long-term level of vt and σ the instantaneous volatility of the variance. Intuitively, the parameter θ, in (1), is such that if the variance is become greater, then the variance is attracted toward θ. Regarding the estimation of the parameters, we have v0 as an unobservable initial state variable (^{Bakshi et al., 1997}). In this paper, v0 is considered as a parameter to be estimated. To estimate the price of a European-type call option it is used:

Here Πjy,v,τ, stands for the probability that the call option expires in the money XT>K. In this case, yt=ln⁡St (natural log of the asset price), vt is variance at time t, τ=T-1 is the term, and Πjy,v,τ=Pln⁡XT>ln⁡K,j=1,2. If the characteristic functions φ1ϕ;y,v and φ2ϕ;y,v related to Πj are known, it is possible to obtain each Πj from its characteristic function by using the inversion theorem (see ^{Waller, Turnbull, and Hardin, 1995}):

At maturity, when t=T, i.e., τ=0, the probabilities are constraints on the terminal condition Πjy,v,0=1y>ln⁡K, which implies that at maturity, when XT>K, the probability that the call option is in the money is equal to 1. In summary, according to ^{Heston (1993)}, the call option price is:

where Π1 and Π2 are defined in (2), and the characteristic functions are given by:

where i=-1 is the imaginary unity, u 1=1/2, u2 =-1/2, b1= κ+λ-ρσ, and b2= κ+λ. Here, λ is the parameter that represents the volatility premium as a function of the asset price, time and volatility. To obtain the option price, it is necessary to estimate the integrals for Πj in (2), and this leads to a numerical integration method. To estimate the put option price, the put-call parity is used:

Another relevant property of this model is that implied volatilities are estimated by option prices generated according to the model, and these have volatility smile skew. Stock options regularly have a negative slope on the implied volatilities. It is remarkable that when calibrating the parameters of Heston’s model by using market data, the correlation parameter value is negative most of the time.

The correlation parameter ρ determines the skew trend. If ρ>0, there is a positive slope, otherwise there is a negative slope. This effect is shown in Figure 1 by generating implied volatilities using the following values: X =100, r= 0.05, τ = T-t = 0.25, κ=2, θ=0.01, λ=0 and v0=0.01, with a strike range oscillating from 90 to 110. If the volatility values increase the variance, the smile skew increases, i.e. a higher σ parameter increases the implied volatility curvature, this can be seen in Figure 2 for the same values but using ρ = 0.

Source: own elaboration.

Figure 1 Effect of the parameters on implied volatility with different ρ values. 

Source: own elaboration.

Figure 2 Effect of the parameters on implied volatility using different values of σ. 

The parameters κ,θ and v0 manage the smile level. The adjustment speed of the mean reversion, κ, also controls the curve up to some point, and using bigger values of κ, the implied volatility curve flattens, as seen in Figure 3.

Source: own elaboration.

Figure 3 Effect of the parameters on implied volatility using different values of κ. 

When the values of θ and v0 either increase or decrease, a parallel movement can be observed in the volatility smile, as seen in Figures 4 and 5. To graph the plots, the bisection algorithm is used with the BSM model to calculate the set of implied volatilities provided that prices are estimated using the Heston model.

Source: own elaboration.

Figure 4 Effect of the parameters on implied volatility with different values of θ. 

Source: own elaboration.

Figure 5 Effect of the parameters on implied volatility using different values of v0 . 

Loss functions

One alternative method for calibrating the parameters of the Heston model is to use loss functions. This method minimizes the difference amongst market prices and estimated prices, or between the marketimplied volatilities and the Heston model. The set of the estimated parameters Θ corresponds to the values that minimize the value of the corresponding loss function in such a way that the model prices or the implied volatilities are close enough to the potential market prices. To do this, it is necessary to use a minimization algorithm with the following constraints: κ>0,θ>0,σ>0,v0>0,ρ∈-0.9,0.9. In this regard, the loss functions include trading asset prices at a point in time (or its implied volatility) leading to risk neutral parameters in the Heston model. Let us assume that there is a set of stock maturities Mt with τt(t=1,…,T) and a set of strike prices Kk(k=1,…,MK). For each combination maturity-strike (τt,Kk), there is a market price Pτt,Kk=Ptk and the corresponding price obtained by Heston’s model Pτt,Kk;Θ=PtkΘ. Moreover, it is possible to link each option to a weight ωtk. There are several methods of defining a loss function, generally categorized as base prices and base-implied volatilities. The former minimizes the error between market quotes and the estimated parameters by a given model. The error is defined as the squared difference between market quotes and the model estimations, or the difference absolute value. It is also possible to use the relative error, i.e., estimating the parameters by using the loss function of the mean squared error (MSE) and minimizing:

With respect to Θ, where M is the number of quotes. The other commonly used loss function is the so-called relative loss of the mean squared error (RLMSE) given by:

Another way of defining the error is by taking the absolute value in (7) or (8).

A shortcoming of the loss function L1 is that options with short maturity, out of the money, and with low value contribute little to the sum in (7); the optimization process adjusts options in the money to the long term, impacting other options. A solution to this issue is to use only options in the money in such way that in (7) call options are used for strike prices smaller than the underlying price and put options for strike prices greater than the underlying price. Another solution is the RLMSE in (8), but the problem with L2 is the opposite. In fact, provided that Ptk is in the denominator, options on the low-price market saturate the sum in (8), although this excess is repaired by allocating the weights ωtk to the terms in the objective function. It turns out that the selection weight process is arbitrary.

The second category of loss functions is associated with minimizing the error between the implied volatilities and the parameters estimated by a theoretical model. This can be also specified as the relative difference, absolute difference or squared difference between the implied volatilities and those given by a theoretical model. This kind of loss function is sensitive because in the derivatives market, options quotes are made in terms of the implied volatility. The adjustment precision of the model is obtained by comparing the implied volatilities and the model volatility, i.e., the calibration parameters, using the loss function of the MSE of the implied volatility:

where σ^tk=σ^(τt,Kk) and σ^tkΘ=σ^(τt,Kk;Θ) are respectively the market-implied volatilities and the model volatilities. It is also possible to use the relative and absolute versions of L3. The principal drawback of (9) is that it needs exhaustive numeric calculation, i.e., for each iteration, the optimization process estimates all the prices with the Heston model PtkΘ, and then a root algorithm is used as a bisection algorithm to calculate the implied volatility σ^tkΘ in PtkΘ and thus finally is obtained σ^tk,σ^tkΘ2. A solution to this problem is to use the loss function of the weighted MSE times squared vega (^{Christoffersen, Heston, and Jacobs, 2009}), which is an approximation to the function in (9). In this case:

where vtk2 is the option price sensitivity in the BSM model in relation to the market-implied volatility vtk evaluated at maturity τt at strike Kt, defined as:

Here, r is the continuously compounded interest rate, q is the continuously compounded dividend yield, and Φ⋅ is the cumulative distribution function of a standard normal random variable.

The main improvement of the loss function L4 is the time reduction in its execution, although there is a loss of precision in the calibration of parameters. Another advantage is that for each evaluation of the loss function, it is necessary to recalculate the option price. The analytic expression of the characteristic functions guarantees that these calculations are done efficiently.

Differential evolution: an algorithm based on evolutionary computation

While traditional optimization methods (i.e., mathematical programming) are limited tools for finding optimal solutions, heuristic algorithms, most of the time, offer a better alternative. In recent years, the blend of computer science and biology has given rise to a new paradigm for solving highly complex problems, the so-called evolutionary computation. There are three main paradigms that make up evolutionary computation: evolutionary strategies, evolutionary programing, and genetic algorithms. It is worth mentioning that genetic algorithms are the most common method of evolutionary computation used in finance ^{Ponsich, López and Coello (2013)}.²

Genetic algorithms were developed by ^{Holland (1975)} and they use the process of natural selection to find the optimal solution through a series of populations that combine and mutate solutions to reach the optimum. Different optimization problems require different methods of mutation and combination of solutions (^{Bradshaw, Walshaw, Ierotheou, and Parrott, 2009, p. 28}). Genetic algorithms are also known as evolutionary algorithms because they are inspired by the theory of evolution of Charles Darwin, and generally work with operators such as reproduction, crossover, and mutation.

The differential evolution (DE) paradigm was developed by Storn and Price (1995). Like many other optimization algorithms, the DE was motivated by the determination of coefficients using the polynomials of Tchebyshev³ and the optimization of the coefficients in the digital filter. The first publication of DE dates back to ^{Storn and Price (1997)}. The DE algorithm is described in four distinct stages (^{Vollrath and Wendland, 2009}): initial selection of the population, mutation, recombination, and selection. These stages are typical in most evolutionary algorithms; however, the DE is unique in the way the parameter sets are mutated. The stages of DE applied to calibration of Heston’s model are listed as follows:

a) Initial population. Generate an n×MP array of randomly chosen parameter values. In this case, n is the number of parameters and MP is the number of members of the population. The rows of the matrix must satisfy the restrictions of the parameters. Lower and upper limits must be chosen in such a way that the random values for the parameters are feasible, v.g., -1<ρ<1,σ<3 and 0<κ<10,θ>0,v0>0. For Step 2, it is necessary that MP≥4. This implies that with n= 5:

and every member of P is a parameter vector zi=κi,θi,σi,v0i,ρi'.

b) Mutation. For each member zi of the population i=1,…,Mp, randomly select another three members zε1,zε2, and zε3, different from each other and from zi. This requires to select indexes ε1,ε2, and ε3 such that i≠ε1≠ε2≠ε3. For every member, create a new giver member:

where F is a constant mutation factor that controls the size of the mutations and takes values in [0,2]. The giver members and components, xi, are used to produce a candidate member that either replaces zi or not, depending on the outcomes of next steps.

c) Recombination. The candidate element ci is generated from the xij elements of the giver and components of the member zij(j=1,…,5) according to the following: for each element j, a random number of a uniform variable Uij is generated and the candidate element is defined as follows:

where E is a randomly selected integer and R is the denominated crossover ratio. E ensures that the candidate is not the same as the member, i.e., ci≠zi. This is because even if all Uij are small and the inequality Uij≤R is not fulfilled at all, the corresponding element xij will turn into the element cij of the candidate ci. This step is continued for each candidate parameter.

d) Candidate selection. The member zi and the candidate ci are both introduced into the objective function f(z). If the candidate’s value is less than that of the member, the candidate replaces zi as the ith term member of the population P in (11). Thus, the instruction is:

When all the elements of the population have admitted the procedure, the population is renovated, and the next generation starts over again in step 2. Therefore, the algorithm is executed for all generations, G . In this way, the set of parameters is chosen for the member of the population in the final generation with the smallest value of f(z).

To further explain the above procedure, it is important to note that when evaluating a loss function given by equations (7)-(10) prices corresponding to a set of different strikes and maturities are calculated. In this case, the characteristic function φ only depends on the time to maturity, not on the strike price. Therefore, we can previously process the variables that are constant for a given maturity and then calculate the prices for all strike prices for this maturity. The following algorithm summarizes this.

Algorithm 5.1 Computing option prices for a given surface.

1. Set parameters, set C= maturities, set K strikes

2. for τ∈C do

3. compute characteristic function φ

4. for X∈K do

5. compute price for strike X, maturity τ

6. end for

7. end for

8. compute objective function L

Finally, the following algorithm summarizes the proposed methodology.

Algorithm 5.2 Differential Evolution.

1. set parameters M_p, N_g, F, R

2. initialize population Pj,i(z),j=1,…,p,i=1,…,Mp

P=KLo+KU-KLo×rand1,Mp;…θLo+θU-θLo×rand1,Mp;…σLo+σU-σLo×rand1,Mp;…v0Lo+v0U-v0Lo×rand1,Mp;…ρ0Lo+ρ0U-ρ0Lo×rand1,Mp;…

3. for k=1 to N_G do

4. P(x)=P(z)

5. for i=1 to M_P do

6. generate ε1,ε2,ε3∈1,…,MP,ε1≠ε2≠ε3≠i

7. Compute Pi(x)=Pε1(z)+F×Pε2(z)-Pε3(z)

8. for j = 1 to p do

9. If U<R then Pj,i(c)=Pj,i(c) else Pj,i(c)=Pj,i(x)

10. end for

11. If fPi(c)<fPi(x) then Pi(z)=Pi(c) else Pi(z)=Pi(x)

12. end for

13. end for

14. find member with lowest value h=findf=minif(Pi(z))

15. solution = Ph(z)

Empirical Results

This section considers a set of option prices for the Mexican Stock Exchange (BMV). The information was obtained from the Bulletin of the Summary of the Options Market since 06/26/2015, published by the Mexican Derivatives Market (MexDer). Four deadlines are considered: 84, 175, 266, and 357 days. The implied volatilities of the options traded were obtained from the price vendor, and the risk-free rate was acquired from the Mexican Central Bank (Banxico). Table 1 shows the relevant data for the calibration of the parameters.

Table 2 shows the expiration dates, strike prices (K), and implied volatilities of the IPC call options obtained from the price vector. There are 31 strike prices, of which 22 are options that are in the money and the remainder is out of the money. This is important since the loss functions previously presented give different results in the calibration under this assumption.

Figure 6 shows the implied market volatility surface and the market prices of IPC options,computed from the information in the bulletin. Notice the reduction in the differences in prices, excluding options with little premium value; for the calibration the choice of the options is a settlement price greater than $100, i.e. prices from K = 49000 to K= 50000 are excluded. Therefore, there are 28 strike prices and 4 expiration dates, a total of 124 implied volatilities. An implied volatility bias is observed for each maturity level. Implied volatility is greater for small strike prices corresponding to call options in the money and put options out of the money, and implied volatility decreases as the strike price increases. Another feature is that this is more common for options with shorter maturities, and is flattened for options with greater expiration. The monotonic relationship between price and volatility shows that put options out of the money have a higher premium than the call option out of the money.

Source: own elaboration with data from BMV, MexDer and Banxico.

Figure 6 Surface of the implied market volatility and IPC market quotes. 

Table 3 shows the parameters of the Heston model calibrated with several loss functions (MSE, LMSE, and MSE of implicit volatility), the time in seconds that the calibration takes, and the value of each objective function. The initial values of the parameters are as follows: κ=9,θ=0.05,σ=0.3,v0=0.05,ρ=-0.8. In this case, the minimum of a nonlinear function of several variables with constraints⁴ is determined, results of table 3 are obtained with the algorithm developed in ^{Nelder and Mead (1965)}. The set of parameters that is chosen is the one that has the lowest value of MSE of implicit volatility, that is, by using the function L4. In all three cases, it is observed that the value of ρ is negative, which means that the slope of the implied volatility curve is negative, as it can be seen in Figure 7.

Source: own elaboration with data from BMV, MexDer and Banxico.

Figure 7 Implied market and Heston volatilities estimated by function L4. 

The above results show that the adjustment of implied volatilities is better for options with long-

term maturities than for short-term maturities, and the time it takes to adjust is much lower compared to DE, as confirmed in Table 4.

Figure 8 shows the surface of implied market volatility and model volatility generated with the parameters obtained from L4. In this case, 31 maturity terms are considered, varying from 84 to 357 days, and there are 271 strike prices, which range from every 50 points from K=35000 up to K=48500. The importance of the surface lies in the fact that while bonds are characterized by maturity, options are distinguished by their term at maturity and their strike price, for this reason requiring a surface instead of a curve. In addition, the yield curve on any given day is a concise description of bond market prices. Another usefulness of the surface is that it is feasible to value an exotic option with the information generated by the volatility surface; in the case that there are deadlines at maturity and strike prices that do not match the deadlines of the surface, an interpolation bilinear process is carried out.

Source: own elaboration with data from BMV, MexDer and Banxico.

Figure 8 Surface of implied market volatility and model volatility using the parameters obtained from function L4. 

Table 4 shows the parameters of the Heston model calibrated with the three previous loss functions but using the DE algorithm. Unlike ^{Vollrath and Wendland (2009)}, we chose the number of population members equal to 45 times the number of parameters for calibration, in such a way that M _P = 225, the probability of the crossover ratio R = 0.5, the mutation factor F = 0.8, and the number of generations G = 600. The initial values of the parameters for the calibration are as follows: κ=9,θ=0.05,σ=0.3,v0=0.05,ρ=-0.8 The set of parameters that is chosen is the one with the lowest value of the quadratic error of implied volatility which is the function L4. In all three cases, it can be observed that the value is negative, which means that the slope of the implied volatility curve is negative, as it can be seen in Figure 9.

Source: own elaboration with data from BMV, MexDer and MexDer.

Figure 9 Implied market and Heston volatilities estimated by function L4 using differential evolution. 

The above results indicate that the adjustment of implied volatilities is somewhat better for

options with long-term maturities than for short-term maturities, and the time it takes to adjust is much greater in comparison with loss functions.⁵

Finally, Figure 10 shows the adjustment for the area of implied market volatilities and the surface implied market volatility, generated with the parameters of L4 using DE. As it can be seen, the implied volatilities for short maturities are more sensitive than those obtained by the first method, which greatly affects options prices within the money as their premiums are greater than options out of the money.

Source: own elaboration with data from BMV, MexDer and Banxico.

Figure 10 Surface of implied market volatility obtained by the parameters of function using differential evolution. 

Conclusions

A relevant feature of the differential evolution algorithm is that it does not require initial parameters since it generates many tests of a population in parallel. Differential evolution performs sampling over the entire parameter space, and if it is carried out enough, it practically guarantees convergence to a global minimum. However, the disadvantage is that it requires a considerable number of repetitions, involving a high cost in time and computational resources.

Another possibility for calibrating the parameters of a stochastic volatility model is by using loss functions. This approach minimizes the distance that measures the error between market prices and model prices, or between implicit market volatilities and those form the theoretical model. In this framework, the Heston model generates consistent market volatility smiles for relatively distant maturities (six months and more). However, for short maturities, it presents problems as it does not generate consistent volatility smiles with the market; this depends to a large extent on the dynamics of the underlying asset, which affects the market quotations at a given time.

The effect of the variation of the parameters of the Heston model on the implied volatility curve has been examined in this research, and it is observed that the correlation parameter ρ determines the direction of skew; ρ>0 corresponds to a positive slope and ρ<0 to a negative slope; while the volatility of the variance σ determines the curvature of the smile. The parameters κ, θ, and v0 control the level of the smile.

The contribution of this study can be summarized as follows, the calibration of parameters of the Heston model with loss functions and differential evolution algorithm, in particular with a loss function of the weighted MSE times squared vega which employs option price sensitivity in the Black and Scholes model in relation to the market-implied volatility. An advantage of this procedure is that it is able to produce precise estimates, but at the expense of increased computation time.

In summary, in this investigation, the parameters of the Heston stochastic volatility model with loss functions have been calibrated using two different methods: the first determines the minimum of a nonlinear function of several variables with constraints; the second employs the algorithm of differential evolution. The two methods are applied to a set of IPC call options with four maturities and twenty-eight strike prices. The criterion for choosing the parameters is the one that has the minimum value of the MSE of implied volatility. The results for the two proposed methods indicate that the function of loss of the MSE weighted by the square of the vega is the one that has the minimum value, but with a minor error and time cost using the first method. In addition, in both methods the adjustment of implied volatilities is better for options with long-term maturities than for short-term maturities, and the ρ value is negative, which means that the slope of the implied volatility curve is negative which reflects a skewed distribution of the underlying asset's return. This implies that the volatility modeled as a stochastic function of the underlying asset price is in line with assumption of stochastic volatility.