A TWO-DIMENSIONAL , TWO-SIDED EULER INVERSION ALGORITHM WITH COMPUTABLE ERROR BOUNDS AND ITS FINANCIAL APPLICATIONS By

In this paper we propose an inversion algorithm with computable error bounds for two-dimensional, two-sided Laplace transforms. The algorithm consists of two discretization parameters and two truncation parameters. Based on the computable error bounds, we can select these parameters appropriately to achieve any desired accuracy. Hence this algorithm is particularly useful to provide benchmarks. In many cases, the error bounds decay quickly (e.g., exponentially), making the algorithm very efficient. We apply this algorithm to price exotic options such as spread options and barrier options under various asset pricing models as well as to evaluate the joint cumulative distribution functions of related state variables. The numerical examples indicate that the inversion algorithm is accurate, fast and easy to implement.

1. Introduction It is well-known that European option prices under the Black-Scholes model (BSM) can be computed directly via the celebrated Black-Scholes formula.However, closed-form solutions for European option prices are usually unavailable in more sophisticated asset pricing models.Carr and Madan [9] applied the fast Fourier transform (FFT) method to value European options under a wide range of models.Specifically, if the characteristic function of the asset return is known analytically, one can derive a closed-form, one-dimensional Fourier transform for the European option value with respect to (w.r.t.) the logarithm of the strike price.Then the FFT can be used to invert this one-dimensional Fourier transform to produce a numerical European option price.For more related literature, we refer to Carr and Madan [9,10] and references therein.
However, as pointed out in [10], the plain FFT breaks down for deep out of the money options and even generates negative values.This problem becomes more serious for model calibration which requires us to value European options for a wide spectrum of strikes repeatedly.Indeed, even in the cases when Carr and Madan's FFT breaks down, its accuracy can be improved by making different choices of the algorithm parameters.For example, Lee [21] studied how to improve the accuracy of one-dimensional Fourier transform method by conducting an error analysis.Alternatively, in the framework of Euler inversion algorithms, Cai et al. [7] generalized the elegant Euler inversion algorithm for one-sided Laplace transforms (Abate and Whitt [1]) to the two-sided case, and the computable error bounds of the resulting algorithm enable us to select the parameters suitably to achieve desired accuracy.
It is worth noting that these inversion algorithms are focused on onedimensional transforms, which are available in closed form for European options in a wide range of models.Nonetheless, for many exotic options even under simple models, there exist only closed-form solutions to twodimensional (two-sided) Laplace transforms rather than one-dimensional (one-sided) ones; e.g., barrier options and step options under the doubleand mixed-exponential jump diffusion models, and spread options under the variance gamma model and even under the BSM.
In this paper, we propose an inversion algorithm with computable error bounds for two-dimensional, two-sided Laplace transforms.Although an important motivation of ours is to price exotic options in financial engineering, the developed algorithm is much more widely applicable and can be potentially used to compute quantities of interest in general contexts with known two-dimensional Laplace transforms, e.g., the time-dependent probability distributions in queueing models (see Choudhury et al. [11]), the operational solutions of fractional diffusion equations (see Valkó and Abate [26]), and the first passage time distributions for spectrally one-sided Lévy processes (see Rogers [24]).The algorithm has three appealing features.
(i) The Laplace inversion formula is very simple to implement (see (5)), involving four parameters-two discretization parameters C 1 and C 2 and two truncation parameters N 1 and N 2 .(ii) Both discretization and truncation errors have explicit expressions.As a result, under mild conditions we can derive computable bounds for both errors in terms of parameters C 1 , C 2 , N 1 and N 2 .By choosing these parameters appropriately based on computable error bounds, the algorithm can achieve any desired accuracy, and therefore is especially useful to provide benchmarks.(iii) In many cases, the discretization error bound decays exponentially, while the truncation error bound decays exponentially or in a power law, leading to fast computation.
Our algorithm essentially extends the one-dimensional, two-sided Laplace inversion algorithm in [7] to the two-dimensional case.Another key difference between these two papers is that we provide a unified approach to finding the discretization error bounds based on the "known" Laplace transforms rather than the "unknown" original functions; see Section 3.2.1 and Proposition 3.2.This technique can also be applied to the one-dimensional case.
Our algorithm can also be viewed as an extension of the multidimensional, one-sided Euler inversion algorithm in [11] to the two-sided case.Petrella [23] proposed a similar algorithm with a scaling factor but without rigorous justification.As pointed out in [7], Petrella's method imposed a constraint on the scaling parameter, which may cause large errors even in the onedimensional case; see Section 4.3 in [7].Furthermore, Petrella did not analyze the truncation errors.Lee [21] conducted an elegant error analysis for the one-dimensional Fourier transform method.By contrast, our paper deals with the so-called Euler inversion algorithm in the two-dimensional case, where the inversion formula and hence the expressions of discretization and truncation errors are different.Therefore, the error analysis especially for discretization errors is also different.
In fact, our algorithm essentially belongs to the category of Euler inversion algorithms which have enjoyed great popularity in the areas of operations research and applied probability (see, e.g., [1,11,23,7]).Motivated by applications in financial engineering, e.g., path-dependent option pricing, we extend the existing Euler inversion algorithms to the two-dimensional, twosided case.An attractive feature of our algorithm is that in many cases, we can derive computable bounds for both discretization and truncation errors (please see .As a result, the numerical outputs of our algorithm consist of not only the inversion result but also the two error bounds (please see the numerical results in Section 8).These error bounds can tell us how accurate our inversion results are.This is useful especially in the absence of reliable benchmarks.Furthermore, by choosing the algorithm parameters appropriately based on the computable error bounds, we can achieve any desired accuracy.Therefore, our algorithm is particularly useful to provide benchmarks.
Besides Fourier and Laplace transforms, other transform-based methods applied in financial engineering include the Hilbert transform in Feng and Linetsky [13,14], the sinc expansion in Feng and Lin [12], and the fast Gaussian transform in Broadie and Yamamoto [2], to name just a few.These algorithms serve different purposes.For example, the Hilbert transform method is a powerful instrument to calculate expectations involving indicator functions.The sinc expansion approach along with the Hilbert transform is very accurate, efficient, and robust in evaluating the cumulative distribution function (cdf) of the asset return and in pricing European options under exponential Lévy models.The fast Gaussian transform algorithm is extremely fast and especially useful when pricing discrete barrier and lookback options under models with Gaussian or mixed-Gaussian returns.
The remainder of this paper is organized as follows.Section 2 presents the two-dimensional, two-sided Laplace inversion formula.We investigate the bounds for discretization and truncation errors in Sections 3 and 4, respectively.The financial applications of our inversion algorithm are discussed in Sections 5-7, and numerical results are given in Section 8.All the proofs are deferred to the appendices.

2.
The main result Before presenting our main result, we introduce several notations used throughout the paper.x denotes a two-dimensional row vector (x 1 , x 2 ) ∈ R 2 or C 2 .For any x and y, their inner product and Hadamard product are expressed by "•" and "•", respectively, i.e., For ease of exposition, we divide the plane R 2 into the following four disjoint regions w.r.t.M = (M 1 , M 2 ); see Figure 1 for an illustration.
2.1.Two-dimensional, two-sided Laplace transforms Consider a function f ( t) defined for t ≡ (t 1 , t 2 ) ∈ R 2 .Its two-sided Laplace transform or bilateral Laplace transform L f ( s) is defined as , where v ≡ ℜ( s) and ω ≡ ℑ( s) are real and imaginary parts of s, respectively.We define the region of absolute convergence (ROAC) of the two-sided Laplace transform to be the interior of the following set It is easy to see that L f ( s) in ( 1) is well defined for any s ∈ROAC.
Note that the ROAC does not depend upon ℑ( s).Therefore, we can simply use the range of ℜ( s) to represent the ROAC throughout the paper.
We point out that both two-dimensional, one-sided Laplace transforms and two-dimensional Fourier transforms are special cases of two-dimensional, two-sided Laplace transforms.Indeed, if f ( t) = 0 for any t ∈ [0, +∞) × [0, +∞), then its two-sided Laplace transform is reduced to the one-sided one If we confine our attention to L f ( s) only for s with ℜ( s) = 0, this results in the Fourier transform F f ( ω) of the function f ( t) It is worth noting that the ROACs of certain two-sided Laplace transforms may not include the imaginary axis.Thus the corresponding Fourier transforms might not be well defined.

2.2.
The two-dimensional, two-sided Laplace inversion formula Under mild conditions, we shall derive a two-dimensional, two-sided Laplace inversion formula that involves the parameters C = (C 1 , C 2 ) and N = (N 1 , N 2 ) for the purpose of controlling discretization and truncation errors, respectively.
Assumption 2.1.The function f ( t) is continuous on R 2 , and for any fixed v ∈ ROAC of L f (•), the function g( t) := e − v• t f ( t) has the following upper bound: where κ, c 1 , and c 2 are positive constants independent of t.
(II) e T ( t, v, C, N ) and e D ( t, v, C) denote the truncation error and the discretization error of the approximation f A ( t, v, C, N ) to f ( t), respectively.
Proof.See Appendix A.
Remark 2.1.Assumption 2.1 and 2.2 are sufficient but not necessary for the inversion formula (4) to hold.Indeed, they are used to justify a twodimensional Poisson summation formula in the proof.Other conditions that can validate the Poisson summation formula will also lead to (4).We adopt Assumption 2.1 and 2.2 mainly because they are mild and easy to check in many financial applications.
Remark 2.2.The inversion algorithm is easy to implement in that the inversion formula ( 5) is very simple.Furthermore, both discretization and truncation errors have closed-form expressions which facilitate respective error controls.As one shall see later, we can usually control them to achieve any desired accuracy by choosing sufficiently large discretization parameter C and truncation parameter N .

Discretization errors
3.1.Necessity of introducing the discretization parameter C First of all, we point out that to invert the two-sided Laplace transforms, the introduction of the discretization parameter C is necessary.Without introducing C (i.e., if C = 0), the discretization error can be quite large no matter what v ∈ ROAC is selected; see the example below.

Exponential decay of discretization errors
Introducing the discretization parameter C is not only necessary but can also make the discretization error achieve an exponential decay, leading to a fast computation for the inversion algorithm.When implementing the numerical inversion in practice, we usually first choose a closed rectangle Without loss of generality, l * j and u * j for j = 1, 2 are assumed to be nonzero.The following theorem shows that for any v ∈ (l * 1 , u * 1 ) × (l * 2 , u * 2 ), the discretization error decays exponentially as C increases.
, the discretization error has the following bound where d j = min{v j − l * j , u * j − v j } for j=1,2, Proof.See Appendix B.
Remark 3.1.As interpreted in Section 3.2.1 and illustrated in Sections 5-7, in many cases δ( y) in ( 8) can be specified explicitly.As a result, the upper bound (9) of the discretization error is computable, and hence we can control the discretization error simply by choosing sufficiently large C 1 and C 2 .

The choice of [l
affects the selection of C and N in a complicated way.Indeed, from (9), we can see that it affects the selection of C through ρ( v, t), ρ 1 (v 1 , t 1 ), ρ 2 (v 2 , t 2 ) (where the values of the function δ( y) at the four corner points of are involved) and d j for j = 1, 2. However, δ( y) is specified case by case and could be selected in different ways even in the same case.Therefore, the effect of on the selection of C also depends on the selection of δ( y) and might need to be analyzed case by case.Besides, since affects the selection of C, (14) and (16) in the case of (14) 2 in the case of (16).Nonetheless, the functions ζ j ( v, ω j ) and the values of α j and ξ j are different if the quantities to compute and/or the models are different.Accordingly, the effect of on the selection of N might also need to be analyzed case by case.
However, it is worth pointing out that although the choice of affects the selection of C and N , as long as it is selected from the ROAC, we can always use the discretization and truncation error bounds to achieve the desired accuracy by choosing sufficiently large C and N .Moreover, numerical examples suggest that the resulting algorithms are accurate and fast.
3.2.1.Specification of the function δ( y) To guarantee that Theorem 3.1 holds, we need to specify the function δ( y) in (8).However, since this depends on the information about the unknown original function f ( t), it does not seem easy.Sometimes we can specify δ( y) based on probability meanings of the unknown function f ( t) in financial applications.This is exactly what is used in Cai et al. [7] in the one-dimensional case.Nonetheless, this problem seems challenging in general cases.
To overcome this difficulty, in Proposition 3.2 we propose a general approach to specifying δ( y) based on the known Laplace transform then it also satisfies (8).

Proof. See Appendix B.
According to Proposition 3.2, as long as we can find a function δ( y) that satisfies (10), then it also satisfies (8).Since L f ( y + i ω) is assumed to be given in closed form, it is usually easy to explicitly specify δ( y) that satisfies (10) and hence (8).Indeed, if we can show that L f ( y + i ω) satisfies certain asymptotic conditions (13) or (15) and moreover, if we can find L( y), L 1 ( y) and L 2 ( y) satisfying (these two conditions are satisfied in many applications; see Sections 5-7) for any y ∈ ROAC, then some simple algebra yields that for any y ∈ ROAC, where U ( y) is defined as ( 12) Then applying Proposition 3.2, we can specify δ( y) := π 2 U ( y), which satisfies (8).As a by-product of Proposition 3.2, we have the following proposition to verify the condition (3) in Assumption 2.1.Proposition 3.3.If there exists a function δ( v) that satisfies (8) for any v in ROAC of L f ( s), then there exist κ, c 1 and c 2 such that (3) holds.
Proof.See Appendix B.

Truncation errors
In addition to the discretization error e D ( t, v, C), our inversion formula (4) has another error source-the truncation error e T ( t, v, C, N ).In Theorem 4.1, we derive truncation error bounds when the Laplace transform satisfies certain asymptotic conditions.It turns out that in many applications, these asymptotic conditions are satisfied and the involved parameters can be specified explicitly; see Sections 5-7.This then leads to computable truncation error bounds, based on which we can control the truncation error simply by choosing sufficiently large N 1 and N 2 .

Proof. See Appendix C.
Since in many cases, both discretization error bounds in Theorem 3.1 and truncation error bounds in Theorem 4.1 can be computed explicitly, we can control these two errors to achieve any desired accuracy, say 10 −n with n ∈ N, by selecting the discretization parameter C and the truncation parameter N in the following two steps.
such that the discretization error is no greater than 0.5 × 10 −n ; • Step 2. For fixed C chosen in Step 1, based on Theorem 4.1, select sufficiently large N = (N 1 , N 2 ) such that the truncation error is no greater than 0.5 × 10 −n .

Application I: Spread options
The spread option, whose payoff depends on the difference between two market variables, has been traded actively in various financial markets, including the equity derivatives market (index spread options), the commodity market (crack spread and crush spread options), and the energy market (spark spread options).However, no closed-form solutions are available for the spread option prices even under the BSM.But their two-dimensional Laplace transforms have explicit expressions under quite general models (Hurd and Zhou [17]).This section Table 1 Connections among some Laplace transforms.Here LSpr( s) is derived based on Hurd and Zhou [17] and B(x, y) denotes the beta function for x and y with ℜ(x) > 0 and ℜ(y) > 0

Original functions
Laplace transforms ROACs The joint pdf f (•) demonstrates that our two-dimensional, two-sided Laplace inversion algorithm can generate highly accurate numerical prices for spread options very fast as well as provide the related error bounds.

Laplace transforms
The payoff function of the spread option is typically defined as (S 1 (T ) − S 2 (T ) − K) + , where K > 0 is the strike price, T is the maturity, and for j = 1, 2, are prices of two correlated assets with respective return processes {X j (t)} and initial prices e x j (0) .Assume ES j (T ) < +∞ for j = 1, 2, and that the two-dimensional, two-sided Laplace transform of the joint pdf f ( u) of X(T ) := (X 1 (T ), X 2 (T )) for u ∈ R 2 is given by Then Table 1 provides the two-dimensional Laplace transform of its joint cdf as well as that of the following spread option price w.r.t.u := (x 1 (0) − log K, x 2 (0) − log K) under the given risk neutral measure P The following proposition shows that under mild conditions our inversion formula (4) applies to evaluate the spread option prices and the joint cdf of the asset returns.
Proof.See Appendix D.

Discretization errors
To compute the discretization error bound in (9), it is required to specify the function δ( y) in ( 8) explicitly.As interpreted in Section 3.2.1, if for any fixed y ∈ ROAC, we can show L f ( y + i ω) satisfies ( 13) or (15) and moreover, if we can find L( y), L 1 ( y) and L 2 ( y) satisfying (11) (it is usually easy to find them because L f ( y + i ω) is assumed to have a closed-form expression; see Section 5.4), then applying Proposition 3.2 and according to Table 1, we can specify respective δ( y) for F (•) and Spr(•) immediately; see Table 2, where U ( y) is defined in (12).

Two examples
The BSM.Consider a two-dimensional BSM with risk-free interest rate r, volatilities σ 1 , σ 2 , dividends q 1 , q 2 , and correlation ).The asset returns X 1 (T ) and X 2 (T ) have a joint normal distribution and the Laplace transform of its joint pdf f ( u) is given by Then the Laplace transforms of the joint cdf of X 1 (T ) and X 2 (T ) and the spread option price can be obtained immediately from Table 1.According to Proposition 5.1 and the following Proposition 5.2, our inversion formula (4) can be used to evaluate them by inverting their transforms.
The truncation error bounds related to the computation of the joint cdf and the spread option price can be derived by Theorem 4.1 with the parameters given in the following proposition.
Proposition 5.2.The two-sided Laplace transform L f ( s) under the two-dimensional BSM satisfies (15) with the parameters The Laplace transforms L F ( s) and L Spr ( s) also satisfy (15) with the parameters given by Table 3.
Proof.The proof follows immediately by noting that The discretization error bounds related to the computation of the joint cdf and the spread option price can be obtained by Theorem 3.1 and Table 2, where The Variance Gamma (VG) Model.Under the two-dimensional VG Model, the Laplace transform of the joint pdf f ( u) of the asset returns X 1 (T ) and X 2 (T ) has the following closed-form expression (see [17]) Like the BSM, from Table 1 we can obtain the Laplace transforms of the joint cdf of X 1 (T ) and X 2 (T ) and the spread option price immediately.Proposition 5.1 and the following Proposition 5.3 indicate that our inversion formula (4) can be applied to evaluate them by inverting their transforms.Furthermore, the truncation error bounds for the joint cdf of X 1 (T ) and X 2 (T ) and the spread option price can also be derived by Theorem 4.1 with the parameters given in Proposition 5.3.
Proof.See Appendix D.
Besides, the discretization error bounds for the joint cdf and the spread option price can be obtained by Theorem 3.1 and Table 2, where we can simply select

Application II: Barrier options
The pricing of barrier options relies heavily upon the joint distribution of the terminal asset return and the first passage time.Accordingly, analytical solutions are usually unavailable under general asset pricing models.Under the double-exponential jump diffusion model (DEM, see Kou [18]), Kou et al. [20] derived a closed-form two-dimensional Laplace transform for the barrier option price w.r.t. the transformed strike and the maturity.This section will apply our inversion formula (4) to compute the barrier option price and the joint cdf of the terminal asset return and the first passage time under the DEM.The computable bounds for both the discretization and truncation errors will also be provided.
6.1.Laplace transforms In the DEM, the asset price S t under the risk neutral measure P is given by S t = S 0 e Xt with the asset return The Lévy exponent of this double exponential jump diffusion process X t is Denote the first passage time of {X t : t ≥ 0} to a flat barrier b by: where b > 0 and X τ b := lim sup t→+∞ X t on the set {τ b = +∞}.Kou and Wang [19] derived the Laplace transform of the joint pdf f τ b ,Xτ b of τ b and X τ b : for ℜ( s) ∈ (0, +∞) × (−η, +∞), where β 1,s 1 and β 2,s 1 are the two roots of the equation G(z) = s 1 with positive real parts.Consider an up-and-in call (UIC) barrier option with the barrier H > S 0 , the strike K, and the maturity T .Its price under the risk neutral measure P is given by where k := − log K and b := log(H/S 0 ) > 0. Besides, we are interested in the joint cdf of τ b and X T , or equivalently the following probability Kou et al. [20] derived closed-form two-dimensional Laplace transforms of the joint cdf F τ b ,X T ( t) w.r.t.t 1 and t 2 as well as that of the UIC barrier option price UIC(T, k) w.r.t.T and k; see Table 4.We point out that it is not easy to determine the ROAC of L f ( s).Kou and Wang [19] showed (0, +∞) × (−η, +∞) ⊂ ROAC of L f ( s), whereas Cai and Sun [8] proved that L f ( s) exists for a larger set, i.e., for The following proposition shows that our inversion formula (4) is applicable to the evaluation of the joint cdf F τ b ,X T (•) and the UIC barrier option price UIC(•).Proposition 6.1.Under the DEM, we have (i) both F τ b ,X T (•) and UIC(•) satisfy Assumption 2.1, and (ii) both L F ( s) and L UIC ( s) satisfy Assumption 2.2.
Proof.See Appendix F.

Discretization errors
To compute the discretization error bounds in (9) via Theorem 3.1, we are required to specify the function δ( y) in (8).
To this end, we can still apply Proposition 3.2 and the method in Section 3.2.1 by using the result in the following Proposition 6.2 in Section 6.3.Here we provide a more straightforward approach alternatively.For any v in the respective ROACs, Consequently, we can simply select δ( y) = e G(y 2 )T for F τ b ,X T (•) and δ( y) = S v 2 +1 0 for UIC(•).
Proof.See Appendix G.

Application III: Computing sensitivities
In Sections 5 and 6, we apply our inversion formula to evaluate certain joint cdfs and exotic option prices.This section will illustrate that it can also be used to compute sensitivities or greeks of options.

Table 5
Laplace transforms of two deltas ∆1( u) and ∆2( u) of the spread option

Functions
Laplace transforms ROACs ∆1( u) L∆ 1 ( s) = e −rT K B(−s 2 ,s 1 +s 2 ) We take two deltas of spread options as an example and use the same notations as in Section 5. Assume L f ( s) is the two-dimensional Laplace transform of the joint pdf of two underlying asset returns X 1 (T ) and X 2 (T ); see (17).Then the Laplace transform L Spr ( s) of the spread option price w.r.t.u ≡ (X 1 (0) − log K, X 2 (0) − log K) is given by (see Table 1) Differentiating ( 25) w.r.t.S 1 (0) and S 2 (0) respectively and interchanging derivatives and integrals based on Theorem A.12 in Schiff [25] yields Laplace transforms of the two deltas; see Table 5.
Under similar mild conditions as in Proposition 5.1, we can show that our inversion algorithm is valid for computing ∆ 1 ( u) and ∆ 2 ( u).
Proposition 7.1.(i) If ES j (T ) < +∞ for j = 1, 2, X(T ) ≡ (X 1 (T ), X 2 (T )) has a continuous distribution under the risk neutral measure P, and 0 belongs to the ROAC of L f ( s), then both of the functions ∆ 1 ( u) and ∆ 2 ( u) Proof.The argument of (ii) is straightforward by Table 7, and (i) can be shown in the same way as for Proposition 5.1.

Discretization errors
To compute discretization error bounds via Theorem 3.1, we derive the required functions δ( y) in (8) in a similar manner as in Section 5.2; see Table 6.

Table 6 Specification of δ( y) for computing discretization error bounds of two deltas. Here U ( y)
is defined in (12) δ( y) in ( 8) Table 7 Connections among the parameters in ( 13) and (15)

Truncation errors
To calculate truncation error bounds via Theorem 4.1, we need to specify the asymptotic behavior of L ∆ 1 ( s) and L ∆ 2 ( s).According to their explicit expressions, we obtain . Consequently, if the asymptotic behavior of L f ( s) satisfies ( 13) or (15) with explicit parameters, then L ∆ 1 ( s) and L ∆ 2 ( s) also satisfy (13) or (15) with parameters given in Table 7.

Numerical examples
In this section we provide numerical results for the quantities of interest discussed in Sections 5-7, by inverting their respective two-dimensional Laplace transforms via our inversion formula (4).The numerical results indicate that our algorithm is very accurate and fast.The major difference from most of other numerical methods is that our algorithm also generates discretization and truncation error bounds, which tell us how accurate the numerical results are.Indeed, we can control these two errors to any desired accuracy by selecting sufficiently large parameters C and N .Therefore, our algorithm is especially suitable to provide benchmarks.All numerical experiments here are conducted via MATLAB R2011a on a desktop computer with 2.85GB of RAM and an Intel Core i5-2500 Table 8 Pricing spread call options and evaluating the joint cdf P (X1(T ) ≤ x1, X2(T ) ≤ x2) under the two-dimensional BSM.The column "Cai & Shi" denotes numerical results obtained via our algorithm.The column "Hurd & Zhou" is taken from Table 1 of [17]."True values" are calculated from the analytical formula of the joint normal cdf.The model parameters are the same as in Table 1 of [17] 8 gives numerical results (denoted by "Cai & Shi") for the spread call options and the joint cdf of the asset returns under the two-dimensional BSM via our inversion formula (4) as well as the related discretization and truncation error bounds.We can see that our algorithm is highly accurate and efficient.(i) Our numerical results for the spread option prices agree with those obtained by Hurd and Zhou [17] (denoted by "Hurd & Zhou") to 6 decimal points.In fact, our algorithm is much more accurate than reflected from the comparison with Hurd and Zhou [17] because the error bounds indicate that our results have achieved the accuracy of 10 −9 ∼ 10 −14 .(ii) Our numerical results for the joint cdf agree with the true values to 10 decimal points.Besides, our algorithm is robust because even in the extreme case x 1 = −1 where the probability is very small, our algorithm is still highly accurate with reliable error bounds.(iii) It takes only 1 second and 0.03 seconds to generate one numerical result for the spread option price and the joint cdf, respectively, via our algorithm.
Similarly, Table 9 provides numerical results for the spread call option prices and the joint cdf of the asset returns under the two-dimensional VG model as well as the associated error bounds.It is indicated that our algorithm is still highly accurate and efficient.8.1.1.Pricing spread options in the "deep out of the money" case As pointed out by Carr and Madan [10], the plain FFT might break down for deep out of the money European call options and even generate negative values.In this section, we intend to apply our two-dimensional, twosided Laplace inversion method to price the spread call options in a similar "deep out of the money" case, where S 2 (0) + K is much larger than S 1 (0).Specifically, under both the two-dimensional BSM and the VG model, we set S 1 (0) = 100 and let S 2 (0) + K vary from 120 to 200 with increment 10.The numerical spread option prices are given in Table 10.We can see that our algorithm is still highly accurate even in the extreme cases, e.g., S 2 (0) + K = 180, 190 and 200.This implies that our algorithm is quite robust.8.1.2.The special case K = 0 under the BSM When K = 0, the spread option price under the two-dimensional BSM has an analytical formula; see (29) in Appendix E. Furthermore, its Laplace transform also becomes simpler in this special case.In fact, its single Laplace transform w.r.t.X 2 (0) has a closed-form expression given by (30) in Appendix E. Then the onedimensional, two-sided Euler inversion algorithm applies as a special case of our two-dimensional algorithm.Moreover, the related discretization and truncation error bounds can be derived similarly.See Appendix E for more details.Table 11 reports the comparison between our Laplace inversion results and the true values obtained from the analytical formula, and also provides the associated discretization and truncation error bounds.It can be seen that the inversion method is still highly accurate.

Table 9
Pricing spread call options and evaluating the joint cdf P (X1(T ) ≤ x1, X2(T ) ≤ x2) under the two-dimensional VG model.The column"Cai & Shi" denotes numerical results obtained via our algorithm.The column "Hurd & Zhou" is taken from Table 3 of [17]."MC values" and "Std.err." are Monte Carlo simulation estimates and associated standard errors, respectively, obtained using a sample size of 10 7 .The model parameters are the same as n Table 3 of

Table 10
Pricing spread options under the two-dimensional BSM and VG model in the "deep out of the money" case.The column"Cai & Shi" denotes numerical results obtained via our algorithm, while "MC values" and "Std.err." are Monte Carlo estimates and associated standard errors, respectively, obtained by simulating 10 7 samples.All the model parameters are the same as in Table 8 for the BSM and in Table 9 for the VG model.
The parameters for the inversion algorithm are v1 = 7, v2 = −2, C1 = C2 = 10, N1 = 400 and N2 = 600 for the BSM, and v1 = 7, v2 = −2, C1 = 7, C2 = 6, N1 = 300 and N2 = 600 for the VG model.It takes around 1 second and 0.6 seconds to produce one numerical result under the BSM and the VG model, respectively Pricing spread call options under the two-dimensional BSM in the "deep out of the money" case  12, we provide numerical results for the up-and-in call barrier option prices and the quantity P (τ b ≤ x 1 , X T ≥ −x 2 ) under the DEM as well as the related error bounds.Note that in this case, the truncation error is not as easy to control as for the spread options.This is primarily because the truncation error bounds decay more slowly in a power way; see Propo-Table 11 Pricing spread options when K = 0 under the two-dimensional BSM.The column "Cai & Shi" denotes numerical results obtained via Laplace inversion algorithm, while "True values" are computed from the analytical formula (29).The model parameters are the same as in Table 8 except that S2(0) varies from 80 to 120.The inversion algorithm parameters are v = −2, C = 2, and N = 80.It takes approximately 0.0003 seconds to generate one spread option price via our algorithm Pricing spread options under the two-dimensional BSM when K = 0  13 gives numerical results for two deltas of spread options under the BSM and VG model as well as the related error bounds.It can be seen that all the numerical results stay within the 95% confidence intervals of the associated Monte Carlo simulation estimates.Indeed, the error bounds indicate that the numerical results are more accurate than reflected from the comparison with the Monte Carlo estimates.Besides, the algorithm is very efficient in that it takes only 0.7 seconds to generate one numerical result.9. Conclusions This paper is devoted to the development of a twodimensional, two-sided Euler inversion algorithm with computable bounds for both discretization and truncation errors.This algorithm is especially useful to provide benchmarks in that the computable error bounds make it possible for us to select the algorithm parameters properly to achieve any desired accuracy.Since the error bounds decay quickly, e.g., exponentially, in many cases, the algorithm is very efficient.Numerical experiments of applying this algorithm to the valuation of some exotic options and joint cdfs suggest that the algorithm is accurate, fast, and simple to implement.
Z-transforms can be considered as the counterpart of Laplace transforms in the discrete case and also have many applications in operations research
Define g(t) := e − v• t f ( t).Then to prove (26) under the assumption t 1 t 2 = 0, it is equivalent to showing the following two-dimensional Poisson summation formula.
Introduce a new function We point out that g p ( 1 2 , 1 2 ) is exactly the LHS of (27).Now we shall prove that the RHS of ( 27) is exactly the Fourier series of g p ( 1 2 , 1 2 ), and moreover, it converges to g p ( 1 2 , 1 2 ).This can be obtained immediately by Proposition 3.3.2 in Grafakos [16] if we can verify the following three conditions: (a) g p ( x) is continuous at the point ( 12 , 1 2 ); (b) g p ( x) ∈ L 1 ([0, 1] × [0, 1]); and (c) the partial sum of the Fourier series of g p ( x) converges at the point ( 12 , 1 2 ).Let us verify (a), (b), and (c) one by one.
As for (a), we claim that g p ( x) is continuous on R 2 .Indeed, consider a sequence x (n) := (x 2 ) → x and assume | x (n) − x| < 1 without loss of generally.Then by (3) we obtain that for any l ∈ Z 2 , with l∈Z 2 g * ( l) < +∞.Applying the dominated convergence theorem yields g p ( x) = lim n→+∞ g p ( x (n) ).
The condition (b) is also satisfied because To show that (c) also holds, we calculate the Fourier coefficients of g p ( x) as follows.
where the second equality holds due to (28) and Fubini's Theorem.Accordingly, for any fixed x ∈ R 2 , the partial sum S N ( x) of the Fourier series of g p ( x) is given by As N 1 and N 2 tend to +∞, the double series S N ( x) is absolutely convergent because of Assumption 2.2.This implies that (c) also holds.
Similarly, it can be shown that Adding the eight inequalities above together completes the proof.
Proof of Proposition 3.2.By the Bromwich contour integral we obtain which justifies Proposition 3.2 immediately.
Proof of Proposition 3.3.Since the ROAC is defined to be the interior of the set (2), there exists ǫ > 0 such that 13), simple algebra yields Similarly, we can obtain Then ( 14) follows immediately.
(ii) Similarly to (i), if L f ( v + i ω) satisfies ( 15), we have where the second inequality holds because the function x −α j e −ρ j x ξ j is de- Thus ( 16) is proved.
APPENDIX D: PROOF OF PROPOSITION 5.1 AND PROPOSITION 5.3 Proof of Proposition 5.1.The proof of (ii) is straightforward according to Table 1.As for (i), F (•) is continuous because X(T ) has a continuous distribution.Besides, applying the dominated convergence theorem and using the fact that Ee X 1 (T ) < +∞ yields that Spr(•) is also continuous.Finally, it is straightforward and elementary to verify that F (•) and Spr(•) satisfy (3) thanks to the assumption that 0 belongs to the ROAC and the details are thus omitted.
APPENDIX E: THE SPECIAL CASE K = 0 OF SPREAD OPTIONS UNDER THE BSM When K = 0, the spread option under the two-dimensional BSM becomes analytically tractable.Indeed, if we define a new measure P d P dP = e −(r−q 2 )T S 2 (T ) S 2 (0) , then by change of measure we obtain a closed-form spread option price as follows •) denotes the standard normal cdf, and the second equality holds because S 1 (T )/S 2 (T ) is still log-normally distributed under the new measure P. In fact, (29) is essentially the Black-Scholes formula.
We claim that for the fixed x > 0 and any y ∈ R, the real parts of the four complex roots are all non-zero.Otherwise, there exist d, y ∈ R such that G(id) = x + iy.However, this contradicts to Note that the real parts of the four roots are all continuous in y.Moreover, by (35) we know ℜ(β j,x+iy ) > 0 and ℜ(γ j,x+iy ) < 0 for j = 1, 2 when x > 0 and y = 0. Therefore, this still holds for the fixed x > 0 and any y ∈ R, i.e., there exists exactly two roots with positive real parts.Now we begin to prove (34).For the first inequality in (34), define On the other hand, apparently f 1 (z) is analytic on K 1 .Moreover, since ℜ(z) > 0 and θ > 0, we have Thus Rouché's Theorem applies and implies that for any |y| > Y (x), the number of zeros of f 1 (z) + g 1 (z) in K 1 is the same as that of f 1 (z) in K 1 .Hence for any |y| > Y (x), G(z) = x + iy has exactly one root in K 1 because G(z) − (x + iy) ≡

Table 4
Connection among several transforms related to barrier options under the DEM for computing the truncation error bounds of two deltas.Here j = 1, 2 and p = 1, . . ., 4. Without loss of generality, we assume M2 ≥ |v2| + 1 here .03 seconds to generate one numerical result for the spread option price and the joint cdf, respectively, via our algorithm Evaluating spread options and related joint CDF in the BSM and VG Table Evaluating the joint cdf P (X1(T ) ≤ x1, X2(T ) ≤ x2) under the two-dimensional BSM 8 second and 0.05 seconds to generate one numerical result for the spread option price and the joint cdf, respectively, via our algorithm Evaluating barrier options and related joint CDF under the DEM In Table . err.UB of trunc.err.True values Abs.err.Computing two deltas of spread options in the BSM and VG model Table

Table 13
Evaluating two deltas ∆1( u) and ∆2( u) of spread options under the BSM and VG model.