Pricing of Discretely Sampled Asian Options Under Levy Processes

We develop a new method for pricing options on discretely sampled arithmetic average in exponential L'evy models. The main idea is the reduction to a backward induction procedure for the difference <i>W_n</i> between the Asian option with averaging over <i>n</i> sampling periods and the price of the European option with maturity one period. This allows for an efficient truncation of the state space. At each step of backward induction, <i>W_n</i> is calculated accurately and fast using a piece-wise interpolation or splines, fast convolution and either iFFT or the parabolic inverse Fourier transform. Numerical results demonstrate the advantages of the method.


Introduction
An Asian option is an option whose terminal payoff depends on the average values of the underlying asset during some period of the option's lifetime. Thus, an Asian option is path-dependent. Asian-style derivatives constitute an important family of derivative securities with a wide variety of applications in financial markets. There are a number of economic reasons to write a contract with the an averaging feature. For example, in the foreign exchange market, Asian options allow one to reduce risk of exchange risk fluctuations, and such an option is typically cheaper than European options. Another reason is that the price manipulation becomes more difficult, even in the thinly traded asset market near option's maturity.
There are Asian options both of the European-style and the American-style, the averaging can be either arithmetic or geometric, and the sampling can be either continuous or discrete. The payoffs of continuously sampled arithmetic average Asian options are of the form 1 where w t are weights; if the sampling is discrete, the payoffs are of the form where 0 ≤ T 0 < T 1 . . . < T N = T are the sampling dates. It is easy to see that the joint distribution of the arithmetic average is quite complicated to characterize analytically, and derivation of efficient explicit formulas for prices of Asian option is very involved, with some exceptions.
1.1. Pricing in the exponential BM model. In the exponential Brownian motion model (Black-Scholes model), pricing of the continuously sampled geometric average options is easy and quite straightforward (and the same holds in Lévy models and a number of other popular models). However, in the same exponential Brownian motion model, pricing the arithmetic average Asians options is far from trivial. In early studies of continuously sampled Asian option, such options were approximated by the geometric average Asian options. This approach significantly underprices the option (see the discussion in Musiela and Rutkowski [44] and the references therein). To overcome this deficiency, a number of authors have suggested various analytical approximations for the distribution of the arithmetic average. For example, Turnbull and Wakeman [48] use the lognormal approximation with matched first and second moments, Milevsky and Posner [43] use the reciprocal gamma approximation. The problem with approximations of this sort is that no reliable error estimates are available.
One of the most celebrated analytical tools is to use the Laplace transform of the Asian option price, introduced by Geman and Yor [32]. Later, the approach has been extensively studied, for example, by Fu, Madan and Wang [29], Carr and Schröder [19,20], Shaw [46]. Linetsky [37] takes a new direction, and derives analytical formulas using spectral expansion for the value of the continuously sampled arithmetic Asian option.
For the discretely sampled Asian options, another approach is to numerically evaluate the density of the sum of random variables as the convolution of individual densities. The second step of this method involves numerical integration of the option's payoff function with respect to this density function. The method is initiated by Carverhill and Clewlow [21], who rely on the use of the fast Fourier transform to evaluate the joint probability density function.
The pricing of Asian option in the Black-Scholes model has been dealt with by a host of researchers, and the list of papers above is by no means complete.
1.2. Pricing in exponential Lévy models. Lévy models provide a better fit to empirical asset price distributions that typically have fatter tails than Gaussian ones, and can reproduce volatility smile phenomena in option prices. For an introduction to applications of these models applied to finance, we refer the reader to [14] and [25]. Pricing the continuously sampled geometric average options in exponential Lévy models is easy and quite straightforward (see [14]) but pricing of arithmetic Asians presents serious mathematical and computational difficulties.
For continuously sampled Asian options, Bayraktar and Xing [7] derive analytical formulas which can be realized numerically fairly fast and accurately. Cai and Kou [17] obtain a closed-form solution for the double-Laplace transform of Asian options under hyper-exponential jump diffusion model, and suggest a numerical procedure for the realization of this closed form solution. The variant of the PDE-approach due to Večeř [49,50] was extended to processes with jumps by Večeř and Xu [51].
For discretely sampled Asian options, Benhamou [8], Fusai and Meucci [30] enhance the method of approximating the joint probability density function of [21], and price Asians under Lévy processes. Since the transition density function of a Lévy process is generally unknown in the closed form, the accuracy of this approach crucially rely on the accuracy on an approximation of the transition density function.Černý and Kyriakou [22] reduce the pricing problem to a sequence of European options in the one-factor model, and use trapezoidal rule as the numerical realization of the (inverse) Fourier transform to approximate the prices of European options.
More recently, Chen et al. [23,24] developed a Monte Carlo algorithm for simulating Lévy processes, based on calculation of the pdf using the Fourier inversion, and applied this algorithm to pricing discretely sampled Asian option. In Albrecher et al. [1,2,3,4], Lemmens et al. [35], pricing bounds for Asian options under Lévy processes are derived and hedging strategies analyzed.
1.3. General structure of our method. In the present paper, we consider discretely averaged Asian options. AsČerný and Kyriakou [22], we first reduce the pricing problem to a series of pricing of European options. However, the efficiency of an approach of this kind strongly depends on the type of the constructed sequence of European options, and the efficiency of the pricing the European options which inevitably have rather complicated payoffs. Therefore, it is necessary to employ additional tricks to enhance the accuracy and speed of calculations, and an efficient control of several sources of errors is of the paramount importance.
We use the reduction to a series of European options, whose terminal payoffs vanish at −∞ and on [0, ∞). Given the error tolerance, we choose x 1 < x M < 0, and, for each option, replace the payoff on (−∞, x 1 ] with the leading term of asymptotics, which is of the form c n e x . This trick is efficient, and |x 1 | can be chosen moderately large in absolute value (typically, less than 5). On [x M , 0], we set the payoff to 0. On [x 1 , x M ], we approximate the payoff by a piece-wise polynomial function or use splines. We derive the bounds for the modification and interpolation errors. The derivation is fairly nontrivial due to a rather complicated structure of the payoff. These bounds are used to give recommendations for the choice of the parameters of our numerical scheme.
At each step of the induction procedure, after the modification of the payoff and interpolation have been made, we can calculate the price of the European option explicitly using the Fourier inversion. We rewrite the resulting formula as a sum of products of the values of the payoff function at the points of the chosen uniform grid, and values of two auxiliary functions, at points of a grid of the form ∆, −M ≤ ≤ M . These two functions are the same at each step of the backward induction procedure (in the case of higher order interpolation, more than two functions are needed). The sum can be represented as two terms plus a discrete convolution, which can be realized very fast using the fast convolution algorithm (this idea goes back to Eydeland [27]). The values of the auxiliary functions can be calculated only once (and stored); various realizations of the inverse Fourier transform can be used for this purpose. We use the simplified trapezoid rule, which is very efficient if the integrand is analytic in a strip, and allows for an efficient error control. The reason is that the disretization error of the infinite trapezoid rule decays exponentially as the mesh approaches 0 if the integrand is an analytic function in a strip around the line of integration. Crucial to this analysis are the fundamental properties of the Whittaker cardinal series (Sinc expansion) for functions that are analytic in a strip, introduced to finance in [28]. The conformal change of variables introduced in [15] (parabolic iFT method) allows one to greatly increase the rate of the decay of the integrand at infinity, and decrease the number of terms in the simplified trapezoid rule and the CPU time.
1.4. Organization of the paper. In Section 2, we explain the main idea and outline the steps of our method. Error estimates and recommendations for the choices of the parameters of the scheme given the error tolerance are derived in Sections 3 and 4. In Section 5, we present an explicit algorithm, and produce numerical examples to illustrate the relative performance of different methods. Section 6 concludes. Technical details are relegated to Appendix A and Appendix B. Several possible directions in which our method can be developed further are outlined in Appendix C.

