Geometric Brownian motion with affine drift and its time-integral
aa r X i v : . [ q -f i n . M F ] D ec Geometric Brownian motion with affine drift and itstime-integral ⋆ Runhuan Feng a , Pingping Jiang b, ∗ , Hans Volkmer c a Department of Mathematics, University of Illinois at Urbana-Champaign, Illinois, USA b School of Management and Economics,The Chinese University of Hong Kong, Shenzhen,Shenzhen, Guangdong, 518172, P.R. ChinaSchool of Management, University of Science and Technology of China, Hefei, Anhui, 230026,P.R.China c Department of Mathematical Sciences, University of Wisconsin–Milwaukee, Wisconsin, USA
Abstract
The joint distribution of a geometric Brownian motion and its time-integral wasderived in a seminal paper by Yor (1992) using Lamperti’s transformation, leadingto explicit solutions in terms of modified Bessel functions. In this paper, we revisitthis classic result using the simple Laplace transform approach in connection to theHeun differential equation. We extend the methodology to the geometric Brownianmotion with affine drift and show that the joint distribution of this process and itstime-integral can be determined by a doubly-confluent Heun equation. Furthermore,the joint Laplace transform of the process and its time-integral is derived from theasymptotics of the solutions. In addition, we provide an application by using theresults for the asymptotics of the double-confluent Heun equation in pricing Asianoptions. Numerical results show the accuracy and efficiency of this new method.
Keywords:
Doubly-confluent Heun equation, geometric Brownian motion withaffine drift, Lamperti’s transformation, asymptotics, boundary value problem. ⋆ R. Feng is supported by an endowment from the State Farm Companies Foundation. P. Jiangis supported by the Postdoctoral Science Foundation of China (No.2020M671853) and the NationalNatural Science Foundation of China (No. 11631004, 71532001). ∗ Corresponding authorEmail address: [email protected] (R. Feng); [email protected] (P. Jiang);[email protected] (H. Volkmer).
Preprint submitted to Elsevier December 18, 2020 . Introduction
Heun’s differential equation (c.f. Ronveaux [22]) is d ydz + (cid:18) γz + δz − ǫz − a (cid:19) dydz + ηz − qz ( z − z − a ) y = 0 , (1.1)where γ, δ, ǫ, η, q and a = 0 , z = 0 , , a, ∞ . Heun’s equation has four confluent forms, the confluent,doubly-confluent, biconfluent and triconfluent Heun equations. In recent years theHeun equation and its confluent forms have found many applications in natural sci-ences. In the theory of black holes, the perturbation equations of massless fields forthe Kerr-de Sitter geometry can be written in the form of separable equations. Theequations have five definite singularities so that the analysis has been expected to bedifficult. Suzuki et al. [23] showed that these equations can be transformed to Heun’sequations thus the known technique could be used for the analysis of the solutions.As Schr¨odinger equation for harmonium and related models may be transformed tothe biconfluent Heun equation, Karwowski and Witek [16] discussed the solubility ofthis equation and its applications in quantum chemistry. Chugunova and Volkmer[5] also investigated that the set of eigenvalues of a non-self-adjoint differential op-erator arising in fluid dynamics were related to the eigenvalues of Heun’s differentialequation.In this paper we present applications of the doubly-confluent Heun equation instochastic analysis. The doubly confluent Heun equation can be obtained from (1.1)as follows. First, by substitution x = hz , we move the singularity z = 1 to x = h .The equation becomes x ( x − h )( x − A ) d ydx + ( cx ( x − A ) + d ( x − A ) − νAx ( x − h )) dydx + A ( τ x + b ) y = 0 , (1.2)where A = ah, c = γ + δ, d = − γh, ǫ = − νA, η = τ A, hq = − bA. Dividing (1.2) by − A and letting h → A → ∞ (double confluence), we obtain x d ydx + (cid:0) νx + cx + d (cid:1) dydx − ( τ x + b ) y = 0 . (1.3)A study of equation (1.3) can be found in Ronveaux [22, Part C]. In this paper, weconsider this equation for the special case where ν = 0, x d ydx + ( cx + d ) dydx − ( τ x + b ) y = 0 . (1.4)2t should be pointed out that, while this equation appeared in the literature as aspecial case, its solutions were not studied in previous literature.In this paper, we explore the close relationship between the doubly-confluent Heunequation (1.4) and the geometric Brownian motion with affine drift. We mainly focuson this process whose infinitesimal generator is given by( G f )( x ) = 2 x f ′′ ( x ) + [2( ν + 1) x + 1] f ′ ( x ) , where ν ∈ R . This process has many natural applications in finance and insurance.For example, it can be used to model the price dynamics of dividend-paying stockswith constant dividend rate, see Lewis [17] for more details. While it is not an absoluterepresentation of reality, the model can be used to approximate the dynamics of actualstocks with dividends paid at discrete points, as an improvement to the classicalBlack-Scholes-Merton model.In the classical Black-Scholes model for option pricing, stock prices are modeledby geometric Brownian motions. Among the best studied exotic options is the Asianoption, whose payoff depends on the average of stock prices. In a continuous-timemodel, the average price is modeled by the integral of geometric Brownian motionover the life of the option divided by the length of that period. Motivated by sucha pricing problem, Yor [24] developed a series of papers on exponential functionalsof Brownian motion B t . At the center of focus for financial application is the jointdistribution of (exp { B ( ν ) t } , A ( ν ) t ) where B ( ν ) t := νt + B t , A ( ν ) t := Z t exp { B ( ν ) s } d s. (1.5)The derivation of the joint distribution was based on a combination of time change,change of measures and their connections to Bessel processes. However, in moresophisticated diffusion models, such as the geometric Brownian motion with affinedrift, a generalization of such techniques seems to be less fruitful. Hence, in thispaper, we use the connection between the joint distribution of time-homogeneousdiffusion and its time-integral and differential equations to develop computationalmethods. Linetsky [18] also used the connection between the time integral of stockprices and the geometric Brownian motion with affine drift for pricing Asian options.The process is also used in the risk management of an insurer’s net liability in variableannuity guaranteed benefits in Feng and Volkmer [7, 8, 9] and Feng and Jing [6].The rest of the paper is organized as follows. In Section 2, we introduce particularsolutions to the Heun equation (1.4) that are required for applications to stochasticanalysis for various applications in finance. While the Heun equation is alreadyknown in analysis, the particular solutions presented here are not known previouslyin the literature. The paper also proposes a method to compute these solutions to3igh accuracy. It is shown in Section 3 that the joint distribution of a solution X t to a stochastic differential equation and its integral Y t = R t X s ds is determined byappropriate solutions to a linear differential equation of the second order. In Section 4,the special case of X t being a geometric Brownian motion with affine drift is discussedand its connection to the Heun equation (1.4) is revealed. Results from Sections 2are used to develop a numerical method to compute the joint distribution of X t and Y t . Section 5 is dedicated to a new stochastic process Z t resulting from the Lampertitransformation X t = Z Y t . The distribution of Z t can again be computed by solvingequation (1.4). An explicit expression for the joint Laplace transform of the geometricBrownian motion with affine drift and its time integral is provided in Section 6,which relies on asymptotics of a solution to a boundary value problem involving aninhomogeneous equation corresponding to (1.4). By using the results, we providean application in pricing Asian options on divided-paying stocks. Numerical resultsshow the accuracy and efficiency of this new method.
2. A doubly confluent Heun equation
Let us consider differential equation (1.4) (with τ replaced by a ) x y ′′ ( x ) + ( cx + d ) y ′ ( x ) − ( ax + b ) y ( x ) = 0 , < x < ∞ , (2.1)containing four real parameters a, b, c, d . The equation has two singularities: x = 0and x = ∞ . If d = 0 then x = 0 is an irregular singularity. If d = 0 then x = 0 is aregular singularity, and the general solution of (2.1) is y ( x ) = x (1 − c ) (cid:8) C I λ (2 √ ax ) + C K λ (2 √ ax ) (cid:9) , (2.2)where I λ , K λ denote modified Bessel functions, and λ = p ( c − + 4 b. (2.3)If a = 0 then x = ∞ is an irregular singularity. If a = 0 then x = ∞ is a regularsingularity, and the general solution of (2.1) is y ( x ) = x − µ (cid:8) C M (cid:0) µ, λ, dx (cid:1) + C U (cid:0) µ, λ, dx (cid:1)(cid:9) , (2.4)where µ = ( c − λ ) , and M, U denote Kummer functions (c.f. Olver et al. [20, Chapter 13]).We now turn to the general case d = 0 and a = 0. We will apply methods ofOlver [21, Chapter 7] to introduce solutions of (2.1). We are interested in finding arecessive solution y ( x ) at x = 0, and a recessive solution y ( x ) at x = ∞ , that is,solutions with the properties that y ( x ) /y ( x ) → x → + and y ( x ) /y ( x ) → x → ∞ for every solution y which is linearly independent of y , y , respectively.4 heorem 2.1. Suppose a > . The differential equation (2.1) has a unique solution y on (0 , ∞ ) such that, as x → ∞ , y ( x ) ∼ x − c exp( − √ ax ) ∞ X k =0 A k x − k , (2.5) where A = 1 and A k ’s are determined by the following recursion for k = 1 , , · · · . A k = 4 b − ( k + − c )( k − + c )4 √ ak A k − + dk A k − , (2.6) with the understanding that A − = 0 . The asymptotic formula (2.5) may be differen-tiated term-by-term. If b ≥ then the solution y is positive and decreasing. Proof.
Chapter 7 of Olver [21] considers the differential equation w ′′ ( z ) + f ( z ) w ′ ( z ) + g ( z ) w ( z ) = 0 , for which the singularity is located at infinity and f ( z ) = ∞ X k =0 f k z k , g ( z ) = ∞ X k =0 g k z k . The first goal is to construct a formal solution of the form w ( z ) = e λz z µ ∞ X k =0 A k z − k . (2.7)In the case of (2.1), we have f = c, f = d, g = − a, g = − b, and all other f i ’s and g i ’s are zero. This falls into the category where f = 4 g andFabry’s transformation should be used. Let z = x / and y ( x ) = w ( z ). Then (2.1)becomes w ′′ ( z ) + (cid:18) c − z + 2 dz (cid:19) w ′ ( z ) − (cid:18) a + bz (cid:19) w ( z ) = 0 . (2.8)We have f = 2 c − , f = 2 d, g = − a, g = − b, and all other f i ’s and g i ’s are zero. The roots of λ + f λ + g = 0 are λ = ± a / .Since we assumed that a >
0, these solutions are distinct. We are interested in arecessive solution of (2.8) at x = ∞ , so we pick λ = − a / . We compute µ = − c A k arecomputed from A = 1 and recursively from Olver [21, (1.11)], which can be writtenas (2.6). By Olver [21, Theorem 2.1, page 232], the ODE (2.1) has a solution y ( x )with the asymptotic behavior (2.5). Since it is the recessive solution, y ( x ) is uniquelydetermined by (2.5).Actually, (2.5) holds in an open sector of the complex plane containing the positivereal axis. This justifies term-by-term differentiation in (2.5).We now prove that the solution y is positive and decreasing if b ≥
0. For every x > x > x such that y ( x ) > y ′ ( x ) <
0. If y is not decreasingon (0 , x ) there would be a point x < x such that y ′ ( x ) = 0 and y is decreasing on( x , x ). This is impossible because (2.1) under the assumption a > , b ≥ y ′′ ( x ) >
0. As x is arbitrary, y ( x ) is positive and decreasing for all x ∈ (0 , ∞ ). (cid:3) Remark 2.1.
When d = 0 the recursive formula (2.6) can be simplified to A k = 12 √ a a k ( λ ) a k − ( λ ) A k − , k = 1 , , · · · , where λ is given by (2.3) and a k ( λ ) = (4 λ − )(4 λ − ) · · · (4 λ − (2 k − ) k !8 k . Hence, we get A k = 1(2 √ a ) k a k ( λ ) . Substituting the expression for A k in (2.5) yields, as x → ∞ , y ( x ) ∼ x / − c/ exp( − √ ax ) ∞ X k =1 a k ( λ )(2 √ ax ) − k . (2.9) Obviously, apart from a constant multiple, y ( x ) is equal to x (1 − c ) K λ (2 √ ax ) . and (2.9) agrees with Hankel’s expansion (Olver et al. [20, (10.40.2)]), K λ ( z ) ∼ (cid:16) π z (cid:17) − e − z ∞ X k =1 a k ( λ ) z k . Theorem 2.2. If d > , the differential equation (2.1) has a unique solution ˆ y on (0 , ∞ ) such that, as x → + , ˆ y ( x ) ∼ ∞ X k =0 A k x k , (2.10)6 here A = 1 and A k ’s are determined by the following recursion for k = 1 , , · · · . A k = b + ( k + c − − k ) kd A k − + akd A k − , (2.11) with the understanding that A − = 0 . If d < , the equation (2.1) has a uniquesolution y on (0 , ∞ ) such that, as x → + , y ( x ) ∼ exp (cid:18) dx (cid:19) x − c ∞ X k =0 A k x k , (2.12) where A = 1 and A s are determined by the following recursion for k = 1 , , · · · . A k = ( k − c + 1) k − bkd A k − − akd A k − , (2.13) with the understanding that A − = 0 . The asymptotic formulas (2.10) , (2.12) may bedifferentiated term-by-term. If a > , b ≥ then y is positive and increasing. Proof.
We substitute z = 1 /x , y ( x ) = w ( z ) in (2.1) and obtain w ′′ ( z ) + (cid:18) − cz − d (cid:19) w ′ ( z ) − (cid:18) az + bz (cid:19) w ( z ) = 0 , < z < ∞ , for which the singularity at z = ∞ corresponds to the singularity at x = 0 in (2.1).We now argue as in the proof of Theorem 2.1. Now we have f = 2 − c, f = − d, g = − b, g = − a, and all other f i ’s and g i ’s are zero.If d >
0, we consider a formal solution of the form (2.7) where λ = 0, µ = 0 and A k ’s satisfy the recursion (2.11). We apply Olver [21, Theorem 2.1, page 232] andobtain a (recessive) solution ˆ y satisfying (2.10). It follows that if b >
0, as x → + ,ˆ y ( x ) ∼ bd x. Or if b = 0, we have as x → + , ˆ y ( x ) ∼ a d x . For any x >
0, there must exist x < x such that ˆ y ( x ) > y ′ ( x ) >
0. If ˆ y isnot increasing on ( x , ∞ ), then there would be a point x > x such that ˆ y ′ ( x ) = 0and ˆ y is increasing on ( x , x ). This is impossible because (2.1) under the assumption7hat a, d > b ≥ y ′′ ( x ) >
0. Hence ˆ y ( x ) is positive and increasingon ( x , ∞ ) for some x < x , given any x >
0. As x is arbitrary, ˆ y ( x ) is increasingand positive for all x ∈ (0 , ∞ ).If d <
0, we consider the formal solution of the form (2.7) with λ = d and µ = c −
2. It follows immediately that A k ’s satisfy the recursion (2.13). Using thesame arguments as in the previous case, we can show that y determined by (2.12) ispositive and increasing on (0 , ∞ ) provided that a > , b ≥ (cid:3) The formal adjoint of equation (2.1) is given by( x w ) ′′ − (( cx + d ) w ) ′ − ( ax + b ) w = 0 . (2.14)This equation is of the form as (2.1) with a new set of parameter values ˜ a , ˜ b , ˜ c , ˜ d ,where ˜ a = a, ˜ b = b + c − , ˜ c = 4 − c, ˜ d = − d. Equation (2.14) can be transformed to (2.1) by setting y ( x ) = x − c exp (cid:18) dx (cid:19) w ( x ) . In particular, the case d < d > y , y introduced in Theorems 2.2, 2.1, respectively. Wecompute y for d > m and, according to(2.10), approximate ˆ y ( x ) ≈ m − X k =0 A k x k . This approximation will be good only for sufficiently large x . By differentiating termby term we also obtain an approximation for y ′ ( x ). Using these approximations forˆ y ( x ), ˆ y ′ ( x ) we solve the initial value problem for equation (2.1) using a suitablenumerical procedure. If we need ˆ y ( x ) to high accuracy it is preferable to use aprocedure like gear implemented in Maple.The computation of y for d < y is similar.
3. Joint distributions
Consider the stochastic differential equation (SDE) on an open interval (0 , ∞ ) dX t = β ( X t ) d t + γ ( X t ) d B t , (3.1)8ith initial condition X = x > Y = { Y t , t ≥ } where Y t = Z t X s d s. Then the pair (
X, Y ) satisfies the SDE systemd X t = β ( X t ) d t + γ ( X t ) d B t , (3.2)d Y t = X t d t, (3.3)under the probability measure P x we have P x ( X = x , Y = 0) = 1.The work of H¨ormander [13] provides a condition for the smoothness of fundamen-tal solutions to hypoelliptic second order differential equations, known as H¨ormander’s“brackets condition” in the literature. Hairer [12] provides an extension to parabolicequations in the context of transition probabilities for diffusions. Consider a secondorder differential operator P = r X j =1 X j + X , where X , X , · · · , X r denote first order homogeneous differential operators in an openset Ω ⊂ R n with C ∞ coefficients and c ∈ C ∞ (Ω). Define a collection of vector fields V k by V = {X i : i > } , V k +1 = V k ∪ { [ U , X j ] : U ∈ V k and j ≥ } . Also define the vector space V k ( x ) = span { V ( x ) : V ∈ V k } . The parabolic H´’ormander’scondition holds if ∪ k ≥ V k ( x ) = R for every x ∈ R . The Lie bracket [ · , · ] is given by[ X j i , X j k ] = X j i X j k −X j k X j i . The infinitesimal generator satisfies the H¨ormander condi-tion, then the Markov process admits a smooth density with respect to Lebesgue mea-sure and the corresponding Markov semigroup maps bounded functions into smoothfunctions.Note that the infinitesimal generator of the two-dimensional process ( X t , Y t ) in(3.2) and (3.3) is given by L = X + X , where the two differential operators are defined by X = β ∗ ( x ) ∂∂x + x ∂∂y , β ∗ ( x ) = β ( x ) − γ ′ ( x ) γ ( x ) X = γ ∗ ( x ) ∂∂x , γ ∗ ( x ) = 1 √ γ ( x ) .
9t is easy to show that[ X , X ] = W ( β ∗ , γ ∗ ) ∂∂x − γ ∗ ( x ) ∂∂y ;[[ X , X ] , X ] = W ( W ( β ∗ , γ ∗ ) , β ∗ ) ∂∂x + [2 β ∗ ( x ) γ ′∗ ( x ) − γ ∗ ( x ) β ′∗ ( x )] ∂∂y where W is the Wronskian. For example, [ X , X ] and [[ X , X ] , X ] are linear inde-pendent at any point, unless there are points at which[2 β ∗ ( x ) γ ′∗ ( x ) − γ ∗ ( x ) β ′∗ ( x )] W ( β ∗ , γ ∗ ) + γ ∗ ( x ) W ( W ( β ∗ , γ ∗ ) , β ∗ ) = 0 . (3.4)Suppose that (3.4) is not true for any point x >
0. Then H¨ormander condi-tion implies that there exists a smooth joint density function of ( X t , Y t ), denoted by p ( t, x, y ) = p ( t, x, y ; x ) , i.e. P x ( X t ∈ d x, Y t ∈ d y ) = p ( t, x, y ) d x d y. The corresponding Kolmogorov forward PDE for p ( t, x, y ) is given by ∂p∂t = − ∂∂x ( β ( x ) p ) − ∂∂y ( xp ) + 12 ∂ ∂x ( γ ( x ) p ) , where the assumption P x ( X = x , Y = 0) = 1 implies that p (0 , x, y ) = δ ( x − x ) δ ( y ) , where δ is the Dirac delta function for which Z ∞ f ( x ) δ ( x ) d x = f (0) , Z ∞−∞ δ ( x ) d x = 1 . We can solve this equation by using the Laplace transform q ( x ) := q ( s, x, w ; x ) = Z ∞ Z ∞ exp( − st ) exp( − wy ) p ( t, x, y ) dt dy. (3.5)For fixed s, w > q ( x ) −
12 ( γ ( x ) q ) ′′ + ( β ( x ) q ) ′ + ( wx + s ) q = δ ( x − x ) . (3.6)This equation is the same as appears in the definition of a Green’s function at x = x .We choose a fundamental system q ( x ) , q ( x ) of solutions of (3.6) with right-hand side0. Then we write the solution q ( x ) of (3.6) in the form q ( x ) = ( C q ( x ) if 0 < x ≤ x , C q ( x ) if x < x < ∞ . (3.7)10efinition (3.5) shows that q ( x ) ≥ R ∞ q ( x ) d x ≤ /s < ∞ , so q ( x ), q ( x ) haveto be integrable on (0 , x ) and ( x , ∞ ), respectively. The constants C and C haveto be chosen such that q ( x ) is continuous at x = x , i.e. C q ( x ) = C q ( x ) , (3.8)and such that 12 γ ( x ) ( C q ′ ( x ) − C q ′ ( x )) = 1 . (3.9)The latter equation follows from (3.6) by integrating both sides x − ǫ to x + ǫ , thenletting 0 < ǫ → q ( x ) = C ( q ( x ) q ( x ) if 0 < x ≤ x , q ( x ) q ( x ) if x < x < ∞ , (3.10)where C = 2 γ ( x ) q ( x ) q ′ ( x ) − q ′ ( x ) q ( x ) . (3.11)We will consider two examples in this work. In this section we take X to begeometric Brownian motion and reproduce the density p ( t, x, y ), which is known fromthe work of Yor. In the next section, we take X to be geometric Brownian motionwith affine drift where no previous results on the joint density p ( t, x, y ) is known inthe current literature.The solution of the SDE (3.1) with β ( x ) = (2 ν + 2) x , ν ∈ R , γ ( x ) = 2 x , x = 1is the geometric Brownian motion X t = exp(2 νt + 2 B t )with integral Y t = R t X s ds . Then the homogeneous form of the ODE (3.6) is − x q ) ′′ + 2( ν + 1)( xq ) ′ + ( wx + s ) q = 0 . (3.12)This is equation (2.1) with parameter values a = w , b = s ν − , c = 3 − ν, d = 0 . According to (2.2), (3.12) has the fundamental system q ( x ) = x ν/ − I λ ( √ wx ) , q ( x ) = x ν/ − K λ ( √ wx ) , where λ := √ s + ν . q ( x ) and q ( x ) are the only solutions of(3.12) that are integrable on (0 ,
1) and (1 , ∞ ), respectively. Using the WronskianOlver et al. [20, 10.28.2], we obtain that q ( x ) q ′ ( x ) − q ′ ( x ) q ( x ) = − x ν − . Setting x = 1 , we obtain C = 1 in (3.11). Therefore, (3.10) gives q ( s, x, w ) = x ν/ − ( K λ ( √ w ) I λ ( √ wx ) if 0 < x ≤ ,I λ ( √ w ) K λ ( √ wx ) if 1 < x < ∞ . (3.13)Essentially we have solved our problem. It remains to invert the Laplace transforms.We may use the formula Gradshteyn and Ryzhik [11, 6.653] Z ∞ exp (cid:18) − z − a + b z (cid:19) I λ (cid:18) abz (cid:19) dzz = 2 ( I λ ( a ) K λ ( b ) if 0 < a < bK λ ( a ) I λ ( b ) if 0 < b < a (3.14)which holds for λ > −
1. If we set a = √ wx , b = √ w , z = 2 wy we obtain Z ∞ exp( − st ) p ( t, x, y ) dt = 12 y x ν/ − exp (cid:18) − x y (cid:19) I λ (cid:18) √ xy (cid:19) (3.15)which is a known result. We can also invert the Laplace transform with respect to t employing the Hartmann-Watson density. Note that Yor obtained the same result us-ing much more complex probabilistic arguments based on Lamperti’s transformationand Girsanov change of measure in Yor [24]. A detailed account of Hartman-Watsondensity function can be found in Barrieu, Rounault and Yor [3].
4. Joint distribution of Geometric Brownian motion with affine drift andits integral
Let ν ∈ R , x > β ( x ) = (2 ν + 2) x + 1, γ ( x ) = 2 x . It is easy to verify that (3.4)does not hold for x >
0. Then the solution of SDE (3.1) is the geometric Brownianmotion with affine drift X t = exp(2 νt + 2 B t ) (cid:18) x + Z t exp( − (2 νs + 2 B s )) ds (cid:19) . (4.1)Its integral is Y t = Z t X s ds. (4.2)12ur goal is to determine the joint density function p ( t, x, y ) of ( X t , Y t ) by themethod proposed in Section 3. The corresponding ODE (3.6) is − x q ) ′′ + { (2( ν + 1) x + 1) q } ′ + ( wx + s ) q = δ ( x − x ) . (4.3)If we replace the right-hand side of (4.3) by 0, we obtain the doubly-confluent Heunequation (2.1) with parameter values a = w > , b = s ν − , c = 3 − ν, d = − < . We choose q j ( x ) = y j ( x ), j = 1 ,
2, with y , y introduced in Theorems 2.2, 2.1.Note that, apart from constant multiples, q ( x ) is the only solution of (4.3) whichis integrable on ( x , ∞ ). The choice of q ( x ) can be justified as follows. We have0 ≤ q ( s, x, w ) ≤ q ( s, x, q ( s, x,
0) is the Laplace transform of the transitiondensity function of X t . It is known that q ( s, x, → x → + . Therefore, q ( s, x, w ) → x → + . Now, apart from constant multiples, q ( x ) is the onlysolution of (4.3) that tends to 0 as x → + .We can then write q ( x ) defined by (3.5) in the form (3.10), where C = 12 x q ( x ) q ′ ( x ) − q ′ ( x ) q ( x ) . We cannot expect a simple formula for the Wronskian of solutions q and q but wecan compute C numerically.We use the following algorithm to determine the joint density p ( t, x, y ). In thefirst step, we use (2.12) to approximate initial values for q and (2.5) for q . In orderto increase the accuracy of the computation, we define u ( x ) = exp (cid:18) − dx (cid:19) x c − q ( x ) , which satisfies the ODE x u ′′ ( x ) + [(4 − c ) x − d ] u ′ ( x ) − ( ax + b + c − u ( x ) = 0 . (4.4)Thus the asymptotics of u ( x ) is determined by the power series in (2.12). We specifythe number m of terms in the partial sum approximation and use an algorithm todetermine the small initial point x ℓ > − k , roughly A m x ml < − k . Then we use numerical methods to find solutionsto (4.4) with initial conditions at x ℓ , which in turn produces the value of q in theinterval ( x l , x ). 13imilarly, we define w ( z ) = z c − / exp(2 √ az ) q ( z ) . Then w ′′ ( z ) + (cid:18) − √ a + 2 dz (cid:19) w ′ ( z ) + (cid:18) − − c + 2 c − bz − d √ az + d (1 − c ) z (cid:19) w ( z ) = 0 . The asymptotics of w ( z ) is given by the power series in (2.7). For each specifiednumber m of terms in the partial sum approximation, we determine the large ini-tial point x r such that the approximation error of w ( z ) is less than 10 − k , roughly A m x − m/ r < − k . We use numerical methods to solve the initial value problem for w ( z ) with initial conditions determined by approximations, thereby leading to solu-tion q in ( x , x r ). Combining the computations of q and q , we find the Laplacetransform q using (3.10).In the second step, we use two-dimensional Laplace inversion routines to find p atvarious point of ( x, y ). Details on various types of two-dimensional Laplace inversioncan be seen in Abate and Whitt [2]. Here we obtain the results using the Talbot-Gaver-Stehfest algorithm p ( t, x, y ) = 2 ln 25 ty M − X k =0 ℜ ( γ k M X k =1 ζ k q (cid:18) δ k t , x, k ln 2 y (cid:19)) , where δ = 2 M , δ k = 2 kπ (cid:18) cot (cid:18) kπM (cid:19) + i (cid:19) , < k < M, γ = e δ ,γ k = (cid:20) ikπM (cid:18) (cid:18) kπM (cid:19)(cid:19) − i cot (cid:18) kπM (cid:19)(cid:21) e δ k , < k < M.ζ k = ( − M + k k ∧ M X j = ⌊ ( k +1) / ⌋ j M +1 M ! (cid:18) Mj (cid:19)(cid:18) jj (cid:19)(cid:18) jk − j (cid:19) . To test the accuracy of the results, we also use the Euler-Gaver-Stehfest algorithm p ( t, x, y ) = 10 M/ ln 25 t t M X k =0 η k M − X k =0 ℜ (cid:26) γ k q (cid:18) β k t , δ k t (cid:19) + γ k q (cid:18) β k t , δ k t (cid:19)(cid:27) , where γ k is the complex conjugate of γ k and η = 12 , η k = ( − k , ≤ k ≤ M, η k = ( − k M X i = k − M (cid:18) Mi (cid:19) − M , M < j ≤ M.
14n this numerical example, we choose ν = 1 . , M = 7. According to Abate andWhitt [2], both algorithms are expected to be accurate up to 0 . M ≈ p (1 , ,
4) rounded up to 7 digitsin Table 1 using the two inversion algorithms with various choices of precision. Theyall agree up to four decimal places. m k
Talbot-Gaver-Stehfest Euler-Gaver-Stehfest30 15 0 . . . . Table 1: Joint density p (1 , , The method can be easily extended to the geometric Brownian motion with affinedrift X t = exp(2 νt + 2 B t ) (cid:18) x − Z t exp( − (2 νs + 2 B s )) ds (cid:19) , Y t = Z t X s d s, with X = x > Y = 0, in which case the Laplace transform q ( x ) satisfies theODE − x q ) ′′ + { (2( ν + 1) x − q } ′ + ( wx + s ) q = 0 , which is also a special case of (2.1) with a = w > , b = s ν − , c = 3 − ν, d = 12 > . We provide the results on p (1 , ,
2) with ν = 1 . M = 7 in Table 2. m k Talbot-Gaver-Stehfest Euler-Gaver-Stehfest30 15 0 . . . . Table 2: Joint density p (1 , ,
5. A diffusion process from Lamperti’s transformation
We recall from Yor [24] that the geometric Brownian motion and its time-integralare connected through Lamperti’s transformation, i.e.exp { B ( ν ) t } = ρ ( ν ) A ( ν ) t , t ≥ , where B ( ν ) t , A ( ν ) t are defined by (1.5), and ρ is a Bessel process with index ν startingfrom 1. Using the same idea, we can connect the geometric Brownian motion withaffine drift and its time-integral through another diffusion process.15 heorem 5.1. Let X t be the solution of the SDE (3.1) with β ( x ) = µx +1 , γ ( x ) = σx ,where µ ∈ R , σ > and X = x > . Then X t = Z R t X s d s , t ≥ , (5.1) where Z is a diffusion process determined by the following SDE d Z t = (cid:18) µ + 1 Z t (cid:19) d t + σ p Z t d B t . (5.2) Proof.
Let τ ( t ) = inf { u : R u X s d s > t } . Note that R t X s d s is a process withcontinuous sample path. Then R τ ( t )0 X s ds = t . Thus, X τ ( t ) = x + Z τ ( t )0 ( µX s + 1) d s + σ Z τ ( t )0 X s d B s = x + µt + Z τ ( t )0 d s + σ Z τ ( t )0 X s d B s . Using the time change formula for Ito integrals (Oksendal [19, Theorem 8.5.7, p156]),we obtain Z τ ( t )0 X s d B s = Z t X τ ( r ) p τ ′ ( r ) d W r , where W t = R τ ( t )0 √ X s d B s is also a Brownian motion. Note that τ ′ ( r ) = 1 /X τ ( r ) using the derivative of inverse function. Thus, Z τ ( t )0 X s d B s = Z t p X τ ( r ) d W r . We also have Z τ ( t )0 d s = Z t τ ′ ( r ) d r = Z t X τ ( r ) d r. Let Z t := X τ ( t ) for t ≥
0. Then, Z t = Z + µt + Z t Z r d r + σ Z t p Z r d W r . Thus the claim follows immediately after reversing the time change. (cid:3)
To the authors’ best knowledge, this diffusion process { Z t , t ≥ } was not previ-ously studied in the probability literature. When starting from a positive initial valueand µ <
0, this process is always positive and possesses a mean-reverting property,which is a desirable property for modeling many physical phenomenons. If we removethe Brownian perturbation by setting σ = 0, then Z t = − µ (cid:16) W ( e − Cµ − µ t − ) (cid:17) , W is the Lambert W-function.The results from Section 2 make it possible to compute the transition densityof the diffusion process Z t without having to use the often expensive Monte Carlosimulations. An account of the theory of time-homogeneous diffusion processes canbe found in Borodin and Salminen [4, Chapter 2]. In this case, it is easier to computethe transition density function p with respect to speed measure, i.e. P z ( Z t ∈ d z ) = p ( t, z , z ) d z = p ( t, z , z ) m ( z ) d z, (5.3)where Z = z >
0. The speed density and scale density of the diffusion process (5.2)are given by m ( z ) = 2 σ z µσ − exp (cid:26) σ (cid:18) − z (cid:19)(cid:27) , s ( z ) = z − µσ exp (cid:26) σ (cid:18) y − (cid:19)(cid:27) . Then p satisfies the forward Kolmogorov equation ∂p∂t = − ∂∂z ( β ( z ) p ) + 12 ∂ ∂z ( γ ( z ) p ) , while p solves the backward Kolmogorov equation. ∂ p ∂t = β ( z ) ∂ p ∂z + 12 γ ( z ) ∂ p ∂z , where β ( z ) = µ + 1 /z , γ ( z ) = σ √ z . Consider the Laplace transform q ( s, z , z ) = Z ∞ e − st p ( t, z , z ) d t. Then q ( z ) := q ( s, z , z ) satisfies the ODE σ z q ′′ ( z ) + (cid:18) µ + 1 z (cid:19) q ′ ( z ) − s q ( z ) = 0 , z > . (5.4)This ODE is a special case of (2.1) with parameter values a = 2 sσ > , b = 0 , c = 2 µσ , d = 2 σ > . Let q , q be the positive monotone solutions y , y introduced in Theorems 2.2,2.1, respectively. Then q ( s, z , z ) = 1 w ( s, z ) ( q ( z ) q ( z ) if 0 < z < z , q ( z ) q ( z ) if z < z , (5.5)17here w ( z ) denotes the Wronskian w ( z ) = 1 s ( z ) [ q ′ ( z ) q ( z ) − q ( z ) q ′ ( z )] . The following algorithm is used in this numerical example to calculate the transitiondensity p ( t, z , z ).1. Approximate the values of q and q at some initial points by the asymptotics(2.10) and (2.5). Take the approximations as initial conditions and use numer-ical methods for initial value problems to determine q and q .2. Substitute the values of q and q in (5.5) to determine values of the Laplacetransform q . Use numerical methods for inverting the Laplace transform toevaluate p . Figure 1: Probability density function p (1 , , z ) of the process Z In our example, we use the parameters µ = 0 . , σ = 1 and initial position z = 1.In the first step, we find approximation of q (1 /
25) and q (25) using the first 40 termsin the asymptotics (2.10) and (2.5). Then we use Maple’s numerical method gear todetermine solutions to the initial value problems, which yield solutions of q and q
18n the interval (1 / , p ( t, z , z ) ≈ n X k =1 ( − n − k k n k !( n − k )! ˜ f k ( t, z , z ) , where ˜ f k ( t, z , z ) = ln 2 t (2 k )! k !( k − k X h =1 ( − h (cid:18) kh (cid:19) q (cid:18) ( h + k ) ln 2 t , z , z (cid:19) , and the Euler algorithm (cf. Abate and Whitt [1]) p ( t, z , z ) ≈ m X k =0 (cid:18) mk (cid:19) − m s n + k ( t ) , where s n ( t ) = e A/ t ℜ ( q ) (cid:18) A t , z , z (cid:19) + e A/ t n X k =1 ( − k ℜ ( q ) (cid:18) A + 2 kπi t , z , z (cid:19) . Because the two algorithms are sufficiently different, we can verify the accuracy ofour numerical method by observing the results on probability density function p fromboth algorithms which agree up to at least four decimal places in all calculations.As recommended by Abate and Whitt [1, page 38], we choose A = 18 .
4. Figure 1displays the graph of the probability density function p (1 , , z ) of the diffusion process Z .
6. Joint Laplace transform of the geometric Brownian motion with affinedrift and its time-integral
An alternative form of the GBM with affine drift to (4.1) is given byd X t = ( rX t − δ ) d t + σX t d B t , X = x > . (6.1)Since such a model typically arises in the context of option pricing in finance literature,we use financial interpretation of model parameters. Assume that X represents thedynamics of stock prices. Dividends are paid out at a constant rate δ > r per timeunit. It should be pointed out that in this model stock prices may hit zero at onepoint, which is considered the time of ruin, prior to a fixed option maturity date T . In19his case, the process is assumed to be absorbed at zero once it hits zero. One can alsofind that the SDE (6.1) is a special case of (3.1) with γ ( x ) = σx and β ( x ) = rx − δ .We focus on the joint Laplace transform of the process and its time-integral whichplays an important role in an application to Asian option. The Laplace transform isdefined by h ( x , T ; w ) := E x h e − w R T X t d t − λX T i , (6.2)where λ, w >
0. Then it follows immediately from the Feymann-Kac formula that h satisfies the PDE for 0 < t < T and x > ∂h∂t + wxh = ( rx − δ ) ∂h∂x + σ x ∂ h∂x , (6.3)subject to ( h ( x ,
0) = e − λx , x > h (0 , t ) = 1 , < t < T. Consider the Laplace transform ˜ h ( x, s ) := R ∞ e − st h ( x, t ) d t for some s >
0. TakingLaplace transforms with respect to t on both sides of (6.3) we obtain the followingODE of ˜ h ( x ) = ˜ h ( x, s ) for x > σ x ˜ h ′′ ( x ) + ( rx − δ )˜ h ′ ( x ) − ( s + wx )˜ h ( x ) = − e − λx , (6.4)subject to boundary conditions ˜ h (0) = e − λx s , (6.5)lim x →∞ ˜ h ( x ) = 0 , (6.6)where the second condition comes from the interpretation of probabilistic represen-tation (6.2) as x → ∞ .We first prove the existence of a unique solution to the boundary value problemand then present the asymptotics of its solution. Let us consider the linear differentialequation x y ′′ + ( cx + d ) y ′ − ( ax + b ) y = 0 , x > , (6.7)where a > b > d < c ∈ R . Based on the discussion of section 2, we knowthat there is a positive number v such that the Wronskian W of y , y is W ( x ) = − vx − c exp dx < . (6.8)20hen we will further consider the inhomogeneous linear differential equation x y ′′ + ( cx + d ) y ′ − ( ax + b ) y = f ( x ) , x > . (6.9)We have the following result. Theorem 6.1.
Let f be a bounded real-valued continuous function on (0 , ∞ ) , and let θ ∈ R . Then there exists a unique solution y of (6.9) with the properties lim x → + y ( x ) = θ, lim x →∞ y ( x ) = 0 . (6.10) Proof.
We first prove uniqueness. Let y, ˜ y be two solutions of (6.9) satisfying (6.10).Then y − ˜ y is a solution of (6.7) with limit 0 as x → + and x → ∞ . Therefore, y − ˜ y must be a multiple of y and y . But y , y are linearly independent by (6.8).Therefore, y = ˜ y .We now prove existence. It is enough to consider θ = 0 because if y is a solutionof (6.9) satisfying (6.10) with θ = 0 then y ( x ) + θ y ( x ) y (0) is the desired solution of (6.9)satisfying (6.10). It is claimed that the solution (for θ = 0) is given by y ( x ) = (cid:18)Z ∞ x y ( t ) f ( t ) t W ( t ) dt (cid:19) y ( x ) + (cid:18)Z x y ( t ) f ( t ) t W ( t ) dt (cid:19) y ( x ) . (6.11)The first integral in (6.11) exists for x >
0. We use that f is bounded, (2.5) and(6.8) to estimate (cid:12)(cid:12)(cid:12)(cid:12) y ( t ) f ( t ) t W ( t ) (cid:12)(cid:12)(cid:12)(cid:12) ≤ Ct − + c exp( − √ at ) , t ≥ t > . (6.12)In this proof C denotes a constant which may have different values in different in-equalities. The second integral in (6.11) exists for x >
0. This follows from (cid:12)(cid:12)(cid:12)(cid:12) y ( t ) f ( t ) t W ( t ) (cid:12)(cid:12)(cid:12)(cid:12) ≤ C, < t < t , (6.13)where we used (2.12) and (6.8). It now follows that y defined by (6.11) is a solutionof (6.9). We show that the first term on the right-hand side of (6.11) converges to 0as x → ∞ . It is easy to find that, for every fixed α ∈ R , Z ∞ u s α e − s ds ≤ Cu α e − u , u ≥ u > . (6.14)If we substitute s = 2 √ at we obtain Z ∞ x t ( α − exp( − √ at ) dt ≤ Cx α exp( − √ ax ) , x ≥ x > . α = − + c and combining with (6.12), Z ∞ x (cid:12)(cid:12)(cid:12)(cid:12) y ( t ) f ( t ) t W ( t ) (cid:12)(cid:12)(cid:12)(cid:12) dt ≤ Cx − + c exp( − √ ax ) , x ≥ x . Following similar calculation to that in the proof of Theorem 2.1, we obtain twoindependent solutions with known asymptotics near x = ∞ , one of which is y deter-mined by the asymptotics in (2.5) and the other, denoted by ˆ y , determined by theasymptotics ˆ y ( x ) ∼ x / − c/ e √ ax ∞ X k =1 A k x − k/ , as x → ∞ , (6.15)with A = 1 and A k ’s determined by the recursive relation A k = ( k + 1 / − c )( k − / c ) − b √ ak A k − + dk A k − , with the understanding that A − = 0. It follows from (2.5), (6.15) that0 < y ( x ) ≤ Cx − c exp(2 √ ax ) , x ≥ x . (6.16)Therefore, Z ∞ x (cid:12)(cid:12)(cid:12)(cid:12) y ( t ) f ( t ) t W ( t ) (cid:12)(cid:12)(cid:12)(cid:12) dt y ( x ) ≤ Cx − , x ≥ x . We show that the second term in (6.11) converges to 0 as x → ∞ . It is easy tosee that, for fixed α ∈ R , Z u s α e s ds ≤ Cu α e u , u ≥ u > . (6.17)If we substitute s = 2 √ at we obtain Z x t ( α − exp(2 √ at ) dt ≤ Cx α exp(2 √ ax ) , x ≥ x . (6.18)From (6.8) and (6.16), we have, (cid:12)(cid:12)(cid:12)(cid:12) y ( t ) f ( t ) t W ( t ) (cid:12)(cid:12)(cid:12)(cid:12) ≤ Ct − + c exp(2 √ at ) , t ≥ t . Therefore, with α = − + c in (6.18), Z x (cid:12)(cid:12)(cid:12)(cid:12) y ( t ) f ( t ) t W ( t ) (cid:12)(cid:12)(cid:12)(cid:12) dt ≤ Cx − + c exp(2 √ ax ) , x ≥ x . Z x (cid:12)(cid:12)(cid:12)(cid:12) y ( t ) f ( t ) t W ( t ) (cid:12)(cid:12)(cid:12)(cid:12) dt y ( x ) ≤ Cx − , x ≥ x . In what follows, we aim to prove that the first term on the right-hand side of(6.11) converges to 0 as x → + . It follows from (2.10), (2.12) that y ( t ) is boundedas t → + . Therefore, we have (cid:12)(cid:12)(cid:12)(cid:12) y ( t ) f ( t ) t W ( t ) (cid:12)(cid:12)(cid:12)(cid:12) ≤ Ct c − exp( − dt ) , < t < t . When we substitute s = − dt , α = − c in (6.17) we obtain Z ∞ x (cid:12)(cid:12)(cid:12)(cid:12) y ( t ) f ( t ) t W ( t ) (cid:12)(cid:12)(cid:12)(cid:12) dt ≤ Cx c exp( − dx ) , < x < x . Therefore, by (2.12), Z ∞ x (cid:12)(cid:12)(cid:12)(cid:12) y ( t ) f ( t ) t W ( t ) (cid:12)(cid:12)(cid:12)(cid:12) dt y ( x ) ≤ Cx , < x < x . Similarly, the second term on the right-hand side of (6.11) converges to 0 as x → + .since from (6.13), we get Z x (cid:12)(cid:12)(cid:12)(cid:12) y ( t ) f ( t ) t W ( t ) (cid:12)(cid:12)(cid:12)(cid:12) dt y ( x ) ≤ Cx, < x < x . This completes the proof. (cid:3)
We will shall the following lemma, whose proof is largely based on integration byparts and hence omitted from the paper.
Lemma 6.1.
Let α ∈ R . We have the following asymptotics Z ∞ u s α e − s ds ∼ u α e − u ∞ X k =0 α k u − k as u → ∞ , (6.19) Z u s α e s ds ∼ u α e u ∞ X k =0 ( − k α k u − k as u → ∞ , (6.20) where α p := α ( α − . . . ( α − p + 1) . Theorem 6.2.
The unique solution y ( x ) of (6.9) , (6.10) with f ( x ) = − u satisfies y ( x ) ∼ ∞ X k =0 C k x k as x → + , (6.21)23 here C = θ , dC = bθ − u and, for k ≥ , kdC k = ( b − c ( k − − ( k − k − C k − + aC k − . (6.22) Moreover, y ( x ) ∼ ∞ X k =0 D k x − k − as x → ∞ , (6.23) where D = ua − and, for k ≥ , aD k = ( k ( k + 1) − ck − b ) D k − + d ( k − D k − , D − := 0 . (6.24) Proof.
We consider first the case θ = 0, u = 1. Then the solution y ( x ) is given by(6.11) with f ( x ) = − f ( x ) = − y ( t ) = K − X k =0 E k t k + O ( t K ) as t → + . (6.25)Therefore, y ( t )( − t W ( t ) = v − t c − exp( − dt ) K − X k =0 E k t k + O ( t K ) ! as t → + . Now we use (6.20) with the substitution s = − dt . Then we obtain Z ∞ x y ( t )( − t W ( t ) dt = x c exp( − dx ) K − X k =0 F k x k + O ( x K ) ! as x → + . Multiplying this integral by y ( x ) and using (2.12), we obtain the desired asymptoticexpansion.We show that the second term on the right-hand side of (6.11) with f ( x ) = − y ( t )( − t W ( t ) = v − K − X k =0 A k t k + O ( t K ) ! as t → + . Integrating this equation on both sides from t = 0 to t = x , multiplying by y ( x ) andusing (6.25) we arrive at the desired asymptotic expansion.24y differentiating (6.11), we show that y ′ ( x ) and y ′′ ( x ) have asymptotic expansionsas x → + which are obtained by differentiating the asymptotic expansion of y ( x )term-by-term. If y ( x ) is the solution for θ = 0 and u = 1, then the solution ˜ y ( x ) forreal θ, u is ˜ y ( x ) = uy ( x ) + θ y ( x ) y (0) . (6.26)Therefore, ˜ y ( x ) also admits an asymptotic expansions of the form (6.21), and it may bedifferentiated. Substituting the expansions for ˜ y ( x ) , ˜ y ′ ( x ) , ˜ y ′′ ( x ) in (6.9) with f ( x ) = − u and comparing coefficients, we obtain the stated recursion for the coefficients C k .Now we show that the first term g ( t ) on the right-hand side of (6.11) with f ( x ) = − g ( t ) ∼ ∞ X k =2 E k x − k/ as x → ∞ . (6.27)We know from (2.5) that y ( t ) = t − c exp( − √ at ) K − X k =0 A k t − k/ + O ( t − K/ ) ! as t → ∞ .Therefore, y ( t )( − t W ( t ) = t − + c exp( − √ at ) K − X k =0 F k t − k/ + O ( t − K/ ) ! as t → ∞ .We now use (6.19) with the substitution s = 2 √ at . Then we obtain Z ∞ x y ( t )( − t W ( t ) dt = x − + c exp( − √ ax ) K − X k =0 G k x − k/ + O ( x − K/ ) ! as x → ∞ . If we multiply this integral by y ( x ) and use (2.5), (6.15), we obtain an asymptoticexpansion of the desired form (6.27).In a similar way we show that the second term on the right-hand side of (6.11)with f ( x ) = − y ( x ) also admits an asymptotic expansion of the form (6.27), and we obtainasymptotic expansions for ˜ y ′ ( x ) , ˜ y ′′ ( x ) by differentiating the asymptotic expansion for˜ y ( x ) term-by-term. By substituting these expansions in (6.9) with f ( x ) = − u andcompare the coefficients, we notice that the coefficients of the terms x − k/ with odd k must vanish. Thus we obtain the recursion (6.24) and the initial value D . (cid:3) T -period Asian call option with the strike price of K , i.e.AC := E x (cid:20) e − rT (cid:18) T Z T X t d t − K (cid:19) + (cid:21) = 1 T e − rT E x [( Y T − K ∗ ) + ] , where Y T = Z T X t d t, K ∗ = KT > . For notational brevity, we assume that the expectation is taken under the risk-neutralprobability measure for which P [ X = x , Y = 0] = 1 . We can compute the following Laplace transform w.r.t. K ∗ . Since the integrands arenonnegative, we exchange the order of integration and obtain Z ∞ e − wK ∗ E x [( Y T − K ∗ ) + ] d K ∗ = E x (cid:20)Z Y T e − wK ∗ ( Y T − K ∗ ) d K ∗ (cid:21) = 1 w E x ( Y T ) − w + 1 w E x [ e − wY T ] . (6.28)An analytic expression for the first term in (6.28) is already obtained in Feng andVolkmer [9, Proposition 3.4]. Therefore, the only unknown quantity to be determinedis the third term in (6.28). Once efficient algorithms for computing both terms areobtained, we would invert the Laplace transform with respect to K ∗ to determinethe price of Asian option. One can easily find that it is a special case of h definedby (6.2) by taking λ = 0. Since the rest of computation for pricing Asian options isirrelevant to the discussion of asymptotics of Heun equation, we shall only focus onthe computation of ˜ h , i.e. the Laplace transform of h , which can be derived by letting λ = 0 in (6.4). In Feng and Volkmer [9], this quantity was represented as the time-integral of the process X t up to the earlier of the first time it hits zero and a fixed time T , i.e. E hR τ ∧ T X t d t i , where τ := inf { t : X t ≤ } . f ( x ) = − u andboundary conditions (6.10), a = 2 wσ , b = 2 sσ , c = 2 rσ , d = − δσ , u = 2 σ , θ = 1 s . Although a numerical algorithm can be applied directly to solve the boundaryvalue problem (6.4) with (6.5) and (6.6), one would have to truncate the domain(0 , ∞ ) to some finite interval for practical reason. Then the asymptotics of the par-ticular solution in (6.21) and (6.23) become useful to determine values at boundarypoints.Here we provide a numerical example to show the solution to the boundary valueproblem using the asymptotics. The following set of parameters are used for compu-tation. s = 2 , w = 1 , σ = 0 . , r = 0 . , δ = 0 . . To avoid the singularities x = 0 and x = ∞ , we consider the left-end-point to be x = 0 .
01 and the right-end-point to be x = 100. We use 30 terms of the asymp-totic formulas (6.21) and (6.23) to find ˜ h ( x ) = 0 . h ( x ) =0 . . Finally, we use Maple’s own default BVP solver to determinethe numerical solution, which is shown in Figure 2.
Figure 2: Solution to the boundary value problem (6.4). eferencesReferences [1] J. Abate and W. Whitt. Numerical inversion of Laplace transform of probabilityfunctions. OSRA Journal on Computing 7 (1995), 36–43.[2] J. Abate and W. Whitt. A unified framework for numerically inverting Laplacetransforms. INFORMS J. Comput. 18 (2006), 408–421.[3] P. Barrieu, A. Rouault and M. Yor. A study of the Hartman-Watson distributionmotivated by numerical problems related to the pricing of Asian options. J. Appl.Probab. 41 (2004), 1049–1058.[4] A. Borodin and P. Salminen. Handbook of Brownian Motion - Facts and Formu-lae. (2002) Second Edition. Springer, Basel.[5] M. Chugunova and H. Volkmer. Study of a differential operator of Heun typearising in fluid dynamics. Stud. Appl. Math. 123 (2009), 291–309.[6] R. Feng and X. Jing. Analytical valuation and hedging of variable annuity guar-anteed lifetime withdrawal benefits. Insurance: Mathematics and Economics ,(2016), 72, 36–48.[7] R. Feng and H.W. Volkmer. Analytical calculation of risk measures for variableannuity guaranteed benefits, Insurance: Mathematics and Economics, (2012) 51(3), 636–648.[8] R. Feng and H. Volkmer. Spectral methods for the calculation of risk measuresfor variable annuity guaranteed benefits. ASTIN Bulletin (2014), 44(3), 653–681.[9] R. Feng and H. Volkmer. An identity of hitting times and its application to thepricing of guaranteed minimum withdrawal benefit. Mathematics and FinancialEconomics. (2016), 10(2), 127-149.[10] T. Gard. Introduction to Stochastic Differential Equations. Dekker, New York.(1988).[11] I. S. Gradshteyn and I. M. Ryzhik. Table of integrals, series and products. Else-vier/Academic Press, seventh edition, Amsterdam 2007.[12] M. Hairer. On Malliavin’s proof of H´’ormander’s theorem.
Acta Math.119