On the Efficacy of Fourier Series Approximations for Pricing European Options

This paper investigates several competing procedures for computing the prices of vanilla European options, such as puts, calls and binaries, in which the underlying model has a characteristic function that is known in semi-closed form. The algorithms investigated here are the half-range Fourier cosine series, the half-range Fourier sine series and the full-range Fourier series. Their performance is assessed in simulation experiments in which an analytical solution is available and also for a simple affine model of stochastic volatility in which there is no closed-form solution. The results suggest that the half-range sine series approximation is the least effective of the three proposed algorithms. It is rather more difficult to distinguish between the performance of the half-range cosine series and the full-range Fourier series. However there are two clear differences. First, when the interval over which the density is approximated is relatively large, the full-range Fourier series is at least as good as the half-range Fourier cosine series, and outperforms the latter in pricing out-of-the-money call options, in particular with maturities of three months or less. Second, the computational time required by the half-range Fourier cosine series is uniformly longer than that required by the full-range Fourier series for an interval of fixed length. Taken together, these two conclusions make a case for pricing options using a full-range range Fourier series as opposed to a half-range Fourier cosine series if a large number of options are to be priced in as short a time as possible.


Introduction
The need to price vanilla European options in a rapid manner arises in numerous activities at financial institutions. Perhaps the most well-known situation in which this need occurs in practice is model calibration, in which exotic options are priced using models with values for the (risk-neutral) parameters chosen in such a way as to ensure that the model reproduces quoted prices for liquid options. For each set of parameters considered in the search space, it is therefore necessary to evaluate the prices of all options in the calibration set using the model, before comparing these with the quoted prices. In a similar vein, there is a growing literature which uses large panels of options data for estimating model parameters ( [1]- [5]) in which the computational complexity is driven by the evaluation of vast numbers of option prices. Consequently, this paper investigates several competing procedures for computing the prices of vanilla European options, such as puts, calls and binaries and assesses their comparative performance on the basis of accuracy and computational speed.
Various strategies have been proposed for calculating the price of option contracts from knowledge of the conditional characteristic function of the underlying model. It is an important fact that a surprisingly large number of models have a semi-closed expression for their conditional characteristic function. For example, the identification of the conditional characteristic function for multivariate affine models with/without jump processes leads to the solution of a family of ordinary differential equations, albeit in the complex plane. In view of the Levy-Khintchine theorem, the identification of the conditional characteristic function for Levy processes is expressed in terms of various integrals with respect to the Levy measure.
The most commonly used techniques for taking advantage of a known conditional characteristic function have at their core the application of the Fast Fourier Transform (FFT). The best documented approaches is due to Carr and Madan [6] who construct an expression for the price of a European call option in terms of an integral over the characteristic function. This integral, which has an oscillatory kernel, is computed by an application of the FFT. Borak, Detlefsen and Hardle [7] apply the FFT strategy and demonstrate its efficacy by comparison with Monte Carlo simulation for a variety of models. Lord, Fang, Bervoets and Oosterlee [8] and Kwok, Leung and Wong [9] demonstrate how Fourier's convolution theorem in combination with the FFT can be used to price certain exotic options from knowledge of the conditional characteristic function of the price of the underlying asset. A different approach pioneered by Fang and Oosterlee [10] uses the characteristic function to directly approximate the marginal transitional probability density of returns by a Fourier cosine series. Recently Zhang, Grzelak and Oosterlee [11] have demonstrated how this methodology can be extended to the pricing of early-exercise commodity options under the Ornstein-Uhlenbeck process.
Rather than describing in detail the nuances of these various strategies, it is useful to point out what overarching assumptions connect them. Recall that the FFT is simply a clever piece of linear algebra that reduces the arithmetical load in implementing the Discrete Fourier Transform (DFT), namely the pair of equations connecting the coefficients of a finite Fourier series with values of the underlying function and vice versa. Therefore the decision to use the FFT implicitly makes the assumption that the underlying function is periodic over an interval of finite length, in practice determined by the range of frequencies submitted to the characteristic function, and that the function has been approximated over the interval by a finite Fourier series. The values of Fourier coefficients calculated from the characteristic function are in error by the extent to which the Finite Fourier Transform 1 differs from the Fourier transform.
Thus techniques using the FFT and those based on the construction of Fourier series share the same common assumptions and deficiencies. However, an important difference between an implementation using the FFT approach and one using the Fourier series approach is that the latter is parsimonious in its use of arithmetic whereas the former typically performs more arithmetic than necessary, albeit in an efficient way. For example, if the FFT is used to determine the value of a probability density function, what is recovered is the value of the function at each node of the interval, whereas all that might be needed is the value of the probability density function over a sub-interval.
The focus of this work is on the algorithm proposed by Fang and Oosterlee [10] who give a convincing demonstration of the efficacy of the Fourier cosine series. This series is more accurately called the half-range 2 cosine series because the actual function to be expanded is defined only on half the interval of periodicity (or range), the function being extended to the full range as an even-valued function. Half-range cosine series usually fail to represent derivatives whereas half-range Fourier sine series usually fail to represent function values. While the use of the half-range Fourier cosine series is a solid idea, Fang and Oosterlee [10] provide no motivation or explanation as to why this choice of approximating transitional density should be preferred over the half-range Fourier sine series or the full-range Fourier series for that matter. For example, intuition would suggest that the latter might perform better simply because it uses higher frequencies, which in turn translate to a more rapidly converging Fourier expansion. Indeed this intuition is borne out in calculation, but of course speed is not the only criterion of relevance in assessing the efficacy of a numerical procedure.
An important but subtle difference shared by the half-range Fourier cosine and full-range Fourier series approximations of density, but different from representations of density based on the half-range Fourier sine series, is that the former assign unit probability to the interval of support when in reality probability lies outside this interval, whereas the latter imposes zero probability density at the endpoints of the interval of support in contravention of reality, but on the other hand does not assign unit probability to the interval of support. Is one approach always superior to the other or is it a case of horses for courses? Intuition might suggest the latter. For example, when pricing a call option the most important component of the pricing error comes from the exclusion of contributions from asset price exterior to the finite interval of support. Because the half-range cosine and full-range Fourier series necessarily capture unit density, intuition might suggest that these approximations provide potential compensation for this component of pricing error. On the other hand intuition would suggest that the same approximations, when used to price binary options, might have a tendency to exaggerate the probability of exercise and therefore overprice this option in contrast to the half-range sine series approximation of probability density.