Pricing Asian calls and puts
2.1. Option specification and the general scheme of calculation. We assume that the riskless rate r ≥ 0 is constant, the market is arbitrage free, the underlying asset pays no dividends, and, under an EMM Q chosen for pricing, the spot price process for the underlying is an exponential of a Lévy process X: S t = e Xt . Let be the arithmetic average of the asset price, where the dates (0 =)T 0 < T 1 < · · · < T N (= T ) are specified in the contract. The terminal payoffs of the Asian put and call options with strike K and maturity date T are respectively; and the time-0 price of the corresponding option is given by where E is the expectation operation under Q, and S 0 is the current spot price of the underlying. We will consider the pricing of the Asian put option; pricing the Asian call option can be reduced to pricing the Asian put option via the put-call parity: , where V is the price of the Asian put option. The first term on the RHS is the weighted sum of S 0 and the moment generating functions of X T j |X 0 , j = 1, 2, . . . , N , hence, it can be easily calculated. It remains to calculate the price of the Asian put option.
Assuming further that the sampling dates are equally spaced: T j = j∆, where∆ = T /N , and taking into account that X is a Lévy process, we rewrite (2.1) as and Using the law of iterated expectations, we obtain To reduce the truncation error in the state space, we introduceW n (x) = V n (x) − V 1 (x). For n = 1, W 1 (x) ≡ 0; for n = 2, and for n = 2, 3, . . . , N − 1, and, for n = 2, 3, . . . , N − 1, define Formally, the calculation of the option price using (2.6) and (2.7) is quite straightforward: and calculate W n+1 ; 4. calculate V N (y) = W N (y) + V 1 (y) at y = ln S 0 − ln((N + 1)S 0 − K), and then the option value using (2.3).
Since V 1 is the price of the put option with strike 1, maturity∆, in the model with zero riskless rate, V 1 can be easily calculated using the Fourier transform technique (see [15] for the review). Since the Fourier transform of the piece-wise polynomial interpolant is a rational function, which can be easily calculated, W n can be calculated almost as easily as V 1 , but it is non-trivial to take the several errors involved into account accurately, and derive sufficiently simple and accurate recommendations for the choice of the parameters of the numerical scheme.
In the remaining part of this section, we list the main formulas of the numerical scheme; the derivation of error bounds, recommendations for the choice of the parameters of the numerical scheme, and proofs are relegated to the following sections.
If, in addition, the values V 1 (x k ) are needed at points of an equally spaced grid x j = x 1 + (j − 1)∆, j = 1, 2, . . . , M , then the iFFT can be used to increase the speed of calculations.
If λ + ≤ min{1, −λ − − 1}, it is advantageous to shift the line of integration below the real line so that the strip of the analyticity of the integrand becomes wider. A shift of this sort is especially useful if x := x + µ∆ is negative and large in absolute value as it is quite often the case in the lower part of the x-grid. For the detailed analysis of various realizations of (2.16), see Section 4.
If ν and/or∆ are small then the integrand decays very slowly at infinity, and too many terms may be needed to satisfy the desired error tolerance. In these cases, the preliminary conformal deformation of the contour of integration in (2.16) with the subsequent change of the variables (parabolic iFT) can be used to greatly decrease the number of terms in the simplified trapezoid rule (see Section 4 and [15]).
An alternative is the refined version of iFFT, which reduces the calculation of the sum in (2.17) to application of several copies of iFFT (see [10,11]), but if the integrand decays very slowly then parabolic iFT becomes indispensable.
2.4.1. Modification of f n . Let f n , n = 1, 2, . . . , N − 1 be as in (2.8) and (2.9). The intervals of integration in (2.6) and (2.7) will be truncated from above at a point x M < 0, equivalently, f n will be set to 0. On (−∞, x 1 ], we replace each integrand in (2.6) and (2.7) with the leading term of the asymptotics as x → −∞. The bounds of the modification errors and recommendations for the choice of x 1 and x M (given the desired error tolerance at each step) will be given in Section 3.1. The starting point is the following lemma.
Lemma 3.2. Let X be a Lévy process, whose characteristic exponent is analytic in a strip Im ξ ∈ (λ − , λ + ) around the real axis, and f n be as in (2.9). Then, for n = 2, 3, . . . , N , and any ω ∈ [0, λ + ), we have Proof. W n are expectations of functions that are bounded by 1 and vanish above 0. Hence, It remains to apply the Taylor formula to (1 − e x ) around x = 0.
To calculate W n+1 using (2.7), we need values of W n and W 2 . When we use an interpolation procedure to approximate f n (x) in (2.8)-(2.9), we need values of the latter function on (−∞, 0). We truncate f n (x) at x = x M < 0, where x M is found from the condition |f n (x)| ≤ , ∀ x ∈ (x M , 0), and > 0 is the error tolerance. Let 1 be the truncation error at each step, and set = 1 /2. Using (3.2), we find an approximate recommendation where ω + ∈ (0, λ + ). If λ + is not very large, and ψ(iλ + ) < ∞ as it is the case with KoBoL and NIG, then we can use (3.4) with ω + = λ + . If λ + is not large and X is VG, or if λ + is very large, we can take ω + so that∆ψ(iω + ) = 0.1 ln , and set On the strength of Lemma 2.1, we choose the lowest point where C is the constant in the O-terms in (2.18) and (2.19). Clearly, we can find x 1 from The constant C is unknown. If is very small, then ln( /C) ≈ ln , and we can use an approximate recommendation (3.6) If X is KoBoL or NIG and λ − is not very large, then we may use (3.6) with ω − = λ − /2. Otherwise, we can take ω − so that∆ψ(iω − ) = 0.1 ln .
Since |E x [u(X∆)]| ≤ u L∞ , the total truncation error is bounded by the sum of truncation errors at each step of backward induction. If tr is the error tolerance for the total truncation error, and N is the number of steps, we choose 1 = tr /(N − 1). When we derive estimates for the other types of errors at the current step, we may assume that there is no truncation at all. Similarly, we may assume that the calculation of the expectation that serves as an input at the next induction step contains no truncation error. The same argument will be used below to account for the impact of the other sources of errors. In the nutshell, each error estimate can be derived as if there have been no errors before.

