1 Introduction
Laplace transforms are powerful tools employed to derive the analytical solutions of differential equations. However, it is too difficult to find or evaluate the analytical inverse transform in closed form. Researchers try to use numerical Laplace transform methods that convert the timedependent equations to parameterdependent problems, which could be parallelized easily in frequency domain. There exist a great deal of approaches that could be used to obtain the inverse transform over the past five decades. Davies and Martin in 1979 wrote a survey [13] investigating many promising methods, in which 14 specific algorithms were tested and compared. Duffy [15] concentrated on three methods developed later, in which two methods are based on existing techniques but considerably improved by other researchers. Cohen wrote a comprehensive review book [10] that shows all aspects of Laplace transform inversion. Finally, Kuhlman surveyed five different methods in [23] and their implementations.
There exist many timeparallel methods, including shootingtype methods, spacetime domain decomposition methods, spacetime multigrid methods, and direct methods [20]. Obviously, Laplace transform belongs to the direct methods where [12] and [36] are cited as pioneering works in this category. Our paper is based on the algorithm described in [12] that is called GaverStehfest algorithm, which has a good accuracy on a fairly wide range of functions as illustrated in [13], thereafter they developed an iterative scheme in [26, 25] through linearization of implied volatility, leading to a series of iterations throughout the timestepping process.
Our purpose in this paper is to employ asynchronous iterative scheme in modeling a Laplacetype solver, which might be extremely flexible in data transmission and exploitation throughout the processing phase. Such iterative scheme was first theoretically proposed in [9] for relaxation algorithms, formalized further in [34, 4, 17] with normbased contraction model and in [5, 6] with nested set model. Another wellknown timeparallel method, called Parareal, has been successfully generalized to the asynchronous iterative scheme [31, 29]. In the next section, we recall and present the classical GaverStehfest algorithm. In Section 3, we introduce the iterative Laplace transform method for an option pricing model with variational volatility. In Section 4, we propose a new Laplace method based on asynchronous iterative scheme with some remarks about the its convergence. Finally, some numerical experiments are given in Section 5 and concluding remarks are presented in Section 6.
2 The GaverStehfest algorithm
Given a secondorder linear elliptic operator , consider the initial value problem
(1) 
where the boundary conditions are supposed defined to have a wellposed problem. In the following we note if the space variable does not hamper the formulation and computation. To solve this equation, we employ the Laplace transform defined as
(2) 
A general contour integral formula for the inverse Laplace transform is given as
(3) 
where is the Bromwich contour that must be further determined, thus yielding a lot of numerical applications in recovering the original function.
We are interested in the GaverStehfest algorithm, which aims to approximate Equation (3) by a sequence of functions
(4) 
where are defined as follows
(5) 
where is an even number denoting the number of processors. This numerical solver could be deduced from Equation (3) by Cauchy integral formula with specific parameters. Notice that
will change sign during iterations. Such behavior is remarkable and indeed affects the numerical performance of GaverStehfest algorithm in some cases, which will be shown in the following sections, thus ignored for the moment.
The implementation of direct Laplace transform method is quite simple. Given the timedependent problem (1), making Laplace transform yields
where is the Laplace transform of depending on , and obviously
Finally, making inverse Laplace transform yields
which leads to Algorithm 1.
We notice that one should operate a reduction operation after Laplace transform, which is a commonly used term borrowed from messagepassing specification, where a designated root node receives data from all nodes and executes an arithmetic or logical operation. Here the specific operation is summation as indicated above. Furthermore, we illustrate the corresponding state diagram in Figure 1.
Note that Laplace transform method operates in the frequency domain, whereas most of other methods are proceeding in the original space, in which timedependent terms must be solved by appropriate temporal methods, such as Euler methods, RungeKutta methods, and multistep methods.
3 Iterative Laplace transform method for quasilinear problem
It is not desirable to solve quasilinear problems directly and thus we deal with them in a different direction. Throughout this paper, we take BlackScholes equation as an example to investigate the behavior of iterative methods (see, e.g., [26]).
We note that option pricing is a crucial target in financial decisionmaking. The breakthrough came with the appearance of the BlackScholes equation [7], which has a huge influence to the financial market and drives an unexpected prosperity in the trading of derivatives. This equation is often called the BlackScholesMerton (BSM) equation since it was further generalized by Merton in important ways [32, 33]. Consider the following European call option pricing problem
where is the option price, depending on stock price and time . As before, we note in the following if possible. Volatility and riskfree interest rate are the constant parameters, with the final and boundary conditions given by
(6) 
where is the maturity of the option and is the strike price. If there exist transaction costs, then the price prediction may become much more complex. The volatility can be treated in different ways, for instance using a modified volatility function . Here we assume that the price of option is a parameter of , thus leading to the BSM equation with implied volatility
(7) 
There are many ways to choose an appropriate . Here we adopt the function described in [26] as following
(8) 
where denotes constant historical volatility.
We should apply Laplace transform method to Equation (7), but this one is too complex to be employed directly. Moreover, from (6) we notice that (7) is a backward equation. Hence, we perform the following variable transformation
substituting into Equation (7) gives
(9) 
with the corresponding conditions
(10) 
Implied volatility defined in (8) becomes
(11) 
Notice that (9) is forward parabolic and the fractional term is flexible for both linear and quasilinear model. If we import a constant volatility, the fractional term will vanish; if we use implied volatility like the one defined (11), the constant coefficient therein will also disappear.
This equation, however, depends on a quasilinear coefficient that should be tackled before getting forward. We share the idea of [26] by introducing a onestep retard into , which linearizes the quasilinear equation (9) in the form
(12) 
We note that can be seen as a “frozen variable” which is assigned the previous value of and thus does not depend on . Now we prepare to apply the Laplace transform method to Equation (12). Setting
yields
(13) 
More specifically, this yields an iterative scheme illustrated in Algorithm 2. The same algorithm has been investigated in [26] under the name of “iterative coefficient–inverse Laplace transform”.
Since all processors need the values of at the beginning of each iteration, we should broadcast such buffer to all neighbors. Note that “broadcast” is a collective communication operation which involves a specific root node communicating its data to all other nodes. We observe that is computed outside because it depends only on the rank and the number of processors . Similarly, the state diagram is shown in Figure 2.
We mention here that there are large numbers of models dealing with the transaction cost, such as those in [8] and [3], some of which are nonlinear differential equations that have different forms and properties. There exist meanwhile many other types of equations related to the option pricing problem, such as the martingale pricing theory [21], the binomial options pricing model [11], and the stochastic volatility option models [22, 14]. They are completely different from the BlackScholes model and thus should be addressed separately. In general, one might employ Algorithm 2 to other types of equations generalized from Equation (1), such as
Substituting the frozen variable and applying Laplace transform yields
where linear operator can be solved by any appropriate discretization methods.
4 Asynchronous Laplace transform method
4.1 Formalization
Asynchronous iterative method releases the restriction of strict data dependency whereby iterations are executed by several processors in arbitrary order without any synchronization during computation, first formalized in [9] for linear systems. Obviously, it must be subject to some constraints that each processor still needs to continuously obtain the latest information. Since there is no synchronization imposed to such scheme, asynchronous iterative algorithm may exhibit incredible flexibility and efficiency even with the expansion of computational nodes, which contributes to overcome the fault tolerance and the load balancing problems.
Consider a vector space
withLet
Practically, denotes the number of processors. Hence
leading to the classical parallel iterative scheme
(14) 
which is presented in Figure 3.
In order to compute th iteration’s in processor , one needs to collect all data of iteration from other processors, which imposes a synchronization point on the computational framework at the end of each iteration.
Now we define the integer subsets such that
and let be nonnegative integervalued functions satisfying
(15) 
with . We introduce here an iterative scheme that generates the following sequence
(16) 
which is actually the mathematical model of asynchronous iterative scheme and illustrated in Figure 4.
Generally, for the sake of convergence analysis, Assumption 1 and Assumption 2 are armed therewith, which ensure that processors read eventually the latest information for each element from essential neighbors and no processor stops updating during iterations.
Assumption 1.
.
Assumption 2.
.
Now we apply asynchronous iterative scheme to the iterative Laplace transform method depicted in Algorithm 2. Since newer value depends on the previous value , thereafter we use and to denote these two values during iterations. let
(17) 
and for a fixed time span , we define
(18) 
Notice that
which leads to a handy operator
(19) 
Combining (17), (18), and (19) gives
Setting
(20) 
yields
which follows exactly the classical parallel iterative scheme (14). Hence, iterative Laplace method can be generalized naturally to the asynchronous iterative scheme subject to (15) and assumed to satisfy Assumption 1 and Assumption 2
(21) 
then leading to Algorithm 3.
The most remarkable notation therein is that exhibits the chaotic behavior of asynchronous iterative scheme, where intrinsically retard term and execution set play important roles during iterations. Meanwhile, another hinge consists in the sending and receiving operations without synchronization point, which yields the asynchronous performance, as opposed to the aforementioned blocking reduction and broadcast operations. Note that one may also employ nonblocking reduction and broadcast operations to implement Algorithm 3. The state diagram is depicted in Figure 5.
4.2 Remarks on convergence
To the best of our knowledge, unfortunately, there is no available theory that can be used to prove the convergence of asynchronous Laplace algorithm. According to the normbased convergence theory (see [17], which is based on the results in [34, 4]), one should prove a relationship in the form
where and denotes the weighted maximum norm with , in order to establish the convergence result. This imposes a contraction condition to the operator . We recall that in our case there exists a summation operator , as seen in (20), and thus is not an assembly but a summation of subvectors. It is clear that the summation operator could not lead to desired partial order under contraction. A more general but similar idea proposed in [5] (see also [6]) is based on the nested set theory, in which a technique known as “box condition” allows to obtain convergence result without using any norm. This can be written as
where for all . This condition involves a superposition of infinite subsets and thus suffers from the same problem as the previous approach. In addition, twostage and nonstationary models asynchronous models were discussed in [18, 19] and some papers made use of the partial ordering such as [16, 35]. These models have gained much popularity in some problems but do not correspond to our case.
The expansion of summation operator is indeed compensated by the frequency decomposition, which may contribute to the convergence. In what follows we will not pursue the convergence issue further; instead, we will focus on the numerical behavior of our method.
5 Numerical experiments
5.1 Environmental setting
We have conducted several experiments for the direct Laplace transform method, using the GaverStehfest algorithm as illustrated in Algorithm 1, as well as the iterative variants based on Algorithms 2 and 3. All tests are executed using (13) with in the linear case, whereas (11) is adopted in the quasilinear case.
All tests are launched on an SGI ICE X cluster connected with InfiniBand, involving MPI library to run parallel applications, which is supported by SGIMPT 2.14. Mathematical operations and linear systems solvers are implemented by the Alinea programming library [27] and iterative algorithms are programmed by JACK [28, 30] for both synchronous and asynchronous variants.
5.2 Direct scheme for linear equation
According to [37], number of processors must be even. We first illustrate the numerical results of Algorithm 1 applied to linear equation, shown in Figure 6.
Here we run programs with different numbers of processes and observe that error behaves like a convex function for each specific time span, where denoting the maturity of the specific option contract. We choose as convergence threshold, which is adequate and conspicuous for the financial data. Table 1 gives the convergence interval of direct method.
convergence interval of  remarks  

