{'fraction_non_alphanumeric': 0.036807009941985026, 'fraction_numerical': 0.012541826147661347, 'mean_word_length': 4.367876986690787, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Fast computation of vanilla prices in time-changed models and implied volatilities using rational approximations

We present a new numerical method to price vanilla options quickly in time-changed Brownian motion models. The method is based on rational function approximations of the Black-Scholes formula. Detailed numerical results are given for a number of widely used models. In particular, we use the variance-gamma model, the CGMY model and the Heston model without correlation to illustrate our results. Comparison to the standard fast Fourier transform method with respect to accuracy and speed appears to favour the newly developed method in the cases considered. We present error estimates for the option prices. Additionally, we use this method to derive a procedure to compute, for a given set of arbitrage-free European call option prices, the corresponding Black-Scholes implied volatility surface. To achieve this, rational function approximations of the inverse of the Black-Scholes formula are used. We are thus able to work out implied volatilities more efficiently than one can by the use of other common methods. Error estimates are presented for a wide range of parameters.


Introduction
The Black-Scholes model introduced in 1973 is arguably the best known result in quantitative finance. Four decades after its introduction, the model continues to be widely used, especially as a universal benchmark model, in part due to its tractability. In a frictionless market in which the asset price is modelled as a geometric Brownian motion (GBM) with constant drift and constant volatility, the price of a European put or call option has a closed form: the celebrated Black-Scholes formula.
The Black-Scholes formula is employed by traders to convert prices into units of implied volatility.
However, the presence of an implied volatility smile in option markets contradicts the assumptions of the Black-Scholes model and demonstrates that the returns are asymmetric and leptokurtic. At shorter maturities the volatility smile becomes more pronounced, showing an increasing deviation from the GBM model. It is therefore well established that the GBM model does not provide a satisfactory description of real market price dynamics. Real market asset prices typically exhibit jumps and volatility clustering, while price paths in the Black-Scholes model are continuous and have constant volatility. These observations are of course well known: we refer the reader to Cont & Tankov [13] and Gatheral [22] for relevant background and further references.
To address the shortcomings of the GBM model, a number of model classes with jumps or stochastic volatility in the asset dynamics have been introduced. These models can be broadly divided into two categories: (a) jump-diffusion or finite activity models; and (b) infinite activity Date: 5th May 2014 models. In the first category, the evolution of prices is modelled by a jump-diffusion process, with a finite number of jumps occurring at random times before a finite time-horizon. The jumps can be interpreted as rare events, crashes or large drawdowns. For models in this category, the jump part is typically given by a compound Poisson process. Examples of such models are the Merton jumpdiffusion model [39], the Kou model [28] and the Bates model with stochastic volatility [4]. The second category comprises models with infinitely many jumps in every time interval. It has been argued that such models give a more realistic description of the price process at various time scales (see Madan [34], Carr et al. [9] and Geman [23]). These models are capable of generating non-trivial small-time behaviour and, as shown in [9], the Brownian component is no longer needed. Many of these models are time-changed models that can be written as Brownian motions time-changed by an increasing stochastic process. The idea of stochastic time-changes was introduced into finance by Clark [12], who modelled the observed price process as Brownian motion run on an independent second process called the clock.
It is important for trading and risk management to be able (a) to price vanilla options in models that provide an accurate description of market prices, and (b) to generate the corresponding arbitrage-free volatility surface quickly and efficiently. Fast calculation of vanilla option prices in the specified model is also essential for carrying out the calibration of the prices to market quotes in a timely fashion. A financial institution active in the derivatives market will seek to evaluate large portfolios of vanilla options at speeds close to real time, in order to control its risk positions. This paper is devoted the development of a new method for the pricing of vanilla options and the implied volatility surface in the class of time-changed models, using rational approximations.
Such approximations are known to offer an efficient and accurate method for computation of the cumulative normal distribution function (see Abramowitz & Stegun [1] for details), and could therefore be expected to perform well for approximations of the Black-Scholes formula, which is given in terms of the cumulative normal distribution function. Rational function approximations generally outperform polynomial approximations in terms of computational efficiency, that is, for a given computational effort smaller (global) errors can be achieved with rational function approximations (see Morris [40]). For some functions, the optimal rational function approximation is known further to be able to achieve substantially higher accuracy than the optimal polynomial approximation with the same number of coefficients. Rational function approximations have been extensively used in theoretical physics, engineering, statistics and economics. For background on rational approximations and their use in economics, refer to Ralston & Rabinowitz [43] (Chapter 7).
The approach proposed here is based on rational function approximation theory: drawing on this theory, an approximation for the value of a vanilla option in a time-changed model is constructed, taking the form of a linear combination of a number of negative exponential moments of the clock.
This yields an explicit approximation for the value of a vanilla option in those time-changed models for which the Laplace transform (and hence any negative exponential moment) of the clock is available in tractable form. For many of the popular time-changed models the Laplace transform of the clock is known in closed-form.
To provide a numerical illustration, the method was implemented for the variance-gamma (VG) model, the CGMY model and the zero-correlation Heston model (the version of this model does not fall into the class of models studied here, but could be analysed by a two-dimensional extension, which is left for future research). The rational function approximation method was found to yield stable, fast and accurate results. Comparison of the method to the standard fast Fourier transform (FFT) method (see e.g. Carr & Madan [11]) with respect to speed and accuracy, appears to be favourable for the rational function approximation method in the cases considered. At this point it is worth mentioning that several alternative and refined Fourier methods have recently been developed for the pricing of European options, such as, for instance, the COS method of Fang & Oosterlee [19] -see Boyarchenko & Levendorskii [6] for an overview. A detailed comparative study of all the different pricing methods is left for future research.
The rational function approximation method is then applied to compute Black-Scholes implied volatilities. The method uses the fact that the computation of the implied volatilities that correspond to a given set of arbitrage-free call option prices can be reduced to the computation of a single function of two variables. Contrary to Li [32], who proposes to approximate the entire resulting surface by a two-dimensional approximation, we perform a number of one-dimensional approximations and interpolate between these for any given combination of input parameters. For a wide range of parameters, the maximal approximation error of this approach is bounded by 5.5 × 10 −5 , refining to our required level of accuracy the approximation of Li [32], which achieves a maximal error of 3.3 × 10 −3 over the same region. It should be noted that this procedure can be used independently of the results derived for option pricing and is not limited to time-changed Brownian motion models. We refer the reader who is primarily interested in the implied volatility computation to Section 6.
The remainder of the paper is organized as follows. In Section 2 the theory for time-changed models is briefly reviewed. In Section 3 the price of an option in a time-changed Brownian motion model is expressed in terms of a normalised Black-Scholes formula, and the approximation to the option price is presented in Section 4. Section 5 is devoted to a study of the error of the method, the application to computation of the implied volatility is given in Section 6, numerical results are reported in Section 7, and concluding comments are given in Section 8.

Time-Changed Brownian motion models for the asset dynamics
We assume frictionless markets and no arbitrage, and take as given an equivalent martingale measure (EMM) Q chosen by the market. All stochastic processes defined in the following are assumed to live on the complete filtered probability space (Ω, F, {F t } t≥0 , Q).
We concentrate on the equity market and consider models in which the stock price process (S t ) t≥0 under Q can be represented as follows: where (W t ) t≥0 is a Brownian motion and (Z t ) t≥0 is an independent stochastic process called the clock.
Throughout the paper, r ≥ 0 and q ≥ 0 denote the constant risk-free interest rate and the constant dividend yield respectively, S 0 denotes the known stock price at time zero and S t denotes the random stock price at time t. The condition that (S t e −(r−q)t ) t≥0 is a martingale will be guaranteed by an appropriate choice of the mean-correcting compensator ω as follows: where E[e Xt ] is assumed to be finite for all t ≥ 0. As can be seen, (X t ) t≥0 arises by subordinating a Brownian motion, with drift θ and volatility σ, by the subordinator (Z t ) t≥0 . Standard references on Lévy subordination include Bertoin [5] (Chapter 3) and Sato [44]. See Cont & Tankov [13] and Geman et al. [24] for background on the role of subordination in financial applications.
The clock (Z t ) t≥0 is required to be an increasing process for which the Laplace transform at time t should be available in tractable form. Typically, (Z t ) t≥0 is modelled by either a Lévy subordinator, as discussed in Section 2.1, or as a time integral of a positive diffusion, as discussed in Section 2.2.
Note that suitable models for the clock can also be constructed by combining these two ingredients and therefore include time-changes of the form can be (i) a Lévy subordinator, (ii) an integral of a non-negative process, or (iii) any other increasing process. It is worth noting that the three stochastic volatility Lévy models discussed by Carr et al. [10] are of this form. This follows since the models considered in [10] are three Lévy A Lévy process (X t ) t≥0 has marginal distributions that are infinitely divisible. Well-known examples of infinitely divisible laws are the Gaussian and the gamma distribution, which will both be utilized in the remainder of this paper. Given a Lévy process (X t ) t≥0 , define the corresponding characteristic function as follows: where φ : R → R is called the characteristic exponent. The characteristic exponent can be shown to exist and shown to be a continuous function. By the Lévy-Khintchine representation, it is also known that the characteristic exponent can be written as with characteristic triplet (A, ν, γ).
A Lévy subordinator (Z t ) t≥0 is a Lévy process that takes values in R + with characteristic triplet (0, ρ, β) satisfying ρ((−∞, 0]) = 0, ∞ 0 (x∧1)ρ(dx) < ∞, and β ≥ 0; that is, (Z t ) t≥0 has no diffusion component, only non-negative jumps of finite variation, and a non-negative drift. As a consequence, the trajectories of Z are almost surely increasing. Since Z t is a positive random variable for any positive t, it is natural to describe its distribution using the Laplace transform: where ψ(u) is called the Laplace exponent of Z. Since (Z t ) t≥0 is an increasing process, it can be interpreted as a "time deformation" and be used as a stochastic clock. Subordinating any Lévy process (in particular, a Brownian motion) by another independent Lévy process will yield a new Lévy process (see Theorem 30.1 in Sato [44]).
We next present some explicit models in this class that we will employ later.
2.1.1. The variance-gamma model. The VG process was first introduced to finance by Madan & Seneta [36] in its symmetric version, and later extended to its asymmetric version by Madan, Carr & Chang [35]. In the asymmetric version discussed here, the variance gamma process has three parameters: θ ∈ R, σ > 0, ν > 0. It is defined by evaluating a Brownian motion with drift θ and volatility σ, as given in Equation (2.2), at an independent gamma time. Specifically, the timechange (Z t ) t≥0 is now given by a gamma process independent of (W t ) t≥0 with marginal distribution at time t given by a gamma distribution G( t ν , ν) with shape parameter t ν and scale parameter ν. The probability density function of Z t , conditional on Z 0 = 0, is given by  [9] and can be seen as a generalization of the VG model discussed above. In particular, the class of tempered stable processes that can be represented as time-changed Brownian motion coincides with this model. The characteristic function of X t in this model can be shown to equal The Laplace transform of Z t , the value of the corresponding clock (Z t ) t≥0 at time t, is given by (see Madan & Yor [37] for details) For this model we set the drift equal to θ = G−M 2 and the volatility to σ = 1.
where (V t ) t≥0 is a mean-reverting non-negative Markov process. The mean-reversion is required to guarantee that the random time-change persists.
2.2.1. Heston model without correlation. The Heston model [26] without correlation is obtained by defining (V t ) t≥0 in Equation (2.11) to be equal to the positive CIR rate process introduced by Cox, Ingersoll & Ross [15]. The stochastic differential equation for this square root process is given by where (B t ) t≥0 is a standard Brownian motion that is independent of Brownian motion (W t ) t≥0 in Equation (2.2). Additionally, the CIR activity rate process starts at V 0 > 0, κ > 0 is the rate of mean reversion, δ > 0 is the long-run activity rate level and ξ > 0 is the active rate volatility. It is further assumed that the Feller condition (2κδ ≥ ξ 2 ) is satisfied to ensure that the process never hits zero.
The characteristic function of Z t is well known from the work of Cox, Ingersoll & Ross [15], and is closely associated with Lévy's stochastic area formula. We recall that the Laplace transform of (Z t ) in this setting is given by where η = 2ξ 2 u + κ 2 . In the setting of (2.1)-(2.2), the Heston model without correlation corresponds to the parameter values σ = 1 and θ = − σ 2 2 = − 1 2 , so that X t = − 1 2 Z t + W Zt and ω = 0.

2.2.2.
Other examples. By taking (V t ) t≥0 to be an affine jump-diffusion process for example, the so-called quadratic models (see Leippold & Wu [30]) arise.

Vanilla Options in Time-Changed Models
A European call option on a stock with maturity date T and strike price K is a contingent claim that gives its holder the right, but not the obligation, to buy the stock at date T for a fixed price K.
The arbitrage-free value of a call option at time zero can be expressed as the discounted conditional expectation of the payoff: (3.1) where the price of the stock S T at time T is given in Equation (2.1) and r ≥ 0 denotes the constant risk-free interest rate as before. In the Black-Scholes model, the risk-neutral dynamics of S are described by the exponential of a Brownian motion with drift, where σ ≥ 0 is the volatility of the asset. The call price (3.1) in this model is given by the Black-Scholes formula: where N represents the standard normal cumulative distribution function. The following parameterisation of the Black-Scholes formula will be considered in what follows: This normalisation is similar to the non-dimensional form of the Black-Scholes formula proposed by Lipton [33] (Chapter 9.3). The function c BS is related to the original Black-Scholes formula as follows: where (W t ) and (Z t ) are independent, the vanilla call option price with underlying S can be expressed as The parameters S 0 , K, T, r, q, θ and σ are all constant and defined above. From Equation (2.3) it follows that the mean-correction ω is available in closed-form: so that ω is specified as follows: is the Laplace exponent of (Z t ) t≥0 . In particular, it should be noted that ω only changes with θ and σ and stays constant once these market parameters have been specified.
In the Heston model with zero correlation, ω is equal to zero, which follows from the fact that e Xt is a true martingale in the Heston model.
Remark. Note that µ T C only depends on the market parameters θ and σ 2 while x T C depends on r, q, S 0 , T , K, θ, σ and ω. In particular, x T C is contract-dependent as it varies with maturity T and strike K. We will refer to µ T C and x T C as the market-dependent parameter and the adjusted log-moneyness, respectively.
Proof. By conditioning on Z T and given that (Z t ) t≥0 and (W t ) t≥0 are independent, we find the following: where we used the definitions (3.3) and (3.5) of C BS and c BS . 2 In the case that (e Xt ) t≥0 is a martingale (or equivalently, θ = − σ 2 2 ), the normalised Black-Scholes formula c BS admits a symmetry: Corollary 1. If µ = 0, then the following holds true for all v ∈ R + and x ∈ R: This identity thus applies, for example, to the Heston model with zero correlation.
Proof. By straightforward algebra we get the following equalities: The discussion so far has only considered call options. However, with the help of the put-call parity, one can easily obtain the value of a put option if the value of the corresponding call option is known. For convenience, we state the put-call parity for the time-changed setting considered here: On account of Proposition 1, the problem of pricing a call option in a time-changed model reduces to the problem of evaluating E[c BS (σ 2 Z T ; x, µ)] for known values of x and µ, which we address in the following section.

Rational Approximations
To evaluate C T C given in Equation To obtain the parameters a x,µ j and b x,µ j , for given x and µ, we use a rational Chebyshev approximation (see Appendix 9.2 for a summary).
Next, using a partial fraction decomposition, we can rewrite the right-hand side of Equation (4.1) as follows: The parameters A x,µ j and B x,µ j will usually be complex numbers. If Re(B j ) < 0, which will typically be the case (see also the remark below), we may rewrite any of the denominators (v − B x,µ j ) −1 in integral form as Interchanging sum and integration yields the following: To approximate the integral in Equation (4.3) efficiently, we truncate the integral and use a Gaussian quadrature: where 0 ≤ c < d < ∞, w k are the quadrature weights, x k the abscissas and f denotes the integrand in (4.3). We compared a number of different Gaussian quadrature methods and found that Gauss-Legendre performs particularly well; it results in small errors for all parameter values and models considered here. Abscissas and weights are calculated by the standard procedure (see Press et al. [41] for details) and do not depend on the integrand f in the case of Gauss-Legendre quadrature. We thus obtain the following approximation: To summarize, given constant values for x and µ, E Q [c BS (σ 2 Z T ; x, µ)] can be approximated by a rational function of the form where A x,µ j and B x,µ j are the parameters of a rational function approximation with degree m (after applying a partial fraction decomposition), all depending on x and µ. Indeed, if Re(B x,µ j ) < 0 for all j, it follows from Equation (4.1) and Equation (4.5) that the following approximation is valid: Substituting v = σ 2 Z T and taking expectations on both sides yields the following expression: which is equal to the right-hand side of (4.6), by the linearity of the expectation.
Combining Equation (4.6) with Proposition 1 leads to an explicit approximation of the call option prices in the time-changed model in terms of the Laplace transform of the stochastic clock Z T .
Remark. Typically, the values B x,µ j in the rational approximation have negative real part, but sometimes a coefficient B x,µ j has non-negative real part. If this is the case, the easiest way to proceed is to increase the order of the rational function approximation. In the cases where switching to a higher order approximation does not resolve the problem, one can deal with the approximation of terms for which Re(B x,µ j ) ≥ 0 separately. For those terms, we define a new variable v * = e −v and approximate the resulting function by a Chebyshev approximation of order n to get (for a review of Chebyshev approximations see Appendix 9.1). Note that if more terms with Re(B x,µ j ) ≥ 0 occur, only a single Chebyshev approximation is needed for all these terms simultaneously. Thus, in Equation (4.6) the terms in the sum corresponding to B x,µ j with non-negative real part should be replaced by an approximation of the form given in Equation 4.7.
Remark. In the case of the Heston model with correlation, a similar conditioning approach would lead to a family of functions of two variables that needs to be approximated. This follows since log (S t /S 0 ) is normally distributed conditional on V t and t 0 V s ds (a fact which can be found in Glasserman & Kim [25] or Broadie & Kaya [8]). For given values µ and x, one would therefore approximate a two-dimensional function, which can then be rewritten in terms of negative exponential moments of the clock. We did so utilizing the downhill simplex search method and results were reasonable. However, a detailed analysis of this case is left for future research, where one should investigate faster and more powerful two-dimensional approximation methods. We assume that S 0 = 1, K = 1.1, r = 0.03, q = 0.01 and T = 1, so that we are pricing an out of the money call option with one-year maturity. It then follows from Equation (3.7) that µ T C and x T C take the following rounded values: A rational approximation of c BS (v; x, µ) with degree m = 5 is then given as follows: for v ∈ [0.0027, 0.0405] after applying the partial fraction decomposition (it should be noted that we rounded these values only to present them here, but use more accurate values to derive the results stated below). Then, if we denote each ratio by A j v−B j , using the gamma distribution, the above can be rewritten as where all five B j s have negative real parts.
In Figure 1  for values in between. As a consequence, one is able to speed up the method dramatically, since any option price with these market parameters can now be calculated only by using interpolation.
Although the rational Chebyshev approximations can be performed very quickly, it should only be done once for a given set of market parameters.
Obviously, when approximating only a few values of x offline, which are then used for the interpolation, one should make sure that the errors of these approximations are sufficiently small. As illustrated in the next section, for some of the x values, errors are much bigger than for others, due to numerical instability. However, each approximation can be carefully checked by evaluating the resulting errors (see Section 5.2), as the computational time for these approximations should not matter as much (given that it is done only once offline). One can therefore search for a rational function approximation that performs well for the given value of x.
Remark. When using this method for fast calibration, one could pre-calculate approximations of c BS (v; x, µ) for a number of parameter combinations of x and µ. The online calculation of an option price would then involve a two-dimensional interpolation between the values of x and µ calculated offline. This should be particularly feasible for calibrations where parameters are not allowed to vary too much from parameters obtained in the previous calibration, as one then has a well-defined range for µ. Investigation of this idea is left for future research.

Numerical Error Estimates
Observe that the error from our method arises from the approximation in Equation  Table 1. Values of the market-dependent parameter (µ T C ), and the maximum and minimum adjusted logmoneyness (x max T C and x min T C ) for five sets of parameters found in the literature (* we have set ρ equal to zero, as our framework only considers the zero-correlation Heston model, even though it was not zero in the original papers).
One can see from the table that the values for µ T C , x max T C and x min T C vary substantially between different models and parameters. Note that µ T C and ω are both zero in the case of the zero-correlation Heston model, implying that µ T C and x T C do not depend on the model parameters at all. From here on we refer to these sets of parameters as Case I-V, respectively. the truncation values of v, we analyse the cdf of Z T for Cases I -IV of Table 1. As stated in Section 2.1.1, Z T is gamma-distributed (Z T ∼ G( T ν , ν)) for Cases I and II. By the scaling property of the gamma distribution, we obtain that σ 2 Z T ∼ G( T ν , σ 2 ν), with probability density function .
where χ 2 γ (λ) denotes a noncentral chi-square random variable with γ degrees of freedom and noncentrality parameter λ. Then, conditional on V 0 and V T , the distribution of Z T = T 0 V s ds is given by Theorem 2.1 in Glasserman & Kim [25] and can therefore be simulated easily. We implemented this procedure and estimated the cdf numerically for Cases III and IV. Results are represented in equal to the 0.1% and 99.9% quantiles of σ 2 Z T , that is, we solve for x setting p equal to 0.001 and 0.999 in the following equation: For the gamma distribution, for example, this inverse is then given by Clearly, the more that is known about the possible range of σ 2 Z T and the smaller we are able to make the range [a, b], the easier it will be to find an appropriate rational approximation for this range. There are also extreme parameter combinations for which the range of [a, b] is rather large and the method does not perform as well (that is if the time change (Z t ) is too wide). Also, for maturities smaller than T = 0.25 we found that the resulting a might be too small for some of the parameter combinations and models, and the rational approximation may result in larger maximum errors.  Table 1 separately. The errors can be calculated as follows:

Error
Calculating these errors for each of the 246 options of Cases I -IV separately results in Table 3, where we represent the maximum absolute error for each case.
As can be seen, maximum logarithmic errors are of order −5 to −6. Even though these errors are already relatively small, we observed that most of the 246 errors are much smaller than these   should yield a great improvement. We will detail this procedure further in the next section and give error estimates in Section 7.

Implied Volatility
The Black-Scholes formula is often used to express option quotes in terms of the implied volatility: Given the option price and all parameters except σ in the Black-Scholes formula (3.3), one searches for the value of σ that solves the equation. This value of σ is referred to as implied volatility and we denote it by σ IV . Since there is no closed-form formula for σ IV , it is usually solved for by using a solver method, which is typically rather slow. and especially Jäckel [27], are faster and more accurate. Adopting the method described in the current paper, we construct an approximation of the inverse function by a rational function, which results in a fast and accurate alternative to the standard solver methods.
We restate the normalisation of the Black-Scholes formula (3.4) with a slightly modified parameterisation: Any call option price can then be written in terms of this normalised formula as follows: Since the normalised Black-Scholes formula c IV is a function of only two variables, finding a suitable rational approximation for any parameter combination becomes easier. Note that the question of finding σ IV reduces to finding v IV in Equation (6.1), where it follows from (6.2) that which are both given constants when searching for the implied volatility. We therefore need to solve for the value of v that solves Equation (6.1) for known values c IV and x. Once we obtain T . It should also be noted that one can restrict, without loss of generality, to the case of call options, as the case of put options follows on account of the put-call parity as in Section 3. In addition, reasoning as in the proof of Corollary 1, the following relation can be seen to hold true: so that the following relation holds for the implied volatility: As a consequence, one can concentrate on the negative half-line x ∈ R − . All results presented in this section so far can be found in a recent article by Li [32], who concentrates on approximating the two-dimensional Black-Scholes implied volatility formula v IV by a single two-dimensional rational approximation. Li approximates the inverse of Equation (6.1) for −0.5 ≤ x ≤ 0 and c LB (x) ≤ This means that bounds are widest for x = 0 and become narrower as x decreases to −0.5. These bounds are well chosen by Li, as it becomes more and more difficult to approximate the function well for small values of x. Still, it should be noted that the range is rather wide.
Note that the inverse is a function of two variables (c IV ,x). Li therefore uses a two-parameter (6.5) and obtains the approximation parameters using the downhill simplex search method. Without further adjustments (like a Newton-Raphson polishing), Li states a maximum error of 3.3×10 −3 for v IV over this region. It should be noted that the error is of similar magnitude for many parameter combinations and that one still needs to divide by √ T to obtain σ IV (making the error even bigger for short maturities). It should also be stressed that the approximation of Li applies to all combinations of c IV and x (within the bounds) and a total of 31 parameters is sufficient for any calculation. Having worked with the downhill simplex search method ourselves, we note that this is a remarkable result and that it should be used wherever this error bound is acceptable.
Rather than the ambitious attempt to approximate the whole function with one set of parameters, where p varies between 7 and 9 for different values of x 1 . The approximation parameters n x i and m x i , where these all depend on x, can be obtained again using the rational Chebyshev approximation (discussed in Appendix 9.2). However, we found that here it is advantageous to use the iterative algorithm discussed in Section 5.13 of Press et al. [41], using a maximum of six iterations to get closer to the so-called ideal "minimax" solution. Once the approximation has been performed for each value of x, approximation parameters are stored to the computer and do not need to be recalculated, as detailed in Section 4.2. Note that these approximations neither depend on the model nor on the market parameters and can therefore be used for any option. To calculate v IV for any value −0.5 ≤ x ≤ 0, one only needs to interpolate between the previously calculated values (using linear interpolation or a cubic spline, for example). We observe that the maximum absolute error for −0.5 ≤ x ≤ −0.0075 is 8.55 × 10 −7 , whereas the maximum absolute error for We would like to stress that this method, like the method presented by Li, does not require any approximation when calculating the implied volatility online. All approximations are done once and are stored to the hard-drive. The parameters one needs to store are all real and the calculation could easily be performed in Microsoft Excel, for example, without further adjustments. Evaluating the implied volatility for any model and any set of parameters (within the bounds) can then be carried out extremely quickly via interpolation. The region that is well approximated by this method is wide and applies to most parameter combinations and options. The region where this method struggles is the region where most methods do: options very close to maturity with only a few hours or days remaining. In these cases, asymptotic methods are appropriate (see e.g. Gao & Lee [21]). If some of the options being priced do not fall within the defined bounds (the bound on x should hardly ever be an issue, but the bounds on c IV might be), there is no harm in using your 1 For some values of x, the approximation is harder and we choose degree 9 for those. For other values, degree 7 is sufficient to obtain the desired error bound. Note that rational approximations of degree 9 can still be performed easily with double precision. One can improve all the approximations performed in this paper by switching to higher order rational approximations, although, high precision arithmetic would become necessary for degrees of 13 and higher. transform can be efficiently computed using the standard FFT algorithm, and one numerical Fourier transform gives option prices for a whole range of strikes with the same maturity. This method was introduced to finance by Carr & Madan [11] and later extended and generalised by a number of authors (e.g. Eberlein et al. [17], Lee [29], Lewis [31] and Raible [42]). More recently, other efficient versions of the Fourier method have been developed as well. However, as the standard FFT method is well-known and widely used in practice, we compare the rational approximation approach presented here to the standard FFT method, employing the Carr-Madan algorithm with the recommended values of parameters (dampening parameter α = 1.5, step-size η = 0.25 and the number of steps N = 2048). Note that for the FFT method, one usually requires the characteristic function of the log of the stock price at time t (where the stock price at time t is assumed to be as in Equation (2.1)), which is As can be seen, one therefore requires the Fourier transform of X t as defined in Equation (2.4) rather than the Laplace transform of Z t . The general formula to obtain the resulting characteristic for all three models discussed above. From this it follows that, for example, the characteristic function of the log-price X t in the VG model is equal to We next compare prices of vanilla call options computed from both methods for the three models discussed in Section 2. For each of the models, we select one set of market parameters from Table 1 and compare plain vanilla call option prices for strikes K ∈ [0.8, 1.2] and maturities T ∈ [0.25, 2.5].
Tables 5 -7 give absolute differences between the option price obtained by the FFT method and the rational approximation method. In particular, Table 5 represents errors for Case I, Table 6 represents errors for Case V and Table 7 represents errors for Case III. In these three tables, we set S 0 = 1, r = 0.03 and q = 0.01. As can be seen, the maximum values attained in Table 5, 6 and 7 are 1.61 × 10 −5 , 1.14 × 10 −6 and 8.59 × 10 −8 , respectively.
For all three tables, we precalculated 30 values of x T C equally spaced between x max T C and x min T C such that strikes K ∈ [0.8, 1.2] and maturities T ∈ [0.25, 2.5] are included. We iterate over ten different values for a and b and calculate errors for each of the ten cases, rather than calculating the value for a and b with the help of the cdf, as this procedure is much quicker. We take those values of a and b that lead to errors smaller than 1 × 10 −6 and delete x values for which no such combination can be found. In Equation (4.7), n was set to 7. The rational Chebyshev approximation in Equation Remark. Note that also these degrees (that is rational approximations with degree 6 and 8) can be easily performed with double precision. If one is interested in improving on the error bounds given here, it might be worth investigating how switching to higher order rational approximations, which will require high precision arithmetics, influences the errors. This could be particularly relevant for the few approximations performed offline, for which computational effort is not the main concern. We also leave a detailed comparison of different Gaussian quadrature methods for further research.     4) for a given strike, which can then be used for any maturity. Therefore, pricing relatively few strikes and relatively many maturities should be advantageous for the rational approximation method and vice versa. Table 8 gives approximate computational times for the two different methods, for different numbers of strikes and maturities. The first row resembles the example of Table 5, where we can see that the rational function approximation clearly outperforms the FFT method. As can also be seen, the rational function approximation is about five to ten times faster than the FFT method, depending on the number of strikes and maturities that need to be priced.   Matlab's method.