3.2.
Interpolation: preliminary estimates. Denote by p∆ the transition density of X. In Appendix A.2, we prove Then, ∀ s ∈ Z + , the following approximate bound holds Then the following approximate bound holds If (2.10) holds with ν > 0, the option price is of class C ∞ ; however, for small ν and/or ∆ the derivatives can be too large, and therefore, a higher order interpolation will have a larger interpolation error. In the VG model, the option price is rather irregular unless ∆ is fairly large. In the result, the justification of even the piece-wise linear interpolation used in the backward induction procedure is possible only for sufficiently large∆, and a higher order interpolation is, essentially, impossible in the sense that it will lead to large discretization errors.
Let f n be as in (2.8)-(2.9). In Appendix A.3, we prove Lemma 3.4. For any integer n ≥ 1 and s ≥ 2 such that p In the case of the piece-wise linear interpolation, we need the simplest case of bound (3.9): for any n ≥ 1, 3.3. Interpolation error estimate, and choice of the interpolation procedure and grid. A grid x on [x 1 , x M ] is chosen to interpolate functions f n ; the integral over (−∞, x 1 ] is calculated replacing f n (x) with c n e x , where c n is given by (2.20) (this is the leading term of asymptotics in (2.18) (resp., (2.19)) if n = 1 (resp., if n ≥ 2)). We use a uniformly spaced grid Assuming that a polynomial or spline interpolation procedure f n;app of order s is chosen, with the error estimate of the form where C s is a universal constant 1 , we need to derive an estimate for the norm of the derivative on the RHS. Lemma 3.4 provides simple approximate upper bounds (3.9) for f (s+1) n L∞ , hence, for the discretization error. For instance, in the case of the piece-wise linear interpolation, using the error estimate (3.11) and the bounds (3.10), we find that the total sum of interpolation errors at all steps of the backward induction procedure, err int , admits an approximate upper bound via By simplifying the RHS, we have For a process of order ν > 0, applying (3.7), we obtain Similarly, one can obtain an error bound for the case of VG model. Let int be the error tolerance allocated for the total interpolation error. Lemmas 3.4 and 3.3 taken together allow us to choose ∆ as a function of int ,∆, parameters of the process, and s, the order of the interpolation procedure. We choose s with the maximal ∆ = ∆ s ; if several ∆ s are rather close, we choose the interpolation procedure with the smallest s.
If one uses an interpolation procedure to obtain values of V 1 and W n at points x = y k := x k − ln(1 − e x k ) + , k = 1, 2, . . . , M (see Step 4 in Section 2.4.3), one needs to take into account the errors of this interpolation. Suppose that a polynomial or spline interpolation procedure V 1;app and W n;app of order s is chosen, with the error estimates similar to (3.11) We will prove in Lemma A.3, that the constants in the RHSs of (3.13) and (3.14) are smaller than the constant in the bound in (3.11), hence, it suffices to choose the interpolation scheme and ∆ as above, for int /2 instead of int .