1  8, 10, 12, 14, 16, 18  : inaccurate; : inaccurate 
5  10, 12, 14, 16, 18  : inaccurate; : inaccurate 
10  14, 16, 18  : inaccurate; : inaccurate 
20  inaccurate, oscillating 
Here, convergence interval and convergence zone are defined as the set of that leads to convergent result. Notice that the larger we assign for time span, the narrower convergence interval we get. When , there is no convergence interval observed throughout experiments. In this case, we should compute results stepbystep, as described in [26]. Here “oscillating” indicates that we could not observe a contracting behavior when increases monotonously. On the contrary, results fall into an inconsistent solution range and oscillate divergently around initial point. We mention here that large , denoting the maturity with the unit of year, is not a normal test case in option pricing problem, which is shown only for the experimental purpose. In the following tests, we will use primarily small to illustrate the real behaviors.
Davies and Martin [13] mentioned that the Laplace transform method based on the GaverStehfest algorithm gives good accuracy on a fairly wide range of functions. Figure 6 and Table 1 show that a wellposed option pricing problem leads to good accuracy with a fairly wide range of convergence interval (e.g., lead to a convergence interval ). Theoretically, the more processors we use, the more precision we obtain through parallel processing [24]. In practice, however, we must consider the arithmetical precision of the test machine, namely, we decrease the truncation error thanks to the large number of processors and meanwhile increase the rounding error due to digit width limit. Unfortunately, the side effect is dramatically large due to the rapid growth of with sign alternating. Thus we could not exploit the power of large scale parallel computing, which acts as an intrinsic defect of the GaverStehfest algorithm.
5.3 Iterative scheme for quasilinear equation
Now we focus on the numerical behavior of iterative Laplace transform methods illustrated in Algorithms 2 and 3. We continue to use the convergence threshold as the termination condition, which may lead to an infinite loop in our program. Experimental results are shown in Tables 2 and 3.
convergence interval of  remarks  

