SWIFT calibration of the Heston model

In the present work, the European option pricing SWIFT method is extended for Heston model calibration. The computation of the option price gradient is simplified thanks to the knowledge of the characteristic function in closed form. The proposed calibration machinery appears to be extremely fast, in particular for a single expiry and multiples strikes, outperforming the state-of-the-art method we compare with. Further, the a priori knowledge of SWIFT parameters makes possible a reliable and practical implementation of the presented calibration method. A wide set of stress, speed and convergence numerical experiments is carried out, with deep in-the-money, at-the-money and deep out-of-the-money options for very short and very long maturities.


Introduction
The Heston model is a well-known stochastic volatility (SV) model for driving the dynamics of the assets. In order to use the Heston model, we need to calibrate its five parameters to real-market data. The goal of calibrating a model using market data is to estimate the model parameters in such a way that, when it is used for option valuation with an appropriate option valuation method, it yields prices similar to the real market ones. Calibration is an important task that requires efficient numerical methods. It encompasses a machinery for option pricing as well as an optimization procedure, aiming at minimizing the differences between the market option prices and the prices given by the valuation method. In regard to the pricing, we use the SWIFT method presented in [31] for European options. It belongs to the class of Fourier inversion methods and has already been used for pricing Bermudan, barrier and Asian options (see [27,24]). The power of SWIFT method partly relies on the knowledge of the characteristic function (ChF) in analytical form. Since the ChF associated to the Heston model is known in closed form, we can tackle the optimization problem by means of the gradient based Levenberg-Marquardt (LM) algorithm [29]. For sake of comparison, we consider the state-of-the-art calibraton method of [8], which is based on the Fourier inversion pricing method of [20] and it also uses the LM optimization algorithm. The main contributions of this work are the following.
• We extend the SWIFT method to the calibration problem by deriving the option price gradient.
• We implement and test the speeding up techniques mentioned in [31] based on multiple strikes valuation.
• We propose a novel method for calibrating the Heston model with a set of options with certain fixed strikes that can be used later on for arbitrary strikes by interpolation.
• We develop and implement speeding up techniques for the option price gradient.
We carry out a wide variety of stress, speed and convergence tests with at-the-money (ATM), deep in-the-money (ITM) and deep out-the-money (OTM) options, ranging from very short to very long maturities. The results show that SWIFT is extremely fast for sets of options with a single expiry and different strikes, being about ten times faster thant the calibration method of [8]. For options with a fixed number of maturities and a fixed number of strikes per maturity, both methods perform similar in terms of speed, beeing more robust the SWIFT method, thanks to the possibility of a priori selecting the parameters of the pricing machinery. This last feature, makes SWIFT method a reliable and practical methodology for real-time updating of Heston model calibration. The paper is organized as follows. We define the calibration and the valuation problem in Section 1. The Heston model, the valuation method of [8] and the calibration challenges are presented in Section 2. Section 3 is devoted to the SWIFT method and its speeding up features. In Section 4, we put forward the calibraton problem with all the mathematical details, and we test our proposed method through the numerical experiments of Section 5. Section 6 concludes and gives some pointers for future research.