Calculation and errors in the dual space
In this section, we consider the case of the piece-wise linear interpolation, and analyze the impact of errors of calculation of values of V 1 and V 2 on the error of W n calculated using (2.23). (The interpolation procedures of higher order can be analyzed similarlysee Appendix B.) We need to satisfy a rather large error tolerance for two sequences (of values of V 1 and V 2 /∆) in L ∞ -norm when a very small factor u 1 is present and a fairly small error tolerance in L 1 -norm of the sequence of values of V 2 /∆. In both cases, it is unnecessary to calculate the values with high accuracy near x = 0, where it is especially difficult to satisfy a small error tolerance (see [15]).
Let c > 0 be the desired error tolerance for the error induced by calculation of V 1 and V 2 at each step. Denote by (V, y) the absolute errors of the calculation of function V = V 1 , V s using flat iFT or parabolic iFT. Then the absolute error of the calculation of E x [u(X∆)] using (2.23) is bounded by Since |u 1 | is small, we can make the first two term in (4.1) very small using a dual grid of a fairly large mesh ζ and rather small number of points M c . For small ∆, the third term admits a fairly accurate bound via the L 1 -norm of d (V 2 , ·) and t (V 2 , ·) times the L ∞ -norm of u. The latter being bounded by 1, it suffices to satisfy the bounds  [15]. We first consider the case of V 1 . Denote by g(ξ) the integrand in (2.16) multiplied by (2π) −1 . Function g is analytic in the strips Im ξ ∈ (0, λ + ), Im ξ ∈ (−1, 0), Im ξ ∈ (λ − , −1). With the choice ω ∈ (0, λ + ), we need to consider the first strip. (In some cases, it is advantageous to push the line of integration down, either into the second strip or the third one. We consider these cases later.) and the Hardy norm is finite: Denote by d (V 1 ; ·; ω, ζ) the error of the replacement of the exact option price (2.16) with the infinite sum (2.17), and, for ω ∈ (µ − , µ + ), set d(ω) = min{ω − µ − , µ + − ω}. Our starting point is Corollary 4.2. [15, Corollary 2.2] Let the error tolerance > 0 for the discretization error be small so that g : To derive fairly accurate prescriptions for an almost optimal choice of ζ for a given error tolerance, the following estimate of the Hardy norm is useful.
, the error of the replacement of the exact integral (2.22) with the infinite sum (4.7) admits an upper bound via given by (4.5) admits an estimate via if s = 4, 6, . . .
The change of variables η → ω tan θ transforms the integral into and the statement of the proposition follows.
The integrand in the Fourier inversion formula which defines V 1 decays at infinity as the integrand in the formula for V 2 , therefore, it suffices to consider the truncation error of the infinite sum in the formula for V s , s ≥ 2. Making the change of variable ξ → η +iω, we rewrite (2.22) as If ζ is small, the truncation error of the replacement of the infinite sum (4.7) with the finite sum Proposition 4.5. In the case of KoBoL of order ν ∈ (0, 2), ν = 1, Therefore, we have an approximate upper bound Making the change of variable:∆C ∞ η ν → t, we obtain the bound (4.12).
Note that the integral ∞ x e −t t a−1 dt is the upper incomplete gamma function (if a > 0) or an exponential integral (if a = 0) or reducible to the upper incomplete gamma function via integration by parts (if a < 0).