Fourier Series and Transform
Suppose that The use of the term "half-range" in describing Expressions (a) and (b) simply refers to the fact that the function a b , has for the construction of the Fourier series been extended into the in- , a b a − as an even-valued function in the case of the half-range cosine series (so that sine contributions vanish) and as an odd-valued function in the case of the half-range sine series (so that the constant and cosine contributions vanish). Thus both half-range series are conventional Fourier series taken over the interval − such that the function represented by the half-range Fourier cosine series is usually not differentiable at y a = , whereas that represented by the half-range Fourier sine series is usually discontinuous at y a = .
In the case of the half-range cosine and sine series in Expressions (1a) and (1b) respectively, the coefficients In terms of the exponential function the real coefficients k a and k b are calculated respectively as the real and imaginary parts of the single complex Equation In the case of the full-range Fourier series in Expression (1c) the real coefficients k ≥ are calculated respectively as the real and imaginary parts of the complex Equation Suppose now that ( ) f y is a transitional probability density function with known characteristic function defined formally by the Equation where ω ∈  and ( ) it can be asserted that ( ) f y ε < for any arbitrary small positive ε . The implication of this observation is that the Fourier coefficients in Equations (3) and (4) can be approximated from knowledge of the characteristic function via the respective formulae k while the coefficients of the full-range Fourier series can be approximated from knowledge of the characteristic function via the formula k The accuracy of approximations (6) and (7) is investigated in Section 6, where it is demonstrated that the error can be made arbitrarily small by choosing a suitably large interval.

Approximating Probability Density Functions
The quality of this practical idea is now explored for three trial probability density functions with known closed-form expressions for their characteristic functions. The first choice is the Gaussian density which may be regarded as representative of distributions with super-exponentially decaying tail density. The second and third choices are the Gamma density and the Cauchy density which are treated as representative examples of distributions with exponentially decaying and algebraically decaying tail densities respectively.