Conclusion
In this paper, we illustrate how to use rational function approximations to derive an alternative approach for calculating vanilla option prices and Black-Scholes implied volatilities. In particular, by approximating a normalisation of the Black-Scholes formula, a two-parameter function of one variable, we derive a pricing formula for any time-changed Brownian motion model. The method was found to be stable and fast, since it is possible to rewrite the rational approximation in terms of negative exponential moments (the Laplace transform) of the time change. The method can therefore be used for any model for which the Laplace transform of the time change is available in tractable form. Although rational Chebyshev approximations can be performed with little computational effort, we show that one can significantly speed up the method by pre-calculating the approximations and interpolating between them for the purpose of option pricing.
By way of illustration, we implemented the method for three popular models from the literature, and found that the method compares favourably to the standard FFT with respect to speed and accuracy.
Additionally, we use the methodology to approximate the inverse of the Black-Scholes formula in order to compute implied volatilities, and found this method to be significantly faster than existing methods. We are able to approximate the resulting two-dimensional function by a total of 105 rational approximations. Implied volatilities for a wide range of input parameters can then be obtained, for any given set of arbitrage-free vanilla option prices, by interpolating between the previously calculated approximations. The absolute error for this method is shown to be bounded by 5.5 × 10 −5 , although this can be refined if necessary by choosing rational approximations with higher order.