4.1.2.
The choice of line of integration. As it was mentioned in Section 2.3, in order to exploit the wider strip of analyticity, a shift of the line of integration might be necessary.
For V 1 in (2.16), the integrand has simple poles at −i and 0, and is analytic in the strips Im ξ ∈ (0, λ + ), Im ξ ∈ (−1, 0) and Im ξ ∈ (0, λ + ). It follows from Theorem 4.1 that the wider the strip of analyticity is, the smaller the discretization error can be made by an appropriate choice of the line of integration inside the strip. Therefore, if the three strips above are of sizably different widths, it is advantageous to choose the widest strip. We can move from one strip to another using the residue theorem.
If we push the line of integration in (2.16) down and cross the pole at zero but remain above the second pole, we obtain , we push the line of integration further down, thereby re-deriving the put-call parity relation (4.14) with ω ∈ (λ − , −1). Similarly, the integrand in (2.22) has a pole of order s at zero. If we push the line of integration in (2.22) down and cross the pole, we obtain for any ω ∈ (λ − , 0).
For the choice of M c , one can apply (4.12) to find the positive integer M c such that M c ζ ≥ Λ.

4.1.4.
Error estimate of (V 2 , ·) L 1 and choices of ω, ζ, M c . Consider the case when flat iFT method is used to calculate V 2 (x ) (x = x + µ∆). Denote by d (V 2 , x ) and t (V 2 , x ) the discretization and truncation errors of the calculation of V 2 (x ). Since it suffices to consider each term on the RHS above separately. Let where ω, d(ω) and D 2 are as in Proposition 4.4, and L t denote the RHS in the upper bounds (4.12) when x = x − µ∆ = −µ∆.
Therefore, (4.20) and (4.21) can be applied for any α ∈ [1, α 0 ). The error estimates and recommendations below are valid for these α. The method is labeled parabolic iFT of order α.

4.2.2.
Errors: preliminary. The upper bounds for the truncation error and a procedure for the choice of Λ = M c ζ are given in Propositions 5.2 and 5.8 in [15]. Below, we consider the error control of the numerical scheme of V 2 . The upper bounds for the truncation error and a procedure for the choice of Λ = M c ζ are given in the following propositions, which are analogues of Proposition 5.2 and 5.8 in [15]. Proposition 4.7. Let ν ∈ (0, 2), ν = 1, and assume that either x > 0 and α ∈ [1, min{1 + 1/ν, 3}), or x = 0 and α ∈ [1, min{1 + 1/ν, 4}). Then a) As η = Re η → ±∞, the real part of the integrand in (4.20) admits an upper bound via (This recommendation is applicable only in the region of sufficiently large Λ, where the LHS is monotone or very close to a monotone function.) Proposition 4.8. Let ν ∈ (0, 2), ν = 1, and assume that either x < 0 and α ∈ [1, min{1 + 1/ν, 3}), or x = 0 and α ∈ [1, min{1 + 1/ν, 4}). Then a) As η = Re η → ±∞, the real part of the integrand in (4.20) admits an upper bound via b) Given an error tolerance for the truncation error, Λ = M c ζ can be chosen as a number, which satisfies (This recommendation is applicable only in the region of sufficiently large Λ, where the LHS is monotone or very close to a monotone function.)

4.2.3.
Choice of ω, ζ, M c . One may use the same length and mesh of the η-grid in order to apply the vectorization procedure in MATLAB. Since the discretization error decays exponentially, we expect that the C++ implementation with the x -dependent truncation parameter will be more efficient.

4.2.4.
Error estimate of (V 2 , ·) L 1 and choices of ω, ζ, M c . Consider the case when parabolic iFT method is used to calculate V 2 (x ) (x = x+µ∆). By analogy with [15], we use the same bound for d (V 2 , ·) L 1 as for flat iFT. One can derive the Hardy norm of the integrand in (4.24) and (4.25), hence, derive accurate error estimate for the discretization and truncation error. However, in our case, we only need simple recommendation for the choice of ω, ζ and Λ. Using (4.18), together with (4.4), a trivial modification of the algorithm in Section 4.2.3 gives the choices of ω, ζ. Set 1 = 0.9 · ∆ 2 · c /4 · π/(2D 2 ). Let where B + and B − are as in (4.28) and (4.32), respectively. Lemma 4.9.
it suffices to consider each term on the RHS above separately. We consider only the case x ≥ 0, the case x ≤ 0 are proved similarly. By Proposition 4.7, the truncation error admits an bound via where A + and L + t are as in (4.27) and (4.38). Integrating e x A + (Λ + )Λ + α w.r.t. x , we obtain Using the above lemma, we find Λ ± such that e x A + (Λ + )Λ α where 1 = 0.9 · ∆ 2 · c /4 · π/(2D 2 ). Positive integer M c can be found from M ± c ζ ± ≥ Λ ± .