Gaussian Density
It is well known that the Gaussian density with mean value µ and variance 2 σ has characteristic function where π 4 n k n = . Figure 1 illustrates the quality of this approximation using 40 frequencies (solid line, 80 N = ) and using 4 frequencies (dashed line, 8 N = ). With as few as 4 frequencies it is clear that the approximating density still provides a good representation of the true density; with 40 frequencies the approximating density function is indistinguishable from the true density function, at least in terms of the resolution in Figure 1. It will be seen that this excellent performance is explained by the fact that the cumulative distribution function of the Gaussian probability density function converges to zero super-exponentially as x → −∞ and to unity as x → ∞ .

Gamma Density
The Gamma density with shape parameter α and scale parameter β has probability density function and characteristic function give by the respective formulae The approximating density is identical to Expression (8) with 8 b a − =. The task is now to reconstruct the Gamma probability density function from ( ) is also Gamma distributed, in this case with shape parameter α and unit scale parameter, then the appropriate comparison for the Gamma distribution is between the approximating Fourier series and a Gamma distribution with unit scale parameter. Figure 2 illustrates the quality of this approximation for 3 ,1 2 ) and 20 frequencies (dashed line, 40 N = ). With the exception of the region very close to the origin, the true density and approximated density are not significantly different for even 20 frequencies, and with 50 frequencies the difference between the true density and the approximating density is not discernible with the exception of the origin which does not present a difficulty since the density to known to be zero there.
The quality of the approximation is again due to the fact that the cumulative distribution function of the Gamma density converges exponentially to zero as x → ∞ . Typical values for the coefficient of mean reversion The important observation from both of these experiments is that distributions with exponentially decaying tail density can also be well described by a relatively small range of frequencies.

Cauchy Density
The Cauchy density with median µ and scale parameter α has probability density function and characteristic function give by the respective formulae The approximating density is again Expression (8) with 20 b a − = . The task is now to reconstruct the Cauchy probability density function from . The cumulative distribution function of the Cauchy density is ( ) which in this case converges algebraically to zero as x → −∞ and algebraically to unity as x → ∞ . Figure 4 illustrates the impact of this inferior rate of decay for the approximation based on 10 frequencies. While some erratic behaviour is evident in the tails of the approximating density, nevertheless the quality of the approximation is remarkably good considering the small number of frequencies in use.
In general, the approximate representations of the Gaussian, Gamma and Cauchy probability density functions all share the common fact that provides an upper bound for the error in the representation of density, where ( ) F x is the conventional cumulative distribution function of the underlying probability density function ( ) f x . In reality the error is significantly less than this upper bound due to cancellations resulting from the fact that the integrand in each of the Gaussian, Gamma and Cauchy densities is the product of a slowly varying probability density function decaying to zero and an independent rapidly oscillating function (in this case trigonometric functions represented in exponential form). However, the presence of this cancellation cannot change the character of the error which is entirely determined by the properties of ( ) F x as x → −∞ and as x → ∞ .