Appendix
The appendix is divided into two parts. We briefly review the general theory of Chebyshev approximations and of rational Chebyshev approximations in Section 9.1 and 9.2, respectively. For further background on rational approximations we refer to Ralston & Rabinowitz [43] (Chapter 7).
9.1. Chebyshev Approximations. The Chebyshev polynomial of degree n is defined as T n (x) = cos(n arc cos x) for − 1 ≤ x ≤ 1, (9.1) where the T has no relation to the T denoting maturity in Section 3. Even though these terms may look trigonometric, it should be noted that explicit expressions for T n (x) can be found by the following recursion T 0 (x) = 1 is exactly equal to f (x) at all N zeros of T n (x). Therefore, if we increase N , the two functions will be exact in more and more points.
If one inserts the expressions given for T n (x) in (9.2) into Equation (9.6), one obtains a polynomial in x, which approximates the function f (x) on the interval [−1, 1]. Even though this Chebyshev approximating polynomial is not equal to the minimax polynomial, which among all polynomials of the same degree has the smallest maximum deviation from the true function, it comes very close.
The minimax polynomial is very difficult and time consuming to calculate, whereas the Chebyshev approximating polynomial is almost identical and is extremely easy to compute.
Finally, if one is interested to approximate a variable y ∈ [a, b], one should use the following change of variable to map it on [−1, 1] and proceed as before b j (c i+r + c |i−r| ) r = 1, 2, ..., 2m + 1, (9.14) where a r = 0 if r > m and b 0 is chosen to be 1 (c j is defined as in Equation (9.5)). This leads to a linear system of equations for any m, which can be solved easily for the parameters a j and b j of Equation (9.9). To illustrate, the following is the resulting system of equations for a 4 × 4 rational This system of equations can be written in matrix form as follows: Substituting the resulting coefficients of a j and b j into Equation (9.9) gives a rational Chebyshev approximation of function f (x).