Numerical algorithm and examples
5.1. Algorithm. We present an explicit algorithm for computing the prices of Asian put and call option. The parabolic iFT method with the choice of numerical parameters described in Section 4.2.1 and 4.2.3 can be used to replace the flat iFT and (refined) FFT method. The inputs are the spot price S 0 , the strike K, the maturity date T , the number of equally spaced sampling dates N + 1, the interest rate r, the parameters of the process, and the error tolerance . If (N + 1)K − S 0 ≤ 0, then the price of the Asian put option is 0. In the algorithm below, we assume (N + 1)K − S 0 > 0, and let V 1 and V 2 be as in (2.16) and (2.22), respectively.
Step I. Choose the error tolerances tr , int and c that will control, respectively, the truncation error, interpolation error, and the impact of the errors of calculation of values of functions V 2 and V 1 .
Let h = e −rT ((N + 1)K − S 0 )/(N + 1), and set tr = tr /h, int = int /h and c = c /h as the error tolerance for the calculation of V N . Set∆ = T /N .
Step II. Set 1 = tr /(N − 1). Using (3.5) and (3.6), find x 1 < x M < 0 so that the errors of the truncation above x M and partial truncation below x 1 do not exceed 1 .
Step III. Choose the order of the interpolation procedure, which (approximately) maximizes mesh ∆ of the grid in the state space given int /2 (recall that we do the interpolation at each step twice),∆ and parameters of the process. Below, we assume that the piece-wise linear interpolation is used. One can choose ∆ so that the RHS of (3.12) does not exceed int /2.
Step IV. Choose a smallest integer M such that (M −1)∆ > x M −x 1 , and re-define , and construct grids: , and z = ∆.
Step VII. Calculate U ≈ W 2 ( x) using (2.23) and fast convolution algorithm. The inputs are u, V 1 , and V 2 . Let W = U .
Step VIII. In the cycle w.r.t. n = 2, 3, . . . , N − 1, • calculate and store array V y ≈ W n ( y) using W as the input and applying the piece-wise linear interpolation procedure; • set u := (1 − e x ) V y and re-set u 1 = c n e x 1 and u M = 0, where c n is given by (2.20); • calculate W ≈ W n+1 ( x) using (2.23), inputs: u, V 2 and V 1 , and fast convolution algorithm; Step IX.
piecewise linear interpolation procedure, and then the Asian put option value using (2.3) and the Asian call option value using (2.2).

Numerical examples.
The results presented below were performed in MATLAB R 7.11.0 (R2010b), on a laptop with characteristics Intel R Celeron R Processor T1600 (1MB Cache, 1.66GHz, 667MHz FSB), under the Genuine Windows Vista TM Home Basic with Service Pack 2 (32-bit) operating system.
In all the examples, the benchmark prices are calculated using our method with very long and fine grid, both in the state space and the frequency domain. The set of numerical parameters guarantee a small error not larger than 10 −10 . [30]. In this subsection , we compare the performance of our algorithm with the performance of MC method and the method developed by Fusai and Meucci [30], for calculating the prices of discretely monitored Asian call options. As in [30], we assume that under a chosen EMM, the log-spot price, X t = ln S t , of the underlying follows a KoBoL process (see (2.12)) with parameters ν = 1.2945, c = 0.0244, λ + = 0.0765, and λ − = −7.5515. Assume the interest rate is r = 0.0367, which allows us to find the remaining parameter µ ≈ 0.138736 from the EMM condition r + ψ(−i) = 0 (where ψ(ξ) is the characteristic exponent of {X t }). The process is of infinite variation since order ν = 1.2945 ≥ 1.

KoBoL: example taken from
We calculate the prices of a discretely sampled Asian call option on the stock S t = e Xt , with spot price S 0 = 100, maturity date T = 1 year and the number of sampling dates N = 12, 50 and 250, respectively.
The results of our calculation are summarized in Table 1. The numerical parameters of our algorithm are specified in the caption to the table. (We use the acronym "MC", "LX(f)" and "FM" to label the results obtained using our method (implemented with flat (refined) iFFT) and the method in [30].) It is well known that the convergence of the MC estimator is very slow, therefore, we used 500, 000 paths. For simulating trajectories of KoBoL processes, we implemented the code 2 of Poirot and Tankov [39]. It is also well known that for processes of infinite variation, the accuracy of the results significantly affected by the simulation bias arise from truncating the small jumps, therefore, we truncate the jump with size less than 1 × 10 −9 . From the table, we observe that MC produce results with relative error of order up to 10 −3 .
For our results in column "LX(f)", we fix a moderately small mesh in the frequency domain, and change the mesh in the state space to see the change of the relative error. Due to the Nyquist relation M ∆ζ = 2π, halving ∆ increases the length of the grid in the frequency domain. One can see that our method produce more accurate results than the method of "MC" and "FM", and faster. The initial length of the grid in the frequency domain depends on the number N of the sampling dates, therefore, we have to choose different M 2 in the refined FFT algorithm for different N . The reason is that as N increases, a longer grid is needed to ensure that the truncation error in the frequency domain is small. We observe that as ∆ is halved, the error decreases by a factor of 10 and more. The exception is the case K = 110,∆ = 1/50 and ∆ ≈ 0.027781, when the error decreases by a factor of 2. This is due to the fact that the truncation error in the frequency domain is rather large. [22]. In this subsection , we compare the performance of our algorithm with the performance of the method developed byČerný and Kyriakou [22]. We use three parameter sets for KoBoL model (see (2.12)) considered in [22]: For these parameter sets, we calculate the prices of a discretely sampled Asian call option on the stock S t = e Xt , with spot price S 0 = 100, maturity date T = 1 year and the number of sampling dates N = 50; the riskless rate r = 0.04. The results are summarized in Table 2. The numerical parameters of our algorithm are specified in the caption to the table. (We use the acronym "LX" and "CK" to label the results obtained using our method and the method in [22], respectively.) We observe that the prices increase as m 2 increases as it is to be expected.