Comparison of True and Approximating Densities
While the illustrations of Figures 1-4 suggest that the Fourier series approximation of density is effective, it would be useful to quantify just how well density is approximated by the various Fourier methods. Instinctively it would seem that one possible way to achieve this objective is to use the Kullback-Leibler (KL) divergence criterion to measure the "distance" or departure of the probability density function ( ) q x from that of ( ) p x , which in the formulation of Equation (11) is taken to be the true probability density function. The measure appeals to the fact that 0 D ≥ with equality if and only if the density functions ( ) p x and ( ) q x are identical. However, the KL criterion is simply an application of Jensen's inequality for the concave function log x , and the properties of the criterion rely critically on the fact that ( ) p x and ( ) q x are probability densities, i.e. they are non-negative functions which integrate to unity.
The central idea of each Fourier approximation is to replace the true probability density function by a function of compact support, that is, ( ) 0 q x ≡ outside its interval of support. Consequently using measure (11) to compute the departure of the approximating density ( ) q x from the true density ( ) p x will typically require division by zero thereby rendering inappropriate the application of the KL criterion when used with ( ) p x as the true probability density. However, inverting the problem, one can use the KL criterion to measure the distance between the approximating density ( ) q x and the true density ( ) p x , that is, one may compute the value of D in Expression (11) taking ( ) p x to be the approximating density and ( ) q x to be the true density. In this case the integrand of expression (11) is well defined for all values of x , and therefore in order to proceed in this direction it remains to ensure that the approximating density integrates to unity. Approximations of density based on either the half range Fourier cosine series or the full range Fourier series are guaranteed to have this property, whereas an approximation based on the half range Fourier sine series need not share this property and therefore is inappropriate for analysis using the KL criterion. Tables 1-3 show values for the KL criterion for the Gaussian distribution with mean value one and unit standard deviation, the Gamma distribution with shape parameter 3 and scale parameter one, and the Cauchy distribution with median two and scale parameter one.
As might be anticipated, the Gaussian distribution is most efficiently approximated by Fourier methods followed by the Gamma distribution and finally the Cauchy distribution. However, it is clear that all of these distributions are well approximated by the half range Fourier cosine and full range Fourier series for sufficiently large intervals of support and an adequate number of frequencies. The results in these tables also reinforce the Table 1. Values of the KL measure (D c for the half range Fourier cosine series and D f for the full Fourier series) are given for the Gaussian density with unit mean and unit standard deviation. Values measure the deviation of the true density deviates from the approximating density for intervals of length 6, 8, 10, 12 and 14 standard deviations centred about the mean of the Gaussian distribution using 5, 10, 25, 50, 100 and 200 frequencies.  idea that the number of frequencies used in the approximation is of secondary importance to the size of the interval of support once sufficient frequencies are in use. This observation accords with intuition in the respect that larger intervals of support capture more of the true density and reduce the distortion in the approximating density, which as has already been commented, will always integrate to one. It would seem that approximations based on 50 frequencies are adequate in all of these examples. The results suggest that using more frequencies provides no meaningful improvement in accuracy. Moreover, for practical purposes there is little to choose between approximation based on the half range Fourier cosine and the full range Fourier series, although the latter has a slight edge for sufficiently large intervals and an adequate number of frequencies.

Pricing European Options
The successfully pricing of European option contracts for affine models of stochastic volatility requires knowledge of the marginal density of the asset price under the risk-neutral measure. The difficulty, however, is that no closed-form expression for this density is available for even the simplest of the multivariate affine models used in finance, although it is well known that such models have characteristic function, ( ) where T is the maturity of the option, T t τ = − is the backward variable and ω is the characteristic variable associated with the non-dimensional variable Y , here defined to be the logarithm of the ratio of asset price to strike price. The state variables (12) denote the values of the latent states of the system at time t . Typically the functions 0 , , M β β  are the solution of a system of ( ) . In overview, there is a well trodden procedure that starts with the specification of the multivariate affine model and ends with the construction of ( ) , t χ ω .

European Call Option
In the case of a European call option with strike price K and maturity T on an asset with spot price 0 S , the price of the option is where ( )

Binary Option
The price of a binary (call) option with strike price K and maturity T on an asset with spot price 0 S is As is the case with the pricing of a call option, the frequencies used in approximation (c) are double those used in approximations (a) and (b) suggesting that calculations based on Expression (c) can be expected to converge more rapidly than those based on Expressions (a) and (b).

The Heston Model of Stochastic Volatility
Heston's [12] risk-neutral model of stochastic volatility has expression where log Y S K = , V is the diffusion of asset price, r is the risk-free rate of interest, q is the dividend rate, κ  is the risk-neutral rate of mean reversion of volatility to the risk-neutral long run value of γ  , σ scales the volatility of diffusion, ρ is the local correlation between returns and volatility and be the characteristic function of ( ) , , , , f Y V t y v T with respect to the forward variables, that is, , , , , , , , , e d d .
with terminal condition Thereafter, it is straightforward to show that the anzatz On a practical note, the fact that the characteristic function of ( ) log S K embeds the difference between the risk-free rate of interest and the dividend rate suggest at first sight that the coefficients 0 , , M β β  must be computed whenever the risk-free rate or dividend rate changes, potentially each day, and not just when the parameters of the model are changed. The key to avoiding this difficulty is to note that the difference ( ) r q − enters the calculation of 0 β alone, and that the behaviours of 1 , , M β β  are independent of the behaviour of 0 β . Consequently, the resolution of this practical dilemma is to divide the calculation of 0 β into two stages with the first stage treating only those calculations that involve the difference ( ) r q − and the second stage treating everything that is independent of the value of ( ) r q − . Both stages are brought together each day in the calculation of the expected payoff, but only the first stage of calculation needs be done on a day-to-day basis. In particular, this calculation is very easy as it amounts to adding ( ) n r q Tω − to the second stage calculation of 0 β in the computation of ( ) n χ ω .