0.01  6, 8, 10, 12  : inaccurate; : divergent 
0.1  4, 6, 8, 10, 12  : inaccurate; : divergent 
1  4, 6, 8, 10, 12  : inaccurate; : divergent 
convergence interval of  remarks  
0.01  6  : inaccurate; : divergent 
0.1  6  : inaccurate; : divergent 
1  4, 6  : inaccurate; : divergent 
We notice that the convergence interval of the asynchronous Laplace transform method is much narrower than the synchronous method. This is a reasonable compensation because asynchronous methods require less waiting time in communication but more computational iterations than classical parallel scheme, which amplify the unstable issue of with rapid growth and sign alternating.
On the other hand, for a large , both synchronous and asynchronous algorithms become divergent. This means that we could not obtain a result underneath the given threshold since the uncertainty propagates throughout iterations. Hence, such behavior is caused by an intrinsic issue of the GaverStehfest algorithm, which is the same as that in the direct method.
We illustrate in Table 4 the results of direct Laplace method with the same interval of as in Tables 2 and 3.
convergence interval of  remarks  

0.01  4, 6, 8, 10, 12, 14, 16, 18, 20, 22  : inaccurate; : inaccurate 
0.1  6, 8, 10, 12, 14, 16, 18, 20  : inaccurate; : inaccurate 
1  6, 8, 10, 12, 14, 16, 18  : inaccurate; : inaccurate 
One should keep in mind that the iterative results are obtained from the BSM equation with implied volatility, whereas the direct method, as seen in Algorithm 2, is applied to the linear equation where is a constant. Thus, comparing the former with the latter could not lead to an obvious conclusion for the performance of our method. We show these results only for clarity and completeness.
Figure 7 illustrates another experiment that divides temporal interval into several smaller parts, as described in [26].
As mentioned before, we address primarily the real situation in option pricing problem where is small. For a small we need no more subdivision in the time domain. Thus, we conduct this test only for evaluating the performance of the algorithm illustrated in [26], where we execute asynchronous Laplace transform method stepbystep along time domain. Similar to the classical temporal integration schemes like Euler methods, we perform computations in each step using always the latest data, but employ a rather wide subinterval to fit the characteristic of the Laplace transform method. The number of processors, the step size, and the number of steps are chosen as , , , respectively. We can see that the execution process in the first time interval is highly structured, except for some unordered data due to the asynchronous scheme. Nonetheless, the last time interval gives different behavior that seems somewhat affected by the unfinished chores due to the asynchronous iterative scheme. This produces larger retards and converges more quickly when receiving the latest data in some steps.
Finally, we test the accuracy of the asynchronous Laplace method in comparison with the synchronous version. We fix the volatility , riskfree rate , maturity , and number of processor . Recall that
The constant volatility and the strike price in the right hand side are fixed, whereas the volatility changes with . Then we change the stoke price and strike price to compare the two methods. Results are shown in Table 5.
Here, and are the option prices obtained from synchronous and asynchronous iterations, respectively. and act as absolute and relative errors such that
We notice that the asynchronous option prices are close enough to the synchronous cases. This observation achieves our expectations, as we speculated that the chaotic iterative process does not affect dramatically the final result. Since the synchronous Laplace transform method has proved to be accurate in predicting the option price values using the BSM equation with implied volatility, we can draw a conclusion that our method yields accurate results as well.
5.4 Further discussion
We have seen that Laplace transform methods are highly influenced by the number of digits of machine precision that are prone to roundoff error propagation, which leads to an unstable behavior. Notice that (5) increases quickly with commutative plusminus sign. The accuracy mounts in the beginning as the number of processes increases and then decreases rapidly.
One could avoid such effect for Algorithm 1 by simply adopting the recommended precision in a multiprecision computational environment as discussed in [1]. Meanwhile, the fixed Talbot method [38] requires the precision , which performs a less restrictive behavior and seems a better alternative to the GaverStehfest algorithm according to the experiments therein. Further investigation gives a unified framework for various types of inverse Laplace transform algorithms [2]. Recall that the contour integral is
Let , we have
Then taking
yields
Finally, using Cauchy integral formula, this becomes
(22) 
Choosing
where is an even integer, and notice that , we get again the GaverStehfest formula (4). We could further rewrite (22) in the form
where
Combining (17) and (19), we could obtain the asynchronous scheme (21) for the unified Laplace transform model, whereby the methods discussed in [2], such as Fourier series methods and Talbot method, might also be executed in asynchronous mode, although no rigorous proofs are available. The contour integral approach with quadrature is another choice and has proved to be very efficient for parabolic problems with timeindependent coefficients [36]. One may expect to combine these strategies to tackle the quasilinear equations in a more effective way.
These issues are beyond the scope of this paper. We leave this work, as well as a rigorous proof of the asynchronous convergence, as a direction of future investigation.
6 Concluding remarks
In this paper, we investigated the Laplace transform and proposed a new method based on the asynchronous iterative scheme. We formalized the mathematical model and gave some remarks on its convergence. Numerical experiments showed that asynchronous Laplace transform method converges chaotically to the correct solution sets. The intrinsic numerical approximation problem affects the accuracy of GaverStehfest algorithm and led to a narrow convergence interval. There exists much remains to be done to prove the convergence of this method, which might require a complete extension of the existing asynchronous theory. Further investigation could focus on the asynchronous formalization for other efficient inversion formulas and nonlinear solvers.
References
 [1] J. Abate and P. P. Valkó. Multiprecision Laplace transform inversion. Int. J. Numer. Methods Eng., 60(5):979–993, 2004.
 [2] J. Abate and W. Whitt. A unified framework for numerically inverting Laplace transforms. INFORMS J. Comput., 18(4):408–421, 2006.
 [3] G. Barles and H. M. Soner. Option pricing with transaction costs and a nonlinear BlackScholes equation. Financ. Stoch., 2(4):369–397, 1998.
 [4] G. M. Baudet. Asynchronous iterative methods for multiprocessors. J. ACM, 25(2):226–244, 1978.
 [5] D. P. Bertsekas. Distributed asynchronous computation of fixed points. Math. Program., 27(1):107–120, 1983.
 [6] D. P. Bertsekas and J. N. Tsitsiklis. Parallel and Distributed Computation: Numerical Methods. PrenticeHall, 1989.
 [7] F. Black and M. Scholes. The pricing of options and corporate liabilities. J. Political Econ., 81(3):637–654, 1973.
 [8] P. P. Boyle and T. Vorst. Option replication in discrete time with transaction costs. J. Financ., 47(1):271–293, 1992.
 [9] D. Chazan and W. Miranker. Chaotic relaxation. Linear Algebra Appl., 2(2):199–222, 1969.
 [10] A. M. Cohen. Numerical Methods for Laplace Transform Inversion. Springer, 2007.
 [11] J. C. Cox, S. A. Ross, and M. Rubinstein. Option pricing: A simplified approach. J. Financ. Econ., 7(3):229–263, 1979.
 [12] D. Crann, A. J. Davies, C.H. Lai, and S. H. Leong. Time domain decomposition for European options in financial modelling. Contemp. Math., 218:486–491, 1998.
 [13] B. Davies and B. Martin. Numerical inversion of the Laplace transform: a survey and comparison of methods. J. Comput. Phys., 33(1):1–32, 1979.
 [14] J.C. Duan. The GARCH option pricing model. Math. Financ., 5(1):13–32, 1995.
 [15] D. G. Duffy. On the numerical inversion of Laplace transforms: Comparison of three new methods on characteristic problems from applications. ACM Trans. Math. Softw., 19(3):333–359, 1993.
 [16] D. El Baz, P. Spitéri, J.C. Miellou, and D. Gazen. Asynchronous iterative algorithms with flexible communication for nonlinear network flow problems. J. Parallel Distrib. Comput., 38(1):1–15, 1996.
 [17] M. N. El Tarazi. Some convergence results for asynchronous algorithms. Numer. Math., 39(3):325–340, 1982. (in French).
 [18] A. Frommer and D. B. Szyld. Asynchronous twostage iterative methods. Numer. Math., 69(2):141–153, 1994.
 [19] A. Frommer and D. B. Szyld. On asynchronous iterations. J. Comput. Appl. Math., 123(12):201–216, 2000.
 [20] M. J. Gander. 50 years of time parallel time integration. In T. Carraro, M. Geiger, S. Körkel, and R. Rannacher, editors, Multiple Shooting and Time Domain Decomposition Methods, pages 69–113. Springer, Cham, 2015.
 [21] J. M. Harrison and D. M. Kreps. Martingales and arbitrage in multiperiod securities markets. J. Econ. Theory, 20(3):381–408, 1979.
 [22] J. Hull and A. White. The pricing of options on assets with stochastic volatilities. J. Financ., 42(2):281–300, 1987.
 [23] K. L. Kuhlman. Review of inverse Laplace transform algorithms for Laplacespace numerical approaches. Numer. Algorithms, 63(2):339–355, 2013.
 [24] A. Kuznetsov. On the convergence of the Gaver–Stehfest algorithm. SIAM J. Numer. Anal., 51(6):2984–2998, 2013.
 [25] C.H. Lai. On transformation methods and the induced parallel properties for the temporal domain. In F. Magoulès, editor, Substructuring Techniques and Domain Decomposition Methods, chapter 3, pages 45–70. SaxeCoburg Publications, Stirlingshire, UK, 2010.
 [26] C.H. Lai, A. K. Parrott, S. Rout, and M. E. Honnor. A distributed algorithm for European options with nonlinear volatility. Comput. Math. Appl., 49(5):885–894, 2005.
 [27] F. Magoulès and A.K. Cheik Ahamed. Alinea: An advanced linear algebra library for massively parallel computations on graphics processing units. Int. J. High Perform. Comput. Appl., 29(3):284–310, 2015.
 [28] F. Magoulès and G. GbikpiBenissan. JACK: An asynchronous communication kernel library for iterative algorithms. J. Supercomput., 73(8):3468–3487, 2017.