KoBoL: the examples taken from
Our method takes much less CPU time than CK method to achieve the same level of accuracy. The CPU time in [22] was recorded on a relatively better computer than ours: Dell Latitude 620 Intel R Core TM 2 Duo Processor T7200 (4MB Cache, 2.00 GHz, 667 MHz FSB) and 2 GB RAM with MATLAB c R15.
For these sets of model and option parameters, the implementation of our approach with flat (refined) iFFT is faster than that with the parabolic iFT. However, for processes with small ν and/or small∆, the advantages of the use of parabolic iFT become significant.

5.2.3.
KoBoL: example with small ν. In this subsection, we compare the performance of flat (refined) iFFT with the performance of parabolic iFT. We take a KoBoL process with parameters ν = 0.2, c = 1.1136, λ + = 3, λ − = −10 from [15]. Assume the interest rate is r = 0.04, which allows us to find the remaining parameter µ ≈ 0.30403 from the EMM condition r + ψ(−i) = 0. The process is of finite variation since the order ν = 0.2 < 1.
For these parameter sets, we calculate the prices of a discretely sampled Asian call option on the stock S t = e Xt , with spot price S 0 = 100, maturity date T = 1 year and the number of sampling dates N = 12. The results are summarized in Table 3. The numerical parameters of our algorithm are specified in the caption to the table. (We use the acronym "LX(f)" and "LX(p)" to label the results obtained using flat (refined) iFFT and parabolic iFT, respectively.) For this process, even when∆ is not small, the parabolic iFT perform much better than flat (refined) iFFT.

BM.
In this subsection, we compare the performance of our method with the performance of the method developed in [30] and [22], for pricing discretely monitored Asian call option under the Brownian motion. The results are summarized in Table 4 and 5.
Since the probability density function of the increment behaves fairly regularly and has very thin tails, one does not need the FFT algorithm to enhance the efficiency and use a uniform grid. We first reduce the expectation (2.4) to the integration on the real line: where p is the normal probability density of the increments X∆ with the mean 3 µ∆ and variance σ 2∆ , and then, use Gaussian quadrature to evaluate the integral at y = x − ln(1 − e x ).
From our results, one can observe that Gaussian quadrature converge to the benchmark price very fast.

Conclusion
We introduced a new method for pricing discretely sample Asian options, and suggest efficient numerical realization. We calculated prices of call options for several sets of parameters in the KoBoL and Brownian motion models, and demonstrated that our method are both more accurate and efficient than the method developed by Fusai and Meucci [30], and the method developed byČerný and Kyriakou [22].
An additional advantage of our approach is that we can independently control errors of calculations in the state space and discrete space, when the inverse Fourier transform is calculated. We derive bounds for all sources of errors, with the prescriptions for the choice of the parameters of the numerical scheme (given error tolerance).
Our method has two main variants. In the first one, the standard simplified trapezoid rule is used, which allows to use the inverse fast Fourier transform (iFFT) to increase the speed of calculations. In the second one, the parabolic iFT method based on an appropriate contour deformations is used; the number of terms in the simplified trapezoid rule can be greatly decreased but iFFT is no longer applicable. For the sets of model parameters in [30] and [22], our approach allows one to achieve the absolute error of 10 −8 within 1 seconds using falt iFT; the implementation with the parabolic iFT needs 2 seconds. However, if the order of the process and/or the interval between two sampling dates are small, the parabolic iFT is significant faster than flat (refined) iFFT.
The bounds for the derivatives in Lemma A.1 are simple, convenient, and, if the number of the sampling dates is not too large, the constants in the bounds are moderate. However, even if we use the piece-wise linear interpolation, we need estimates for the second derivatives of V 1 and W n . If we use the piece-wise cubic interpolation or cubic splines, we need bounds for the derivatives of order 4. We can derive bounds for the derivatives of V 1 and W n of order s + 1 > 1 in terms of the L 1 -norm of the derivatives p (s) ∆ of the transition kernel. Lemma A.2. Let g be continuous with a measurable bounded derivative. Then Proof. For processes under consideration, the transition kernel is infinitely differentiable on R, with the exception of 0 for VG, therefore, For all s such that p (s) ∆ L 1 < ∞, and all n ≥ 1, Proof. Let x < 0. Change the variable y = x − ln(1 − e x ). Then we have the following sequence of equivalent equalities and therefore, where u(y) = −e y /(1 + e y ). Direct calculations show that By applying Lemma A.4 to g = V 1 and g = W n and taking into account Lemma A.3, we obtain Lemma 3.4.