Error Analysis
Let ( )

Truncation Error
Truncation error occurs when the semi-infinite interval in Expression (13) or (16) where ( ) The inference from this observation is that approximations of ( ) f y  which place unit density in [ ] , a b may potentially compensate for some of the error incurred in using the truncation procedure to price a call option. The half-range Fourier cosine and full-range Fourier series approximations of transitional probability density both fall into this category, which in turn suggests that the use of the truncation procedure with these approximations of marginal density may well provide rather smaller pricing errors than might casually be anticipated. By contrast, the half-range Fourier sine series approximation of marginal density does not capture unit density in the interval [ ] , a b , and the corresponding conjecture is that this approximation of marginal density will perform less well than the half-range Fourier cosine and full-range Fourier series approximations with regard to the pricing of a call option. Table 4 indicates that this conjecture is borne out in practice.

Approximation Errors
Suppose now that the transitional density  Table 4. Percentage errors recorded by the half-range cosine series, half-range sine series and full-range Fourier series of pricing a European call option and a binary option with strike K = 1200 and maturities one month (panel (a)), three months (panel (b)) and six months (panel (c)). The factor refers to the multiple of T σ used to establish the interval over which to compute the numerical approximation.
Consequently a misspecification error arises in the Fourier coefficients. Second, the truncation of the Fourier series (27) to a finite number of terms, say N terms, introduces another approximation error. Therefore the price of a call option based on this strategy is The misspecification error in the coefficients k a and k b due to the use of the coefficients k A and k B respectively is determined from the identity ( ) ( It therefore follows directly that 1 1 , , for all values of k . Thus the misspecification error in the Fourier coefficients is driven entirely by the choice of = . In practice this error will be significantly smaller than the maximum bound given in Equation (33) as a result of arithmetical cancellation due to the oscillatory nature of the integrand, but as commented previously, the decay of this error as [ ] , a b →  will be determined by the behaviour of ( ) F y  as y → −∞ and as y → ∞ .
In conclusion, the total error in pricing a call option, namely The contributions to the error from the first, second, third and fourth terms on the right hand side of inequality (36) are dominated by the behaviour of ( ) The well known result that if In conclusion, the error in pricing a call option can be made arbitrarily small by restricting the marginal probability density to a finite interval [ ] , a b and approximating the density in that interval by a misspecified Fourier series constructed from the characteristic function at the appropriate frequencies.

Performance under Simulation
A series of simulation experiments was undertaken in order to examine the efficacy of the half-range Fourier cosine series, half-range Fourier sine series and full-range Fourier series in respect of how accurately these approximations price European call options and binary options. The first experiment prices options in a Black-Scholes world so that a closed-form solution may be used to assess the pricing error, whereas the second experiment prices the same options using Heston's model of stochastic volatility.

Black-Scholes Pricing
Assume that the asset price, S , follows a geometric Brownian motion in which dW is the increment of a Wiener process. The advantage of using this specification is that, following the seminal work of Black and Scholes [13], exact prices are known for each type of option for all combinations of spot price 0 S , strike price K , and maturity T . It follows, therefore, that the relative percentage pricing error incurred using numerical methods based on Fourier series for a variety of different strike prices and maturities can be identified exactly.  Tables 4-6 with each table reporting the results for maturities of one month (panel (a)), three months (panel (b)) and six months (panel (c)).
Two very clear general conclusions emerge from these results. 1) For options that are deep in-the-money ( Table 6) in which 0 1000 S = and 800 K = , or at-the-money ( Table 5) in which 0 1000 S = and K = 1000, it matters little (if at all) which of these three pricing algorithms are chosen irrespective of the maturity of the option. Indeed the relative pricing error is zero for all practical purposes.
2) For options that are deep out-of-the-money (Table 4) in which 0 1000 S = and 1200 K = , the choice of algorithm is more important. The results may be summarized succinctly as follows.
a) The half-range Fourier sine series does not perform as well as the other two approximations and its use is therefore not recommended. This finding accords well with our previous intuition based on a consideration of the contribution to the price of a call option lost as a result of truncating marginal density.
b) The half-range Fourier cosine series and the full-range Fourier series both perform relatively well. When the size of the interval of approximation is a relatively small multiple of T σ , namely either 5 or 6, then the half-range Fourier cosine series performs better that the full-range Fourier series. As the multiple of T σ increases and the size of the interval of approximation becomes larger, the full-range Fourier series begins to dominate. When the interval of approximation has size 10 T σ , then the full-range Fourier series is unambiguously superior to the half-range Fourier cosine series, particularly for options of short maturity.
On the basis of this analysis and on accuracy grounds, it is hard to ignore the claim that the full-range Fourier series is the algorithm of choice when using Fourier methods to price options. Moreover, the full-range Fourier series converges faster than either the half-range sine or cosine series and is therefore likely to price options more rapidly. These themes are explored in more detail in the pricing of call options under Heston's model of stochastic volatility.