[29]
F. Magoulès and G. GbikpiBenissan.
Asynchronous Parareal time discretization for partial differential equations.
SIAM J. Sci. Comput., 40(6):C704–C725, 2018.  [30] F. Magoulès and G. GbikpiBenissan. JACK2: An MPIbased communication library with nonblocking synchronization for asynchronous iterations. Adv. Eng. Softw., 119:116–133, 2018.
 [31] F. Magoulès, G. GbikpiBenissan, and Q. Zou. Asynchronous iterations of Parareal algorithm for option pricing models. Mathematics, 6(4):1–18, 2018.
 [32] R. C. Merton. Theory of rational option pricing. Bell J. Econ., 4(1):141–183, 1973.
 [33] R. C. Merton. On the pricing of corporate debt: the risk structure of interest rates. J. Financ., 29(2):449–470, 1974.
 [34] J.C. Miellou. Algorithmes de relaxation chaotique à retards. ESAIM: M2AN, 9(R1):55–82, 1975. (in French).
 [35] J.C. Miellou, D. El Baz, and P. Spitéri. A new class of asynchronous iterative algorithms with order intervals. Math. Comput., 67(221):237–255, 1998.
 [36] D. Sheen, I. H. Sloan, and V. Thomée. A parallel method for timediscretization of parabolic problems based on contour integral representation and quadrature. Math. Comput., 69(229):177–195, 2000.
 [37] H. Stehfest. Numerical inversion of Laplace transforms. Commun. ACM, 13(1):47–49, 1970.
 [38] A. Talbot. The accurate numerical inversion of Laplace transforms. IMA J. Appl. Math., 23(1):97–120, 1979.