Appendix B. Interpolation of higher order: cubic spline
We want to evaluate the action of Fourier transform to a measurable function f on R: Given the function e Im ξ·x f (x) ∈ L 1 (R), the right hand side of (B.1) converge absolutely.
Typically, the function f has a kink or point of discontinuity, in such cases, it is advantageous to take into account this kink, and consider two sets of piecewise smooth polynomial approximation.
In the evaluation of (B.1), both for a smoother approximation and for a more efficient approximation, one has to go to piecewise polynomial approximations with higher order pieces. In this section, we describe the scheme for the piecewise cubic spline interpolation.
We assume that we are given the values of f on a uniformly spaced grid The coefficients a k , b k , c k and d k are determined from the following continuity conditions: , and S k (x k+1 ) = S k+1 (x k+1 ); and two initial conditions: S 1 (x 1 ) = s 1 and S M (x M ) = s M . Let us write S k (x k ) = s k , k = 1, . . . , M . Next, we replace the function f (x) that appears in the integrand in (B.1) with the cubic spline functions we just described. This replaces (F x→ξ f )(ξ) with a sum of integrals of the form A direct calculation ultimately leads to the following approximation: By providing the slops at the end points, the slopes s k , k = 2, . . . , M − 1 can be obtained by solving a tridiagonal linear system of equations: where x 0 = x 1 − ∆ and x M +1 = x 1 + ∆, V 1 is as in (2.16), and V s is as in (2.22).

Appendix C. Alternative calculations
We outline several possible directions in which the method of Section 2 can be developed further. In the following setup, since (i)FFT and fast convolution algorithm are not applicable, and parabolic iFT is more efficient and accurate than flat iFT for point-wise calculation, we would recommend to use parabolic iFT when the calculations in the dual space are needed. For the choice of parameters of parabolic iFT, see Section 4.2.3. C.1. Approach 1. In (2.6)-(2.7), we have to evaluate V 1 in (2.16) and W n in (2.23) at points of the form y k := x k − ln(1 − e x k ) + , where x k are points of an equally spaced grid. By taking into account this feature, we re-arrange the algorithm in Section 2.4.3 as follows: 1. calculate V 1 ( ∆ − ln(1 − e x +1 )) for = 0, 1, . . . , M −1, and V 2 on a M ×M matrix, with the entries (i, j) = (1 − i)∆ − ln(1 − e x j ), using parabolic iFT; 2. calculate E y k [u(X∆)], k = 1, 2, . . . , M 1 , using (2.23) and matrix multiplication.
Moreover, one can work with non-uniformly spaced grid of where V 2 is as in (2.22).
C.2. Approach 2. First, write the expectations (2.6) and (2.7) as the integrations on the real line: and, for n ≥ 2, where f 1 and f n are as in (2.8) and (2.9), respectively, and p is the the probability density of the increments X∆. Then, any standard quadrature can be used to evaluate the truncated integral at y = z − ln(1 − e z ) + , provided one can calculate f n and p on chosen grids. Explicitly, the integrals are calculated as follows. First, as in Section 2, we truncate the intervals of integrations in (C.2) and (C.3) from above at a point z M < 0, and use a partial truncation from below, that is, on (−∞, z 1 ], we replace f n (z) with c n e z , where c 1 = −e −∆ψ(−i) , and c n is given by (2.20) if n ≥ 2. Next, we choose a quadrature, construct a grid of points z in [z 1 , z M ], and evaluate f n on z. Set y = z − ln(1 − e z ), and use parabolic iFT to accurately evaluate p(z i − y j ). Finally, apply the quadrature to approximate the integral: where w are the weight of the chosen quadrature, and This scheme is especially useful in cases when the probability density function of the increment behave relatively regular and its tail is relatively thin ( Column 2 contains the benchmark prices obtained using our method. Column labeled "MC" contain the relative difference between the benchmark prices and the results obtained using Monte Carlo method. Column labeled "LX(f)" contain the relative difference between the benchmark prices and the results obtained using our method. Column labeled "FM" contains the relative difference between the benchmark prices and the results obtained using the method of Fusai and Meucci [30]. The results are taken from the tables in op. cit.. The example is taken from [30].  Column 2 contains the benchmark prices obtained using our method. Columns labeled "LX(f)" and "LX(p)" contain the relative difference between the benchmark prices and the results obtained using our method, implemented with flat (refined) iFFT and parabolic iFT, respectively. Column labeled "CK" contains the relative difference between the benchmark prices and the results obtained using the method ofČerný and Kyriakou [22]. The results are the same with our benchmark prices, they are taken from the tables in op. cit..  Column 2 contains the benchmark prices obtained using our method. Columns labeled "LX(f)" and "LX(p)" contain the relative difference between the benchmark prices and the results obtained using our method, implemented with flat (refined) iFFT and parabolic iFT, respectively. The example is taken from [15].  Column 2 contains the benchmark prices obtained using our method with the numerical parameters as follows. Column labeled "LX" contain the relative difference between the benchmark prices and the results obtained using our method. Column labeled "FM" contains the relative difference between the benchmark prices and the results obtained using the method of Fusai and Meucci [30]. The results are taken from the tables in op. cit.. The example is taken from [30].   Column 2 contains the benchmark prices obtained using our method with the numerical parameters as follows. Column labeled "LX" contain the relative difference between the benchmark prices and the results obtained using our method. For the method ofČerný and Kyriakou [22], the results are taken from the tables in op. cit. and are the same as our benchmark prices. The calculations of "LX" presented were performed in MATLAB c 7.11.0 (R2010b), on a laptop with characteristics Intel R Celeron R Processor T1600 (1MB Cache, 1.66GHz, 667MHz FSB) and 1 GB RAM, under the Genuine Windows Vista TM Home Basic with Service Pack 2 (32-bit) operating system. The calculation of "CK" presented were performed in MATLAB c R15, on a Dell Latitude 620 Intel R Core