Heston's Model
A total of approximately 40,000 options over ten years were generated by simulation of Heston's model. Sixteen Table 5. Percentage errors recorded by the half-range cosine series, half-range sine series and full-range Fourier series of pricing a European call option and a binary option with strike K = 1000 and maturities one month (panel (a)), three months (panel (b)) and six months (panel (c)). The factor refers to the multiple of T σ used to establish the interval over which to compute the numerical approximation. options were generated each day spread over 4 maturities ranging from 92 to 5 days and 4 strikes two of which are initialised at 20% out of and into the money and two of which are initialised at 7% out of and into the money. The half-range Fourier cosine series and the full-range Fourier series are then used to price call options and binary options with these strikes. In this instance no exact solutions are available, and so the accuracy of each method in respect of each type of option is gauged by comparison against values calculated using a large interval. The left hand and middle columns of Table 7 show respectively the 2 L and 1 L relative pricing errors for out-of-the-money call options calculated using the half-range Fourier cosine series and the full-range Fourier series. The right hand column of Table 7  1 252 DT = . A clear finding from Table 7 is that the half-range Fourier cosine series performs well in respect of the 2 L measure for shorter intervals and in the 1 L measure for almost all choices of interval length. By contrast the Table 6. Percentage errors recorded by the half-range cosine series, half-range sine series and full-range Fourier series of pricing a European call option and a binary option with strike K = 800 and maturities one month (panel (a)), three months (panel (b)) and six months (panel (c)). The factor refers to the multiple of T σ used to establish the interval over which the numerical approximation is computed. ( performance of the full-range Fourier series is poor with regard to both measures for shorter intervals. However, the quality of approximation provided by the half-range cosine series is erratic as the size of the Fourier window increases whereas that provided by the full-range Fourier series improves systematically to the extent that its performance surpasses that of the half-range Fourier cosine series for intervals of length 24 standard deviations or more. Furthermore, this level of accuracy is achieved by the full-range Fourier series in approximately 25% quicker than that required by the half-range Fourier cosine series. Fang and Oosterlee [10] suggest choosing intervals of length 20 standard deviations. In this problem, it is evident that the half-range Fourier cosine series still enjoys an advantage over the full-range Fourier series for intervals of this length. The suggestion of this investigation is that 20 standard deviations should be regarded as a minimum length of interval. Table 8 shows the 2 L and 1 L relative pricing errors calculated from the halfrange Fourier cosine series and the full-range Fourier series in respect of in-the-money call options based on Heston's model. Table 7. The L 2 and L 1 relative pricing errors calculated from the half-range Fourier cosine series and the full-range Fourier series for out-of-the-money call options based on Heston's model. The timings refer to the times (seconds) needed to price 12 daily call options over 10 years with maturities ranging from 5 to 92 days. The factor refers to the multiple of σ Y used to establish the interval over which the numerical approximation is computed.  The results reported in Table 8 exhibit a similar pattern of behaviour to those reported in Table 7, namely that for intervals of short length the half-range Fourier cosine series outperforms the full-range Fourier series with respect to both the 2 L and 1 L measures of accuracy, but that this dominance vanishes when the length of the interval of approximation is 24 standard deviations or more. In particular, both approaches generate noticeably smaller relative errors for in-the-money options than for out-of-the-money options. Because in-the-money call options carry a higher price, one inference of this observed reduction in relative error is that the absolute pricing error in using each Fourier representation is insensitive to the moneyness of the option. Table 8 and Table  9 report the performance of the half-range Fourier cosine series and the full-range Fourier series when used to price binary options.
The results presented in Table 9 and Table 10 demonstrate that both the half-range Fourier cosine series and the full-range Fourier series both generate good estimates of the price of binary options with the former performing better than the latter for intervals of short length, but with this advantage disappearing when the length of the interval is increased. Table 9. The L 2 and L 1 relative pricing errors calculated from the half-range cosine series and the full-range Fourier series for out-of-the-money binary options based on Heston's model. The timings refer to the times (seconds) needed to price 12 daily call options over 10 years with maturities ranging from 5 to 92 days. The factor refers to the multiple of σ Y used to establish the interval over which the numerical approximation is computed.  Table 10. The L 2 and L 1 relative pricing errors calculated from the half-range cosine series and the full-range Fourier series for in-the-money binary options based on Heston's model. The timings refer to the times (seconds) needed to price 12 daily call options over 10 years with maturities ranging from 5 to 92 days. The factor refers to the multiple of σ Y used to establish the interval over which the numerical approximation is computed. The results reported in Tables 7-10 are characterised by two common denominators. First, the computational time required by the half-range Fourier cosine series is uniformly longer than that required by the full-range Fourier series for an interval of fixed length. The simple explanation for this observation is that the full-range Fourier series uses larger frequencies for a given length of interval, and therefore the full-range Fourier series converges more rapidly. Second, the pricing of call options and binary options using the half-range Fourier cosine series representation of transitional density is noticeably better than the corresponding pricing using the full-range Fourier series for short intervals [ ] , a b , but this advantage vanishes in the case of the binary option when the interval becomes suitably large and is reversed in the case of a call option. The primary explanation for this behaviour lies in the realisation that tail density is characterised by long wavelengths, or equivalently by the presence of low frequency terms, whereas peaky density is characterised by short wavelengths, or equivalently high frequency terms.
For an interval of given length, the frequencies present in the half-range Fourier cosine expansion are smaller that the frequencies present in the full-range Fourier series, and therefore when the length of the interval [ ] , a b is small, the half-range Fourier cosine series better captures tail behaviour. Increasing the length of the interval lowers the range of frequencies admitted to the full-range Fourier series expansion, which in turn improves its ability to handle tail density. However, the same increase in length will benefit the half-range Fourier cosine series representation of transitional density only provided there remains significant tail density to be characterised by these new frequencies. On the other hand, the half-range Fourier cosine series will always suffer from the drawback that the (periodic) function it represents has gradient zero at x a = and x b = (because the series is not generally differentiable at these points). Consequently, the half-range Fourier cosine series approximation will always misrepresent the gradient of the true density function at x a = and x b = in contrast to the behaviour of the approximation provided by the full-range Fourier series.