Option valuation
The option valuation will be tackled under the framework of the expected discounted payoff pricing formula that we next recall. Consider a European option contract with strike price K, expiring at time T , with τ := T − t the time to maturity, and S t the price of its underlying asset at time t. Then, if one considers state variables x and y that fully describe S t and the random variable S T respectively, the general pricing formula becomes, where v(x, t) denotes the option value at time t, v(y, T ) is the payoff, r the risk-neutral interest rate, E Q the expectation under the risk-neutral measure, and f (y|x) is the probability density of y given x. A common choice for the state variables x and y is to define, x = ln S t K . ( More precisely, given the states variables choice, we can express the payoff of a European option in log-asset price as, v(y, T ) = [α · K (e y − 1)] + := K · g(y), with, α = 1, for a call, −1, for a put, (4) where g(y) := max (α (e y − 1) , 0) denotes the strike-free payoff.

Heston dynamics and calibration issues
The widely used Black and Scholes (BS) model fails to capture essential well-known properties of the real world market dynamics of the underlying return distributions, as its high kurtosis, its negative skew, the correlation between the underlying price and its volatility, the risk premium investors give deep (ITM) or (OTM) options, etc. All these properties result in the well-known BS implied volatility surface. SV models, treat both the underlying price and its volatility as (potentially correlated) stochastic processes, which allows to better capture some of these properties. One of the first and most well-known SV models is the Heston model [20], defined by the following system of stochastic differential equations.
Definition 1. The Heston model price-volatility dynamics is defined by, where W (1) t and W (2) t are two correlated Wiener processes dW (1) t dW (2) t = ρdt and ν t is the variance of the underlying asset price at time t. Then, if one specifies the initial value of the variance ν 0 , the model is properly defined. From now on, θ := (ν 0 , ν, ρ, κ, σ) T will refer to the vector of model parameters.
Several works have shown the relationship between the Heston model parameters and the shape of the implied volatility surface [7,15,17,21] necessary to obtain the same prices with a BS model. It can be summarized as follows, • ν 0 controls the position of the volatility surface, • ρ controls its skewness, • κ and σ control the convexity of the surface, • κ(ν 0 − ν) controls the term structure of implied volatility.
A method for calibrating the Heston model is presented in [8]. This method starts from the expression of the price of a European call option presented in the original work by Heston [20], and it is adapted here in terms of the state variables x and y defined in Section 1.
where P 1 and P 2 are defined as, andf (u) is the initial state independent ChF of the process. It is worth remarking that, as we will see in Section 3, the SWIFT method will benefit from the fact that the ChFf (u) does not depend on the initial state variable.
Remark 1. The dependence of the ChF in time and the model parameters is omitted for readability.
Remark 2. The expression (5) omits the dividend yield term q, which appears in [8] but is assumed to be 0 here for readability. The results presented here are valid for any constant value q.
Remark 3. Typically, the ChF of a random variable with density function f is defined asf (u) = R f (x)e iux dx. However, to be consistent with [31], it is defined in this work asf (u) We can see that there is a sign difference in all the u-dependent equations and expressions in [8].
We write expression (5) in a more compact form in the following lemma, Lemma 1 (Heston's pricing method). Let V (θ; x, τ ) be the price of a European call option with strike K and Heston dynamics, given by expression (5). If we definef (u; x) = e −iuxf (u), then, Proof. Having into account that S t /K = e log(St/K) = e x then, from expression (5), (6) and (7) we have, Finally, sincef (u; x) = e −iuxf (u), then the result follows.

Calibration challenges
As opposed to simpler 1-dimensional models, Heston model calibration is a multidimensional optimization problem with 5 degrees of freedom given by θ := (ν 0 , ν, ρ, κ, σ) T . Furthermore, the structure of this optimization problem is not known. According to [8] no consensus exists among researchers regarding whether the objective function of this optimization problem is convex or not. Some results point to a non-convex function, as the calibration methods proposed in [4,28] (which yielded different results for different initial points) and one must use long or short term approximations and rules to provide a convenient initial guess. Recent research [16] claims to provide methods that reach a unique solution independently of the initial point which, according to that study, indicates some structure that, even if not necessarily convex, tends to lead an initial guess to a stationary result. There is also no consensus on whether the problem always has a single optimum. In particular, it is known that there exist dependencies between the parameters that yield similar results. For example, lim t→∞ Var(ν t ) = σ 2 ν κ (where Var(·) refers to the variance operator), so large values of κ and σ can provide a model that prices options similarly to one with proportionally smaller values of these two parameters. The work by [8] claims that this results in the objective function of the optimization problem being flat close to the optimum. As said above, there is no guarantee that a gradient-based method converges to the global optimum of the model parameters, but even obtaining a local optimum has been traditionally difficult. Many papers in the literature uses numerical gradients (see [16]) for these methods when trying to solve the Heston calibration problem (which are less accurate and more computationally consuming), because no simple analytical gradients were available and the ones obtained with symbolic algebra packages from the expressions of the ChF were intractable. Prior to [8], the existing methods could be summarized as follows.
• Heuristic based models. Using the relationships outlined above, some works reduce the dimension of the optimization problem by assuming some values or relationships between the parameters from the observation of a specific volatility surface. For example [15] sets ν 0 to the short-term ATM implied variance obtained by using a BS model, a heuristic further justified by [4], where the linearity between ν 0 and the BS implied volatility was verified for short maturities (less than 2 months). Other heuristics used in the industry are κ = 2.75 τ and setting ν to the BS short term volatility [7]. These assumptions may restrict the optimization problem domain and exclude the optimum.
• Stochastic methods. They are typically used in combination with deterministic search methods, as the Nelder-Mead simplex method [23] and would avoid the pitfalls of the gradient-based methods if the optimization problem is not convex. Some examples are used in [4] and differential evolution and particle swar are used in [18]. These methods are too computationally expensive for real-time use as [12] that employs GPU computations to calibrate options using a SV model called SABR, and it took 421.72 seconds to calibrate 12 instruments with a tolerance of 10 −2 using 2 NVIDIA Geforce GTX470 GPUs.
In this work, we consider the analytical expression for the ChF provided in [8].

The characteristic function
For long-term maturities, [22] shows that the original ChF provided in [20] has discontinuities as u increases, which can lead to numerical problems. The source of these discontinuities was discussed in [1], and an alternative expression which was continuous in the full parameters space was presented in [33]. A more compact version of the ChF was later derived in [10] from the moment generating function of the process. This expression also has the benefit of having simpler analytical expressions of the gradient of the ChF than in previous expressions, but it also presents discontinuities as u increases. Finally, an expression with both, simple derivatives and continuity in the full parameters domain, was provided in [8], where, ξ := κ + σρiu, d := ξ 2 + σ 2 (u 2 − iu), Further, its gradient with respect to the Heston model parameters θ = (ν 0 , ν, σ, κ, ρ) T is given by, where the partial derivatives of A, A 2 , B, and d are given in [8] and can be seen in Appendix A.

European option valuation and calibration with SWIFT
In this section we give a brief overview on the SWIFT method, originally developed for pricing European options in [31]. In this work the method will be extended to European options calibration. For sake of completeness a section is devoted to the basic theory on Shannon wavelets.

Multi-resolution analysis and Shannon wavelets
Consider the space of square-integrable functions, denoted by L 2 (R), where, Then we can define a useful structure for function approximation called multi-resolution analysis (MRA). Let us start with a family of closed nested subspaces in L 2 (R), where, , If these conditions are met, then there exists a function φ ∈ V 0 , known as scaling function that generates a family of orthonormal bases of V m , denoted {φ m,k } k∈Z , These families allow to approximate any f ∈ L 2 (R) with increasing resolution by means of the projection map P m : dx, and · is the complex conjugate operator.
Increasing the considered number of elements of the finite family will increase the resolution of the approximation, converging to a perfect representation when all the functions of the original family are used (see [26]). Apart from increasing the resolution by means of m, wavelets can be moved by means of k and stretched or compressed by means of m to represent local properties of a function. A basic reference on wavelets is [9].
In this work, we employ Shannon wavelets [5], since they are regular and orthogonal functions with compact support in the Fourier domain. The regularity, as opposed to the Haar family used in [30], gives us much better accuracy. The rapid decay of the scaling function in the Fourier domain replicates the behaviour of the ChF of the heavy-tail processes that we encounter in finance. A set of Shannon scaling functions φ m,k in the subspace V m is defined as, where, is the scaling function, also called cardinal sine function.
The following lemma about the bound on the error of the orthogonal projection ǫ m (x) := f (x) − P m f (x) is provided in [27]. where, is the normailzed mass of the two-side tails off .

SWIFT method
The SWIFT method belongs to the class of Fourier inversion methods for pricing European options within the discounted expected payoff formula (1). The density function f of (1) is replaced by a finite combination of Shannon wavelets, and the coefficients of the approximation are computed from its ChF.
The overall process provides an efficient way to obtain the value of an option, and it can be summarized, as in [31,27], by a set of consecutive approximation steps, which are described below.
• Wavelet projection: the function f is replaced by its Shannon wavelet projection at scale m ∈ Z, with • Series truncation: the set of values of k involved in the sum of expression (14) is reduced to a finite interval [k 1 , k 2 ], It is important to notice that the first approximation lets us express, which quickly justifies that, for any given x, the density coefficients vanish as |k| increases, because the mass at the tails of f tends to 0. It is also worth noting that increasing m will result in this mapping being less favorable. That is, for each k, D m,k will be bounded by a point closer to the center of the density function, potentially requiring to increase the range of values for k in interval Remark 4. From this point onward a symmetric interval [1 − η, η] will be considered both for convenience and for consistency with the code implementation.
• Density coefficients approximation: the integral required to compute D m,k is replaced by an approximation D * m,k as it will be shown in Section 3.2.1, We can then define V m,k : For European options, one can instead express it in terms of the strike-free payoff by defining • Payoff coefficients approximation: the integral required to compute U m,k is approximated in an analogous way as the integral to compute the density coefficients, and U m,k is replaced by an approximation U * m,k as it will be shown in Section 3.
These coefficients can be precomputed when initializing the SWIFT procedure and shared across different strikes and maturities, saving computation time.

Density and payoff coefficients approximation
Density and payoff coefficients calculation rely on the approximation of sinc(x) by a finite combination of cosines (all the details can be found in [31] and [27] and the references therein). As a result of this approximation, a new parameter called J appears. This parameter will be labeled J d and J p to denote density and payoff coefficients, respectively. Following the notation of [27], the density coefficients are given by, where ℜ denotes de real part, u d j = π 2J d (2j − 1), and the payoff coefficients are given by, where, and, We can see that the value of I j (c) is periodic on c. In general, all the approximations to sinc(x) used in the SWIFT method are periodic, which can give rise to boundary issues and undervaluation of option prices when the option strikes approach to the boundary of (−c, c). This issue also appears in the COS option pricing method [11], another Fourier-transform-based option pricing method closely rleated to the SWIFT method, and is discussed in [27]. In that work, the authors use the independence between the parameter c regulating the payoff integral domain and the parameter η regulating the wavelet series truncation, to carefully choose a value for c to avoid this problem. An iterative method to choose appropriate values for m, η, J d , and J p is provided in [31,27].
As most operators used in the SWIFT method are linear, one can easily obtain an expression for the option price gradient that will be used for calibration. In particular, the only dependence the European option price formula has on the model parameters vector θ appears in the term D * m,k , so we can simply define, and replace it into expression (19) to obtain the expression of the gradient.

Speeding up the SWIFT method
We elaborate on several enhancements of the SWIFT method in terms of efficiency, either on the pricing or during the calibration process. In Section 3.3.1, it is explained how the density and payoff coefficients can be computed by means of the Fast Fourier Transform 1 (FFT). Section 3.3.2 is devoted to the adaptation of the SWIFT pricing machinery for multiple strikes valuation. Those two transformations where already pointed out in [31]. Moreover, when the vector of strikes meets a certain property, then the calibration can be carried out in a two stage procedure detailed in Section 3.3.3. Finally, Section 3.3.4 describes another improvement in regard to the option price gradient calculation.

Fast computation of the density and payoff coefficients
We start with a general expression of the summation term that appears in both the density and payoff coefficients approximation, We can extend it by defining g j = 0 for j = 0 and J < j < 2J and take the j-independent terms outside of the summation, obtaining, This last summation expression coincides with the Discrete Fourier Transform (DFT) of length 2J of {g j }, and the computation of all the values f k can then be speeded up by the FFT implementation.
Remark 5. Note that computing the density or payoff coefficients imposes a restriction on the wavelet series truncation parameter 2η < J.

Valuation with multiple strikes
A key property of the SWIFT method is that, given a scale of approximation m, the payoff and density coefficient associated to each wavelet φ m,k can be computed through two FFTs (one for all the density coefficients, and one for all payoff coefficients). Without this property, the SWIFT computation speed would not be competitive with other numerical option pricing methods [13].
In the option calibration problem, one usually needs to consider the option prices of several options at different strikes. In this specific case, if one were to compute the option prike of M options at strikes K := (K 1 , . . . , K M ) T , then the formulation proposed in expression (19) would need to recompute the density and payoffs coefficients for every strike K i . This involves evaluating the ChF η · J p · M times, an operation which, for the Heston model is more costly than evaluating the strike-free payoff function, or its integral. As stated in [31] one can improve the computation time of option pricing for multiple strikes whenf (u; x) =f (u)e −iux , a property present in both Lévy and Heston models.
As stated in [27], let us start from expression (19), and considering the previously mentioned vector of strikes K, with its associate vector of initial states x := (log(S 0 /K 1 ), . . . , log(S 0 /K M )) T . We can then substitute the density coefficients approximation (20) into the option price expression (19) and interchange the two resulting summations, obtaining, where,Ũ The original formulation from expression (19) requires the following computations, • For each of the M strikes: -1 FFT of length 2J d to compute 2η density coefficients.
-J d evaluations of the ChFf (u d j 2 m ; x).
• 1 FFT of length 2J p to compute 2η payoff coefficients.
• J p evaluations of the strike-free payoff integral I j (c) defined in expression (21).
The payoff computations are independent of the strike price and can be computed only once, and be reused for all strikes.
On the other hand, the alternative formulation provided in expression (29) requires, • For each of the M strikes: • J d evaluations of the x-independent ChFf (u d j 2 m ).
• 2 FFT of lengths 2η and 2J p to compute the J d values ofŨ j (−c, c).
• J p evaluations of the strike-free payoff integral I j (c) defined in expression (21).
The required values off (u d j 2 m ),Ũ j (−c, c), and the values of I j (c) required to compute the latter, can be computed only once and be reused for all strikes.
Remark 6. In general, whenever the dependency on x inf (u; x) can be easily isolated and is cheap to compute, one can benefit from the alternative formulation proposed in this section.
The computation of the ChF tends to be more expensive than the computation of the payoff integral, so a SWIFT implementation through expression (29) tends to outperform the one through expression (19) when several strike prices are involved. A discussion on the benefits of using a formulation equivalent to the one provided in (29) for multiple strikes appears in [31] where it is shown that is possible to define F j :=f (u d j 2 m ) and compute in advance its J d required values once and reuse them for all strikes.

Fixed set of strikes
Let us consider, Then expression (29) can be rearranged as, If one chooses selectively the values of the strikes K l , so that 2 m x l is an integer number, this computation can be speeded up by the use of an FFT algorithm. If one chooses x k := 2k − J d 2 m+1 , then expression (31) becomes, where,G j := G j e πij = G j (−1) j .
Remark 7. Note that, as with other FFT-based computations presented in this work, this approach imposes a bound M ≤ J d on the number of different strikes that can be computed with the FFT.
If one considers the domain D of x = log(S t /K), this approach allows pricing options in a symmetrical boundary (− J d 2 m+1 , J d 2 m+1 ) ∈ D at J d uniformly distributed points at distance 2 m . We can not usually choose the strike prices at which to price the options, particularly not when calibrating a model with real market data, as only a limited set of strike values are listed on any exchange market, but this method could be used to quickly compute the option prices of an already calibrated model at a grid of points that could be tuned by the choice of m and J d . Then the option prices at any intermediate strike could be interpolated with the help of a derivative-free spline (or, if the derivative with respect to K in expression (32) preserves the same speed properties, with the help of any spline method that uses derivatives).

Option price gradient
The option price gradient must be computed during the calibration process that will be presented in Section 4. All the aforementioned techniques can be applied to the option price gradient. We can also enumerate three more speed up properties, • The value of e −iu d j 2 m x l can be reused for the price as well as for the gradient computations.
• If the parameters of the SWIFT method, are not changed during the gradient descent used in the calibration problem, then the values of bothŨ j (−c, c) and e −iu d j 2 m x l can be reused throughout all the calibration steps.
• We can reuse the values off (u d j 2 m ) from the price computation to compute the gradient.
So combining all the speed properties above, when solving a gradient-based calibration problem, we only need to first computeŨ j (−c, c) and e −iu d j 2 m x l , and then in each gradient-descent step, one can simultaneously calculate both the price and the gradient of all strikes by computing once for each j ∈ [1, J d ] the values off (u d j 2 m ) and h(u d j 2 m ).

Calibration
Calibrating a model for the asset underlying the option, is a sophisticated procedure that requires highly efficient numerical methods. In particular, the pricing of the options used for calibration should be carried out by means of an accurate, fast and robust valuation method. In this work, we calibrate the Heston model by means of the SWIFT method, and compare it with the Heston's pricing method of [8] that we have written in a more compact form in Lemma 1. It is worth noting that the choice of a specific objective function will have an impact on how accurately the model of the underlying asset that we will obtain through callibration will describe different real market scenarios [6]. For comparison sake, the same one as in [8] will be used. Let V * (x i , τ i ) be the market price of a European call option and V (θ; x i , τ i ) be the price at the same strikes and maturities obtained by using either the SWIFT method or the Heston pricing formula in expression (8) of Lemma 1. Let us also assume that we use a set of n different options to calibrate the model, so that i ∈ [1, n] ⊂ Z. Then, the calibration of the model is defined as the minimization problem, where r(θ) is the n-dimensional vector of the residuals obtained when pricing the options considered for calibration using the model parameters. That is, If we calculate the Jacobian of r, gives us, where, The Hessian matrix of the residual element r i reads, where, Then, the gradient and Hessian matrix of the objective function defined by expressions (34) and (35) are, We solve the optimization problem (34) by means of the LM 2 method. This method is as a blend of gradient descent (GD) and Gauss-Newton (GN) iteration, depending on whether the current guess is close or far from the optimum. The exact expression of the step ∆θ is, where I is the identity matrix and JJ T + µI substitutes the Hessian matrix used in the Newton method. When the current guess is far from the optimum, a large value is given to µ so that, and a small step of a steepest-descent method is taken. When the current guess is close to the optimum, a small value is given to µ so that, and the Hessian usually used in the Newton method is replaced by its GN approximation. This approximation is reliable if either r i or H(r i ) are small, and [8] justifies its usage by conjecturing that f is nearly linear, a condition that guarantees the latter. We should note that even if f were not linear, then LM should only use small values of µ when |r| is small at the current step of the optimization problem. The iterative algorithm implemented for the LM method stops when, at a certain iteration n, any of the following criteria is fulfilled, The first stopping criteria (45) is fulfilled when the objective function defined by expressions (34) and (35) has reached a value closer to zero than the prescribed tolerance. It is only when the method stops due to this criteria that we will consider that the model has been properly calibrated. The second criteria (46) corresponds to a flat gradient, and the third (47) corresponds to a stagnating update (this last one has never happened while testing the convergence during the Heston model calibration).

Numerical results
In this word, the SWIFT 3 method is used to calibrate a Heston model with European call options price data at different strikes and maturities, and it will be compared to the pricing and calibration method based on expression (8) proposed by [8], which for the sake of readability will be called Cui pricer (CP). CP will be implemented using a Gauss-Legendre quadrature with 64 nodes for its numerical integration step. The upper limit of the integral will be truncated, whenever possible, at u = 200, but will be adjusted if necessary. The calibration process will consist of applying an LM method to the objective function defined in expression (34). The SWIFT method will be implemented using the ChF expression and its derivatives provided also in [8]. We perform a wide variety of tests that can be summarized as follows, • Stress tests: the CP and SWIFT methods will be tested with several combinations of extreme strikes (ATM and deep ITM and OTM) as well as with long-term and short-term maturities to detect any possible limitation or numerical issue in a wide usage range.
• Speed 4 tests: the option calibration speeds for the regular SWIFT method (defined by expression (19)) and the one devised to quickly compute several option prices with different strike and same maturity (defined by expression (29), which will be denoted KSWIFT), will be compared against CP for three different strike and expiry sets to check whether the multiple-strike alternative formulation is necessary to obtain a competitive option calibration method. These scenarios will represent, -A single expiry and multiple strikes.
-A fixed number of maturities and a fixed number of strikes per maturity.
-Different expiries for each strike.
When computing options with more than one different strike, a combination of OTM and ITM options will be used to provide an heterogeneous sample of contracts. Similarly, when more than one maturity is considered, a sample of long-and short-term expiries will be used.
• Realistic convergence tests: as in [8], convergence of the method will be tested for realistic model parameters representative of long-dated Foreign Exchange (FX), interest rate, and equity options, as they are relevant and, according to [19],challenging for simulations of the Heson model.
Several sets of Heston parameters will be used for the different numerical tests and are presented in Table 1. The last three sets of parameters are representative of long-term FX, interest rate, and equity options respectively [3].  Table 1: Set of Heston parameters used in the numerical tests.
Remark 8. θ (1) is obtained from [8] and θ (2) is a plausible set of parameters proposed by us. It may not be representative of any real world market, but it is only used as our objective value in our speed tests. In Section 5.2, we will use an initial guess of θ

Stress tests
Deep ITM and OTM call options are priced together with ATM call options for long-and short-term maturities using the SWIFT method at two different scales of approximation (m = 3 and m = 7) and the CP method. The time until maturity τ is given in years. Thus, the expiries of 0.04 attempt to simulate a situation of around two weeks until expiration of the option contract.
• CP and V 7 SW produced not a number (nan) results when evaluating very long expiries. Looking at the option price execution with the integrated debugger of the GDB compiler [34] showed that expression (9) runs into numerical overflow when the exponent dτ 2 of its hyperbolic functions is big enough (the same error happened when using the original expression provided in [10]). In most of the tests above, the overflow could be avoided when carefully setting an appropriate value for the upper bound u of integral in expression (8), and by using a smaller value of the scale m. The error can also be avoided by selecting the ChF expression provided in [33] (we use it later on, and we denote the obtained prices by V SH and present the results in Table 3).
• The SWIFT method at scale m = 3 tends to underprice short expiry options. After checking the SWIFT parameters obtained through the parameter choice method defined in Section 3, it was observed that the initial value for η, obtained by simply using the cumulant expression proposed in [31], resulted in a truncated Shannon wavelet expansion that did not cover a sufficient domain of the density function f (y|x). A dynamical choice of the parameter η based on the calculation of the area underneath the curve of the density function, as described in [31], can avoid this issue.
Increasing the value of m also fixes the problem.
• None of the methods can handle the deep OTM option with a short expiry. The expected value should be close to but bigger than 0, as there are only 10 trading days to expiry and the price of the underlying should increase 50% so that the option contract would not expire worthless. CP value seems too high and, in fact, when moving the value of u in the interval [100, 400] the price never clearly converges to a certain value, and it can give higher estimates for u > 200 than 1.079e − 3, or even negative values. Changing the ChF expression does not fix this issue. SWIFT consistently gives it a price of 0. The contribution that makes the price different than zero probably lies on the tails of the distribution function, and one would require a really big value of c so that a point with a positive payoff is even considered in expression (29). Table 3 shows some results either selecting an appropriate value u or keeping u equal to 200 and using the ChF from [33]. The column u indicates the value at which CP integral is truncated, and the option price obtained is shown in column V CP . Column V SH shows the price obtained when keeping u = 200 but implementing CP using the ChF provided in [33]. The last row is an example of the negative values obtained in the deep OTM short-term call.  (1) 100 200 0.04 300 -1.174e-5 1.079e-3 Table 3: Results for different u and/or using the ChF from [33].
As it has been shown so far, a crucial step of the calibracion process is the selection of u for the CP method and m for the SWIFT method. A method to set an optimal value u is not provided in [8], and it is therefore a matter of trial and error, since it must be manually determined when changing the time to expiry of the options one wants to price. As for the SWIFT method, we can use the iterative parameter choice provided in [27]. In particular, a suitable scale of approximation m can be selected by means of Lemma 2. Once the level of approximation m is fixed, the parameter η can be adaptively calculated in order to determine the wavelet series truncation accurately.

Speed tests
The calibration speed has been tested for three different sets of strike prices and maturities, which are available in Appendix B. In order to be sure that the calibration problem was properly converging, we chose θ (2) as the objective value for the Heston parameters. When testing each set we perform the following actions, • Generate option price values for each strike-expiry pair using θ (2) as input. For this step, we used SWIFT method (we also generated them using CP method and checked that the difference between both results stayed under 10 −7 ).
• Solve the calibration problem with the desired method using θ (2) 0 as initial guess and use as inputs the strike-expiry pairs and the prices obtained in the first step. Set 2 has been obtained from the code provided in [8], and represents a total of 40 options, distributed in 8 different maturities, each maturity with 5 different strikes. Set 1 and set 3 are extreme cases derived from set 2, in order to test the best and worst calibration points distribution for KSWIFT. Set 1 has the same 40 different strike prices than set 2, but only a single maturity. What we denote as set 3 is not really a different data set, since it consists of running the calibration problem with set 2 and preventing the KSWIFT algorithm from applying the speed up techniques discussed in Section 3.3.2. For this reason, in Table 4   The values for KSWIFT and CP have been averaged over 100 executions of the calibration to provide a good estimate of the required calibration time. It can be seen that regular SWIFT is orders of magnitude slower without averaging the required time over several executions. Hence, the multiple-strike alternative formulation presented in Section 3.3.2 and all the speed-up techniques discussed through Section 3 are necessary to provide a competitive method that can be used for real-time model updating.
KSWIFT performance is comparable to CP for set 2, an order of magnitude faster for set 1, and an order of magnitude slower for set 3. It can be argued that both set 1 and set 3 are extreme cases that are not really relevant for real option trading situations. One would rarely use a single strike per expiry to calibrate an option pricing model, and using data from a single expiry only seems reasonable when trading a single option expiry (in this case, one could benefit from the speed properties of KSWIFT on scenarios like set 1). According to [8], a reasonable calibration scenario consists of using option prices from strikes at 0%, ±25%, and ±50% BS delta (derivative of the option price with respect to the underlying price value. It has a closed analytical expression for European BS options).
The calibration time of KSWIFT and CP is about 0.05 seconds, which seems sufficient for real-time model updating to provide market information to a human trader. In a more computationally demanding trading environments, like high-frequency trading neither KSWIFT nor CP would be competitive enough.
Remark 9. All the single expiry tests (the first three tests on Table 4) converged to an approximated value different than θ (2) but approximated all the option prices properly. Using different initial guesses lead to different approximated values which minimized the objective function. It would be interesting to see whether this is a property of the Heston distribution (that is, it has at least a degree of freedom when defined from option prices in a single expiry) or it is due to the specific scenario being tested.

Realistic convergence tests
We use the same procedure as in previous section to solve several different calibration problems. For each case, we use one of the proposed realistic parameter sets as objective value and generate option prices for each strike-expiry pair. Then, for each objective parameter set, 100 different initial guesses are generated. Each component of the initial parameters guess is drawn uniformly and randomly within ±10% distance of the optimal value. According to [8] this is representative of real option calibration as, usually, the initial guess used for a certain calibration problem is the last available parameter estimation. If the calibration is updated fast enough, it is expected that the initial guess will be this close to the optimum. The maturities used in [8] are not available, so for these tests the strike-expiry set 2 will be used. As can be seen in Table 5, even under challenging parameters setups representative of real option  trading, KSWIFT is able to provide accurate estimations of the Heston parameters, taking on average a computation time of hundreds of milliseconds. These results, both in terms of speed as in terms of accuracy, are comparable to the tests in [8], so it is concluded that KSWIFT is as efficient as CP for real market scenarios. Further, if we take into account the robustness of KSWIFT in terms of the a priori knowledge of its parameters, as stated in Section 3.2, we conclude that KSWIFT is a very competitive method for calibration.

Conclusions and future research
We have investigated the problem of calibrating the Heston model, which belongs to the class of stochastic volatility models. An extension of the SWIFT method has been provided in this work for European options calibration, along with novel speed-up techniques, which can radically improve the performance when several of the priced and calibrated options have the same time to maturity. Some numerical issues arise with the ChF for very long-term expiries. Following the a priori knowledge of parameters selection for the SWIFT method seems to be enough to avoid these problems, while the parameters of CP need to be adjusted manually. The proposed speed-up techniques are deemed necessary in order to make SWIFT a competitive calibration method, as it has been seen in the numerical speed tests. In particular, it has been shown that the only situation where the proposed calibration is significantly slower than CP is when one calibrates the model with many different maturities with no more than one or two strikes per maturity. As the number of strikes per expiry increases, the relative speed of the SWIFT method increases, and it is about ten times faster than CP when calibrating 40 options with a single maturity. Both extreme situations are not representative of real option trading needs, and for a reasonably real situation of 5 strikes per expiry, the SWIFT technique is slightly faster than CP. A SWIFT implementation without the previously discussed speed-up techniques has also been tested and deemed non-competitive, with calibration times that reached the dozens of seconds. Further, the proposed calibration strategy passes the realistic calibration tests for challenging Heston model parameters setups presented.
In summary, the proposed SWIFT method is a robust and efficient machinery for real-time updating of option models used in human-supervised trading schemes. Neither SWIFT nor CP method are suitable for the most demanding algorithmic trading situations, like high-frequency trading. Several future work may encompass the following topics, • Most of the calibration tests with a single expiry have run into an optimal value different than the original one. It is to be seen if this is a property of the Heston model or if this was due instead to the specific parameter or strike/maturity values being used.
• It would be interesting to study the properties of the SWIFT implementation proposed for a chosen set of strikes in expression (32). We could interpolate the values at all strikes with splines methods that require derivatives, and not only derivative-free ones.
• Options with very long maturities may hamper the calibration process due to numerical overflows during the pricing step. The problem of long maturities has been tackled with Haar wavelets in [30]. It might be worth investigating whether we can do the same with Shannon wavelets.
• Deep OTM options with very short maturity are challenging to price. The problem seems to be the lack of accuracy of the approximation on the tails of the density function.
• Comparison with other calibration methods based on approximation formulae, like for instance the work by [2].
is only be applied to KSWIFT, and consists of the same strikes and maturities as set 2. The code implementation of KSWIFT generated for this work receives as inputs a vector of expiry-defined-data (EDD). Each element of the vector of EDD contains a single expiry and a vector of strikes. Thus, set 2 will have all the strikes with the same expiry grouped in a single EDD, and set 3 will have EDD consisting on a single strike. This will enforce full recomputation of the density and payoff coefficients for each strike.   Table 7: Set 2 of strikes and expiries.