Conclusion
One clear conclusion from these calculations is that the half-range Fourier cosine series and the full-range Fourier series perform uniformly better than the half-range Fourier sine series. The half-range Fourier cosine series and the full-range Fourier series both perform with credit. When the length of the interval [ ] , a b is relatively small, say ten or so standard deviations, it is clear that the half-range Fourier cosine series outperforms the full-range Fourier series over the same interval. On the other hand for intervals of larger length the full-range Fourier series is at least as good as the half-range Fourier cosine series, and outperforms the latter in pricing out-of-the-money call options, in particular, with maturities of three months or less. In general the full-range Fourier series outperforms the half-range Fourier cosine series in all circumstances provided that the interval [ ] , a b is suitably large, although the effect is so small as not to be significant in practice. The explanation for this behaviour lies in the fact that the half-range Fourier cosine series, although not representing a differentiable function, always uses a lower spectrum of frequencies than the full-range Fourier series and therefore enjoys an initial advantage in describing tail density. As the interval length is increased, this advantage vanishes. On the other hand the larger spectrum of the full-range Fourier series guarantees more rapid convergence. Of course, timing issues may not be important if small numbers of model call prices are to be calculated, but when facing calibration problems or estimation problems as described in the introductory remarks, timing issues are indeed significant once adequate numerical accuracy is assured.